LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-07-11 22:03:50 +02:00
Code Releases Activity

add HF module

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-07-17 13:35:24 +02:00
parent ae540679cd
commit 2c7653b7bc
17 changed files with 370 additions and 528 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
input/methods
View File

@ -1,5 +1,5 @@
# RHF UHF KS MOM
T F F F
# RHF UHF RMOM UMOM KS
T F F F F
# MP2* MP3
F F
# CCD pCCD DCD CCSD CCSD(T)
@ -11,9 +11,9 @@
# phRPA* phRPAx* crRPA ppRPA
F F F F
# G0F2* evGF2* qsGF2* G0F3 evGF3
T F F F F
F F F F F
# G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
F F F F F F
# G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
F F F T F F
F F F F F F
# * unrestricted version available

2
input/options
View File

@ -15,4 +15,4 @@
# ACFDT: AC Kx XBS
F T T
# BSE: phBSE phBSE2 ppBSE dBSE dTDA evDyn
T F F T F F
F F F T F F

11
src/GT/G0T0eh.f90
View File

@ -1,6 +1,6 @@
subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,ppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, &
ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eGT)
ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc)
! Perform ehG0T0 calculation
@ -58,15 +58,14 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rhoL_RPA(:,:,:)
double precision,allocatable :: rhoR_RPA(:,:,:)
double precision,allocatable :: rhoL_RPA(:,:,:,:)
double precision,allocatable :: rhoR_RPA(:,:,:,:)
double precision,allocatable :: eGT(:)
double precision,allocatable :: eGTlin(:)
! Output variables
double precision :: eGT(nBas)
! Hello world
write(*,*)
@ -100,7 +99,7 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evD
! Memory allocation
allocate(Sig(nBas),SigX(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), &
rhoL_RPA(nBas,nBas,nS),rhoR_RPA(nBas,nBas,nS),eGTlin(nBas))
rhoL_RPA(nBas,nBas,nS,2),rhoR_RPA(nBas,nBas,nS,2),eGT(nBas),eGTlin(nBas))
!-------------------!
! Compute screening !

11
src/GT/G0T0pp.f90
View File

@ -1,5 +1,5 @@
subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eGT)
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc)
! Perform one-shot calculation with a T-matrix self-energy (G0T0)
@ -62,12 +62,11 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp
double precision,allocatable :: SigX(:)
double precision,allocatable :: Sig(:)
double precision,allocatable :: Z(:)
double precision,allocatable :: eGT(:)
double precision,allocatable :: eGTlin(:)
! Output variables
double precision,intent(out) :: eGT(nBas)
! Hello world
write(*,*)
@ -88,11 +87,11 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp
allocate(Om1ab(nVVab),X1ab(nVVab,nVVab),Y1ab(nOOab,nVVab), &
Om2ab(nOOab),X2ab(nVVab,nOOab),Y2ab(nOOab,nOOab), &
rho1ab(nBas,nBas,nVVab),rho2ab(nBas,nBas,nOOab), &
rho1ab(nBas,nBas,nVVab),rho2ab(nBas,nBas,nOOab), &
Om1aa(nVVaa),X1aa(nVVaa,nVVaa),Y1aa(nOOaa,nVVaa), &
Om2aa(nOOaa),X2aa(nVVaa,nOOaa),Y2aa(nOOaa,nOOaa), &
rho1aa(nBas,nBas,nVVaa),rho2aa(nBas,nBas,nOOaa), &
SigX(nBas),Sig(nBas),Z(nBas),eGTlin(nBas))
rho1aa(nBas,nBas,nVVaa),rho2aa(nBas,nBas,nOOaa), &
SigX(nBas),Sig(nBas),Z(nBas),eGT(nBas),eGTlin(nBas))
!----------------------------------------------
! alpha-beta block

20
src/GT/GTeh_excitation_density.f90
View File

