working on qsUGF2

This commit is contained in:
Pierre-Francois Loos 2021-03-07 23:04:47 +01:00
parent 29ed74cba0
commit 2d826d0f1e
6 changed files with 70 additions and 30 deletions

View File

@ -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

View File

@ -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)') &

View File

@ -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(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'

View File

@ -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
!------------------------------------------------------------------------

View File

@ -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
!------------------------------------------------------------------------

View File

@ -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)