From 2d826d0f1e3f8f65b6ff38035c733e4d7f5137b4 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sun, 7 Mar 2021 23:04:47 +0100 Subject: [PATCH] working on qsUGF2 --- input/methods | 10 +++++----- src/GF/print_qsGF2.f90 | 24 ++++++++++-------------- src/GF/print_qsUGF2.f90 | 4 ++-- src/GF/qsGF2.f90 | 40 +++++++++++++++++++++++++++++++++++++--- src/GF/qsUGF2.f90 | 8 ++++---- src/QuAcK/QuAcK.f90 | 14 ++++++++++++-- 6 files changed, 70 insertions(+), 30 deletions(-) diff --git a/input/methods b/input/methods index d790e69..f719fe3 100644 --- a/input/methods +++ b/input/methods @@ -1,7 +1,7 @@ # RHF UHF KS MOM - T F F F + F T F F # MP2* MP3 MP2-F12 - T F F + F F F # CCD DCD CCSD CCSD(T) F F F F # drCCD rCCD lCCD pCCD @@ -10,10 +10,10 @@ F F F F # RPA* RPAx* ppRPA F F F -# G0F2 evGF2 qsGF2 G0F3 evGF3 - F F F F F +# G0F2 evGF2 qsGF2* G0F3 evGF3 + F F T F F # G0W0* evGW* qsGW* - F F T + F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/src/GF/print_qsGF2.f90 b/src/GF/print_qsGF2.f90 index 02c694d..5e3faf9 100644 --- a/src/GF/print_qsGF2.f90 +++ b/src/GF/print_qsGF2.f90 @@ -1,4 +1,5 @@ -subroutine print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,ENuc,P,T,V,J,K,F,SigC,Z,EqsGF2,Ec,dipole) +subroutine print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,P,T,V,J,K,F,SigC,Z, & + ENuc,ET,EV,EJ,Ex,Ec,EqsGF2,dipole) ! Print one-electron energies and other stuff for qsGF2 @@ -20,40 +21,35 @@ subroutine print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,ENuc,P,T,V,J,K,F,SigC double precision,intent(in) :: T(nBas,nBas),V(nBas,nBas) double precision,intent(in) :: J(nBas,nBas),K(nBas,nBas),F(nBas,nBas) double precision,intent(in) :: Z(nBas),SigC(nBas,nBas) + double precision,intent(in) :: ET + double precision,intent(in) :: EV + double precision,intent(in) :: EJ + double precision,intent(in) :: Ex double precision,intent(in) :: Ec + double precision,intent(in) :: EqsGF2 double precision,intent(in) :: dipole(ncart) ! Local variables integer :: q,ixyz,HOMO,LUMO - double precision :: Gap,ET,EV,EJ,Ex + double precision :: Gap double precision,external :: trace_matrix ! Output variables - double precision,intent(out) :: EqsGF2 - ! HOMO and LUMO HOMO = nO LUMO = HOMO + 1 Gap = eGF2(LUMO)-eGF2(HOMO) -! Compute energies - - ET = trace_matrix(nBas,matmul(P,T)) - EV = trace_matrix(nBas,matmul(P,V)) - EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) - Ex = 0.25d0*trace_matrix(nBas,matmul(P,K)) - EqsGF2 = ET + EV + EJ + Ex - ! Dump results write(*,*)'-------------------------------------------------------------------------------' if(nSCF < 10) then - write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent qsG',nSCF,'W',nSCF,' calculation' + write(*,'(1X,A21,I1,A2,A12)')' Self-consistent qsG',nSCF,'F2',' calculation' else - write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent qsG',nSCF,'W',nSCF,' calculation' + write(*,'(1X,A21,I2,A2,A12)')' Self-consistent qsG',nSCF,'F2',' calculation' endif 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)') & diff --git a/src/GF/print_qsUGF2.f90 b/src/GF/print_qsUGF2.f90 index 5b825e9..4db8f23 100644 --- a/src/GF/print_qsUGF2.f90 +++ b/src/GF/print_qsUGF2.f90 @@ -69,9 +69,9 @@ subroutine print_qsUGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,cGF2,PGF2,Ov,T,V,J,K, write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' if(nSCF < 10) then - write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent qsG',nSCF,'W',nSCF,' calculation' + write(*,'(1X,A21,I1,A2,A12)')' Self-consistent qsG',nSCF,'F2',' calculation' else - write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent qsG',nSCF,'W',nSCF,' calculation' + write(*,'(1X,A21,I2,A2,A12)')' Self-consistent qsG',nSCF,'F2',' calculation' endif write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' diff --git a/src/GF/qsGF2.f90 b/src/GF/qsGF2.f90 index 5f810cf..2f39f70 100644 --- a/src/GF/qsGF2.f90 +++ b/src/GF/qsGF2.f90 @@ -52,6 +52,10 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, double precision :: rcond double precision,external :: trace_matrix double precision :: dipole(ncart) + double precision :: ET + double precision :: EV + double precision :: EJ + double precision :: Ex double precision :: Ec double precision :: EcBSE(nspin) @@ -175,11 +179,41 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) - ! Print results + !------------------------------------------------------------------------ + ! Compute total energy + !------------------------------------------------------------------------ + + ! Kinetic energy + + ET = trace_matrix(nBas,matmul(P,T)) + + ! Potential energy + + EV = trace_matrix(nBas,matmul(P,V)) + + ! Coulomb energy + + EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) + + ! Exchange energy + + Ex = -0.25d0*trace_matrix(nBas,matmul(P,K)) + + ! Total energy + + EqsGF2 = ET + EV + EJ + Ex + + ! Correlation energy + + call MP2(nBas,nC,nO,nV,nR,ERI_MO,ENuc,EqsGF2,eGF2,Ec) + + !------------------------------------------------------------------------ + ! Print results + !------------------------------------------------------------------------ - call MP2(nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,eGF2,Ec) call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,ENuc,P,T,V,J,K,F,SigCp,Z,EqsGF2,Ec,dipole) + call print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,P,T,V,J,K,F,SigCp,Z, & + ENuc,ET,EV,EJ,Ex,Ec,EqsGF2,dipole) enddo !------------------------------------------------------------------------ diff --git a/src/GF/qsUGF2.f90 b/src/GF/qsUGF2.f90 index 93d2833..eaa0dbb 100644 --- a/src/GF/qsUGF2.f90 +++ b/src/GF/qsUGF2.f90 @@ -268,14 +268,14 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved, Ex(is) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,is),K(:,:,is))) end do - ! Correlation energy - - call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,eGF2,Ec) - ! Total energy EqsGF2 = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) + ! Correlation energy + + call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EqsGF2,eGF2,Ec) + !------------------------------------------------------------------------ ! Print results !------------------------------------------------------------------------ diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 9647ef1..64e4819 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -848,8 +848,18 @@ program QuAcK call cpu_time(start_GF2) - call qsGF2(maxSCF_GF,thresh_GF,n_diis_GF,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GF,nNuc,ZNuc,rNuc,ENuc, & - nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + if(unrestricted) then + + call qsUGF2(maxSCF_GF,thresh_GF,n_diis_GF,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GF, & + nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO, & + ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) + + else + + call qsGF2(maxSCF_GF,thresh_GF,n_diis_GF,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GF,nNuc,ZNuc,rNuc,ENuc, & + nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + + end if call cpu_time(end_GF2)