@ -14,17 +14,17 @@ subroutine GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR)
! Local variables
integer :: m,jb,p,q,j,b
double precision :: X,Y
double precision :: X,Y,Xt,Yt
! Output variables
double precision,intent(out) :: rhoL(nBas,nBas,nS)
double precision,intent(out) :: rhoR(nBas,nBas,nS)
double precision,intent(out) :: rhoL(nBas,nBas,nS,2)
double precision,intent(out) :: rhoR(nBas,nBas,nS,2)
! Initialization
rhoL(:,:,:) = 0d0
rhoR(:,:,:) = 0d0
rhoL(:,:,:,:) = 0d0
rhoR(:,:,:,:) = 0d0
!$OMP PARALLEL &
!$OMP SHARED(nC,nBas,nR,nO,nS,rhoL,rhoR,ERI,XpY,XmY) &
@ -42,8 +42,14 @@ subroutine GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR)
X = 0.5d0*(XpY(m,jb) + XmY(m,jb))
Y = 0.5d0*(XpY(m,jb) - XmY(m,jb))
rhoL(p,q,m) = rhoL(p,q,m) + ERI(p,j,b,q)*X + ERI(p,b,j,q)*Y
rhoR(p,q,m) = rhoR(p,q,m) + (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*X + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*Y
Xt = 0.5d0*(XpY(jb,m) + XmY(jb,m))
Yt = 0.5d0*(XpY(jb,m) - XmY(jb,m))
rhoL(p,q,m,1) = rhoL(p,q,m,1) + ERI(p,j,b,q)*X + ERI(p,b,j,q)*Y
rhoL(p,q,m,2) = rhoL(p,q,m,2) + ERI(p,b,j,q)*X + ERI(p,j,b,q)*Y
rhoR(p,q,m,1) = rhoR(p,q,m,1) + (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*X + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*Y
rhoR(p,q,m,2) = rhoR(p,q,m,2) + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*X + (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*Y
enddo
enddo

28
src/GT/GTeh_self_energy_diag.f90
View File

@ -16,12 +16,12 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig
integer,intent(in) :: nS
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rhoL(nBas,nBas,nS)
double precision,intent(in) :: rhoR(nBas,nBas,nS)
double precision,intent(in) :: rhoL(nBas,nBas,nS,2)
double precision,intent(in) :: rhoR(nBas,nBas,nS,2)
! Local variables
integer :: i,a,p,q,m
integer :: i,a,p,m
double precision :: num,eps
! Output variables
@ -46,7 +46,7 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig
do m=1,nS
eps = e(p) - e(i) + Om(m)
num = rhoL(i,p,m)*rhoR(i,p,m)
num = rhoL(p,i,m,2)*rhoR(i,p,m,1)
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
@ -61,7 +61,7 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig
do m=1,nS
eps = e(p) - e(a) - Om(m)
num = rhoL(p,a,m)*rhoR(p,a,m)
num = rhoL(p,a,m,1)*rhoR(a,p,m,2)
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
@ -72,17 +72,17 @@ subroutine GTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,Sig
! GM correlation energy
EcGM = 0d0
do i=nC+1,nO
do a=nO+1,nBas-nR
do m=1,nS
! do i=nC+1,nO
! do a=nO+1,nBas-nR
! do m=1,nS
eps = e(a) - e(i) + Om(m)
num = rhoL(i,a,m)*rhoR(i,a,m)
EcGM = EcGM - num*eps/(eps**2 + eta**2)
! eps = e(a) - e(i) + Om(m)
! num = rhoL(i,a,m)*rhoR(i,a,m)
! EcGM = EcGM - num*eps/(eps**2 + eta**2)
end do
end do
end do
! end do
! end do
! end do
! Compute renormalization factor from derivative

16
src/GT/dynamic_Tmatrix_A.f90
View File

@ -1,4 +1,4 @@
subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,Omega2,rho1,rho2,OmBSE,TA,ZA)
subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Om1,Om2,rho1,rho2,OmBSE,TA,ZA)
! Compute the dynamic part of the Bethe-Salpeter equation matrices for GT
@ -23,8 +23,8 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
double precision,intent(in) :: OmBSE
double precision,intent(in) :: Omega1(nVV)
double precision,intent(in) :: Omega2(nOO)
double precision,intent(in) :: Om1(nVV)
double precision,intent(in) :: Om2(nOO)
double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: rho2(nBas,nBas,nOO)
@ -47,7 +47,7 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
! Build dynamic A matrix
jb = 0
!$omp parallel do default(private) shared(TA,ZA,Omega1,Omega2,OmBSE,eGT,rho1,rho2,nO,nBas,nVV,nOO,chi,eps,eta,nC,nR,lambda)
!$omp parallel do default(private) shared(TA,ZA,Om1,Om2,OmBSE,eGT,rho1,rho2,nO,nBas,nVV,nOO,chi,eps,eta,nC,nR,lambda)
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = (b-nO) + (j-1)*(nBas-nO)
@ -60,12 +60,12 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
chi = 0d0
do cd=1,nVV
eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j))
eps = + OmBSE - Om1(cd) + (eGT(i) + eGT(j))
chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*eps/(eps**2 + eta**2)
end do
do kl=1,nOO
eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b))
eps = + OmBSE + Om2(kl) - (eGT(a) + eGT(b))
chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*eps/(eps**2 + eta**2)
end do
@ -74,12 +74,12 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
chi = 0d0
do cd=1,nVV
eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j))
eps = + OmBSE - Om1(cd) + (eGT(i) + eGT(j))
chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
do kl=1,nOO
eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b))
eps = + OmBSE + Om2(kl) - (eGT(a) + eGT(b))
chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do

5
src/GT/evGTeh.f90
View File

