From a8b51ab800a1c2663c926f5b1025fd5f8ba6fb3b Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Thu, 12 Dec 2024 19:25:33 +0100 Subject: [PATCH] HF on HPC branch: OK --- src/HF/RHF_hpc.f90 | 169 ++++++++++++++++++++++++++------------------- 1 file changed, 98 insertions(+), 71 deletions(-) diff --git a/src/HF/RHF_hpc.f90 b/src/HF/RHF_hpc.f90 index 740bc94..a4d3325 100644 --- a/src/HF/RHF_hpc.f90 +++ b/src/HF/RHF_hpc.f90 @@ -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 @@ -93,6 +94,9 @@ 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) @@ -100,70 +104,60 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh 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 @@ -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 @@ -219,7 +234,7 @@ 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 @@ -227,10 +242,20 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh 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, & @@ -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 @@ -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