Skip to content

Commit

Permalink
HF on HPC branch: OK
Browse files Browse the repository at this point in the history
  • Loading branch information
@AbdAmmar
AbdAmmar committed Dec 12, 2024
1 parent 84a989b commit a8b51ab
Showing 1 changed file with 98 additions and 71 deletions.
169 changes: 98 additions & 71 deletions src/HF/RHF_hpc.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,7 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
double precision,allocatable :: Fp(:,:)
double precision,allocatable :: ERI_chem(:)
double precision,allocatable :: ERI_phys(:,:,:,:), J_deb(:,:), K_deb(:,:)
double precision,allocatable :: tmp1(:,:), FX(:,:)


! Output variables
Expand Down Expand Up @@ -93,77 +94,70 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
allocate(err_diis(nBas_Sq,max_diis))
allocate(F_diis(nBas_Sq,max_diis))

allocate(tmp1(nBas,nBas))
allocate(FX(nBas,nOrb))

! Guess coefficients and density matrix
call mo_guess(nBas,nOrb,guess_type,S,Hc,X,c)

call dgemm('N', 'T', nBas, nBas, nO, 2.d0, &
c(1,1), nBas, c(1,1), nBas, &
0.d0, P(1,1), nBas)


ERI_size = (nBas * (nBas + 1)) / 2
ERI_size = (ERI_size * (ERI_size + 1)) / 2
ERI_size = shiftr(nBas * (nBas + 1), 1)
ERI_size = shiftr(ERI_size * (ERI_size + 1), 1)
allocate(ERI_chem(ERI_size))
call read_2e_integrals_hpc(working_dir, ERI_size, ERI_chem)

call wall_time(t1)
call Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, J)
call wall_time(t2)
print*, " J built in (sec):", (t2-t1)

call wall_time(t1)
call exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K)
call wall_time(t2)
print*, " K built in (sec):", (t2-t1)


allocate(ERI_phys(nBas,nBas,nBas,nBas))
allocate(J_deb(nBas,nBas))
allocate(K_deb(nBas,nBas))

call read_2e_integrals(working_dir, nBas, ERI_phys)

call wall_time(t1)
call Hartree_matrix_AO_basis(nBas, P, ERI_phys, J_deb)
call wall_time(t2)
print*, " J_deb built in (sec):", (t2-t1)

call wall_time(t1)
call exchange_matrix_AO_basis(nBas, P, ERI_phys, K_deb)
call wall_time(t2)
print*, " K_deb built in (sec):", (t2-t1)

print*, "max error on J = ", maxval(dabs(J - J_deb))
diff = 0.d0
do ii = 1, nBas
do jj = 1, nBas
diff_loc = dabs(J(jj,ii) - J_deb(jj,ii))
if(diff_loc .gt. 1d-10) then
print*, 'error in J on: ', jj, ii
print*, J(jj,ii), J_deb(jj,ii)
stop
endif
diff = diff + diff_loc
enddo
enddo
print*, 'total diff on J = ', diff

print*, "max error on K = ", maxval(dabs(K - K_deb))
diff = 0.d0
do ii = 1, nBas
do jj = 1, nBas
diff_loc = dabs(K(jj,ii) - K_deb(jj,ii))
if(diff_loc .gt. 1d-10) then
print*, 'error in K on: ', jj, ii
print*, K(jj,ii), K_deb(jj,ii)
stop
endif
diff = diff + diff_loc
enddo
enddo
print*, 'total diff on K = ', diff

stop
!call wall_time(t1)
!call Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, J)
!call wall_time(t2)
!print*, " J built in (sec):", (t2-t1)
!call wall_time(t1)
!call exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K)
!call wall_time(t2)
!print*, " K built in (sec):", (t2-t1)
!allocate(ERI_phys(nBas,nBas,nBas,nBas))
!allocate(J_deb(nBas,nBas))
!allocate(K_deb(nBas,nBas))
!call read_2e_integrals(working_dir, nBas, ERI_phys)
!call wall_time(t1)
!call Hartree_matrix_AO_basis(nBas, P, ERI_phys, J_deb)
!call wall_time(t2)
!print*, " J_deb built in (sec):", (t2-t1)
!call wall_time(t1)
!call exchange_matrix_AO_basis(nBas, P, ERI_phys, K_deb)
!call wall_time(t2)
!print*, " K_deb built in (sec):", (t2-t1)
!print*, "max error on J = ", maxval(dabs(J - J_deb))
!diff = 0.d0
!do ii = 1, nBas
! do jj = 1, nBas
! diff_loc = dabs(J(jj,ii) - J_deb(jj,ii))
! if(diff_loc .gt. 1d-10) then
! print*, 'error in J on: ', jj, ii
! print*, J(jj,ii), J_deb(jj,ii)
! stop
! endif
! diff = diff + diff_loc
! enddo
!enddo
!print*, 'total diff on J = ', diff
!print*, "max error on K = ", maxval(dabs(K - K_deb))
!diff = 0.d0
!do ii = 1, nBas
! do jj = 1, nBas
! diff_loc = dabs(K(jj,ii) - K_deb(jj,ii))
! if(diff_loc .gt. 1d-10) then
! print*, 'error in K on: ', jj, ii
! print*, K(jj,ii), K_deb(jj,ii)
! stop
! endif
! diff = diff + diff_loc
! enddo
!enddo
!print*, 'total diff on K = ', diff
!stop