@ -1,6 +1,6 @@
subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,ppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF, &
cHF,eHF,Vxc,eG0T0)
cHF,eHF,Vxc)
! Perform self-consistent eigenvalue-only ehGT calculation
@ -41,7 +41,6 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: Vxc(nBas)
double precision,intent(in) :: eG0T0(nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
@ -121,7 +120,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,
Conv = 1d0
e_diis(:,:) = 0d0
error_diis(:,:) = 0d0
eGT(:) = eG0T0(:)
eGT(:) = eHF(:)
eOld(:) = eGT(:)
Z(:) = 1d0
rcond = 0d0

8
src/GT/evGTpp.f90
View File

@ -1,6 +1,5 @@
subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nBas, &
nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0T0)
subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, &
singlet,triplet,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc)
! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT)
@ -41,7 +40,6 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: eG0T0(nBas)
! Local variables
@ -114,7 +112,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
Conv = 1d0
e_diis(:,:) = 0d0
error_diis(:,:) = 0d0
eGT(:) = eG0T0(:)
eGT(:) = eHF(:)
eOld(:) = eGT(:)
Z(:) = 1d0
rcond = 0d0

9
src/GW/G0W0.f90
View File

@ -1,6 +1,6 @@
subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, &
ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eGW)
ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc)
! Perform G0W0 calculation
@ -61,6 +61,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: eGW(:)
double precision,allocatable :: eGWlin(:)
integer :: nBas2
@ -70,12 +71,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD
integer :: nR2
integer :: nS2
double precision,allocatable :: seHF(:),seGW(:),sERI(:,:,:,:)
! Output variables
double precision :: eGW(nBas)
! Hello world
write(*,*)
@ -115,7 +112,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD
! Memory allocation
allocate(SigC(nBas),SigX(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS),eGWlin(nBas))
allocate(SigC(nBas),SigX(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS),eGW(nBas),eGWlin(nBas))
!-------------------!
! Compute screening !

13
src/GW/evGW.f90
View File

@ -1,6 +1,6 @@
subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF, &
cHF,eHF,Vxc,eG0W0)
cHF,eHF,Vxc)
! Perform self-consistent eigenvalue-only GW calculation
@ -42,7 +42,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: Vxc(nBas)
double precision,intent(in) :: eG0W0(nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
@ -137,7 +136,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
Conv = 1d0
e_diis(:,:) = 0d0
error_diis(:,:) = 0d0
eGW(:) = eG0W0(:)
eGW(:) = eHF(:)
eOld(:) = eGW(:)
Z(:) = 1d0
rcond = 0d0
@ -312,10 +311,10 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcppBSE(1)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcppBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy =',EcppBSE(1) + 3d0*EcppBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 total energy =',ENuc + ERHF + EcppBSE(1) + 3d0*EcppBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (singlet) =',EcppBSE(1)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (triplet) =',3d0*EcppBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy =',EcppBSE(1) + 3d0*EcppBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW total energy =',ENuc + ERHF + EcppBSE(1) + 3d0*EcppBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)

6
src/GW/qsGW.f90
View File

@ -1,6 +1,6 @@
subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,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)
subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn, &
singlet,triplet,eta,regularize,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)
! Perform a quasiparticle self-consistent GW calculation

136
src/HF/HF.f90 Normal file
View File

