mirror of https://github.com/pfloos/quack
GMP3 and GG0F2
This commit is contained in:
Pierre-Francois Loos
parent
c7093e0c2c
commit
a60aa8d8ee
|
|
@ -1,7 +1,7 @@
|
|||
# RHF UHF GHF ROHF
|
||||
F F T F
|
||||
# MP2 MP3
|
||||
F F
|
||||
T T
|
||||
# CCD pCCD DCD CCSD CCSD(T)
|
||||
F F F F F
|
||||
# drCCD rCCD crCCD lCCD
|
||||
|
|
@ -13,6 +13,6 @@
|
|||
# G0F2 evGF2 qsGF2 G0F3 evGF3
|
||||
F F F F F
|
||||
# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
|
||||
T F F F F F
|
||||
F F F F F F
|
||||
# G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
|
||||
F F F F F F
|
||||
|
|
|
|||
|
|
@ -0,0 +1,120 @@
|
|||
subroutine GG0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, &
|
||||
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
|
||||
|
||||
! Perform a one-shot second-order Green function calculation
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: dophBSE
|
||||
logical,intent(in) :: doppBSE
|
||||
logical,intent(in) :: TDA
|
||||
logical,intent(in) :: dBSE
|
||||
logical,intent(in) :: dTDA
|
||||
logical,intent(in) :: linearize
|
||||
double precision,intent(in) :: eta
|
||||
logical,intent(in) :: regularize
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nS
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: ERHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: Ec
|
||||
double precision :: EcBSE(nspin)
|
||||
double precision,allocatable :: eGFlin(:)
|
||||
double precision,allocatable :: eGF(:)
|
||||
double precision,allocatable :: SigC(:)
|
||||
double precision,allocatable :: Z(:)
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'************************************************'
|
||||
write(*,*)'| One-shot second-order Green function |'
|
||||
write(*,*)'************************************************'
|
||||
write(*,*)
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(SigC(nBas),Z(nBas),eGFlin(nBas),eGF(nBas))
|
||||
|
||||
! Frequency-dependent second-order contribution
|
||||
|
||||
if(regularize) then
|
||||
|
||||
! call GF2_reg_self_energy_diag(eta,nBas,nC,nO,nV,nR,eHF,ERI,SigC,Z)
|
||||
|
||||
else
|
||||
|
||||
call GGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,eHF,ERI,SigC,Z)
|
||||
|
||||
end if
|
||||
|
||||
eGFlin(:) = eHF(:) + Z(:)*SigC(:)
|
||||
|
||||
if(linearize) then
|
||||
|
||||
write(*,*) '*** Quasiparticle energies obtained by linearization ***'
|
||||
|
||||
eGF(:) = eGFlin(:)
|
||||
|
||||
else
|
||||
|
||||
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
|
||||
write(*,*)
|
||||
|
||||
call GGF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,ERI,eGFlin,eHF,eGF,Z)
|
||||
|
||||
end if
|
||||
|
||||
! Print results
|
||||
|
||||
call GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF,Ec)
|
||||
call print_G0F2(nBas,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec)
|
||||
|
||||
! Perform BSE2 calculation
|
||||
|
||||
! if(dophBSE) then
|
||||
!
|
||||
! call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
|
||||
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (singlet) =',EcBSE(1)
|
||||
! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (triplet) =',EcBSE(2)
|
||||
! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy =',sum(EcBSE(:))
|
||||
! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 total energy =',ENuc + EHF + sum(EcBSE(:))
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
|
||||
! end if
|
||||
|
||||
! Perform ppBSE2 calculation
|
||||
|
||||
! if(doppBSE) then
|
||||
!
|
||||
! call GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
|
||||
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (singlet) =',EcBSE(1),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (triplet) =',3d0*EcBSE(2),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy =',EcBSE(1) + 3d0*EcBSE(2),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2),' au'
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
|
||||
! end if
|
||||
|
||||
end subroutine
|
||||
|
|
@ -0,0 +1,109 @@
|
|||
subroutine GGF(doG0F2,doevGF2,doqsGF2,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, &
|
||||
nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF)
|
||||
|
||||
! Green's function module
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
logical :: doG0F2
|
||||
logical :: doevGF2
|
||||
logical :: doqsGF2
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
logical,intent(in) :: dophBSE
|
||||
logical,intent(in) :: doppBSE
|
||||
logical,intent(in) :: TDA
|
||||
logical,intent(in) :: dBSE
|
||||
logical,intent(in) :: dTDA
|
||||
logical,intent(in) :: linearize
|
||||
double precision,intent(in) :: eta
|
||||
logical,intent(in) :: regularize
|
||||
|
||||
integer,intent(in) :: nNuc
|
||||
double precision,intent(in) :: ZNuc(nNuc)
|
||||
double precision,intent(in) :: rNuc(nNuc,ncart)
|
||||
double precision,intent(in) :: ENuc
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC(nspin)
|
||||
integer,intent(in) :: nO(nspin)
|
||||
integer,intent(in) :: nV(nspin)
|
||||
integer,intent(in) :: nR(nspin)
|
||||
integer,intent(in) :: nS(nspin)
|
||||
|
||||
double precision,intent(in) :: EHF
|
||||
double precision,intent(in) :: epsHF(nBas)
|
||||
double precision,intent(in) :: cHF(nBas,nBas)
|
||||
double precision,intent(in) :: PHF(nBas,nBas)
|
||||
double precision,intent(in) :: S(nBas,nBas)
|
||||
double precision,intent(in) :: T(nBas,nBas)
|
||||
double precision,intent(in) :: V(nBas,nBas)
|
||||
double precision,intent(in) :: Hc(nBas,nBas)
|
||||
double precision,intent(in) :: X(nBas,nBas)
|
||||
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
|
||||
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: start_GF ,end_GF ,t_GF
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute G0F2 electronic binding energies
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(doG0F2) then
|
||||
|
||||
call wall_time(start_GF)
|
||||
call GG0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, &
|
||||
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
|
||||
call wall_time(end_GF)
|
||||
|
||||
t_GF = end_GF - start_GF
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute evGF2 electronic binding energies
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(doevGF2) then
|
||||
|
||||
call wall_time(start_GF)
|
||||
! call evGGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
|
||||
! linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
|
||||
! ERI,dipole_int,epsHF)
|
||||
call wall_time(end_GF)
|
||||
|
||||
t_GF = end_GF - start_GF
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform qsGF2 calculation
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(doqsGF2) then
|
||||
|
||||
call wall_time(start_GF)
|
||||
! call qsGGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,eta,regularize,nNuc,ZNuc,rNuc,ENuc, &
|
||||
! nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF)
|
||||
call wall_time(end_GF)
|
||||
|
||||
t_GF = end_GF - start_GF
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGF2 = ',t_GF,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
|
|
@ -0,0 +1,78 @@
|
|||
subroutine GGF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,ERI,eGFlin,eOld,eGF,Z)
|
||||
|
||||
! Compute the graphical solution of the GF2 QP equation
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: eGFlin(nBas)
|
||||
double precision,intent(in) :: eOld(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: p
|
||||
integer :: nIt
|
||||
integer,parameter :: maxIt = 64
|
||||
double precision,parameter :: thresh = 1d-6
|
||||
double precision,external :: GGF2_SigC,GGF2_dSigC
|
||||
double precision :: SigC,dSigC
|
||||
double precision :: f,df
|
||||
double precision :: w
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: eGF(nBas)
|
||||
double precision,intent(out) :: Z(nBas)
|
||||
|
||||
! Run Newton's algorithm to find the root
|
||||
|
||||
write(*,*)'-----------------------------------------------------'
|
||||
write(*,'(A5,1X,A3,1X,A15,1X,A15,1X,A10)') 'Orb.','It.','e_GFlin (eV)','e_GF (eV)','Z'
|
||||
write(*,*)'-----------------------------------------------------'
|
||||
|
||||
do p=nC+1,nBas-nR
|
||||
|
||||
w = eGFlin(p)
|
||||
nIt = 0
|
||||
f = 1d0
|
||||
|
||||
do while (abs(f) > thresh .and. nIt < maxIt)
|
||||
|
||||
nIt = nIt + 1
|
||||
|
||||
SigC = GGF2_SigC(p,w,eta,nBas,nC,nO,nV,nR,eOld,ERI)
|
||||
dSigC = GGF2_dSigC(p,w,eta,nBas,nC,nO,nV,nR,eOld,ERI)
|
||||
f = w - eHF(p) - SigC
|
||||
df = 1d0/(1d0 - dSigC)
|
||||
|
||||
w = w - df*f
|
||||
|
||||
end do
|
||||
|
||||
if(nIt == maxIt) then
|
||||
|
||||
eGF(p) = eGFlin(p)
|
||||
write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,eGFlin(p)*HaToeV,eGF(p)*HaToeV,Z(p),'Cvg Failed!'
|
||||
|
||||
else
|
||||
|
||||
eGF(p) = w
|
||||
Z(p) = df
|
||||
|
||||
write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,eGFlin(p)*HaToeV,eGF(p)*HaToeV,Z(p)
|
||||
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
end subroutine
|
||||
|
|
@ -0,0 +1,50 @@
|
|||
double precision function GGF2_SigC(p,w,eta,nBas,nC,nO,nV,nR,eHF,ERI)
|
||||
|
||||
! Compute diagonal of the correlation part of the self-energy
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: p
|
||||
double precision,intent(in) :: w
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas,nC,nO,nV,nR
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,j,a,b
|
||||
double precision :: eps
|
||||
|
||||
GGF2_SigC = 0d0
|
||||
|
||||
! Occupied part of the correlation self-energy
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=nO+1,nBas-nR
|
||||
|
||||
eps = w + eHF(a) - eHF(i) - eHF(j)
|
||||
GGF2_SigC = GGF2_SigC + 0.5d0*(ERI(p,a,i,j) - ERI(p,a,j,i))**2*eps/(eps**2 + eta**2)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do i=nC+1,nO
|
||||
do a=nO+1,nBas-nR
|
||||
do b=nO+1,nBas-nR
|
||||
|
||||
eps = w + eHF(i) - eHF(a) - eHF(b)
|
||||
GGF2_SigC = GGF2_SigC + 0.5d0*(ERI(p,i,a,b) - ERI(p,i,b,a))**2*eps/(eps**2 + eta**2)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
end function
|
||||
|
|
@ -0,0 +1,52 @@
|
|||
double precision function GGF2_dSigC(p,w,eta,nBas,nC,nO,nV,nR,eHF,ERI)
|
||||
|
||||
! Compute diagonal of the correlation part of the self-energy
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: p
|
||||
double precision,intent(in) :: w
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas,nC,nO,nV,nR
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,j,a,b
|
||||
double precision :: eps
|
||||
|
||||
! Initialize
|
||||
|
||||
GGF2_dSigC = 0d0
|
||||
|
||||
! Occupied part of the correlation self-energy
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=nO+1,nBas-nR
|
||||
|
||||
eps = w + eHF(a) - eHF(i) - eHF(j)
|
||||
GGF2_dSigC = GGF2_dSigC - 0.5d0*(ERI(p,a,i,j) - ERI(p,a,j,i))**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do i=nC+1,nO
|
||||
do a=nO+1,nBas-nR
|
||||
do b=nO+1,nBas-nR
|
||||
|
||||
eps = w + eHF(i) - eHF(a) - eHF(b)
|
||||
GGF2_dSigC = GGF2_dSigC - 0.5d0*(ERI(p,i,a,b) - ERI(p,i,b,a))**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
end function
|
||||
|
|
@ -0,0 +1,72 @@
|
|||
subroutine GGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
|
||||
|
||||
! Compute diagonal part of the GF2 self-energy and its renormalization factor
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
double precision,intent(in) :: e(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,j,a,b
|
||||
integer :: p
|
||||
double precision :: eps
|
||||
double precision :: num
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: SigC(nBas)
|
||||
double precision,intent(out) :: Z(nBas)
|
||||
|
||||
! Initialize
|
||||
|
||||
SigC(:) = 0d0
|
||||
Z(:) = 0d0
|
||||
|
||||
! Compute GF2 self-energy
|
||||
|
||||
do p=nC+1,nBas-nR
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=nO+1,nBas-nR
|
||||
|
||||
eps = e(p) + e(a) - e(i) - e(j)
|
||||
num = 0.5d0*(ERI(p,a,i,j) - ERI(p,a,j,i))**2
|
||||
|
||||
SigC(p) = SigC(p) + num*eps/(eps**2 + eta**2)
|
||||
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
do p=nC+1,nBas-nR
|
||||
do i=nC+1,nO
|
||||
do a=nO+1,nBas-nR
|
||||
do b=nO+1,nBas-nR
|
||||
|
||||
eps = e(p) + e(i) - e(a) - e(b)
|
||||
num = 0.5d0*(ERI(p,i,a,b) - ERI(p,i,b,a))**2
|
||||
|
||||
SigC(p) = SigC(p) + num*eps/(eps**2 + eta**2)
|
||||
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
Z(:) = 1d0/(1d0 - Z(:))
|
||||
|
||||
end subroutine
|
||||
|
|
@ -24,7 +24,7 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole)
|
|||
! Local variables
|
||||
|
||||
integer :: ixyz
|
||||
integer :: i,j
|
||||
integer :: mu,nu
|
||||
integer :: HOMO
|
||||
integer :: LUMO
|
||||
double precision :: Gap
|
||||
|
|
@ -49,42 +49,37 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole)
|
|||
|
||||
allocate(Paa(nBas2,nBas2),Pab(nBas2,nBas2),Pba(nBas2,nBas2),Pbb(nBas2,nBas2))
|
||||
|
||||
! Paa(:,:) = P( 1:nBas , 1:nBas )
|
||||
! Pab(:,:) = P( 1:nBas ,nBas+1:nBas2)
|
||||
! Pba(:,:) = P(nBas+1:nBas2, 1:nBas )
|
||||
! Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2)
|
||||
Paa(:,:) = P( 1:nBas , 1:nBas )
|
||||
Pab(:,:) = P( 1:nBas ,nBas+1:nBas2)
|
||||
Pba(:,:) = P(nBas+1:nBas2, 1:nBas )
|
||||
Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2)
|
||||
|
||||
allocate(Ca(nBas,nBas2),Cb(nBas,nBas2))
|
||||
|
||||
Ca(:,:) = C( 1:nBas ,1:nBas2)
|
||||
Cb(:,:) = C(nBas+1:nBas2,1:nBas2)
|
||||
|
||||
Paa = matmul(transpose(Ca),matmul(P( 1:nBas , 1:nBas ),Ca))
|
||||
Pab = matmul(transpose(Ca),matmul(P( 1:nBas ,nBas+1:nBas2),Cb))
|
||||
Pba = matmul(transpose(Cb),matmul(P(nBas+1:nBas2, 1:nBas ),Ca))
|
||||
Pbb = matmul(transpose(Cb),matmul(P(nBas+1:nBas2,nBas+1:nBas2),Cb))
|
||||
|
||||
! Compute expectation values of S^2
|
||||
|
||||
Sx2 = 0.25d0*trace_matrix(nBas2,Paa+Pbb) + 0.25d0*trace_matrix(nBas2,Pab+Pba)**2
|
||||
do i=1,nBas2
|
||||
do j=1,nBas2
|
||||
Sx2 = Sx2 - 0.5d0*(Paa(i,j)*Pbb(j,i) + Pab(i,j)*Pab(j,i))
|
||||
Sx2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab+Pba)**2
|
||||
do mu=1,nBas
|
||||
do nu=1,nBas
|
||||
Sx2 = Sx2 - 0.5d0*(Paa(mu,nu)*Pbb(nu,mu) + Pab(mu,nu)*Pab(nu,mu))
|
||||
end do
|
||||
end do
|
||||
|
||||
Sy2 = 0.25d0*trace_matrix(nBas2,Paa+Pbb) - 0.25d0*trace_matrix(nBas2,Pab+Pba)**2
|
||||
do i=1,nBas2
|
||||
do j=1,nBas2
|
||||
Sy2 = Sy2 - 0.5d0*(Paa(i,j)*Pbb(j,i) - Pab(i,j)*Pab(j,i))
|
||||
Sy2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) - 0.25d0*trace_matrix(nBas,Pab+Pba)**2
|
||||
do mu=1,nBas
|
||||
do nu=1,nBas
|
||||
Sy2 = Sy2 - 0.5d0*(Paa(mu,nu)*Pbb(nu,mu) - Pab(mu,nu)*Pab(nu,mu))
|
||||
end do
|
||||
end do
|
||||
|
||||
Sz2 = 0.25d0*trace_matrix(nBas2,Paa+Pbb) + 0.25d0*trace_matrix(nBas2,Pab-Pba)**2
|
||||
do i=1,nBas2
|
||||
do j=1,nBas2
|
||||
Sz2 = Sz2 - 0.25d0*(Paa(i,j)*Pbb(j,i) - Pab(i,j)*Pab(j,i))
|
||||
Sz2 = Sz2 + 0.25d0*(Pab(i,j)*Pba(j,i) - Pba(i,j)*Pab(j,i))
|
||||
Sz2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab-Pba)**2
|
||||
do mu=1,nBas
|
||||
do nu=1,nBas
|
||||
Sz2 = Sz2 - 0.25d0*(Paa(mu,nu)*Pbb(nu,mu) - Pab(mu,nu)*Pab(nu,mu))
|
||||
Sz2 = Sz2 + 0.25d0*(Pab(mu,nu)*Pba(nu,mu) - Pba(mu,nu)*Pab(nu,mu))
|
||||
end do
|
||||
end do
|
||||
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
subroutine GMP(doMP2,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
|
||||
subroutine GMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
|
||||
|
||||
! Moller-Plesset module
|
||||
|
||||
|
|
@ -8,6 +8,7 @@ subroutine GMP(doMP2,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
|
|||
! Input variables
|
||||
|
||||
logical,intent(in) :: doMP2
|
||||
logical,intent(in) :: doMP3
|
||||
|
||||
logical,intent(in) :: regularize
|
||||
integer,intent(in) :: nBas
|
||||
|
|
@ -23,6 +24,7 @@ subroutine GMP(doMP2,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
|
|||
! Local variables
|
||||
|
||||
double precision :: start_MP ,end_MP ,t_MP
|
||||
double precision :: Ec
|
||||
|
||||
! Output variables
|
||||
|
||||
|
|
@ -33,7 +35,23 @@ subroutine GMP(doMP2,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
|
|||
if(doMP2) then
|
||||
|
||||
call wall_time(start_MP)
|
||||
call GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
|
||||
call GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF,Ec)
|
||||
call wall_time(end_MP)
|
||||
|
||||
t_MP = end_MP - start_MP
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute MP3 energy
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(doMP3) then
|
||||
|
||||
call wall_time(start_MP)
|
||||
call GMP3(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
|
||||
call wall_time(end_MP)
|
||||
|
||||
t_MP = end_MP - start_MP
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e)
|
||||
subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
|
||||
|
||||
! Perform second-order Moller-Plesset calculation with and without regularizers
|
||||
|
||||
|
|
@ -30,10 +30,11 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e)
|
|||
double precision :: E2x,E2xs,E2xs2,E2xk
|
||||
|
||||
double precision :: EcsMP2,Ecs2MP2,EckMP2
|
||||
double precision :: EcMP2
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: EcMP2
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
|
|
|
|||
|
|
@ -0,0 +1,170 @@
|
|||
subroutine GMP3(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e)
|
||||
|
||||
! Perform third-order Moller-Plesset calculation
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
double precision,intent(in) :: ENuc,EHF
|
||||
double precision,intent(in) :: e(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: eps,E2,EcMP2
|
||||
double precision :: eps1,eps2,E3a,E3b,E3c
|
||||
double precision :: EcMP3
|
||||
|
||||
integer :: i,j,k,l,a,b,c,d
|
||||
|
||||
double precision,allocatable :: dbERI(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: OOOO(:,:,:,:)
|
||||
double precision,allocatable :: OOVV(:,:,:,:)
|
||||
double precision,allocatable :: OVVO(:,:,:,:)
|
||||
double precision,allocatable :: VVOO(:,:,:,:)
|
||||
double precision,allocatable :: VVVV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: eO(:)
|
||||
double precision,allocatable :: eV(:)
|
||||
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'************************************************'
|
||||
write(*,*)'| Moller-Plesset third-order calculation |'
|
||||
write(*,*)'************************************************'
|
||||
write(*,*)
|
||||
|
||||
! Antysymmetrize ERIs
|
||||
|
||||
allocate(dbERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call antisymmetrize_ERI(2,nBas,ERI,dbERI)
|
||||
|
||||
! Form energy denominator
|
||||
|
||||
! Create integral batches
|
||||
|
||||
allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VVOO(nV,nV,nO,nO),VVVV(nV,nV,nV,nV))
|
||||
|
||||
OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO )
|
||||
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
|
||||
OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO )
|
||||
VVOO(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas, 1:nO , 1:nO )
|
||||
VVVV(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas,nO+1:nBas,nO+1:nBas)
|
||||
|
||||
deallocate(dbERI)
|
||||
|
||||
allocate(eO(nO),eV(nV))
|
||||
|
||||
eO(:) = e(1:nO)
|
||||
eV(:) = e(nO+1:nBas)
|
||||
|
||||
! Compute MP2 energy
|
||||
|
||||
E2 = 0d0
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
|
||||
eps = eO(i) + eO(j) - eV(a) - eV(b)
|
||||
|
||||
E2 = E2 + OOVV(i,j,a,b)*OOVV(i,j,a,b)/eps
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
EcMP2 = 0.25d0*E2
|
||||
|
||||
! Compute MP3 energy
|
||||
|
||||
E3a = 0d0
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do l=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
|
||||
eps1 = eO(i) + eO(j) - eV(a) - eV(b)
|
||||
eps2 = eO(k) + eO(l) - eV(a) - eV(b)
|
||||
|
||||
E3a = E3a + OOVV(i,j,a,b)*OOOO(k,l,i,j)*VVOO(a,b,k,l)/(eps1*eps2)
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
E3b = 0d0
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
do c=1,nV-nR
|
||||
do d=1,nV-nR
|
||||
|
||||
eps1 = eO(i) + eO(j) - eV(a) - eV(b)
|
||||
eps2 = eO(i) + eO(j) - eV(c) - eV(d)
|
||||
|
||||
E3b = E3b + OOVV(i,j,a,b)*VVVV(a,b,c,d)*VVOO(c,d,i,j)/(eps1*eps2)
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
E3c = 0d0
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
do c=1,nV-nR
|
||||
|
||||
eps1 = eO(i) + eO(j) - eV(a) - eV(b)
|
||||
eps2 = eO(i) + eO(k) - eV(a) - eV(c)
|
||||
|
||||
E3c = E3c + OOVV(i,j,a,b)*OVVO(k,b,c,j)*VVOO(a,c,i,k)/(eps1*eps2)
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
EcMP3 = 0.125d0*E3a + 0.125d0*E3b + E3c
|
||||
|
||||
write(*,*)
|
||||
write(*,'(A32)') '-----------------------'
|
||||
write(*,'(A32)') ' MP3 calculation '
|
||||
write(*,'(A32)') '-----------------------'
|
||||
write(*,'(A32,1X,F16.10)') ' MP2 contribution ',EcMP2
|
||||
write(*,'(A32,1X,F16.10)') ' MP3 contribution ',EcMP3
|
||||
write(*,'(A32)') '-----------------------'
|
||||
write(*,'(A32,1X,F16.10)') ' MP3 correlation energy', EcMP2 + EcMP3
|
||||
write(*,'(A32,1X,F16.10)') ' MP3 total energy',ENuc + EHF + EcMP2 + EcMP3
|
||||
write(*,'(A32)') '-----------------------'
|
||||
write(*,*)
|
||||
|
||||
end subroutine
|
||||
|
|
@ -1,4 +1,4 @@
|
|||
subroutine MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e)
|
||||
subroutine MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
|
||||
|
||||
! Perform second-order Moller-Plesset calculation with and without regularizers
|
||||
|
||||
|
|
@ -29,10 +29,11 @@ subroutine MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e)
|
|||
double precision :: E2x,E2xs,E2xs2,E2xk
|
||||
|
||||
double precision :: EcsMP2,Ecs2MP2,EckMP2
|
||||
double precision :: EcMP2
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: EcMP2
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
|
|
|
|||
|
|
@ -24,6 +24,7 @@ subroutine RMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
|
|||
! Local variables
|
||||
|
||||
double precision :: start_MP ,end_MP ,t_MP
|
||||
double precision :: Ec
|
||||
|
||||
! Output variables
|
||||
|
||||
|
|
@ -34,7 +35,7 @@ subroutine RMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
|
|||
if(doMP2) then
|
||||
|
||||
call wall_time(start_MP)
|
||||
call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
|
||||
call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF,Ec)
|
||||
call wall_time(end_MP)
|
||||
|
||||
t_MP = end_MP - start_MP
|
||||
|
|
|
|||
|
|
@ -26,6 +26,7 @@ subroutine UMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbb
|
|||
! Local variables
|
||||
|
||||
double precision :: start_MP ,end_MP ,t_MP
|
||||
double precision :: Ec(nsp)
|
||||
|
||||
! Output variables
|
||||
|
||||
|
|
@ -36,7 +37,7 @@ subroutine UMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbb
|
|||
if(doMP2) then
|
||||
|
||||
call wall_time(start_MP)
|
||||
call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
|
||||
call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF,Ec)
|
||||
call wall_time(end_MP)
|
||||
|
||||
t_MP = end_MP - start_MP
|
||||
|
|
|
|||
|
|
@ -29,10 +29,11 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec)
|
|||
double precision :: Edab,Exab,Ecab
|
||||
double precision :: Edbb,Exbb,Ecbb
|
||||
double precision :: Ed,Ex
|
||||
double precision :: Ec(nsp)
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: Ec(nsp)
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, &
|
||||
subroutine GQuAcK(doGHF,dostab,doMP2,doMP3,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, &
|
||||
nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, &
|
||||
maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, &
|
||||
spin_conserved,spin_flip,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, &
|
||||
|
|
@ -11,6 +11,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
|
|||
logical,intent(in) :: doGHF
|
||||
logical,intent(in) :: dostab
|
||||
logical,intent(in) :: doMP2
|
||||
logical,intent(in) :: doMP3
|
||||
logical,intent(in) :: dophRPA,dophRPAx,doppRPA
|
||||
logical,intent(in) :: doG0F2,doevGF2,doqsGF2
|
||||
logical,intent(in) :: doG0W0,doevGW,doqsGW
|
||||
|
|
@ -180,7 +181,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
|
|||
if(doMP) then
|
||||
|
||||
call wall_time(start_MP)
|
||||
call GMP(doMP2,reg_MP,nBas2,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
|
||||
call GMP(doMP2,doMP3,reg_MP,nBas2,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
|
||||
call wall_time(end_MP)
|
||||
|
||||
t_MP = end_MP - start_MP
|
||||
|
|
@ -198,8 +199,8 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
|
|||
if(doRPA) then
|
||||
|
||||
call wall_time(start_RPA)
|
||||
call GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel, &
|
||||
nBas2,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF,cHF,S)
|
||||
call GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas2,nC,nO,nV,nR,nS,ENuc,EHF, &
|
||||
ERI_MO,dipole_int_MO,epsHF,cHF,S)
|
||||
call wall_time(end_RPA)
|
||||
|
||||
t_RPA = end_RPA - start_RPA
|
||||
|
|
@ -217,10 +218,8 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
|
|||
if(doGF) then
|
||||
|
||||
call wall_time(start_GF)
|
||||
! call GGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, &
|
||||
! dophBSE,doppBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,lin_GF,eta_GF,reg_GF, &
|
||||
! nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO, &
|
||||
! dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
|
||||
call GGF(doG0F2,doevGF2,doqsGF2,maxSCF_GF,thresh_GF,max_diis_GF,dophBSE,doppBSE,TDA,dBSE,dTDA,lin_GF,eta_GF,reg_GF, &
|
||||
nNuc,ZNuc,rNuc,ENuc,nBas2,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
|
||||
call wall_time(end_GF)
|
||||
|
||||
t_GF = end_GF - start_GF
|
||||
|
|
@ -238,10 +237,9 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
|
|||
if(doGW) then
|
||||
|
||||
call wall_time(start_GW)
|
||||
call GGW(doG0W0,doevGW,doqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, &
|
||||
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA, &
|
||||
lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas2,nC,nO,nV,nR,nS, &
|
||||
EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
|
||||
call GGW(doG0W0,doevGW,doqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT,exchange_kernel,doXBS, &
|
||||
dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc, &
|
||||
nBas2,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
|
||||
call wall_time(end_GW)
|
||||
|
||||
t_GW = end_GW - start_GW
|
||||
|
|
|
|||
|
|
@ -215,7 +215,7 @@ program QuAcK
|
|||
!--------------------------!
|
||||
|
||||
if(doGQuAcK) &
|
||||
call GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, &
|
||||
call GQuAcK(doGHF,dostab,doMP2,doMP3,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, &
|
||||
nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, &
|
||||
maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, &
|
||||
spin_conserved,spin_flip,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, &
|
||||
|
|
|
|||
|
|
@ -70,7 +70,7 @@ subroutine GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO,
|
|||
if(doppRPA) then
|
||||
|
||||
call wall_time(start_RPA)
|
||||
! call ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,epsHF)
|
||||
call ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,epsHF)
|
||||
call wall_time(end_RPA)
|
||||
|
||||
t_RPA = end_RPA - start_RPA
|
||||
|
|
|
|||
|
|
@ -0,0 +1,100 @@
|
|||
subroutine ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
|
||||
|
||||
! Perform ppGRPA calculation
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: TDA
|
||||
logical,intent(in) :: doACFDT
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: EHF
|
||||
double precision,intent(in) :: e(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: ispin
|
||||
integer :: nOO
|
||||
integer :: nVV
|
||||
double precision,allocatable :: Bpp(:,:)
|
||||
double precision,allocatable :: Cpp(:,:)
|
||||
double precision,allocatable :: Dpp(:,:)
|
||||
double precision,allocatable :: Om1(:)
|
||||
double precision,allocatable :: X1(:,:)
|
||||
double precision,allocatable :: Y1(:,:)
|
||||
double precision,allocatable :: Om2(:)
|
||||
double precision,allocatable :: X2(:,:)
|
||||
double precision,allocatable :: Y2(:,:)
|
||||
|
||||
double precision :: EcRPA
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'****************************************'
|
||||
write(*,*)'| particle-particle GRPA calculation |'
|
||||
write(*,*)'****************************************'
|
||||
write(*,*)
|
||||
|
||||
! Initialization
|
||||
|
||||
EcRPA = 0d0
|
||||
|
||||
ispin = 4
|
||||
|
||||
nOO = nO*(nO-1)/2
|
||||
nVV = nV*(nV-1)/2
|
||||
|
||||
allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), &
|
||||
Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
|
||||
|
||||
if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
|
||||
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,e,ERI,Cpp)
|
||||
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,e,ERI,Dpp)
|
||||
|
||||
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
|
||||
|
||||
! call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
|
||||
|
||||
call print_excitation_energies('ppRPA@GHF (N+2)',ispin,nVV,Om1)
|
||||
call print_excitation_energies('ppRPA@GHF (N-2)',ispin,nOO,Om2)
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy =',EcRPA
|
||||
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy =',ENuc + EHF + EcRPA
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
! Compute the correlation energy via the adiabatic connection
|
||||
|
||||
! if(doACFDT) then
|
||||
|
||||
! write(*,*) '--------------------------------------------------------'
|
||||
! write(*,*) 'Adiabatic connection version of ppRPA correlation energy'
|
||||
! write(*,*) '--------------------------------------------------------'
|
||||
! write(*,*)
|
||||
|
||||
! call ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ERI,e,EcRPA)
|
||||
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (singlet) =',EcRPA(1),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (triplet) =',EcRPA(2),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy =',EcRPA(1) + EcRPA(2),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy =',ENuc + EHF + EcRPA(1) + EcRPA(2),' au'
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
|
||||
! end if
|
||||
|
||||
end subroutine
|
||||
Loading…
Reference in New Issue