From 0e432942113e5d21f8497d16d077cfaee206953a Mon Sep 17 00:00:00 2001 From: pfloos Date: Tue, 17 Sep 2024 20:13:18 +0200 Subject: [PATCH] clean up in ccG0W0 --- src/AOtoMO/Hartree_matrix_AO_basis.f90 | 4 +- src/GW/RGW.f90 | 1 - src/GW/ccRG0W0.f90 | 80 ++++++++++++-------------- src/HF/RHF.f90 | 38 ++++++------ 4 files changed, 59 insertions(+), 64 deletions(-) diff --git a/src/AOtoMO/Hartree_matrix_AO_basis.f90 b/src/AOtoMO/Hartree_matrix_AO_basis.f90 index e53efe4..03000a0 100644 --- a/src/AOtoMO/Hartree_matrix_AO_basis.f90 +++ b/src/AOtoMO/Hartree_matrix_AO_basis.f90 @@ -21,10 +21,10 @@ subroutine Hartree_matrix_AO_basis(nBas,P,G,H) H(:,:) = 0d0 - do mu=1,nBas + do si=1,nBas do nu=1,nBas do la=1,nBas - do si=1,nBas + do mu=1,nBas H(mu,nu) = H(mu,nu) + P(la,si)*G(mu,la,nu,si) end do end do diff --git a/src/GW/RGW.f90 b/src/GW/RGW.f90 index da23c82..c038cad 100644 --- a/src/GW/RGW.f90 +++ b/src/GW/RGW.f90 @@ -167,7 +167,6 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii call wall_time(start_GW) call ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) -! call ccRG0W0_mat(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW diff --git a/src/GW/ccRG0W0.f90 b/src/GW/ccRG0W0.f90 index d8a3a5a..a599719 100644 --- a/src/GW/ccRG0W0.f90 +++ b/src/GW/ccRG0W0.f90 @@ -31,9 +31,7 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH integer :: nSCF double precision :: Conv - double precision :: x - - double precision,allocatable :: eGW(:) + double precision,allocatable :: Sig(:) double precision,allocatable :: Z(:) double precision,allocatable :: del(:,:,:) @@ -55,50 +53,29 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH write(*,*)'*****************************' write(*,*) -! Form energy denominator and guess amplitudes +! Memory allocation allocate(del(nO,nV,nOrb)) allocate(vec(nO,nV,nOrb)) allocate(res(nO,nV,nOrb)) allocate(amp(nO,nV,nOrb)) - allocate(eGW(nOrb),Z(nOrb)) + allocate(Sig(nOrb)) + allocate(Z(nOrb)) - allocate(r_diis(nO*nV*nOrb,max_diis)) - allocate(t_diis(nO*nV*nOrb,max_diis)) + allocate(r_diis(nO*nV*nOrb,max_diis)) + allocate(t_diis(nO*nV*nOrb,max_diis)) ! Initialization - eGW(:) = eHF(:) + Sig(:) = 0d0 Z(:) = 1d0 - ! Compute energy differences - - do i=nC+1,nO - do j=nC+1,nO - do a=1,nV-nR - - del(i,a,j) = eHF(i) + eHF(j) - eHF(nO+a) - - end do - end do - end do - - do i=nC+1,nO - do a=1,nV-nR - do b=1,nV-nR - - del(i,a,nO+b) = eHF(nO+a) + eHF(nO+b) - eHF(i) - - end do - end do - end do - !-------------------------! ! Main loop over orbitals ! !-------------------------! - do p=nC+1,nO+1 + do p=nO,nO+1 ! Initialization @@ -114,6 +91,26 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH res(:,:,:) = 0d0 ! Compute energy differences + + do i=nC+1,nO + do j=nC+1,nO + do a=1,nV-nR + + del(i,a,j) = eHF(i) + eHF(j) - eHF(nO+a) - eHF(p) + + end do + end do + end do + + do i=nC+1,nO + do a=1,nV-nR + do b=1,nV-nR + + del(i,a,nO+b) = eHF(nO+a) + eHF(nO+b) - eHF(i) - eHF(p) + + end do + end do + end do do i=nC+1,nO do a=1,nV-nR @@ -144,7 +141,7 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH write(*,*)'| CC-based G0W0 calculation |' write(*,*)'-------------------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & - '|','#','|','HF','|','G0W0','|','Conv','|' + '|','#','|','Sig_c (eV)','|','e_GW (eV)','|','Conv','|' write(*,*)'-------------------------------------------------------------' do while(Conv > thresh .and. nSCF < maxSCF) @@ -155,7 +152,7 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH ! Compute residual for 2h1p sector - res(:,:,:) = vec(:,:,:) + (del(:,:,:) - eGW(p))*amp(:,:,:) + res(:,:,:) = vec(:,:,:) + (del(:,:,:) - Sig(p))*amp(:,:,:) do i=nC+1,nO do a=1,nV-nR @@ -197,7 +194,7 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH ! Update amplitudes - amp(:,:,:) = amp(:,:,:) - res(:,:,:)/(del(:,:,:) - eHF(p)) + amp(:,:,:) = amp(:,:,:) - res(:,:,:)/del(:,:,:) ! DIIS extrapolation @@ -210,13 +207,13 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH ! Compute quasiparticle energy - eGW(p) = eHF(p) + Sig(p) = 0d0 do i=nC+1,nO do a=1,nV-nR do q=nC+1,nOrb-nR - eGW(p) = eGW(p) + vec(i,a,q)*amp(i,a,q) + Sig(p) = Sig(p) + vec(i,a,q)*amp(i,a,q) end do end do @@ -225,11 +222,12 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH ! Dump results write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X)') & - '|',nSCF,'|',eHF(p)*HaToeV,'|',eGW(p)*HaToeV,'|',Conv,'|' + '|',nSCF,'|',Sig(p)*HaToeV,'|',(eHF(p)+Sig(p))*HaToeV,'|',Conv,'|' end do write(*,*)'-------------------------------------------------------------' + write(*,*) !------------------------------------------------------------------------ ! End of SCF loop !------------------------------------------------------------------------ @@ -244,19 +242,17 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - stop - end if write(*,*)'-------------------------------------------------------------------------------' - write(*,*)' CC-G0W0 calculation ' + write(*,*)'| CC-based G0W0 calculation |' write(*,*)'-------------------------------------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & - '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' + '|','Orb','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' write(*,*)'-------------------------------------------------------------------------------' write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X)') & - '|',p,'|',eHF(p)*HaToeV,'|',(eGW(p)-eHF(p))*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' + '|',p,'|',eHF(p)*HaToeV,'|',Sig(p)*HaToeV,'|',Z(p),'|',(eHF(p)+Sig(p))*HaToeV,'|' write(*,*)'-------------------------------------------------------------------------------' end do diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index e66de73..61510d4 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -90,10 +90,10 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN call mo_guess(nBas,nOrb,guess_type,S,Hc,X,c) - !P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO))) - call dgemm('N', 'T', nBas, nBas, nO, 2.d0, & - c(1,1), nBas, c(1,1), nBas, & - 0.d0, P(1,1), nBas) + P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO))) +! call dgemm('N', 'T', nBas, nBas, nO, 2.d0, & +! c(1,1), nBas, c(1,1), nBas, & +! 0.d0, P(1,1), nBas) ! Initialization @@ -127,10 +127,10 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN call exchange_matrix_AO_basis(nBas,P,ERI,K) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - if(nBas .ne. nOrb) then - call AOtoMO(nBas, nOrb, c(1,1), F(1,1), Fp(1,1)) - call MOtoAO(nBas, nOrb, S(1,1), c(1,1), Fp(1,1), F(1,1)) - endif +! if(nBas .ne. nOrb) then +! call AOtoMO(nBas, nOrb, c(1,1), F(1,1), Fp(1,1)) +! call MOtoAO(nBas, nOrb, S(1,1), c(1,1), Fp(1,1), F(1,1)) +! endif ! Check convergence @@ -174,24 +174,24 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN ! Diagonalize Fock matrix - if(nBas .eq. nOrb) then +! if(nBas .eq. nOrb) then Fp = matmul(transpose(X),matmul(F,X)) cp(:,:) = Fp(:,:) call diagonalize_matrix(nOrb,cp,eHF) c = matmul(X,cp) - else - Fp = matmul(transpose(c),matmul(F,c)) - cp(:,:) = Fp(:,:) - call diagonalize_matrix(nOrb,cp,eHF) - c = matmul(c,cp) - endif +! else +! Fp = matmul(transpose(c),matmul(F,c)) +! cp(:,:) = Fp(:,:) +! call diagonalize_matrix(nOrb,cp,eHF) +! c = matmul(c,cp) +! endif ! Density matrix - !P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO))) - call dgemm('N', 'T', nBas, nBas, nO, 2.d0, & - c(1,1), nBas, c(1,1), nBas, & - 0.d0, P(1,1), nBas) + P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO))) +! call dgemm('N', 'T', nBas, nBas, nO, 2.d0, & +! c(1,1), nBas, c(1,1), nBas, & +! 0.d0, P(1,1), nBas) ! Dump results