@ -0,0 +1,136 @@
subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,guess_type,mix,level_shift, &
nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,EHF,epsHF,cHF,PHF,vHF)
! Hartree-Fock module
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: doRHF
logical,intent(in) :: doUHF
logical,intent(in) :: doRMOM
logical,intent(in) :: doUMOM
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
double precision,intent(in) :: thresh
double precision,intent(in) :: level_shift
logical,intent(in) :: mix
integer,intent(in) :: nBas
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
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(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
integer :: nSCF
integer :: n_diis
double precision :: ET
double precision :: EV
double precision :: EJ
double precision :: EK
double precision :: dipole(ncart)
double precision :: Conv
double precision :: Gap
double precision :: rcond
double precision,external :: trace_matrix
double precision,allocatable :: error(:,:)
double precision,allocatable :: error_diis(:,:)
double precision,allocatable :: F_diis(:,:)
double precision,allocatable :: J(:,:)
double precision,allocatable :: K(:,:)
double precision,allocatable :: cp(:,:)
double precision,allocatable :: Fp(:,:)
! Output variables
logical,intent(out) :: unrestricted
double precision,intent(out) :: EHF
double precision,intent(out) :: epsHF(nBas)
double precision,intent(out) :: cHF(nBas,nBas)
double precision,intent(out) :: PHF(nBas,nBas)
double precision,intent(out) :: vHF(nBas)
double precision,intent(out) :: F(nBas,nBas)
!------------------------------------------------------------------------
! Compute RHF energy
!------------------------------------------------------------------------
if(doRHF) then
! Check that RHF calculation is worth doing...
if(nO(1) /= nO(2)) then
write(*,*) ' !!! The system does not appear to be closed shell !!!'
write(*,*)
stop
end if
call RHF(maxSCF,thresh,n_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,EHF,epsHF,cHF,PHF,vHF)
end if
!------------------------------------------------------------------------
! Compute UHF energy
!------------------------------------------------------------------------
if(doUHF) then
! Switch on the unrestricted flag
unrestricted = .true.
call UHF(maxSCF,thresh,n_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF,vHF)
end if
!------------------------------------------------------------------------
! Restricted maximum overlap method
!------------------------------------------------------------------------
if(doRMOM) then
! Check that RMOM calculation is worth doing...
if(nO(1) /= nO(2)) then
write(*,*) ' !!! The system does not appear to be closed shell !!!'
write(*,*)
stop
end if
! call RMOM(maxSCF,thresh,n_diis,guess_type,nNuc,ZNuc,rNuc,ENuc, &
! nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF,vHF)
end if
!------------------------------------------------------------------------
! Unrestricted maximum overlap method
!------------------------------------------------------------------------
if(doUMOM) then
! Switch on the unrestricted flag
unrestricted = .true.
! call UMOM(maxSCF,thresh,n_diis,guess_type,nNuc,ZNuc,rNuc,ENuc, &
! nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF,vHF)
end if
end subroutine

195
src/HF/MOM.f90
View File

@ -1,195 +0,0 @@
subroutine MOM(maxSCF,thresh,max_diis,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,c,e,P)
! Maximum overlap method
implicit none
! Input variables
integer,intent(in) :: maxSCF,max_diis
double precision,intent(in) :: thresh
integer,intent(in) :: nBas,nO
double precision,intent(in) :: ENuc
double precision,intent(in) :: S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),X(nBas,nBas)
! Local variables
integer :: iBas,jBas
integer :: nSCF,nBasSq,n_diis
double precision :: ET,EV,EJ,EK,Conv,Gap
double precision :: rcond
double precision,external :: trace_matrix
double precision,allocatable :: error(:,:),error_diis(:,:),F_diis(:,:)
double precision,allocatable :: J(:,:),K(:,:),cp(:,:),F(:,:),Fp(:,:)
double precision,allocatable :: cG(:,:),ON(:)
! Output variables
double precision,intent(inout):: ERHF,c(nBas,nBas),e(nBas),P(nBas,nBas)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Maximum overlap method |'
write(*,*)'************************************************'
write(*,*)
! Useful quantities
nBasSq = nBas*nBas
! Memory allocation
allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas), &
cp(nBas,nBas),Fp(nBas,nBas),F(nBas,nBas), &
cG(nBas,nBas),ON(nBas), &
error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Set up guess orbitals
cG(:,:) = c(:,:)
! Set up occupation numbers
ON(1:nO) = 1d0
ON(nO+1:nBas) = 0d0
! HOMO-LUMO transition
ON(nO) = 0d0
ON(nO+1) = 1d0
write(*,*)
write(*,*) ' --- Initial MO occupations --- '
write(*,*)
call matout(nBas,1,ON)
write(*,*)
! Compute density matrix
call density_matrix(nBas,ON,c,P)
! Initialization
n_diis = 0
F_diis(:,:) = 0d0
error_diis(:,:) = 0d0
Conv = 1d0
nSCF = 0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
write(*,*)
write(*,*)'----------------------------------------------------'
write(*,*)'| MOM calculation |'
write(*,*)'----------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
'|','#','|','HF energy','|','Conv','|','HL Gap','|'
write(*,*)'----------------------------------------------------'
do while(Conv > thresh .and. nSCF < maxSCF)
! Increment
nSCF = nSCF + 1
! Build Fock matrix
call Coulomb_matrix_AO_basis(nBas,P,ERI,J)
call exchange_matrix_AO_basis(nBas,P,ERI,K)
F(:,:) = Hc(:,:) + J(:,:) + 0.5*K(:,:)
! Check convergence
error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
Conv = maxval(abs(error))
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
! Reset DIIS if required
if(abs(rcond) < 1d-15) n_diis = 0
! Diagonalize Fock matrix
Fp = matmul(transpose(X),matmul(F,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c = matmul(X,cp)
! MOM overlap
call MOM_overlap(nBas,nO,S,cG,c,ON)
! Density matrix
call density_matrix(nBas,ON,c,P)
! Compute HF energy
ERHF = trace_matrix(nBas,matmul(P,Hc)) &
+ 0.5d0*trace_matrix(nBas,matmul(P,J)) &
+ 0.5d0*trace_matrix(nBas,matmul(P,K))
! Compute HOMO-LUMO gap
if(nBas > nO) then
Gap = e(nO+1) - e(nO)
else
Gap = 0d0
endif
! Dump results
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
'|',nSCF,'|',ERHF+ENuc,'|',Conv,'|',Gap,'|'
enddo
write(*,*)'----------------------------------------------------'
!------------------------------------------------------------------------
! End of SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
stop
endif
write(*,*)
write(*,*) ' --- Final MO occupations --- '
write(*,*)
call matout(nBas,1,ON)
write(*,*)
! Compute HF energy
ET = trace_matrix(nBas,matmul(P,T))
EV = trace_matrix(nBas,matmul(P,V))
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
EK = 0.5d0*trace_matrix(nBas,matmul(P,K))
ERHF = ET + EV + EJ + EK
call print_RHF(nBas,nO,e,c,ENuc,ET,EV,EJ,EK,ERHF)
end subroutine MOM

152
src/HF/RMOM.f90
View File

@ -1,68 +1,40 @@
subroutine RMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,e,c,P,Vx)
subroutine RMOM(maxSCF,thresh,max_diis,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,c,e,P)
! Perform restricted Hartree-Fock calculation with MOM algorithm
! Maximum overlap method
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: maxSCF,max_diis,guess_type
integer,intent(in) :: maxSCF,max_diis
double precision,intent(in) :: thresh
integer,intent(in) :: nBas
integer,intent(in) :: nO
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
integer,intent(in) :: nBas,nO
double precision,intent(in) :: ENuc
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(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),X(nBas,nBas)
! Local variables
integer :: nSCF
integer :: nBasSq
integer :: n_diis
double precision :: ET
double precision :: EV
double precision :: EJ
double precision :: EK
double precision :: dipole(ncart)
double precision :: Conv
double precision :: Gap
integer :: iBas,jBas
integer :: nSCF,nBasSq,n_diis
double precision :: ET,EV,EJ,EK,Conv,Gap
double precision :: rcond
double precision,external :: trace_matrix
double precision,allocatable :: error(:,:)
double precision,allocatable :: error_diis(:,:)
double precision,allocatable :: F_diis(:,:)
double precision,allocatable :: J(:,:)
double precision,allocatable :: K(:,:)
double precision,allocatable :: cp(:,:)
double precision,allocatable :: F(:,:)
double precision,allocatable :: Fp(:,:)
double precision,allocatable :: ON(:)
double precision,allocatable :: error(:,:),error_diis(:,:),F_diis(:,:)
double precision,allocatable :: J(:,:),K(:,:),cp(:,:),F(:,:),Fp(:,:)
double precision,allocatable :: cG(:,:),ON(:)
! Output variables
double precision,intent(out) :: ERHF
double precision,intent(out) :: e(nBas)
double precision,intent(out) :: c(nBas,nBas)
double precision,intent(out) :: P(nBas,nBas)
double precision,intent(out) :: Vx(nBas)
double precision,intent(inout):: ERHF,c(nBas,nBas),e(nBas),P(nBas,nBas)
! Hello world
write(*,*)
write(*,*)'*********************************************'
write(*,*)'| Restricted Maximum Overlap Method |'
write(*,*)'*********************************************'
write(*,*)'************************************************'
write(*,*)'| Maximum overlap method |'
write(*,*)'************************************************'
write(*,*)
! Useful quantities
@ -71,22 +43,35 @@ subroutine RMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,
! Memory allocation
allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas), &
cp(nBas,nBas),Fp(nBas,nBas),F(nBas,nBas),ON(nBas), &
allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas), &
cp(nBas,nBas),Fp(nBas,nBas),F(nBas,nBas), &
cG(nBas,nBas),ON(nBas), &
error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Guess coefficients and eigenvalues
call mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
! ON(:) = 0d0
! do i=1,nO
! ON(i) = 1d0
! ON(i) = dble(2*i-1)
! end do
! call density_matrix(nBas,ON,c,P)
! Set up guess orbitals
cG(:,:) = c(:,:)
! Set up occupation numbers
ON(1:nO) = 1d0
ON(nO+1:nBas) = 0d0
! HOMO-LUMO transition
ON(nO) = 0d0
ON(nO+1) = 1d0
write(*,*)
write(*,*) ' --- Initial MO occupations --- '
write(*,*)
call matout(nBas,1,ON)
write(*,*)
! Compute density matrix
call density_matrix(nBas,ON,c,P)
! Initialization
n_diis = 0
@ -100,7 +85,7 @@ subroutine RMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,
!------------------------------------------------------------------------
write(*,*)
write(*,*)'----------------------------------------------------'
write(*,*)'| RHF calculation |'
write(*,*)'| MOM calculation |'
write(*,*)'----------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
'|','#','|','HF energy','|','Conv','|','HL Gap','|'
@ -117,7 +102,7 @@ subroutine RMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,
call Coulomb_matrix_AO_basis(nBas,P,ERI,J)
call exchange_matrix_AO_basis(nBas,P,ERI,K)
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
F(:,:) = Hc(:,:) + J(:,:) + 0.5*K(:,:)
! Check convergence
@ -127,11 +112,11 @@ subroutine RMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
if(abs(rcond) > 1d-7) then
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
else
n_diis = 0
end if
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
! Reset DIIS if required
if(abs(rcond) < 1d-15) n_diis = 0
! Diagonalize Fock matrix
@ -140,17 +125,19 @@ subroutine RMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,
call diagonalize_matrix(nBas,cp,e)
c = matmul(X,cp)
! MOM overlap
call MOM_overlap(nBas,nO,S,cG,c,ON)
! Density matrix
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
call density_matrix(nBas,ON,c,P)
! call density_matrix(nBas,ON,c,P)
! Compute HF energy
ERHF = trace_matrix(nBas,matmul(P,Hc)) &
+ 0.5d0*trace_matrix(nBas,matmul(P,J)) &
+ 0.25d0*trace_matrix(nBas,matmul(P,K))
ERHF = trace_matrix(nBas,matmul(P,Hc)) &
+ 0.5d0*trace_matrix(nBas,matmul(P,J)) &
+ 0.5d0*trace_matrix(nBas,matmul(P,K))
! Compute HOMO-LUMO gap
@ -189,21 +176,20 @@ subroutine RMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,
endif
write(*,*)
write(*,*) ' --- Final MO occupations --- '
write(*,*)
call matout(nBas,1,ON)
write(*,*)
! Compute HF energy
ET = trace_matrix(nBas,matmul(P,T))
EV = trace_matrix(nBas,matmul(P,V))
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
EK = 0.25d0*trace_matrix(nBas,matmul(P,K))
ET = trace_matrix(nBas,matmul(P,T))
EV = trace_matrix(nBas,matmul(P,V))
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
EK = 0.5d0*trace_matrix(nBas,matmul(P,K))
ERHF = ET + EV + EJ + EK
! Compute dipole moments
call print_RHF(nBas,nO,e,c,ENuc,ET,EV,EJ,EK,ERHF)
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_RHF(nBas,nO,e,C,ENuc,ET,EV,EJ,EK,ERHF,dipole)
! Compute Vx for post-HF calculations
call exchange_potential(nBas,c,K,Vx)
end subroutine RMOM
end subroutine