! Initialization

Expand Down Expand Up @@ -192,25 +186,46 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
nSCF = nSCF + 1

! Build Fock matrix
call Hartree_matrix_AO_basis(nBas,P,ERI_phys,J)
call exchange_matrix_AO_basis(nBas,P,ERI_phys,K)
call Hartree_matrix_AO_basis_hpc (nBas, ERI_size, P(1,1), ERI_chem(1), J(1,1))
call exchange_matrix_AO_basis_hpc(nBas, ERI_size, P(1,1), ERI_chem(1), K(1,1))
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)

! Check convergence
err = matmul(F, matmul(P, S)) - matmul(matmul(S, P), F)
call dgemm("N", "N", nBas, nBas, nBas, 1.d0, &
S(1,1), nBas, P(1,1), nBas, &
0.d0, tmp1(1,1), nBas)
call dgemm("N", "N", nBas, nBas, nBas, 1.d0, &
tmp1(1,1), nBas, F(1,1), nBas, &
0.d0, err(1,1), nBas)
call dgemm("N", "T", nBas, nBas, nBas, 1.d0, &
F(1,1), nBas, tmp1(1,1), nBas, &
-1.d0, err(1,1), nBas)
!err = matmul(F, matmul(P, S)) - matmul(matmul(S, P), F)
if(nSCF > 1) Conv = maxval(abs(err))

! Kinetic energy
ET = trace_matrix(nBas, matmul(P, T))
call dgemm("N", "N", nBas, nBas, nBas, 1.d0, &
P(1,1), nBas, T(1,1), nBas, &
0.d0, tmp1(1,1), nBas)
ET = trace_matrix(nBas, tmp1(1,1))

! Potential energy
EV = trace_matrix(nBas, matmul(P, V))
call dgemm("N", "N", nBas, nBas, nBas, 1.d0, &
P(1,1), nBas, V(1,1), nBas, &
0.d0, tmp1(1,1), nBas)
EV = trace_matrix(nBas, tmp1(1,1))

! Hartree energy
EJ = 0.5d0*trace_matrix(nBas, matmul(P, J))
call dgemm("N", "N", nBas, nBas, nBas, 1.d0, &
P(1,1), nBas, J(1,1), nBas, &
0.d0, tmp1(1,1), nBas)
EJ = 0.5d0*trace_matrix(nBas, tmp1(1,1))

! Exchange energy
EK = 0.25d0*trace_matrix(nBas, matmul(P, K))
call dgemm("N", "N", nBas, nBas, nBas, 1.d0, &
P(1,1), nBas, K(1,1), nBas, &
0.d0, tmp1(1,1), nBas)
EK = 0.25d0*trace_matrix(nBas, tmp1(1,1))

! Total energy
ERHF = ET + EV + EJ + EK
Expand All @@ -219,18 +234,28 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F)
endif
endif

! Level shift
if(level_shift > 0d0 .and. Conv > thresh) then
call level_shifting(level_shift,nBas,nOrb,nO,S,c,F)
endif

! Diagonalize Fock matrix
Fp = matmul(transpose(X), matmul(F, X))
call dgemm("N", "N", nBas, nOrb, nBas, 1.d0, &
F(1,1), nBas, X(1,1), nBas, &
0.d0, FX(1,1), nBas)
call dgemm("T", "N", nOrb, nOrb, nBas, 1.d0, &
X(1,1), nBas, FX(1,1), nBas, &
0.d0, Fp(1,1), nOrb)
!Fp = matmul(transpose(X), matmul(F, X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nOrb,cp,eHF)
c = matmul(X,cp)
!c = matmul(X, cp)
call dgemm("N", "N", nBas, nOrb, nOrb, 1.d0, &
X(1,1), nBas, cp(1,1), nOrb, &
0.d0, c(1,1), nBas)


! Density matrix
call dgemm('N', 'T', nBas, nBas, nO, 2.d0, &
Expand Down Expand Up @@ -258,6 +283,7 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
write(*,*)

deallocate(J,K,err,cp,Fp,err_diis,F_diis)
deallocate(tmp1, FX, ERI_chem)

stop

Expand All @@ -280,5 +306,6 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
end if

deallocate(J,K,err,cp,Fp,err_diis,F_diis)
deallocate(tmp1, FX, ERI_chem)

end subroutine

0 comments on commit a8b51ab

Please sign in to comment.