254
src/QuAcK/QuAcK.f90
View File

@ -3,9 +3,8 @@ program QuAcK
implicit none
include 'parameters.h'
logical :: doSph
logical :: unrestricted = .false.
logical :: doRHF,doUHF,doMOM
logical :: doHF,doRHF,doUHF,doRMOM,doUMOM
logical :: dostab
logical :: doKS
logical :: doMP2,doMP3
@ -32,9 +31,6 @@ program QuAcK
double precision,allocatable :: cHF(:,:,:),eHF(:,:),PHF(:,:,:)
double precision,allocatable :: Vxc(:,:)
double precision,allocatable :: eG0W0(:,:)
double precision,allocatable :: eG0T0(:,:)
logical :: doACFDT
logical :: exchange_kernel
logical :: doXBS
@ -122,14 +118,12 @@ program QuAcK
write(*,*)
! Spherium calculation?
doSph = .false.
call wall_time(start_QuAcK)
! Which calculations do you want to do?
call read_methods(doRHF,doUHF,doKS,doMOM, &
call read_methods(doRHF,doUHF,doRMOM,doUMOM,doKS, &
doMP2,doMP3, &
doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_crCCD,do_lCCD, &
@ -168,7 +162,7 @@ program QuAcK
! nS = number of single excitation
! = nO*nV
call read_molecule(nNuc,nEl(:),nO(:),nC(:),nR(:))
call read_molecule(nNuc,nEl,nO,nC,nR)
allocate(ZNuc(nNuc),rNuc(nNuc,ncart))
! Read geometry
@ -193,24 +187,16 @@ program QuAcK
! Memory allocation for one- and two-electron integrals
allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas,nspin),eG0T0(nBas,nspin),PHF(nBas,nBas,nspin), &
S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas), &
dipole_int_AO(nBas,nBas,ncart),dipole_int_MO(nBas,nBas,ncart),Vxc(nBas,nspin),F_AO(nBas,nBas))
allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),PHF(nBas,nBas,nspin),S(nBas,nBas),T(nBas,nBas), &
V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas),dipole_int_AO(nBas,nBas,ncart), &
dipole_int_MO(nBas,nBas,ncart),Vxc(nBas,nspin),F_AO(nBas,nBas))
! Read integrals
call wall_time(start_int)
if(doSph) then
call read_integrals_sph(nBas,S,T,V,Hc,ERI_AO)
else
call read_integrals(nBas,S,T,V,Hc,ERI_AO)
call read_dipole_integrals(nBas,dipole_int_AO)
end if
call read_integrals(nBas,S,T,V,Hc,ERI_AO)
call read_dipole_integrals(nBas,dipole_int_AO)
call wall_time(end_int)
@ -224,52 +210,27 @@ program QuAcK
call orthogonalization_matrix(ortho_type,nBas,S,X)
!------------------------------------------------------------------------
! Compute RHF energy
! Hartree-Fock module
!------------------------------------------------------------------------
if(doRHF) then
doHF = doRHF .or. doUHF .or. doRMOM .or. doUMOM
! Check that RHF calculation is worth doing...
if(nO(1) /= nO(2)) then
write(*,*) ' !!! The system does not appear to be closed shell !!!'
write(*,*)
stop
end if
if(doHF) then
call wall_time(start_HF)
call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,F_AO,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc)
call HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF_HF,thresh_HF,n_diis_HF, &
guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F_AO, &
ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RHF = ',t_HF,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for HF = ',t_HF,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute UHF energy
!------------------------------------------------------------------------
if(doUHF) then
! Switch on the unrestricted flag
unrestricted = .true.
call cpu_time(start_HF)
call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,eHF,cHF,PHF,Vxc)
call cpu_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for UHF = ',t_HF,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute KS energy
! KS module
!------------------------------------------------------------------------
if(doKS) then
@ -290,33 +251,6 @@ program QuAcK
end if
!------------------------------------------------------------------------
! Maximum overlap method
!------------------------------------------------------------------------
if(doMOM) then
call cpu_time(start_HF)
if(unrestricted) then
! call UMOM()
else
call MOM(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc)
end if
call cpu_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MOM = ',t_HF,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! AO to MO integral transform for post-HF methods
!------------------------------------------------------------------------
@ -327,81 +261,70 @@ program QuAcK
write(*,*) 'AO to MO transformation... Please be patient'
write(*,*)
if(doSph) then
if(unrestricted) then
allocate(ERI_MO(nBas,nBas,nBas,nBas))
ERI_MO(:,:,:,:) = ERI_AO(:,:,:,:)
print*,'!!! MO = AO !!!'
deallocate(ERI_AO)
! Read and transform dipole-related integrals
allocate(dipole_int_aa(nBas,nBas,ncart),dipole_int_bb(nBas,nBas,ncart))
dipole_int_aa(:,:,:) = dipole_int_AO(:,:,:)
dipole_int_bb(:,:,:) = dipole_int_AO(:,:,:)
do ixyz=1,ncart
call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_aa(:,:,ixyz))
call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_bb(:,:,ixyz))
end do
! Memory allocation
allocate(ERI_MO_aaaa(nBas,nBas,nBas,nBas),ERI_MO_aabb(nBas,nBas,nBas,nBas),ERI_MO_bbbb(nBas,nBas,nBas,nBas))
! 4-index transform for (aa|aa) block
bra1 = 1
bra2 = 1
ket1 = 1
ket2 = 1
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_aaaa)
! 4-index transform for (aa|bb) block
bra1 = 1
bra2 = 1
ket1 = 2
ket2 = 2
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_aabb)
! 4-index transform for (bb|bb) block
bra1 = 2
bra2 = 2
ket1 = 2
ket2 = 2
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_bbbb)
else
if(unrestricted) then
! Memory allocation
allocate(ERI_MO(nBas,nBas,nBas,nBas))
allocate(F_MO(nBas,nBas))
! Read and transform dipole-related integrals
allocate(dipole_int_aa(nBas,nBas,ncart),dipole_int_bb(nBas,nBas,ncart))
dipole_int_aa(:,:,:) = dipole_int_AO(:,:,:)
dipole_int_bb(:,:,:) = dipole_int_AO(:,:,:)
do ixyz=1,ncart
call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_aa(:,:,ixyz))
call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_bb(:,:,ixyz))
end do
! Read and transform dipole-related integrals
dipole_int_MO(:,:,:) = dipole_int_AO(:,:,:)
do ixyz=1,ncart
call AOtoMO_transform(nBas,cHF,dipole_int_MO(:,:,ixyz))
end do
! Memory allocation
allocate(ERI_MO_aaaa(nBas,nBas,nBas,nBas),ERI_MO_aabb(nBas,nBas,nBas,nBas),ERI_MO_bbbb(nBas,nBas,nBas,nBas))
! 4-index transform for (aa|aa) block
bra1 = 1
bra2 = 1
ket1 = 1
ket2 = 1
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_aaaa)
! 4-index transform for (aa|bb) block
bra1 = 1
bra2 = 1
ket1 = 2
ket2 = 2
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_aabb)
! 4-index transform for (bb|bb) block
bra1 = 2
bra2 = 2
ket1 = 2
ket2 = 2
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_bbbb)
! 4-index transform
bra1 = 1
bra2 = 1
ket1 = 1
ket2 = 1
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO)
else
! Memory allocation
allocate(ERI_MO(nBas,nBas,nBas,nBas))
allocate(F_MO(nBas,nBas))
! Read and transform dipole-related integrals
dipole_int_MO(:,:,:) = dipole_int_AO(:,:,:)
do ixyz=1,ncart
call AOtoMO_transform(nBas,cHF,dipole_int_MO(:,:,ixyz))
end do
! 4-index transform
bra1 = 1
bra2 = 1
ket1 = 1
ket2 = 1
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO)
F_MO(:,:) = F_AO(:,:)
call AOtoMO_transform(nBas,cHF,F_MO)
end if
F_MO(:,:) = F_AO(:,:)
call AOtoMO_transform(nBas,cHF,F_MO)
end if
@ -917,8 +840,6 @@ program QuAcK
! Perform G0W0 calculatiom
!------------------------------------------------------------------------
eG0W0(:,:) = eHF(:,:)
if(doG0W0) then
call cpu_time(start_GW)
@ -926,11 +847,11 @@ program QuAcK
call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,dophBSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, &
linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, &
dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0W0)
dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc)
else
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE,singlet,triplet, &
linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc)
end if
call cpu_time(end_GW)
@ -953,13 +874,13 @@ program QuAcK
call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,dophBSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc, &
EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb, &
PHF,cHF,eHF,Vxc,eG0W0)
PHF,cHF,eHF,Vxc)
else
call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX, &
dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE,singlet,triplet,linGW,eta_GW,regGW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc)
end if
call cpu_time(end_GW)
@ -1042,7 +963,7 @@ program QuAcK
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0W0 = ',t_GW,' seconds'
write(*,*)
if(dophBSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0)
! if(dophBSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0)
end if
@ -1061,7 +982,7 @@ program QuAcK
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufGW = ',t_GW,' seconds'
write(*,*)
if(dophBSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0)
! if(dophBSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0)
end if
@ -1069,8 +990,6 @@ program QuAcK
! Perform G0T0pp calculatiom
!------------------------------------------------------------------------
eG0T0(:,:) = eHF(:,:)
if(doG0T0pp) then
call cpu_time(start_GT)
@ -1081,14 +1000,13 @@ program QuAcK
call UG0T0(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,evDyn, &
spin_conserved,spin_flip,linGT,eta_GT,regGT,nBas,nC,nO,nV, &
nR,nS,ENuc,EUHF,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, &
dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0T0)
dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc)
else
! call soG0T0(eta_GT,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF)
call G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,evDyn,doppBSE,singlet,triplet, &
linGT,eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO, &
PHF,cHF,eHF,Vxc,eG0T0)
linGT,eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc)
end if
@ -1114,14 +1032,14 @@ program QuAcK
dophBSE,TDA_T,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,&
eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,ERI_AO, &
ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa, &
dipole_int_bb,PHF,cHF,eHF,Vxc,eG0T0)
dipole_int_bb,PHF,cHF,eHF,Vxc)
else
call evGTpp(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, &
dophBSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GT,regGT, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO, &
PHF,cHF,eHF,Vxc,eG0T0)
PHF,cHF,eHF,Vxc)
end if
@ -1168,8 +1086,6 @@ program QuAcK
! Perform G0T0eh calculatiom
!------------------------------------------------------------------------
eG0T0(:,:) = eHF(:,:)
if(doG0T0eh) then
call cpu_time(start_GT)
@ -1181,7 +1097,7 @@ program QuAcK
else
call G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE,singlet,triplet, &
linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0T0)
linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc)
end if
@ -1206,7 +1122,7 @@ program QuAcK
call evGTeh(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, &
dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,evDyn,doppBSE,singlet,triplet,linGT,eta_GT,regGT, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0T0)
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc)
end if
call cpu_time(end_GT)

24
src/QuAcK/read_methods.f90
View File

@ -1,4 +1,4 @@
subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
subroutine read_methods(doRHF,doUHF,doRMOM,doUMOM,doKS, &
doMP2,doMP3, &
doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_crCCD,do_lCCD, &
@ -17,7 +17,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
! Input variables
logical,intent(out) :: doRHF,doUHF,doKS,doMOM
logical,intent(out) :: doRHF,doUHF,doRMOM,doUMOM,doKS
logical,intent(out) :: doMP2,doMP3
logical,intent(out) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
logical,intent(out) :: do_drCCD,do_rCCD,do_crCCD,do_lCCD
@ -38,10 +38,11 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
! Set all the booleans to false
doRHF = .false.
doUHF = .false.
doKS = .false.
doMOM = .false.
doRHF = .false.
doUHF = .false.
doRMOM = .false.
doUMOM = .false.
doKS = .false.
doMP2 = .false.
doMP3 = .false.
@ -91,11 +92,12 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
! Read mean-field methods
read(1,*)
read(1,*) answer1,answer2,answer3,answer4
if(answer1 == 'T') doRHF = .true.
if(answer2 == 'T') doUHF = .true.
if(answer3 == 'T') doKS = .true.
if(answer4 == 'T') doMOM = .true.
read(1,*) answer1,answer2,answer3,answer4,answer5
if(answer1 == 'T') doRHF = .true.
if(answer2 == 'T') doUHF = .true.
if(answer3 == 'T') doRMOM = .true.
if(answer4 == 'T') doUMOM = .true.
if(answer5 == 'T') doKS = .true.
! Read MPn methods