PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:34:46 +01:00
Code Releases Activity

cleanup

Browse Source
This commit is contained in:
Pierre-Francois Loos 2019-03-19 10:13:18 +01:00
parent 1623792de0
commit fdba3ff0ee
124 changed files with 4 additions and 12236 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
GoDuck
View File

@ -13,6 +13,6 @@ then
cp examples/basis."$1"."$2" input/basis
cp examples/basis."$1"."$2" input/weight
./bin/IntPak
./bin/MCQC
./bin/QuAcK
fi

2
GoSph
View File

@ -13,5 +13,5 @@ cp int/Sph_ERI_"$1".dat int/ERI.dat
cp int/Sph_Kin_"$1".dat int/Kin.dat
cp int/Sph_Nuc_"$1".dat int/Nuc.dat
cp int/Sph_Ov_"$1".dat int/Ov.dat
./bin/MCQC
./bin/QuAcK
fi

4
input/molecule
View File

@ -1,4 +1,4 @@
# nAt nEl nCore nRyd
1 2 0 0
# nAt nEl nEla nElb nCore nRyd
1 2 0 0 0 0
# Znuc x y z
He 0.0 0.0 0.0

48
src/MCQC/ADC.f90
View File

@ -1,48 +0,0 @@
subroutine ADC(singlet_manifold,triplet_manifold,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
! ADC main routine
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: singlet_manifold,triplet_manifold
integer,intent(in) :: maxSCF
double precision,intent(in) :: thresh
integer,intent(in) :: max_diis
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: ispin
! Hello world
write(*,*)
write(*,*)'**********************'
write(*,*)'| ADC(n) module |'
write(*,*)'**********************'
write(*,*)
! ADC(2) calculation for singlet manifold
if(singlet_manifold) then
ispin = 1
call ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
endif
! ADC(2) calculation for triplet manifold
if(triplet_manifold) then
ispin = 2
call ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
endif
end subroutine ADC

359
src/MCQC/ADC2.f90
View File

@ -1,359 +0,0 @@
subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
! Compute ADC(2) excitation energies: see Schirmer, Cederbaum & Walter, PRA, 28 (1983) 1237
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: ispin
integer,intent(in) :: maxSCF
double precision,intent(in) :: thresh
integer,intent(in) :: max_diis
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: nH,nP,nHH,nPP,nSCF,n_diis
double precision :: Conv
double precision,external :: Kronecker_delta
double precision,allocatable :: B_ADC(:,:),X_ADC(:,:),e_ADC(:),SigInf(:,:),G_ADC(:,:)
double precision,allocatable :: db_ERI(:,:,:,:),eOld(:),error_diis(:,:),e_diis(:,:)
integer :: i,j,k,l
integer :: a,b,c,d
integer :: p,q,r,s
integer :: nADC,iADC,jADC
! Hello world
write(*,*)
write(*,*)'***********************************'
write(*,*)'| 2nd-order ADC calculation |'
write(*,*)'***********************************'
write(*,*)
! Number of holes
nH = nO - nC
nHH = nH*(nH+1)/2
! Number of particles
nP = nV - nR
nPP = nP*(nP+1)/2
write(*,*) 'Total states: ',nH+nP
write(*,*) 'Hole states: ',nH
write(*,*) 'Particle states: ',nP
! Size of ADC(2) matrices
nADC = nH + nP + nH*nPP + nHH*nP
write(*,'(1X,A25,I3,A6,I6)') 'Size of ADC(2) matrices: ',nADC,' x ',nADC
! Memory allocation
allocate(db_ERI(nBas,nBas,nBas,nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis),eOld(nADC), &
B_ADC(nADC,nADC),X_ADC(nADC,nADC),e_ADC(nADC),G_ADC(nADC,nADC),SigInf(nADC,nADC))
! Create double-bar MO integrals
call antisymmetrize_ERI(ispin,nBas,ERI,db_ERI)
! Initialization
Conv = 1d0
nSCF = 0
n_diis = 0
e_diis(:,:) = 0d0
error_diis(:,:) = 0d0
SigInf(:,:) = 0d0
eOld(:) = 0d0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
!
! | e + SigInf (U^I)^t (U^II)^t |
! | |
! B = | U^I K^I + C^I 0 |
! | |
! | U^II 0 K^II + C^II |
!
!
do while(Conv > thresh .and. nSCF < maxSCF)
!
! Build ADC(2) B matrix -- Eq. (38b) --
!
write(*,'(1X,A7,1X,I4)') 'Cycle: ',nSCF
!
! Diagonal part: static self-energy and epsilon
!
B_ADC(:,:) = 0d0
B_ADC(nC+1:nV,nC+1:nV) = SigInf(nC+1:nV,nC+1:nV)
jADC = 0
do p=nC+1,nV
jADC = jADC + 1
B_ADC(jADC,jADC) = e(p)
enddo
!
! U matrices -- Eq. (40a) --
!
do p=nC+1,nV
iADC = p - nC
jADC = nH + nP
! U^I: 2p-1h -- Eqs. (40a) and (41a) --
do i=nC+1,nO
do a=nO+1,nV-nR
do b=a,nV-nR
jADC = jADC + 1
B_ADC(iADC,jADC) = db_ERI(p,i,a,b)
enddo
enddo
enddo
! U^II: 2h-1p -- Eqs. (40a) and (41b) --
do i=nC+1,nO
do j=i,nO
do a=nO+1,nV-nR
jADC = jADC + 1
B_ADC(iADC,jADC) = db_ERI(p,a,i,j)
enddo
enddo
enddo
enddo
!
! K matrices -- Eq. (40b) --
!
! K^I: 2p-1h -- Eqs. (40b) and (41a) --
jADC = nH + nP
do i=nC+1,nO
do a=nO+1,nV-nR
do b=a,nV-nR
jADC = jADC + 1
B_ADC(jADC,jADC) = e(a) + e(b) - e(i)
enddo
enddo
enddo
! K^II: 2h-1p -- Eqs. (40b) and (41b) --
do i=nC+1,nO
do j=i,nO
do a=nO+1,nV
jADC = jADC + 1
B_ADC(jADC,jADC) = e(i) + e(j) - e(a)
enddo
enddo
enddo
!
! C matrices -- Eq. (42c)
!
! C^I: 2p-1h-TDA -- Eqs. (42a) and (42c) --
iADC = nH + nP
do i=nC+1,nO
do a=nO+1,nV-nR
do b=a,nV-nR
iADC = iADC + 1
jADC = nH + nP
do j=nC+1,nO
do c=nO+1,nV
do d=c,nV-nR
jADC = jADC + 1
B_ADC(iADC,jADC) = B_ADC(iADC,jADC) &
+ db_ERI(a,b,c,d)*Kronecker_delta(i,j) &
- db_ERI(j,b,i,d)*Kronecker_delta(a,c) &
- db_ERI(j,a,i,c)*Kronecker_delta(b,d) &
+ db_ERI(b,a,c,d)*Kronecker_delta(i,j) &
- db_ERI(j,a,i,d)*Kronecker_delta(b,c) &
- db_ERI(j,b,i,c)*Kronecker_delta(a,d)
enddo
enddo
enddo
enddo
enddo
enddo
! C^II: 2p-1h-TDA -- Eqs. (42b) and (42c) --
iADC = nH + nP + nH * nPP
do i=nC+1,nO
do j=i,nO
do a=nO+1,nV-nR
iADC = iADC + 1
jADC = nH + nP + nH*nPP
do k=nC+1,nO
do l=k,nO
do b=nO+1,nV-nR
jADC = jADC + 1
B_ADC(iADC,jADC) = B_ADC(iADC,jADC) &
- db_ERI(i,j,k,l)*Kronecker_delta(a,b) &
+ db_ERI(b,j,a,l)*Kronecker_delta(i,k) &
+ db_ERI(b,i,a,k)*Kronecker_delta(j,l) &
- db_ERI(j,i,k,l)*Kronecker_delta(a,b) &
+ db_ERI(b,i,a,l)*Kronecker_delta(j,k) &
+ db_ERI(b,j,a,k)*Kronecker_delta(i,l)
enddo
enddo
enddo
enddo
enddo
enddo
! Fold B onto itself
do iADC=1,nADC
do jADC=iADC+1,nADC
B_ADC(jADC,iADC) = B_ADC(iADC,jADC)
enddo
enddo
! Diagonalize B to obtain X and E -- Eq. (38a) --
X_ADC(:,:) = B_ADC(:,:)
call diagonalize_matrix(nADC,X_ADC,e_ADC)
! print results
write(*,*) '================================='
write(*,*) 'ADC(2) excitation energies (eV)'
do iADC=1,nADC
if(NORM2(X_ADC(1:nH+nP,iADC)) > 0.1d0 ) &
write(*,'(2(2X,F12.6))') e_ADC(iADC)*HaToeV,NORM2(X_ADC(1:nH+nP,iADC))
enddo
write(*,*) '================================='
! Convergence criteria
Conv = maxval(abs(e_ADC - eOld))
! Store result for next iteration
eOld(:) = e_ADC(:)
! Compute W -- Eq (11) --
SigInf(:,:) = 0d0
do i=nC+1,nO
do p=nC+1,nV-nR
do q=nC+1,nV-nR
SigInf(p,q) = SigInf(p,q) - db_ERI(p,i,q,i)
enddo
enddo
enddo
! Compute the one-particle Greeen function -- Eq. (28) --
G_ADC(:,:) = 0d0
do iADC=1,nADC
if(e_ADC(iADC) > 0d0 ) cycle
do p=nC+1,nV-nR
do q=nC+1,nV-nR
G_ADC(p,q) = G_ADC(p,q) + X_ADC(p,iADC)*X_ADC(q,iADC)
enddo
enddo
enddo
! Compute static self-energy for next iteration -- Eq. (25) --
do p=nC+1,nV-nR
do q=nC+1,nV-nR
do r=nC+1,nV-nR
do s=nC+1,nV-nR
SigInf(p,q) = SigInf(p,q) + db_ERI(p,r,q,s)*G_ADC(r,s)
enddo
enddo
enddo
enddo
! Print results
! call print_ADC2(nBas,nO,nSCF,Conv,e,eADC)
! Increment
nSCF = nSCF + 1
enddo
!------------------------------------------------------------------------
! End main SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF+1) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
endif
end subroutine ADC2

108
src/MCQC/AO_values.f90
View File

@ -1,108 +0,0 @@
subroutine AO_values(doDrift,nBas,nShell,nWalk,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,r,AO,dAO)
! Compute values of the AOs and their derivatives (if required)
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: doDrift
integer,intent(in) :: nBas,nShell,nWalk
double precision,intent(in) :: CenterShell(maxShell,3)
integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell)
double precision,intent(in) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK)
double precision,intent(in) :: r(nWalk,2,3)
! Local variables
integer :: atot,nShellFunction,a(3)
integer,allocatable :: ShellFunction(:,:)
double precision :: rASq,xA,yA,zA,NormCoeff,prim
integer :: iSh,iShF,iK,iW,iEl,iBas,ixyz
! Output variables
double precision,intent(out) :: AO(nWalk,2,nBas),dAO(nWalk,2,3,nBas)
! Initialization
AO = 0d0
if(doDrift) dAO = 0d0
iBas = 0
!------------------------------------------------------------------------
! Loops over shells
!------------------------------------------------------------------------
do iSh=1,nShell
atot = TotAngMomShell(iSh)
nShellFunction = (atot*atot + 3*atot + 2)/2
allocate(ShellFunction(1:nShellFunction,1:3))
call generate_shell(atot,nShellFunction,ShellFunction)
do iShF=1,nShellFunction
iBas = iBas + 1
a(1) = ShellFunction(iShF,1)
a(2) = ShellFunction(iShF,2)
a(3) = ShellFunction(iShF,3)
do iW=1,nWalk
do iEl=1,2
xA = r(iW,iEl,1) - CenterShell(iSh,1)
yA = r(iW,iEl,2) - CenterShell(iSh,2)
zA = r(iW,iEl,3) - CenterShell(iSh,3)
! Calculate distance for exponential
rASq = xA**2 + yA**2 + zA**2
!------------------------------------------------------------------------
! Loops over contraction degrees
!-------------------------------------------------------------------------
do iK=1,KShell(iSh)
! Calculate the exponential part
prim = DShell(iSh,iK)*NormCoeff(ExpShell(iSh,iK),a)*exp(-ExpShell(iSh,iK)*rASq)
AO(iW,iEl,iBas) = AO(iW,iEl,iBas) + prim
if(doDrift) then
prim = -2d0*ExpShell(iSh,iK)*prim
do ixyz=1,3
dAO(iW,iEl,ixyz,iBas) = dAO(iW,iEl,ixyz,iBas) + prim
enddo
endif
enddo
if(doDrift) then
dAO(iW,iEl,1,iBas) = xA**(a(1)+1)*yA**a(2)*zA**a(3)*dAO(iW,iEl,1,iBas)
if(a(1) > 0) dAO(iW,iEl,1,iBas) = dAO(iW,iEl,1,iBas) + dble(a(1))*xA**(a(1)-1)*yA**a(2)*zA**a(3)*AO(iW,iEl,iBas)
dAO(iW,iEl,2,iBas) = xA**a(1)*yA**(a(2)+1)*zA**a(3)*dAO(iW,iEl,2,iBas)
if(a(2) > 0) dAO(iW,iEl,2,iBas) = dAO(iW,iEl,2,iBas) + dble(a(2))*xA**a(1)*yA**(a(2)-1)*zA**a(3)*AO(iW,iEl,iBas)
dAO(iW,iEl,3,iBas) = xA**a(1)*yA**a(2)*zA**(a(3)+1)*dAO(iW,iEl,3,iBas)
if(a(3) > 0) dAO(iW,iEl,3,iBas) = dAO(iW,iEl,3,iBas) + dble(a(3))*xA**a(1)*yA**a(2)*zA**(a(3)-1)*AO(iW,iEl,iBas)
endif
! Calculate polynmial part
AO(iW,iEl,iBas) = xA**a(1)*yA**a(2)*zA**a(3)*AO(iW,iEl,iBas)
enddo
enddo
enddo
deallocate(ShellFunction)
enddo
!------------------------------------------------------------------------
! End loops over shells
!------------------------------------------------------------------------
end subroutine AO_values

81
src/MCQC/AOtoMO_integral_transform.f90
View File

@ -1,81 +0,0 @@
subroutine AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis)
! AO to MO transformation of two-electron integrals
! Semi-direct O(N^5) algorithm
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas),c(nBas,nBas)
! Local variables
double precision,allocatable :: scr(:,:,:,:)
integer :: mu,nu,la,si,i,j,k,l
! Output variables
double precision,intent(out) :: ERI_MO_basis(nBas,nBas,nBas,nBas)
! Memory allocation
allocate(scr(nBas,nBas,nBas,nBas))
scr(:,:,:,:) = 0d0
do l=1,nBas
do si=1,nBas
do la=1,nBas
do nu=1,nBas
do mu=1,nBas
scr(mu,nu,la,l) = scr(mu,nu,la,l) + ERI_AO_basis(mu,nu,la,si)*c(si,l)
enddo
enddo
enddo
enddo
enddo
ERI_MO_basis(:,:,:,:) = 0d0
do l=1,nBas
do la=1,nBas
do nu=1,nBas
do i=1,nBas
do mu=1,nBas
ERI_MO_basis(i,nu,la,l) = ERI_MO_basis(i,nu,la,l) + c(mu,i)*scr(mu,nu,la,l)
enddo
enddo
enddo
enddo
enddo
scr(:,:,:,:) = 0d0
do l=1,nBas
do k=1,nBas
do la=1,nBas
do nu=1,nBas
do i=1,nBas
scr(i,nu,k,l) = scr(i,nu,k,l) + ERI_MO_basis(i,nu,la,l)*c(la,k)
enddo
enddo
enddo
enddo
enddo
ERI_MO_basis(:,:,:,:) = 0d0
do l=1,nBas
do k=1,nBas
do j=1,nBas
do i=1,nBas
do nu=1,nBas
ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j)*scr(i,nu,k,l)
enddo
enddo
enddo
enddo
enddo
end subroutine AOtoMO_integral_transform

85
src/MCQC/AOtoMO_oooa.f90
View File

@ -1,85 +0,0 @@
subroutine AOtoMO_oooa(nBas,nO,nA,cO,cA,O,ooOoa)
! AO to MO transformation of two-electron integrals for the block oooa
! Semi-direct O(N^5) algorithm
implicit none
! Input variables
integer,intent(in) :: nBas,nO,nA
double precision,intent(in) :: cO(nBas,nO),cA(nBas,nA),O(nBas,nBas,nBas,nBas)
! Local variables
double precision,allocatable :: scr1(:,:,:,:),scr2(:,:,:,:)
integer :: mu,nu,la,si,i,j,k,x
! Output variables
double precision,intent(out) :: ooOoa(nO,nO,nO,nA)
! Memory allocation
allocate(scr1(nBas,nBas,nBas,nBas),scr2(nBas,nBas,nBas,nBas))
write(*,*)
write(*,'(A42)') '----------------------------------------'
write(*,'(A42)') ' AO to MO transformation for oooa block '
write(*,'(A42)') '----------------------------------------'
write(*,*)
scr1 = 0d0
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do si=1,nBas
do x=1,nA
scr1(mu,nu,la,x) = scr1(mu,nu,la,x) + O(mu,nu,la,si)*cA(si,x)
enddo
enddo
enddo
enddo
enddo
scr2 = 0d0
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do i=1,nO
do x=1,nA
scr2(i,nu,la,x) = scr2(i,nu,la,x) + cO(mu,i)*scr1(mu,nu,la,x)
enddo
enddo
enddo
enddo
enddo
scr1 = 0d0
do nu=1,nBas
do la=1,nBas
do i=1,nO
do k=1,nO
do x=1,nA
scr1(i,nu,k,x) = scr1(i,nu,k,x) + scr2(i,nu,la,x)*cO(la,k)
enddo
enddo
enddo
enddo
enddo
ooOoa = 0d0
do nu=1,nBas
do i=1,nO
do j=1,nO
do k=1,nO
do x=1,nA
ooOoa(i,j,k,x) = ooOoa(i,j,k,x) + cO(nu,j)*scr1(i,nu,k,x)
enddo
enddo
enddo
enddo
enddo
deallocate(scr1,scr2)
end subroutine AOtoMO_oooa

85
src/MCQC/AOtoMO_oooo.f90
View File

@ -1,85 +0,0 @@
subroutine AOtoMO_oooo(nBas,nO,cO,O,ooOoo)
! AO to MO transformation of two-electron integrals for the block oooo
! Semi-direct O(N^5) algorithm
implicit none
! Input variables
integer,intent(in) :: nBas,nO
double precision,intent(in) :: cO(nBas,nO),O(nBas,nBas,nBas,nBas)
! Local variables
double precision,allocatable :: scr1(:,:,:,:),scr2(:,:,:,:)
integer :: mu,nu,la,si,i,j,k,l
! Output variables
double precision,intent(out) :: ooOoo(nO,nO,nO,nO)
! Memory allocation
allocate(scr1(nBas,nBas,nBas,nBas),scr2(nBas,nBas,nBas,nBas))
write(*,*)
write(*,'(A42)') '----------------------------------------'
write(*,'(A42)') ' AO to MO transformation for oooo block '
write(*,'(A42)') '----------------------------------------'
write(*,*)
scr1 = 0d0
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do si=1,nBas
do l=1,nO
scr1(mu,nu,la,l) = scr1(mu,nu,la,l) + O(mu,nu,la,si)*cO(si,l)
enddo
enddo
enddo
enddo
enddo
scr2 = 0d0
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do i=1,nO
do l=1,nO
scr2(i,nu,la,l) = scr2(i,nu,la,l) + cO(mu,i)*scr1(mu,nu,la,l)
enddo
enddo
enddo
enddo
enddo
scr1 = 0d0
do nu=1,nBas
do la=1,nBas
do i=1,nO
do k=1,nO
do l=1,nO
scr1(i,nu,k,l) = scr1(i,nu,k,l) + scr2(i,nu,la,l)*cO(la,k)
enddo
enddo
enddo
enddo
enddo
ooOoo = 0d0
do nu=1,nBas
do i=1,nO
do j=1,nO
do k=1,nO
do l=1,nO
ooOoo(i,j,k,l) = ooOoo(i,j,k,l) + cO(nu,j)*scr1(i,nu,k,l)
enddo
enddo
enddo
enddo
enddo
deallocate(scr1,scr2)
end subroutine AOtoMO_oooo

135
src/MCQC/AOtoMO_oooooo.f90
View File

@ -1,135 +0,0 @@
subroutine AOtoMO_oooooo(nBas,nO,cO,O,oooOooo)
! AO to MO transformation of three-electron integrals for the block oooooo
! Semi-direct O(N^7) algorithm
implicit none
! Input variables
integer,intent(in) :: nBas,nO
double precision,intent(in) :: cO(nBas,nO),O(nBas,nBas,nBas,nBas,nBas,nBas)
! Local variables
double precision,allocatable :: scr1(:,:,:,:,:,:),scr2(:,:,:,:,:,:)
integer :: mu,nu,la,si,ka,ta,i,j,k,l,m,n
! Output variables
double precision,intent(out) :: oooOooo(nO,nO,nO,nO,nO,nO)
! Memory allocation
allocate(scr1(nBas,nBas,nBas,nBas,nBas,nBas),scr2(nBas,nBas,nBas,nBas,nBas,nBas))
write(*,*)
write(*,'(A42)') '------------------------------------------'
write(*,'(A42)') ' AO to MO transformation for oooooo block '
write(*,'(A42)') '------------------------------------------'
write(*,*)
scr1 = 0d0
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do si=1,nBas
do ka=1,nBas
do ta=1,nBas
do n=1,nO
scr1(mu,nu,la,si,ka,n) = scr1(mu,nu,la,si,ka,n) + O(mu,nu,la,si,ka,ta)*cO(ta,n)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
scr2 = 0d0
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do si=1,nBas
do ka=1,nBas
do i=1,nO
do n=1,nO
scr2(i,nu,la,si,ka,n) = scr2(i,nu,la,si,ka,n) + cO(mu,i)*scr1(mu,nu,la,si,ka,n)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
scr1 = 0d0
do nu=1,nBas
do la=1,nBas
do si=1,nBas
do ka=1,nBas
do i=1,nO
do m=1,nO
do n=1,nO
scr1(i,nu,la,si,m,n) = scr1(i,nu,la,si,m,n) + scr2(i,nu,la,si,m,n)*cO(ka,m)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
scr2 = 0d0
do nu=1,nBas
do la=1,nBas
do si=1,nBas
do i=1,nO
do j=1,nO
do m=1,nO
do n=1,nO
scr2(i,j,la,si,m,n) = scr2(i,j,la,si,m,n) + cO(nu,j)*scr1(i,nu,la,si,m,n)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
scr1 = 0d0
do la=1,nBas
do si=1,nBas
do i=1,nO
do j=1,nO
do l=1,nO
do m=1,nO
do n=1,nO
scr1(i,j,la,l,m,n) = scr1(i,j,la,l,m,n) + scr2(i,j,la,si,m,n)*cO(si,l)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
oooOooo = 0d0
do si=1,nBas
do i=1,nO
do j=1,nO
do k=1,nO
do l=1,nO
do m=1,nO
do n=1,nO
oooOooo(i,j,k,l,m,n) = oooOooo(i,j,k,l,m,n) + cO(la,k)*scr1(i,j,la,k,m,n)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
deallocate(scr1,scr2)
end subroutine AOtoMO_oooooo

77
src/MCQC/AOtoMO_oovv.f90
View File

@ -1,77 +0,0 @@
subroutine AOtoMO_oovv(nBas,nO,nV,cO,cV,O,ooOvv)
! AO to MO transformation of two-electron integrals for the block oovv
! Semi-direct O(N^5) algorithm
implicit none
! Input variables
integer,intent(in) :: nBas,nO,nV
double precision,intent(in) :: cO(nBas,nO),cV(nBas,nV),O(nBas,nBas,nBas,nBas)
! Local variables
double precision,allocatable :: scr1(:,:,:,:),scr2(:,:,:,:)
integer :: mu,nu,la,si,i,j,a,b
! Output variables
double precision,intent(out) :: ooOvv(nO,nO,nV,nV)
! Memory allocation
allocate(scr1(nBas,nBas,nBas,nBas),scr2(nBas,nBas,nBas,nBas))
scr1 = 0d0
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do si=1,nBas
do b=1,nV
scr1(mu,nu,la,b) = scr1(mu,nu,la,b) + O(mu,nu,la,si)*cV(si,b)
enddo
enddo
enddo
enddo
enddo
scr2 = 0d0
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do i=1,nO
do b=1,nV
scr2(i,nu,la,b) = scr2(i,nu,la,b) + cO(mu,i)*scr1(mu,nu,la,b)
enddo
enddo
enddo
enddo
enddo
scr1 = 0d0
do nu=1,nBas
do la=1,nBas
do i=1,nO
do a=1,nV
do b=1,nV
scr1(i,nu,a,b) = scr1(i,nu,a,b) + scr2(i,nu,la,b)*cV(la,a)
enddo
enddo
enddo
enddo
enddo
ooOvv = 0d0
do nu=1,nBas
do i=1,nO
do j=1,nO
do a=1,nV
do b=1,nV
ooOvv(i,j,a,b) = ooOvv(i,j,a,b) + cO(nu,j)*scr1(i,nu,a,b)
enddo
enddo
enddo
enddo
enddo
end subroutine AOtoMO_oovv

18
src/MCQC/AOtoMO_transform.f90
View File

@ -1,18 +0,0 @@
subroutine AOtoMO_transform(nBas,c,A)
! Perform AO to MO transformation of a matrix A for given coefficients c
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: c(nBas,nBas)
! Output variables
double precision,intent(inout):: A(nBas,nBas)
A = matmul(transpose(c),matmul(A,c))
end subroutine AOtoMO_transform

44
src/MCQC/Bethe_Salpeter_A_matrix.f90
View File

@ -1,44 +0,0 @@
subroutine Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,A_lr)
! Compute the extra term for Bethe-Salpeter equation for linear response
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega(nS),rho(nBas,nBas,nS)
! Local variables
double precision :: chi
integer :: i,j,a,b,ia,jb,kc
! Output variables
double precision,intent(out) :: A_lr(nS,nS)
ia = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
ia = ia + 1
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
chi = 0d0
do kc=1,nS
chi = chi + rho(i,j,kc)*rho(a,b,kc)/Omega(kc)
enddo
A_lr(ia,jb) = A_lr(ia,jb) - ERI(i,a,j,b) + 4d0*chi
enddo
enddo
enddo
enddo
end subroutine Bethe_Salpeter_A_matrix

44
src/MCQC/Bethe_Salpeter_B_matrix.f90
View File

@ -1,44 +0,0 @@
subroutine Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,B_lr)
! Compute the extra term for Bethe-Salpeter equation for linear response
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega(nS),rho(nBas,nBas,nS)
! Local variables
double precision :: chi
integer :: i,j,a,b,ia,jb,kc
! Output variables
double precision,intent(out) :: B_lr(nS,nS)
ia = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
ia = ia + 1
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
chi = 0d0
do kc=1,nS
chi = chi + rho(i,b,kc)*rho(a,j,kc)/Omega(kc)
enddo
B_lr(ia,jb) = B_lr(ia,jb) - ERI(i,a,b,j) + 4d0*chi
enddo
enddo
enddo
enddo
end subroutine Bethe_Salpeter_B_matrix

203
src/MCQC/CCD.f90
View File

@ -1,203 +0,0 @@
subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
! CCD module
implicit none
! Input variables
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
double precision,intent(in) :: thresh
integer,intent(in) :: nBas,nEl
double precision,intent(in) :: ENuc,ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: nBas2
integer :: nO
integer :: nV
integer :: nSCF
double precision :: Conv
double precision :: EcMP2,EcMP3,EcMP4
double precision :: ECCD,EcCCD
double precision,allocatable :: seHF(:)
double precision,allocatable :: sERI(:,:,:,:)
double precision,allocatable :: dbERI(:,:,:,:)
double precision,allocatable :: eO(:)
double precision,allocatable :: eV(:)
double precision,allocatable :: delta_OOVV(:,:,:,:)
double precision,allocatable :: OOOO(:,:,:,:)
double precision,allocatable :: OOVV(:,:,:,:)
double precision,allocatable :: OVOV(:,:,:,:)
double precision,allocatable :: VVVV(:,:,:,:)
double precision,allocatable :: X1(:,:,:,:)
double precision,allocatable :: X2(:,:)
double precision,allocatable :: X3(:,:)
double precision,allocatable :: X4(:,:,:,:)
double precision,allocatable :: u(:,:,:,:)
double precision,allocatable :: v(:,:,:,:)
double precision,allocatable :: r2(:,:,:,:)
double precision,allocatable :: t2(:,:,:,:)
! Hello world
write(*,*)
write(*,*)'**************************************'
write(*,*)'| CCD calculation |'
write(*,*)'**************************************'
write(*,*)
! Spatial to spin orbitals
nBas2 = 2*nBas
allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI)
! Antysymmetrize ERIs
allocate(dbERI(nBas2,nBas2,nBas2,nBas2))
call antisymmetrize_ERI(2,nBas2,sERI,dbERI)
deallocate(sERI)
! Define occupied and virtual spaces
nO = nEl
nV = nBas2 - nO
! Form energy denominator
allocate(eO(nO),eV(nV))
allocate(delta_OOVV(nO,nO,nV,nV))
eO(:) = seHF(1:nO)
eV(:) = seHF(nO+1:nBas2)
call form_delta_OOVV(nO,nV,eO,eV,delta_OOVV)
deallocate(seHF)
! Create integral batches
allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVOV(nO,nV,nO,nV),VVVV(nV,nV,nV,nV))
OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO )
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas2,nO+1:nBas2)
OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas2, 1:nO ,nO+1:nBas2)
VVVV(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2,nO+1:nBas2,nO+1:nBas2)
deallocate(dbERI)
! MP2 guess amplitudes
allocate(t2(nO,nO,nV,nV))
t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.))
EcMP4 = 0d0
! Initialization
allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(nO,nO,nV,nV))
allocate(X1(nO,nO,nO,nO),X2(nV,nV),X3(nO,nO),X4(nO,nO,nV,nV))
Conv = 1d0
nSCF = 0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
write(*,*)
write(*,*)'----------------------------------------------------'
write(*,*)'| CCD calculation |'
write(*,*)'----------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
'|','#','|','E(CCD)','|','Ec(CCD)','|','Conv','|'
write(*,*)'----------------------------------------------------'
do while(Conv > thresh .and. nSCF < maxSCF)
! Increment
nSCF = nSCF + 1
! Form linear array
call form_u(nO,nV,OOOO,VVVV,OVOV,t2,u)
! Form interemediate arrays
call form_X(nO,nV,OOVV,t2,X1,X2,X3,X4)
! Form quadratic array
call form_v(nO,nV,X1,X2,X3,X4,t2,v)
! Compute residual
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + u(:,:,:,:) + v(:,:,:,:)
! Check convergence
Conv = maxval(abs(r2(:,:,:,:)))
! Update amplitudes
t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:)
! Compute correlation energy
EcCCD = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.))
if(nSCF == 1) EcMP3 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2 + v/delta_OOVV,.true.))
! Dump results
ECCD = ERHF + EcCCD
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
'|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|'
enddo
write(*,*)'----------------------------------------------------'
!------------------------------------------------------------------------
! End of SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
stop
endif
! Moller-Plesset energies
write(*,*)
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP2) = ',EcMP2
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP3) = ',EcMP3
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP4-SDQ) = ',EcMP4
write(*,*)
end subroutine CCD

259
src/MCQC/CCSD.f90
View File

@ -1,259 +0,0 @@
subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF)
! CCSD module
implicit none
! Input variables
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
double precision,intent(in) :: thresh
logical,intent(in) :: doCCSDT
integer,intent(in) :: nBas,nEl
double precision,intent(in) :: ENuc,ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
double precision :: start_CCSDT,end_CCSDT,t_CCSDT
integer :: nBas2
integer :: nO
integer :: nV
integer :: nSCF
double precision :: Conv
double precision :: EcMP2
double precision :: ECCSD,EcCCSD
double precision :: EcCCT
double precision,allocatable :: seHF(:)
double precision,allocatable :: sERI(:,:,:,:)
double precision,allocatable :: dbERI(:,:,:,:)
double precision,allocatable :: delta_OV(:,:)
double precision,allocatable :: delta_OOVV(:,:,:,:)
double precision,allocatable :: OOOO(:,:,:,:)
double precision,allocatable :: OOOV(:,:,:,:)
double precision,allocatable :: OVOO(:,:,:,:)
double precision,allocatable :: VOOO(:,:,:,:)
double precision,allocatable :: OOVV(:,:,:,:)
double precision,allocatable :: OVVO(:,:,:,:)
double precision,allocatable :: OVVV(:,:,:,:)
double precision,allocatable :: VOVV(:,:,:,:)
double precision,allocatable :: VVVO(:,:,:,:)
double precision,allocatable :: VVVV(:,:,:,:)
double precision,allocatable :: eO(:)
double precision,allocatable :: eV(:)
double precision,allocatable :: hvv(:,:)
double precision,allocatable :: hoo(:,:)
double precision,allocatable :: hvo(:,:)
double precision,allocatable :: gvv(:,:)
double precision,allocatable :: goo(:,:)
double precision,allocatable :: aoooo(:,:,:,:)
double precision,allocatable :: bvvvv(:,:,:,:)
double precision,allocatable :: hovvo(:,:,:,:)
double precision,allocatable :: r1(:,:)
double precision,allocatable :: r2(:,:,:,:)
double precision,allocatable :: t1(:,:)
double precision,allocatable :: t2(:,:,:,:)
double precision,allocatable :: tau(:,:,:,:)
! Hello world
write(*,*)
write(*,*)'**************************************'
write(*,*)'| CCSD calculation |'
write(*,*)'**************************************'
write(*,*)
! Spatial to spin orbitals
nBas2 = 2*nBas
allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI)
! Antysymmetrize ERIs
allocate(dbERI(nBas2,nBas2,nBas2,nBas2))
call antisymmetrize_ERI(2,nBas2,sERI,dbERI)
deallocate(sERI)
! Define occupied and virtual spaces
nO = nEl
nV = nBas2 - nO
! Form energy denominator
allocate(eO(nO),eV(nV))
allocate(delta_OV(nO,nV),delta_OOVV(nO,nO,nV,nV))
eO(:) = seHF(1:nO)
eV(:) = seHF(nO+1:nBas2)
call form_delta_OV(nO,nV,eO,eV,delta_OV)
call form_delta_OOVV(nO,nV,eO,eV,delta_OOVV)
deallocate(seHF)
! Create integral batches
allocate(OOOO(nO,nO,nO,nO), &
OOOV(nO,nO,nO,nV),OVOO(nO,nV,nO,nO),VOOO(nV,nO,nO,nO), &
OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO), &
OVVV(nO,nV,nV,nV),VOVV(nV,nO,nV,nV),VVVO(nV,nV,nV,nO), &
VVVV(nV,nV,nV,nV))
OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO )
OOOV(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO ,nO+1:nBas2)
OVOO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas2, 1:nO , 1:nO )
VOOO(:,:,:,:) = dbERI(nO+1:nBas2, 1:nO , 1:nO , 1:nO )
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas2,nO+1:nBas2)
OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas2,nO+1:nBas2, 1:nO )
OVVV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas2,nO+1:nBas2,nO+1:nBas2)
VOVV(:,:,:,:) = dbERI(nO+1:nBas2, 1:nO ,nO+1:nBas2,nO+1:nBas2)
VVVO(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2,nO+1:nBas2, 1:nO )
VVVV(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2,nO+1:nBas2,nO+1:nBas2)
deallocate(dbERI)
! MP2 guess amplitudes
allocate(t1(nO,nV),t2(nO,nO,nV,nV),tau(nO,nO,nV,nV))
t1(:,:) = 0d0
t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
call form_tau(nO,nV,t1,t2,tau)
EcMP2 = 0.5d0*dot_product(pack(OOVV,.true.),pack(tau,.true.))
write(*,'(1X,A10,1X,F10.6)') 'Ec(MP2) = ',EcMP2
! Initialization
allocate(hvv(nV,nV),hoo(nO,nO),hvo(nV,nO), &
gvv(nV,nV),goo(nO,nO), &
aoooo(nO,nO,nO,nO),bvvvv(nV,nV,nV,nV),hovvo(nO,nV,nV,nO), &
r1(nO,nV),r2(nO,nO,nV,nV))
Conv = 1d0
nSCF = 0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
write(*,*)
write(*,*)'----------------------------------------------------'
write(*,*)'| CCSD calculation |'
write(*,*)'----------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
'|','#','|','E(CCSD)','|','Ec(CCSD)','|','Conv','|'
write(*,*)'----------------------------------------------------'
do while(Conv > thresh .and. nSCF < maxSCF)
! Increment
nSCF = nSCF + 1
! Scuseria Eqs. (5), (6) and (7)
call form_h(nO,nV,eO,eV,OOVV,t1,tau,hvv,hoo,hvo)
! Scuseria Eqs. (9), (10), (11), (12) and (13)
call form_g(nO,nV,hvv,hoo,VOVV,OOOV,t1,gvv,goo)
call form_abh(nO,nV,OOOO,OVOO,OOVV,VVVV,VOVV,OVVO,OVVV,t1,tau,aoooo,bvvvv,hovvo)
! Compute residuals
call form_r1(nO,nV,OVVO,OVVV,OOOV,hvv,hoo,hvo,t1,t2,tau,r1)
call form_r2(nO,nV,OOVV,OVOO,OVVV,OVVO,gvv,goo,aoooo,bvvvv,hovvo,t1,t2,tau,r2)
! Check convergence
Conv = max(maxval(abs(r1(:,:))),maxval(abs(r2(:,:,:,:))))
! Update
t1(:,:) = t1(:,:) - r1(:,:) /delta_OV (:,:)
t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:)
call form_tau(nO,nV,t1,t2,tau)
! Compute correlation energy
EcCCSD = 0.5d0*dot_product(pack(OOVV,.true.),pack(tau,.true.))
! Dump results
ECCSD = ERHF + EcCCSD
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
'|',nSCF,'|',ECCSD+ENuc,'|',EcCCSD,'|',Conv,'|'
end do
write(*,*)'----------------------------------------------------'
!------------------------------------------------------------------------
! End of SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
stop
end if
! Deallocate memory
deallocate(hvv,hoo,hvo, &
delta_OV,delta_OOVV, &
gvv,goo, &
aoooo,bvvvv,hovvo, &
tau, &
r1,r2)
!------------------------------------------------------------------------
! (T) correction
!------------------------------------------------------------------------
if(doCCSDT) then
call cpu_time(start_CCSDT)
call CCSDT(nO,nV,eO,eV,OOVV,VVVO,VOOO,t1,t2,EcCCT)
call cpu_time(end_CCSDT)
t_CCSDT = end_CCSDT - start_CCSDT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for (T) = ',t_CCSDT,' seconds'
write(*,*)
write(*,*)
write(*,*)'----------------------------------------------------'
write(*,*)' CCSDT(T) energy '
write(*,*)'----------------------------------------------------'
write(*,'(1X,A20,1X,F15.10)')' E(CCSD(T)) = ',ECCSD + EcCCT
write(*,'(1X,A20,1X,F10.6)') ' Ec(CCSD(T)) = ',EcCCSD + EcCCT
write(*,*)'----------------------------------------------------'
write(*,*)
end if
end subroutine CCSD

45
src/MCQC/CCSDT.f90
View File

@ -1,45 +0,0 @@
subroutine CCSDT(nO,nV,eO,eV,OOVV,VVVO,VOOO,t1,t2,EcCCT)
! Compute the (T) correction of the CCSD(T) energy
implicit none
! Input variables
integer,intent(in) :: nO,nV
double precision,intent(in) :: eO(nO)
double precision,intent(in) :: eV(nV)
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
double precision,intent(in) :: VVVO(nV,nV,nV,nO)
double precision,intent(in) :: VOOO(nV,nO,nO,nO)
double precision,intent(in) :: t1(nO,nV)
double precision,intent(in) :: t2(nO,nO,nV,nV)
! Local variables
double precision,allocatable :: delta_OOOVVV(:,:,:,:,:,:)
double precision,allocatable :: ub(:,:,:,:,:,:)
double precision,allocatable :: ubb(:,:,:,:,:,:)
! Output variables
double precision,intent(out) :: EcCCT
! Memory allocation
allocate(delta_OOOVVV(nO,nO,nO,nV,nV,nV),ub(nO,nO,nO,nV,nV,nV),ubb(nO,nO,nO,nV,nV,nV))
! Form CCSD(T) quantities
call form_delta_OOOVVV(nO,nV,eO,eV,delta_OOOVVV)
call form_ub(nO,nV,OOVV,t1,ub)
call form_ubb(nO,nV,VVVO,VOOO,t2,ubb)
call form_T(nO,nV,delta_OOOVVV,ub,ubb,EcCCT)
end subroutine CCSDT

85
src/MCQC/CIS.f90
View File

@ -1,85 +0,0 @@
subroutine CIS(singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,eHF)
! Perform configuration interaction single calculation`
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: singlet_manifold,triplet_manifold
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),eHF(nBas)
! Local variables
logical :: dRPA
logical :: dump_matrix = .false.
logical :: dump_trans = .false.
integer :: ispin
double precision,allocatable :: A(:,:),Omega(:)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Configuration Interaction Singles |'
write(*,*)'************************************************'
write(*,*)
! Switch on exchange for CIS
dRPA = .false.
! Memory allocation
allocate(A(nS,nS),Omega(nS))
! Compute CIS matrix
if(singlet_manifold) then
ispin = 1
call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,eHF,ERI,A)
if(dump_matrix) then
print*,'CIS matrix (singlet state)'
call matout(nS,nS,A)
write(*,*)
endif
call diagonalize_matrix(nS,A,Omega)
call print_excitation('CIS ',ispin,nS,Omega)
if(dump_trans) then
print*,'Singlet CIS transition vectors'
call matout(nS,nS,A)
write(*,*)
endif
endif
if(triplet_manifold) then
ispin = 2
call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,eHF,ERI,A)
if(dump_matrix) then
print*,'CIS matrix (triplet state)'
call matout(nS,nS,A)
write(*,*)
endif
call diagonalize_matrix(nS,A,Omega)
call print_excitation('CIS ',ispin,nS,Omega)
if(dump_trans) then
print*,'Triplet CIS transition vectors'
call matout(nS,nS,A)
write(*,*)
endif
endif
end subroutine CIS

34
src/MCQC/Coulomb_matrix_AO_basis.f90
View File

@ -1,34 +0,0 @@
subroutine Coulomb_matrix_AO_basis(nBas,P,G,J)
! Compute Coulomb matrix in the AO basis
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: G(nBas,nBas,nBas,nBas)
! Local variables
integer :: mu,nu,la,si
! Output variables
double precision,intent(out) :: J(nBas,nBas)
J = 0d0
do si=1,nBas
do nu=1,nBas
do la=1,nBas
do mu=1,nBas
J(mu,nu) = J(mu,nu) + P(la,si)*G(mu,la,nu,si)
enddo
enddo
enddo
enddo
end subroutine Coulomb_matrix_AO_basis

26
src/MCQC/Coulomb_matrix_MO_basis.f90
View File

@ -1,26 +0,0 @@
subroutine Coulomb_matrix_MO_basis(nBas,c,P,G,J)
! Compute Coulomb matrix in the MO basis
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: c(nBas,nBas),P(nBas,nBas)
double precision,intent(in) :: G(nBas,nBas,nBas,nBas)
! Output variables
double precision,intent(out) :: J(nBas,nBas)
! Compute Hartree Hamiltonian in the AO basis
call Coulomb_matrix_AO_basis(nBas,P,G,J)
! Transform Coulomb matrix in the MO basis
J = matmul(transpose(c),matmul(J,c))
end subroutine Coulomb_matrix_MO_basis

61
src/MCQC/DIIS_extrapolation.f90
View File

@ -1,61 +0,0 @@
subroutine DIIS_extrapolation(n_err,n_e,n_diis,error,e,error_in,e_inout)
! Perform DIIS extrapolation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: n_err,n_e
double precision,intent(in) :: error_in(n_err),error(n_err,n_diis),e(n_e,n_diis)
! Local variables
double precision :: rcond
double precision,allocatable :: A(:,:),b(:),w(:)
! Output variables
integer,intent(inout) :: n_diis
double precision,intent(inout):: e_inout(n_e)
! Memory allocaiton
allocate(A(n_diis+1,n_diis+1),b(n_diis+1),w(n_diis+1))
! Update DIIS "history"
call prepend(n_err,n_diis,error,error_in)
call prepend(n_e,n_diis,e,e_inout)
! Build A matrix
A(1:n_diis,1:n_diis) = matmul(transpose(error),error)
A(1:n_diis,n_diis+1) = -1d0
A(n_diis+1,1:n_diis) = -1d0
A(n_diis+1,n_diis+1) = +0d0
! Build x matrix
b(1:n_diis) = +0d0
b(n_diis+1) = -1d0
! Solve linear system
call linear_solve(n_diis+1,A,b,w,rcond)
! Extrapolate
if(rcond > 1d-14) then
e_inout(:) = matmul(w(1:n_diis),transpose(e(:,1:n_diis)))
else
n_diis = 0
endif
end subroutine DIIS_extrapolation

132
src/MCQC/G0W0.f90
View File

@ -1,132 +0,0 @@
subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,P,ERI_AO_basis,ERI_MO_basis,cHF,eHF,eG0W0)
! Perform G0W0 calculation
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ENuc,ERHF
double precision,intent(in) :: cHF(nBas,nBas),eHF(nBas),Hc(nBas,nBas),P(nBas,nBas)
double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas),ERI_MO_basis(nBas,nBas,nBas,nBas)
! Local variables
logical :: dRPA
integer :: ispin
double precision :: EcRPA,EcGM
double precision,allocatable :: H(:,:),SigmaC(:),Z(:)
double precision,allocatable :: Omega(:,:),XpY(:,:,:),rho(:,:,:,:),rhox(:,:,:,:)
! Output variables
double precision :: eG0W0(nBas)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| One-shot G0W0 calculation |'
write(*,*)'************************************************'
write(*,*)
! SOSEX correction
if(SOSEX) write(*,*) 'SOSEX correction activated!'
write(*,*)
! Switch off exchange for G0W0
dRPA = .true.
! Spin manifold
ispin = 1
! Memory allocation
allocate(H(nBas,nBas),SigmaC(nBas),Z(nBas), &
Omega(nS,nspin),XpY(nS,nS,nspin), &
rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin))
! Compute Hartree Hamiltonian in the MO basis
call Hartree_matrix_MO_basis(nBas,cHF,P,Hc,ERI_AO_basis,H)
! Compute linear response
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin))
! Compute correlation part of the self-energy
call excitation_density_from_MO(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
if(SOSEX) call excitation_density_SOSEX_from_MO(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rhox(:,:,:,ispin))
call self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,eHF, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigmaC)
! COHSEX static approximation
if(COHSEX) then
Z(:) = 1d0
else
call renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,eHF,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z)
endif
! Solve the quasi-particle equation
eG0W0(:) = eHF(:) + Z(:)*SigmaC(:)
! Dump results
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigmaC,Z,eG0W0,EcRPA,EcGM)
! Plot stuff
call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eG0W0,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin))
! Perform BSE calculation
if(BSE) then
! Singlet manifold
if(singlet_manifold) then
ispin = 1
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0W0,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
endif
! Triplet manifold
if(triplet_manifold) then
ispin = 2
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,cHF,ERI_AO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0W0,ERI_MO_basis, &
rho(:,:,:,1),EcRPA,Omega(:,ispin),XpY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
endif
endif
end subroutine G0W0

131
src/MCQC/GF2.f90
View File

@ -1,131 +0,0 @@
subroutine GF2(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,V,e0)
! Perform second-order Green function calculation in diagonal approximation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: maxSCF
double precision,intent(in) :: thresh
integer,intent(in) :: max_diis
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: e0(nBas),V(nBas,nBas,nBas,nBas)
! Local variables
integer :: nSCF,n_diis
double precision :: eps,Conv
double precision,allocatable :: eGF2(:),eOld(:),Bpp(:,:,:),error_diis(:,:),e_diis(:,:)
integer :: i,j,a,b,p,q
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Second-order Green function calculation |'
write(*,*)'************************************************'
write(*,*)
! Memory allocation
allocate(Bpp(nBas,nBas,2),eGF2(nBas),eOld(nBas), &
error_diis(nBas,max_diis),e_diis(nBas,max_diis))
! Initialization
Conv = 1d0
nSCF = 0
n_diis = 0
e_diis(:,:) = 0d0
error_diis(:,:) = 0d0
eGF2(:) = e0(:)
eOld(:) = e0(:)
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
do while(Conv > thresh .and. nSCF < maxSCF)
! Frequency-dependent second-order contribution
Bpp(:,:,:) = 0d0
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = eGF2(p) + e0(a) - e0(i) - e0(j)
Bpp(p,q,1) = Bpp(p,q,1) &
+ (2d0*V(p,a,i,j) - V(p,a,j,i))*V(q,a,i,j)/eps
enddo
enddo
enddo
enddo
enddo
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = eGF2(p) + e0(i) - e0(a) - e0(b)
Bpp(p,q,2) = Bpp(p,q,2) &
+ (2d0*V(p,i,a,b) - V(p,i,b,a))*V(q,i,a,b)/eps
enddo
enddo
enddo
enddo
enddo
print*,'Sig2 in GF2'
call matout(nBas,nBas,Bpp(:,:,1) + Bpp(:,:,2))
! eGF2(:) = e0(:) &
! + Bpp(:,1) + Bpp(:,2)
Conv = maxval(abs(eGF2 - eOld))
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(nBas,nBas,n_diis,error_diis,e_diis,eGF2-eOld,eGF2)
eOld = eGF2
! Print results
call print_GF2(nBas,nO,nSCF,Conv,e0,eGF2)
! Increment
nSCF = nSCF + 1
enddo
!------------------------------------------------------------------------
! End main SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF+1) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
endif
end subroutine GF2

124
src/MCQC/GF2_diag.f90
View File

@ -1,124 +0,0 @@
subroutine GF2_diag(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,V,e0)
! Perform second-order Green function calculation in diagonal approximation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: maxSCF
double precision,intent(in) :: thresh
integer,intent(in) :: max_diis
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: e0(nBas),V(nBas,nBas,nBas,nBas)
! Local variables
integer :: nSCF,n_diis
double precision :: eps,Conv
double precision,allocatable :: eGF2(:),eOld(:),Bpp(:,:),error_diis(:,:),e_diis(:,:)
integer :: i,j,a,b,p
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Second-order Green function calculation |'
write(*,*)'************************************************'
write(*,*)
! Memory allocation
allocate(Bpp(nBas,2),eGF2(nBas),eOld(nBas), &
error_diis(nBas,max_diis),e_diis(nBas,max_diis))
! Initialization
Conv = 1d0
nSCF = 0
n_diis = 0
e_diis(:,:) = 0d0
error_diis(:,:) = 0d0
eGF2(:) = e0(:)
eOld(:) = e0(:)
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
do while(Conv > thresh .and. nSCF < maxSCF)
! Frequency-dependent second-order contribution
Bpp(:,:) = 0d0
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = eGF2(p) + e0(a) - e0(i) - e0(j)
Bpp(p,1) = Bpp(p,1) &
+ (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)/eps
enddo
enddo
enddo
enddo
do p=nC+1,nBas-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = eGF2(p) + e0(i) - e0(a) - e0(b)
Bpp(p,2) = Bpp(p,2) &
+ (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)/eps
enddo
enddo
enddo
enddo
eGF2(:) = e0(:) &
+ Bpp(:,1) + Bpp(:,2)
Conv = maxval(abs(eGF2 - eOld))
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(nBas,nBas,n_diis,error_diis,e_diis,eGF2-eOld,eGF2)
eOld = eGF2
! Print results
call print_GF2(nBas,nO,nSCF,Conv,e0,eGF2)
! Increment
nSCF = nSCF + 1
enddo
!------------------------------------------------------------------------
! End main SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF+1) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
endif
end subroutine GF2_diag

488
src/MCQC/GF3_diag.f90
View File

@ -1,488 +0,0 @@
subroutine GF3_diag(maxSCF,thresh,max_diis,renormalization,nBas,nC,nO,nV,nR,V,e0)
! Perform third-order Green function calculation in diagonal approximation
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: thresh
integer,intent(in) :: maxSCF,max_diis,renormalization
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: e0(nBas),V(nBas,nBas,nBas,nBas)
! Local variables
integer :: nSCF,n_diis
double precision :: eps,eps1,eps2,Conv
double precision,allocatable :: Sig2(:),SigInf(:),Sig3(:),eGF3(:),eOld(:)
double precision,allocatable :: App(:,:),Bpp(:,:),Cpp(:,:),Dpp(:,:)
double precision,allocatable :: Z(:),X2h1p(:),X1h2p(:),Sig2h1p(:),Sig1h2p(:)
double precision,allocatable :: error_diis(:,:),e_diis(:,:)
integer :: i,j,k,l,a,b,c,d,p
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Third-order Green function calculation |'
write(*,*)'************************************************'
write(*,*)
! Memory allocation
allocate(eGF3(nBas),eOld(nBas), &
Sig2(nBas),SigInf(nBas),Sig3(nBas), &
App(nBas,6),Bpp(nBas,2),Cpp(nBas,6),Dpp(nBas,6), &
Z(nBas),X2h1p(nBas),X1h2p(nBas),Sig2h1p(nBas),Sig1h2p(nBas), &
error_diis(nBas,max_diis),e_diis(nBas,max_diis))
!------------------------------------------------------------------------
! Compute third-order frequency-independent contribution
!------------------------------------------------------------------------
App(:,:) = 0d0
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps1 = e0(j) + e0(i) - e0(a) - e0(b)
eps2 = e0(k) + e0(i) - e0(a) - e0(b)
App(p,1) = App(p,1) &
- (2d0*V(p,k,p,j) - V(p,k,j,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(a,b,k,i)/(eps1*eps2)
enddo
enddo
enddo
enddo
enddo
enddo
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
eps1 = e0(j) + e0(i) - e0(a) - e0(b)
eps2 = e0(j) + e0(i) - e0(a) - e0(c)
App(p,2) = App(p,2) &
+ (2d0*V(p,c,p,b) - V(p,c,b,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(j,i,c,a)/(eps1*eps2)
enddo
enddo
enddo
enddo
enddo
enddo
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
eps1 = e0(j) + e0(i) - e0(a) - e0(b)
eps2 = e0(j) - e0(c)
App(p,3) = App(p,3) &
+ (2d0*V(p,c,p,j) - V(p,c,j,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(a,b,c,i)/(eps1*eps2)
enddo
enddo
enddo
enddo
enddo
enddo
App(:,4) = App(:,3)
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps1 = e0(j) + e0(i) - e0(a) - e0(b)
eps2 = e0(k) - e0(b)
App(p,5) = App(p,5) &
- (2d0*V(p,b,p,k) - V(p,b,k,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(i,j,k,a)/(eps1*eps2)
enddo
enddo
enddo
enddo
enddo
enddo
App(:,6) = App(:,5)
! Frequency-independent part of the third-order self-energy
SigInf(:) = App(:,1) + App(:,2) + App(:,3) + App(:,4) + App(:,5) + App(:,6)
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
nSCF = 0
n_diis = 0
Conv = 1d0
Sig2(:) = 0d0
Sig3(:) = 0d0
e_diis(:,:) = 0d0
error_diis(:,:) = 0d0
eGF3(:) = e0(:)
eOld(:) = e0(:)
do while(Conv > thresh .and. nSCF < maxSCF)
! Frequency-dependent second-order contribution
Bpp(:,:) = 0d0
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = eGF3(p) + e0(a) - e0(i) - e0(j)
Bpp(p,1) = Bpp(p,1) &
+ (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)/eps
enddo
enddo
enddo
enddo
do p=nC+1,nBas-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = eGF3(p) + e0(i) - e0(a) - e0(b)
Bpp(p,2) = Bpp(p,2) &
+ (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)/eps
enddo
enddo
enddo
enddo
! Total second-order Green function
Sig2(:) = Bpp(:,1) + Bpp(:,2)
! Frequency-dependent third-order contribution: "C" terms
Cpp(:,:) = 0d0
do p=nC+1,nBas-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
eps1 = eGF3(p) + e0(i) - e0(a) - e0(b)
eps2 = eGF3(p) + e0(i) - e0(c) - e0(d)
Cpp(p,1) = Cpp(p,1) &
+ (2d0*V(p,i,a,b) - V(p,i,b,a))*V(a,b,c,d)*V(p,i,c,d)/(eps1*eps2)
enddo
enddo
enddo
enddo
enddo
enddo
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps1 = eGF3(p) + e0(i) - e0(a) - e0(b)
eps2 = e0(j) + e0(k) - e0(a) - e0(b)
Cpp(p,2) = Cpp(p,2) &
+ (2d0*V(p,i,a,b) - V(p,i,b,a))*V(a,b,j,k)*V(p,i,j,k)/(eps1*eps2)
enddo
enddo
enddo
enddo
enddo
enddo
Cpp(:,3) = Cpp(:,2)
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
eps1 = eGF3(p) + e0(a) - e0(i) - e0(j)
eps2 = e0(i) + e0(j) - e0(b) - e0(c)
Cpp(p,4) = Cpp(p,4) &
+ (2d0*V(p,a,i,j) - V(p,a,j,i))*V(i,j,b,c)*V(p,a,b,c)/(eps1*eps2)
enddo
enddo
enddo
enddo
enddo
enddo
Cpp(:,5) = Cpp(:,4)
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do l=nC+1,nO
do a=nO+1,nBas-nR
eps1 = eGF3(p) + e0(a) - e0(i) - e0(j)
eps2 = eGF3(p) + e0(a) - e0(k) - e0(l)
Cpp(p,6) = Cpp(p,6) &
- (2d0*V(p,a,k,l) - V(p,a,l,k))*V(k,l,i,j)*V(p,a,i,j)/(eps1*eps2)
enddo
enddo
enddo
enddo
enddo
enddo
! Frequency-dependent third-order contribution: "D" terms
Dpp(:,:) = 0d0
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
eps1 = eGF3(p) + e0(i) - e0(a) - e0(b)
eps2 = eGF3(p) + e0(j) - e0(b) - e0(c)
Dpp(p,1) = Dpp(p,1) &
+ V(p,i,a,b)*(V(a,j,i,c)*( V(p,j,c,b) - 2d0*V(p,j,b,c)) &
+ V(a,j,c,i)*( V(p,j,b,c) - 2d0*V(p,j,c,b)))/(eps1*eps2)
Dpp(p,1) = Dpp(p,1) &
+ V(p,i,b,a)*(V(a,j,i,c)*(4d0*V(p,j,b,c) - 2d0*V(p,j,c,b)) &
+ V(a,j,c,i)*( V(p,j,c,b) - 2d0*V(p,j,b,c)))/(eps1*eps2)
enddo
enddo
enddo
enddo
enddo
enddo
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
eps1 = eGF3(p) + e0(i) - e0(a) - e0(c)
eps2 = e0(i) + e0(j) - e0(a) - e0(b)
Dpp(p,2) = Dpp(p,2) &
+ V(p,i,c,a)*(V(a,b,i,j)*(4d0*V(p,b,c,j) - 2d0*V(p,b,j,c)) &
+ V(a,b,j,i)*( V(p,b,j,c) - 2d0*V(p,b,c,j)))/(eps1*eps2)
Dpp(p,2) = Dpp(p,2) &
+ V(p,i,a,c)*(V(a,b,i,j)*( V(p,b,j,c) - 2d0*V(p,b,c,j)) &
+ V(a,b,j,i)*( V(p,b,c,j) - 2d0*V(p,b,j,c)))/(eps1*eps2)
enddo
enddo
enddo
enddo
enddo
enddo
Dpp(:,3) = Dpp(:,2)
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps1 = eGF3(p) + e0(a) - e0(j) - e0(k)
eps2 = e0(i) + e0(j) - e0(a) - e0(b)
Dpp(p,4) = Dpp(p,4) &
+ V(p,a,k,j)*(V(j,i,a,b)*(4d0*V(p,i,k,b) - 2d0*V(p,i,b,k)) &
+ V(j,i,b,a)*( V(p,i,b,k) - 2d0*V(p,i,k,b)))/(eps1*eps2)
Dpp(p,4) = Dpp(p,4) &
+ V(p,a,j,k)*(V(j,i,a,b)*( V(p,i,b,k) - 2d0*V(p,i,k,b)) &
+ V(j,i,b,a)*( V(p,i,k,b) - 2d0*V(p,i,b,k)))/(eps1*eps2)
enddo
enddo
enddo
enddo
enddo
enddo
Dpp(:,5) = Dpp(:,4)
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps1 = eGF3(p) + e0(a) - e0(i) - e0(k)
eps2 = eGF3(p) + e0(b) - e0(j) - e0(k)
Dpp(p,6) = Dpp(p,6) &
- V(p,a,k,i)*(V(i,b,a,j)*(4d0*V(p,b,k,j) - 2d0*V(p,b,j,k)) &
+ V(i,b,j,a)*( V(p,b,j,k) - 2d0*V(p,b,k,j)))/(eps1*eps2)
Dpp(p,6) = Dpp(p,6) &
- V(p,a,i,k)*(V(i,b,a,j)*( V(p,b,j,k) - 2d0*V(p,b,k,j)) &
+ V(i,b,j,a)*( V(p,b,k,j) - 2d0*V(p,b,j,k)))/(eps1*eps2)
enddo
enddo
enddo
enddo
enddo
enddo
! Compute renormalization factor (if required)
Z(:) = 1d0
if(renormalization == 0) then
Sig3(:) = SigInf(:) &
+ Cpp(:,1) + Cpp(:,2) + Cpp(:,3) + Cpp(:,4) + Cpp(:,5) + Cpp(:,6) &
+ Dpp(:,1) + Dpp(:,2) + Dpp(:,3) + Dpp(:,4) + Dpp(:,5) + Dpp(:,6)
elseif(renormalization == 1) then
Sig3(:) = SigInf(:) &
+ Cpp(:,1) + Cpp(:,2) + Cpp(:,3) + Cpp(:,4) + Cpp(:,5) + Cpp(:,6) &
+ Dpp(:,1) + Dpp(:,2) + Dpp(:,3) + Dpp(:,4) + Dpp(:,5) + Dpp(:,6)
Z(:) = Cpp(:,2) + Cpp(:,3) + Cpp(:,4) + Cpp(:,5) &
+ Dpp(:,2) + Dpp(:,3) + Dpp(:,4) + Dpp(:,5)
Z(nC+1:nBas-nR) = Z(nC+1:nBas-nR)/Sig2(nC+1:nBas-nR)
Z(:) = 1d0/(1d0 - Z(:))
Sig3(:) = Z(:)*Sig3(:)
elseif(renormalization == 2) then
Sig2h1p(:) = Cpp(:,4) + Cpp(:,5) + Cpp(:,6) + Dpp(:,4) + Dpp(:,5) + Dpp(:,6)
Sig1h2p(:) = Cpp(:,1) + Cpp(:,2) + Cpp(:,3) + Dpp(:,1) + Dpp(:,2) + Dpp(:,3)
X2h1p(:) = Cpp(:,4) + Cpp(:,5) + Dpp(:,4) + Dpp(:,5)
X1h2p(:) = Cpp(:,2) + Cpp(:,3) + Dpp(:,2) + Dpp(:,3)
X2h1p(nC+1:nBas-nR) = X2h1p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,1)
X1h2p(nC+1:nBas-nR) = X1h2p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,2)
Sig3(:) = SigInf(:) + &
+ 1d0/(1d0 - X2h1p(:))*Sig2h1p(:) &
+ 1d0/(1d0 - X1h2p(:))*Sig1h2p(:)
elseif(renormalization == 3) then
Sig3(:) = SigInf(:) &
+ Cpp(:,1) + Cpp(:,2) + Cpp(:,3) + Cpp(:,4) + Cpp(:,5) + Cpp(:,6) &
+ Dpp(:,1) + Dpp(:,2) + Dpp(:,3) + Dpp(:,4) + Dpp(:,5) + Dpp(:,6)
Sig2h1p(:) = Cpp(:,4) + Cpp(:,5) + Cpp(:,6) + Dpp(:,4) + Dpp(:,5) + Dpp(:,6)
Sig1h2p(:) = Cpp(:,1) + Cpp(:,2) + Cpp(:,3) + Dpp(:,1) + Dpp(:,2) + Dpp(:,3)
X2h1p(:) = Cpp(:,4) + Cpp(:,5) + Dpp(:,4) + Dpp(:,5)
X1h2p(:) = Cpp(:,2) + Cpp(:,3) + Dpp(:,2) + Dpp(:,3)
X2h1p(nC+1:nBas-nR) = X2h1p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,1)
X1h2p(nC+1:nBas-nR) = X1h2p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,2)
Z(:) = X2h1p(:)*Sig2h1p(:) + X1h2p(:)*Sig1h2p(:)
Z(nC+1:nBas-nR) = Z(nC+1:nBas-nR)/(Sig3(nC+1:nBas-nR) - SigInf(nC+1:nBas-nR))
Z(:) = 1d0/(1d0 - Z(:))
Sig3(:) = Z(:)*Sig3(:)
endif
! Total third-order Green function
eGF3(:) = e0(:) + Sig2(:) + Sig3(:)
! Convergence criteria
Conv = maxval(abs(eGF3 - eOld))
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(nBas,nBas,n_diis,error_diis,e_diis,eGF3-eOld,eGF3)
! Store result for next iteration
eOld(:) = eGF3(:)
! Print results
call print_GF3(nBas,nO,nSCF,Conv,e0,Z,eGF3)
! Increment
nSCF = nSCF + 1
enddo
!------------------------------------------------------------------------
! End main SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF+1) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
endif
end subroutine GF3_diag

65
src/MCQC/Green_function.f90
View File

@ -1,65 +0,0 @@
subroutine Green_function(nBas,nO,nV,nWalk,nWP,cO,cV,eO_Quad,eV_Quad,AO, &
o1MO,o2MO,v1MO,v2MO,o11,o12,o21,o22,v11,v12,v21,v22)
! Calculate the Green functions
implicit none
include 'parameters.h'
include 'quadrature.h'
! Input variables
integer,intent(in) :: nBas,nO,nV,nWalk,nWP
double precision,intent(in) :: AO(nWalk,2,nBas),cO(nBas,nO),cV(nBas,nV), &
eO_Quad(nQuad,nO),eV_Quad(nQuad,nV)
! Local variables
integer :: kW,lW,klW,i,a,q
double precision :: o1MO(nWalk,nO),o2MO(nWalk,nO),v1MO(nWalk,nV),v2MO(nWalk,nV)
! Output variables
double precision,intent(out) :: o11(nQuad,nWP),o12(nQuad,nWP),o21(nQuad,nWP),o22(nQuad,nWP)
double precision,intent(out) :: v11(nQuad,nWP),v12(nQuad,nWP),v21(nQuad,nWP),v22(nQuad,nWP)
! Calculate occupied and virtual MOs
o1MO = matmul(AO(:,1,:),cO)
o2MO = matmul(AO(:,2,:),cO)
v1MO = matmul(AO(:,1,:),cV)
v2MO = matmul(AO(:,2,:),cV)
! Compute occupied Green functions
o11 = 0d0
o12 = 0d0
o21 = 0d0
o22 = 0d0
v11 = 0d0
v12 = 0d0
v21 = 0d0
v22 = 0d0
do q=1,nQuad
klW = 0
do kW=1,nWalk-1
do lW=kW+1,nWalk
klW = klW + 1
do i=1,nO
o11(q,klW) = o11(q,klW) + o1MO(kW,i)*o1MO(lW,i)*eO_Quad(q,i)
o12(q,klW) = o12(q,klW) + o1MO(kW,i)*o2MO(lW,i)*eO_Quad(q,i)
o21(q,klW) = o21(q,klW) + o2MO(kW,i)*o1MO(lW,i)*eO_Quad(q,i)
o22(q,klW) = o22(q,klW) + o2MO(kW,i)*o2MO(lW,i)*eO_Quad(q,i)
enddo
do a=1,nV
v11(q,klW) = v11(q,klW) + v1MO(kW,a)*v1MO(lW,a)*eV_Quad(q,a)
v12(q,klW) = v12(q,klW) + v1MO(kW,a)*v2MO(lW,a)*eV_Quad(q,a)
v21(q,klW) = v21(q,klW) + v2MO(kW,a)*v1MO(lW,a)*eV_Quad(q,a)
v22(q,klW) = v22(q,klW) + v2MO(kW,a)*v2MO(lW,a)*eV_Quad(q,a)
enddo
enddo
enddo
enddo
end subroutine Green_function

33
src/MCQC/Hartree_matrix_AO_basis.f90
View File

@ -1,33 +0,0 @@
subroutine Hartree_matrix_AO_basis(nBas,P,Hc,G,H)
! Compute Hartree matrix in the AO basis
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas),G(nBas,nBas,nBas,nBas)
! Local variables
integer :: mu,nu,la,si
! Output variables
double precision,intent(out) :: H(nBas,nBas)
H = Hc
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do si=1,nBas
H(mu,nu) = H(mu,nu) + P(la,si)*G(mu,la,nu,si)
enddo
enddo
enddo
enddo
end subroutine Hartree_matrix_AO_basis

26
src/MCQC/Hartree_matrix_MO_basis.f90
View File

@ -1,26 +0,0 @@
subroutine Hartree_matrix_MO_basis(nBas,c,P,Hc,G,H)
! Compute Hartree matrix in the MO basis
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: c(nBas,nBas),P(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas),G(nBas,nBas,nBas,nBas)
! Output variables
double precision,intent(out) :: H(nBas,nBas)
! Compute Hartree matrix in the AO basis
call Hartree_matrix_AO_basis(nBas,P,Hc,G,H)
! Transform Hartree matrix in the MO basis
H = matmul(transpose(c),matmul(H,c))
end subroutine Hartree_matrix_MO_basis

344
src/MCQC/MCMP2.f90
View File

@ -1,344 +0,0 @@
subroutine MCMP2(doDrift,nBas,nC,nO,nV,c,e,EcMP2, &
nMC,nEq,nWalk,dt,nPrint, &
nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
Norm, &
EcMCMP2,Err_EcMCMP2,Var_EcMCMP2)
! Perform Monte Carlo MP2 calculation
implicit none
include 'parameters.h'
include 'quadrature.h'
! Input variables
logical,intent(in) :: doDrift
integer,intent(in) :: nBas,nC,nO,nV,nMC,nEq,nWalk,nPrint
double precision,intent(inout):: dt
double precision,intent(in) :: EcMP2(3)
double precision,intent(in) :: c(nBas,nBas),e(nBas)
integer,intent(in) :: nShell
integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell)
double precision,intent(in) :: CenterShell(maxShell,3),DShell(maxShell,maxK),ExpShell(maxShell,maxK)
! Local variables
logical :: AcPh,EqPh,Accept,dump
double precision :: start_Eq,end_Eq,t_Eq,start_Ac,end_Ac,t_Ac
integer :: nWP
double precision :: Norm,NormSq,nData,tau
double precision,allocatable :: chi1(:,:,:),chi2(:,:,:),eta(:)
double precision,allocatable :: cO(:,:),cV(:,:),eO(:),eV(:),P(:,:),eO_Quad(:,:),eV_Quad(:,:)
double precision,allocatable :: r(:,:,:), r12(:), gAO(:,:,:), g(:,:), w(:)
double precision,allocatable :: rp(:,:,:),r12p(:),gAOp(:,:,:), gp(:,:),wp(:)
double precision,allocatable :: o1MO(:,:),o2MO(:,:),v1MO(:,:),v2MO(:,:)
double precision,allocatable :: o11(:,:),o12(:,:),o21(:,:),o22(:,:)
double precision,allocatable :: v11(:,:),v12(:,:),v21(:,:),v22(:,:)
double precision,allocatable :: fd_Quad(:,:),fx_Quad(:,:),fd(:),fx(:),fdx(:)
double precision,allocatable :: dgAO(:,:,:,:),dg(:,:,:),dgAOp(:,:,:,:),dgp(:,:,:)
double precision,allocatable :: F(:,:,:),Fp(:,:,:),T(:),Tp(:)
double precision :: acceptance,D
double precision :: eloc_MP2(3),mean_MP2(3),variance_MP2(3)
integer :: iW,kW,lW,klW,iMC,q
! Output variables
double precision,intent(out) :: EcMCMP2(3),Err_EcMCMP2(3),Var_EcMCMP2(3)
! Number of distinct walker pairs
nWP = nWalk*(nWalk-1)/2
! Diffusion coefficient
D = 0.5d0
! Do diffusion-drift moves?
if(doDrift) then
write(*,*)
write(*,*) '*** Diffusion-drift algorithm ***'
write(*,*)
else
write(*,*)
write(*,*) '*** Diffusion-only algorithm ***'
write(*,*)
endif
! Print results
dump = .true.
if(dump) open(unit=13,file='results/data')
!------------------------------------------------------------------------
! Memory allocation
!------------------------------------------------------------------------
allocate(cO(nBas,nO),cV(nBas,nV),eO(nO),eV(nV), &
eO_Quad(nQuad,nO),eV_Quad(nQuad,nV), &
P(nBas,nBas),r(nWalk,2,3),rp(nWalk,2,3), &
chi1(nWalk,2,3),chi2(nWalk,2,3),eta(nWalk), &
r12(nWalk),r12p(nWalk),w(nWalk),wp(nWalk), &
g(nWalk,2),gp(nWalk,2),gAO(nWalk,2,nBas),gAOp(nWalk,2,nBas), &
dg(nWalk,2,3),dgp(nWalk,2,3),dgAO(nWalk,2,3,nBas),dgAOp(nWalk,2,3,nBas), &
o1MO(nWalk,nO),v1MO(nWalk,nV),o2MO(nWalk,nO),v2MO(nWalk,nV), &
o11(nQuad,nWP),v11(nQuad,nWP),o12(nQuad,nWP),v12(nQuad,nWP), &
o21(nQuad,nWP),v21(nQuad,nWP),o22(nQuad,nWP),v22(nQuad,nWP), &
fd_Quad(nQuad,nWP),fd(nWP),fx_Quad(nQuad,nWP),fx(nWP),fdx(nWP), &
T(nWalk),Tp(nWalk),F(nWalk,2,3),Fp(nWalk,2,3))
! Split MOs into occupied and virtual sets
eO(1:nO) = e(nC+1:nC+nO)
eV(1:nV) = e(nC+nO+1:nBas)
do q=1,nQuad
tau = 1d0/rQuad(q)
eO_Quad(q,1:nO) = exp(+eO(1:nO)*(tau-1d0))*sqrt(tau)
eV_Quad(q,1:nV) = exp(-eV(1:nV)*(tau-1d0))*sqrt(tau)
enddo
cO(1:nBas,1:nO) = c(1:nBas,nC+1:nC+nO)
cV(1:nBas,1:nV) = c(1:nBas,nC+nO+1:nBas)
! Compute norm of the trial wave function
call norm_trial(nBas,nO,cO,P,Norm,NormSq)
!------------------------------------------------------------------------
! Initialize MC-MP2 calculation
!------------------------------------------------------------------------
! Initialize electron coordinates
call random_number(r)
r = 2d0*r - 1d0
! Compute initial interelectronic distances
call rij(nWalk,r,r12)
! Compute initial AO values and their derivatives (if required)
call AO_values(doDrift,nBas,nShell,nWalk,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,r,gAO,dgAO)
! Compute initial weight function
call density(doDrift,nBas,nWalk,P,gAO,dgAO,g,dg)
! Compute initial weights
w(1:nWalk) = g(1:nWalk,1)*g(1:nWalk,2)/r12(1:nWalk)
! Compute initial quantum force
if(doDrift) call drift(nWalk,r,r12,g,dg,F)
! Equilibration or Accumulation?
AcPh = .false.
EqPh = .true.
! Initialization
nData = 0d0
acceptance = 0d0
mean_MP2 = 0d0
variance_MP2 = 0d0
T = 1d0
Tp = 1d0
!------------------------------------------------------------------------
! Start main Monte Carlo loop
!------------------------------------------------------------------------
call cpu_time(start_Eq)
do iMC=1,nEq+nMC
! Timings
if(iMC == nEq + 1) then
AcPh = .true.
EqPh = .false.
write(*,*) 'Time step value at the end of equilibration: dt = ',dt
write(*,*)
call cpu_time(end_Eq)
t_Eq = end_Eq - start_Eq
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for equilibration = ',t_Eq,' seconds'
write(*,*)
call cpu_time(start_Ac)
endif
! Optimize time step to reach 50% acceptance
if(EqPh .and. mod(iMC,100) == 0) call optimize_timestep(nWalk,iMC,acceptance,dt)
! Move electrons
call random_number(chi1)
call random_number(chi2)
! Diffusion
rp(:,:,:) = r(:,:,:) + sqrt(2d0*D*dt)*sqrt(-2d0*log(chi1(:,:,:)))*cos(2d0*pi*chi2(:,:,:))
! Drift
if(doDrift) rp(:,:,:) = rp(:,:,:) + D*dt*F(:,:,:)
! Compute new interelectronic distances
call rij(nWalk,rp,r12p)
! Compute new AO values and their derivatives (if required)
call AO_values(doDrift,nBas,nShell,nWalk,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,rp,gAOp,dgAOp)
call Density(doDrift,nBas,nWalk,P,gAOp,dgAOp,gp,dgp)
! Compute new weights
wp(1:nWalk) = gp(1:nWalk,1)*gp(1:nWalk,2)/r12p(1:nWalk)
! Compute new quantum force and transition probability
if(doDrift) then
call Drift(nWalk,rp,r12p,gp,dgp,Fp)
call transition_probability(nWalk,dt,D,r,rp,F,Fp,T,Tp)
endif
! Move for walkers
call random_number(eta)
do iW=1,nWalk
Accept = (wp(iW)*Tp(iW))/(w(iW)*T(iW)) > eta(iW)
if(Accept) then
acceptance = acceptance + 1d0
r(iW,1:2,1:3) = rp(iW,1:2,1:3)
gAO(iW,1:2,1:nBas) = gAOp(iW,1:2,1:nBas)
r12(iW) = r12p(iW)
w(iW) = wp(iW)
if(doDrift) F(iW,1:2,1:3) = Fp(iW,1:2,1:3)
endif
enddo
! Accumulation phase
if(AcPh) then
nData = nData + 1d0
! Calculate Green functions
call Green_function(nBas,nO,nV,nWalk,nWP,cO,cV,eO_Quad,eV_Quad,gAO, &
o1MO,o2MO,v1MO,v2MO,o11,o12,o21,o22,v11,v12,v21,v22)
! Compute local energy
fd_Quad = o11*o22*v11*v22 + o12*o21*v12*v21
fx_Quad = o11*o22*v12*v21 + o12*o21*v11*v22
fd = matmul(wQuad,fd_Quad)
fx = matmul(wQuad,fx_Quad)
eloc_MP2 = 0d0
klW = 0
do kW=1,nWalk-1
do lW=kW+1,nWalk
klW = klW + 1
eloc_MP2(2) = eloc_MP2(2) + fd(klW)/(r12(kW)*r12(lW)*w(kW)*w(lW))
eloc_MP2(3) = eloc_MP2(3) + fx(klW)/(r12(kW)*r12(lW)*w(kW)*w(lW))
enddo
enddo
eloc_MP2(2) = -2d0*eloc_MP2(2)/dble(2*nWP)
eloc_MP2(3) = eloc_MP2(3)/dble(2*nWP)
fdx = -2d0*fd + fx
eloc_MP2(1) = eloc_MP2(2) + eloc_MP2(3)
! Accumulate results
mean_MP2 = mean_MP2 + eloc_MP2
variance_MP2 = variance_MP2 + eloc_MP2*eloc_MP2
! Print results
if(mod(iMC,nPrint) == 0) then
call compute_error(nData,mean_MP2,variance_MP2,Err_EcMCMP2)
EcMCMP2 = mean_MP2/nData
Var_EcMCMP2 = variance_MP2/nData
EcMCMP2 = Norm*EcMCMP2
Var_EcMCMP2 = Norm*Var_EcMCMP2
Err_EcMCMP2 = Norm*Err_EcMCMP2
write(*,*)
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,I15)') 'Number of data points ','|',int(nData)
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,10I15)') 'acceptance ','|',int(100*acceptance/dble(nWalk*iMC))
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,10F15.8)') 'MP2 correlation energy Total ','|',EcMCMP2(1)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Direct ','|',EcMCMP2(2)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Exchange ','|',EcMCMP2(3)
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,10F15.8)') 'Statistical error Total ','|',Err_EcMCMP2(1)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Direct ','|',Err_EcMCMP2(2)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Exchange ','|',Err_EcMCMP2(3)
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,10F15.8)') 'Variance Total ','|',Var_EcMCMP2(1)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Direct ','|',Var_EcMCMP2(2)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Exchange ','|',Var_EcMCMP2(3)
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,10F15.8)') 'Dev. wrt deterministic Total ','|',EcMCMP2(1) - EcMP2(1)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Direct ','|',EcMCMP2(2) - EcMP2(2)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Exchange ','|',EcMCMP2(3) - EcMP2(3)
write(*,*)'-------------------------------------------------------'
if(dump) write(13,*) int(nData),EcMCMP2(1),Err_EcMCMP2(1)
endif
endif
!------------------------------------------------------------------------
! End main Monte Carlo loop
!------------------------------------------------------------------------
enddo
! Timing
call cpu_time(end_Ac)
t_Ac = end_Ac - start_Ac
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for accumulation = ',t_Ac,' seconds'
write(*,*)
! Close files
if(dump) close(unit=13)
end subroutine MCMP2

446
src/MCQC/MCMP2.f90.x
View File

@ -1,446 +0,0 @@
subroutine MCMP2(varmin,doDrift,nBas,nEl,nC,nO,nV,c,e,EcMP2, &
nMC,nEq,nWalk,dt,nPrint, &
nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
TrialType,Norm,cTrial,gradient,hessian, &
EcMCMP2,Err_EcMCMP2,Var_EcMCMP2)
! Perform Monte Carlo MP2 calculation
implicit none
include 'parameters.h'
include 'quadrature.h'
! Input variables
logical,intent(in) :: varmin,doDrift
integer,intent(in) :: nBas,nEl,nC,nO,nV,nMC,nEq,nWalk,nPrint
double precision,intent(inout):: dt
double precision,intent(in) :: EcMP2(3)
double precision,intent(in) :: c(nBas,nBas),e(nBas)
integer,intent(in) :: nShell
integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell)
double precision,intent(in) :: CenterShell(maxShell,3),DShell(maxShell,maxK),ExpShell(maxShell,maxK)
! Local variables
logical :: AcPh,EqPh,Accept,dump
double precision :: start_Eq,end_Eq,t_Eq,start_Ac,end_Ac,t_Ac
integer :: nWP
double precision :: Norm,NormSq,nData,tau
double precision,allocatable :: chi1(:,:,:),chi2(:,:,:),eta(:)
double precision,allocatable :: cO(:,:),cV(:,:),eO(:),eV(:),P(:,:),eO_Quad(:,:),eV_Quad(:,:)
double precision,allocatable :: r(:,:,:), r12(:), gAO(:,:,:), g(:,:), w(:)
double precision,allocatable :: rp(:,:,:),r12p(:),gAOp(:,:,:), gp(:,:),wp(:)
double precision,allocatable :: o1MO(:,:),o2MO(:,:),v1MO(:,:),v2MO(:,:)
double precision,allocatable :: o11(:,:),o12(:,:),o21(:,:),o22(:,:)
double precision,allocatable :: v11(:,:),v12(:,:),v21(:,:),v22(:,:)
double precision,allocatable :: fd_Quad(:,:),fx_Quad(:,:),fd(:),fx(:),fdx(:)
double precision,allocatable :: dgAO(:,:,:,:),dg(:,:,:),dgAOp(:,:,:,:),dgp(:,:,:)
double precision,allocatable :: F(:,:,:),Fp(:,:,:),T(:),Tp(:)
double precision :: acceptance,D
double precision :: eloc_MP2(3),mean_MP2(3),variance_MP2(3)
integer :: iW,kW,lW,klW,iMC,i,a,q,iTr,jTr
double precision :: el,del
double precision,allocatable :: psiTr(:,:),dw(:,:),deloc(:),edeloc(:),mean_de(:),mean_ede(:),mean_dede(:,:)
double precision,allocatable :: dwdw(:),mean_dw(:)
! Output variables
double precision,intent(out) :: EcMCMP2(3),Err_EcMCMP2(3),Var_EcMCMP2(3)
integer,intent(in) :: TrialType
double precision,intent(inout):: cTrial(nBas),gradient(nBas),hessian(nBas,nBas)
! Number of distinct walker pairs
nWP = nWalk*(nWalk-1)/2
! Diffusion coefficient
D = 0.5d0
! Do diffusion-drift moves?
write(*,*)
if(doDrift) then
write(*,*) '*** Diffusion-drift algorithm ***'
else
write(*,*) '*** Diffusion-only algorithm ***'
endif
write(*,*)
! Print results
dump = .true.
if(dump) open(unit=13,file='results/data')
! Variance minimization
if(varmin) then
open(unit=14,file='results/varmin')
endif
!------------------------------------------------------------------------
! Memory allocation
!------------------------------------------------------------------------
allocate(cO(nBas,nO),cV(nBas,nV),eO(nO),eV(nV), &
eO_Quad(nQuad,nO),eV_Quad(nQuad,nV), &
P(nBas,nBas),r(nWalk,2,3),rp(nWalk,2,3), &
chi1(nWalk,2,3),chi2(nWalk,2,3),eta(nWalk), &
r12(nWalk),r12p(nWalk),w(nWalk),wp(nWalk), &
g(nWalk,2),gp(nWalk,2),gAO(nWalk,2,nBas),gAOp(nWalk,2,nBas), &
dg(nWalk,2,3),dgp(nWalk,2,3),dgAO(nWalk,2,3,nBas),dgAOp(nWalk,2,3,nBas), &
o1MO(nWalk,nO),v1MO(nWalk,nV),o2MO(nWalk,nO),v2MO(nWalk,nV), &
o11(nQuad,nWP),v11(nQuad,nWP),o12(nQuad,nWP),v12(nQuad,nWP), &
o21(nQuad,nWP),v21(nQuad,nWP),o22(nQuad,nWP),v22(nQuad,nWP), &
fd_Quad(nQuad,nWP),fd(nWP),fx_Quad(nQuad,nWP),fx(nWP),fdx(nWP), &
T(nWalk),Tp(nWalk),F(nWalk,2,3),Fp(nWalk,2,3))
allocate(psiTr(nWalk,2),dw(nWalk,nBas),deloc(nBas),edeloc(nBas), &
mean_de(nBas),mean_ede(nBas),mean_dede(nBas,nBas))
allocate(dwdw(nBas),mean_dw(nBas))
! Split MOs into occupied and virtual sets
eO(1:nO) = e(nC+1:nC+nO)
eV(1:nV) = e(nC+nO+1:nBas)
do q=1,nQuad
tau = 1d0/rQuad(q)
do i=1,nO
eO_Quad(q,i) = exp(+eO(i)*(tau-1d0))*sqrt(tau)
enddo
do a=1,nV
eV_Quad(q,a) = exp(-eV(a)*(tau-1d0))*sqrt(tau)
enddo
enddo
cO(1:nBas,1:nO) = c(1:nBas,nC+1:nC+nO)
cV(1:nBas,1:nV) = c(1:nBas,nC+nO+1:nBas)
! Compute norm of the trial wave function
if(TrialType == 0) then
call NormTrial(TrialType,nBas,nO,cO,P,Norm,NormSq)
elseif(TrialType == 1) then
call NormTrial(TrialType,nBas,1,cTrial,P,Norm,NormSq)
endif
!------------------------------------------------------------------------
! Initialize MC-MP2 calculation
!------------------------------------------------------------------------
! Initialize electron coordinates
call random_number(r)
r = 2d0*r - 1d0
! Compute initial interelectronic distances
call rij(nWalk,r,r12)
! Compute initial AO values and their derivatives (if required)
call AO_values(doDrift,nBas,nShell,nWalk,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,r,gAO,dgAO)
! Compute initial weight function
call Density(doDrift,nBas,nWalk,P,gAO,dgAO,g,dg)
! Compute initial weights
w(1:nWalk) = g(1:nWalk,1)*g(1:nWalk,2)/r12(1:nWalk)
! Compute initial quantum force
if(doDrift) call Drift(nWalk,r,r12,g,dg,F)
! Equilibration or Accumulation?
AcPh = .false.
EqPh = .true.
! Initialization
nData = 0d0
acceptance = 0d0
mean_MP2 = 0d0
variance_MP2 = 0d0
if(varmin) then
mean_de = 0d0
mean_ede = 0d0
mean_dede = 0d0
! mean_dw = 0d0
endif
T = 1d0
Tp = 1d0
!------------------------------------------------------------------------
! Start main Monte Carlo loop
!------------------------------------------------------------------------
call cpu_time(start_Eq)
do iMC=1,nEq+nMC
! Timings
if(iMC == nEq + 1) then
AcPh = .true.
EqPh = .false.
write(*,*) 'Time step value at the end of equilibration: dt = ',dt
write(*,*)
call cpu_time(end_Eq)
t_Eq = end_Eq - start_Eq
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for equilibration = ',t_Eq,' seconds'
write(*,*)
call cpu_time(start_Ac)
endif
! Optimize time step to reach 50% acceptance
if(EqPh .and. mod(iMC,100) == 0) call optimize_timestep(nWalk,iMC,acceptance,dt)
! Move electrons
call random_number(chi1)
call random_number(chi2)
! Diffusion
rp = r + sqrt(2d0*D*dt)*sqrt(-2d0*log(chi1))*cos(2d0*pi*chi2)
! Drift
if(doDrift) rp = rp + D*dt*F
! Compute new interelectronic distances
call rij(nWalk,rp,r12p)
! Compute new AO values and their derivatives (if required)
call AO_values(doDrift,nBas,nShell,nWalk,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,rp,gAOp,dgAOp)
call Density(doDrift,nBas,nWalk,P,gAOp,dgAOp,gp,dgp)
! Compute new weights
wp(1:nWalk) = gp(1:nWalk,1)*gp(1:nWalk,2)/r12p(1:nWalk)
! Compute new quantum force and transition probability
if(doDrift) then
call Drift(nWalk,rp,r12p,gp,dgp,Fp)
call transition_probability(nWalk,dt,D,r,rp,F,Fp,T,Tp)
endif
! Move for walkers
call random_number(eta)
do iW=1,nWalk
Accept = (wp(iW)*Tp(iW))/(w(iW)*T(iW)) > eta(iW)
if(Accept) then
acceptance = acceptance + 1d0
r(iW,1:2,1:3) = rp(iW,1:2,1:3)
gAO(iW,1:2,1:nBas) = gAOp(iW,1:2,1:nBas)
r12(iW) = r12p(iW)
w(iW) = wp(iW)
if(doDrift) F(iW,1:2,1:3) = Fp(iW,1:2,1:3)
endif
enddo
! Accumulation phase
if(AcPh) then
nData = nData + 1d0
! Calculate Green functions
call Green_function(nBas,nO,nV,nWalk,nWP,cO,cV,eO_Quad,eV_Quad,gAO, &
o1MO,o2MO,v1MO,v2MO,o11,o12,o21,o22,v11,v12,v21,v22)
! Compute local energy
fd_Quad = o11*o22*v11*v22 + o12*o21*v12*v21
fx_Quad = o11*o22*v12*v21 + o12*o21*v11*v22
fd = matmul(wQuad,fd_Quad)
fx = matmul(wQuad,fx_Quad)
eloc_MP2 = 0d0
klW = 0
do kW=1,nWalk-1
do lW=kW+1,nWalk
klW = klW + 1
eloc_MP2(2) = eloc_MP2(2) + fd(klW)/(r12(kW)*r12(lW)*w(kW)*w(lW))
eloc_MP2(3) = eloc_MP2(3) + fx(klW)/(r12(kW)*r12(lW)*w(kW)*w(lW))
enddo
enddo
eloc_MP2(2) = -2d0*eloc_MP2(2)/dble(2*nWP)
eloc_MP2(3) = eloc_MP2(3)/dble(2*nWP)
fdx = -2d0*fd + fx
eloc_MP2(1) = eloc_MP2(2) + eloc_MP2(3)
! Accumulate results
mean_MP2 = mean_MP2 + eloc_MP2
variance_MP2 = variance_MP2 + eloc_MP2*eloc_MP2
! Accumulation for variane minimization
if(varmin) then
psiTr = 0d0
do iTr=1,nBas
psiTr(:,:) = psiTr(:,:) + cTrial(iTr)*gAO(:,:,iTr)
enddo
do iW=1,nWalk
do iTr=1,nBas
dw(iW,iTr) = gAO(iW,1,iTr)/psiTr(iW,1) + gAO(iW,2,iTr)/psiTr(iW,2)
enddo
enddo
deloc = 0d0
edeloc = 0d0
dwdw = 0d0
do iTr=1,nBas
klW = 0
do kW=1,nWalk-1
do lW=kW+1,nWalk
klW = klW + 1
el = fdx(klW)/(r12(kW)*r12(lW)*w(kW)*w(lW))
del = dw(kW,iTr) + dw(lW,iTr)
deloc(iTr) = deloc(iTr) + del*el
edeloc(iTr) = edeloc(iTr) + del*el*el
dwdw(iTr) = dwdw(iTr) + del
enddo
enddo
enddo
deloc = -2d0*deloc/dble(2*nWP)
edeloc = -2d0*edeloc/dble(2*nWP)
dwdw = 2d0*dwdw/dble(2*nWP)
mean_de(:) = mean_de(:) + deloc(:)
mean_ede(:) = mean_ede(:) + edeloc(:)
mean_dw(:) = mean_dw(:) + dwdw(:)
do iTr=1,nBas
do jTr=1,nBas
mean_dede(iTr,jTr) = mean_dede(iTr,jTr) + deloc(iTr)*deloc(jTr)
enddo
enddo
endif
! Print results
if(mod(iMC,nPrint) == 0) then
ecMCMP2 = mean_MP2/nData
Var_EcMCMP2 = variance_MP2/nData
call CalcError(nData,EcMCMP2,Var_EcMCMP2,Err_EcMCMP2)
EcMCMP2 = Norm*EcMCMP2
Var_EcMCMP2 = Norm*Var_EcMCMP2
Err_EcMCMP2 = Norm*Err_EcMCMP2
write(*,*)
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,I15)') 'Number of data points ','|',int(nData)
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,10I15)') 'acceptance ','|',int(100*acceptance/dble(nWalk*iMC))
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,10F15.8)') 'MP2 correlation energy Total ','|',EcMCMP2(1)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Direct ','|',EcMCMP2(2)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Exchange ','|',EcMCMP2(3)
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,10F15.8)') 'Statistical error Total ','|',Err_EcMCMP2(1)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Direct ','|',Err_EcMCMP2(2)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Exchange ','|',Err_EcMCMP2(3)
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,10F15.8)') 'Variance Total ','|',Var_EcMCMP2(1)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Direct ','|',Var_EcMCMP2(2)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Exchange ','|',Var_EcMCMP2(3)
write(*,*)'-------------------------------------------------------'
write(*,'(1X,A36,1X,A1,1X,10F15.8)') 'Dev. wrt deterministic Total ','|',EcMCMP2(1) - EcMP2(1)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Direct ','|',EcMCMP2(2) - EcMP2(2)
write(*,'(1X,A36,1X,A1,1X,10F15.8)') ' Exchange ','|',EcMCMP2(3) - EcMP2(3)
write(*,*)'-------------------------------------------------------'
if(dump) write(13,*) int(nData),EcMCMP2(1),Err_EcMCMP2(1)
! Compute gradient and hessian for variance minimization
if(varmin) then
gradient = 2d0*(mean_ede - mean_MP2(1)*mean_de/nData)
do iTr=1,nBas
do jTr=1,nBas
hessian(iTr,jTr) = 2d0*(mean_dede(iTr,jTr) - mean_de(iTr)*mean_de(jTr)/nData)
enddo
enddo
gradient = gradient/nData
hessian = hessian/nData
print*,'gradient'
call matout(nBas,1,gradient)
print*,'hessian'
call matout(nBas,nBas,hessian)
endif
endif
endif
!------------------------------------------------------------------------
! End main Monte Carlo loop
!------------------------------------------------------------------------
enddo
! Timing
call cpu_time(end_Ac)
t_Ac = end_Ac - start_Ac
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for accumulation = ',t_Ac,' seconds'
write(*,*)
! Close files
if(dump) close(unit=13)
if(varmin) close(unit=14)
end subroutine MCMP2

190
src/MCQC/MOM.f90
View File

@ -1,190 +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,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(:,:) + 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(nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
! 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

51
src/MCQC/MOM_overlap.f90
View File

@ -1,51 +0,0 @@
subroutine MOM_overlap(nBas,nO,S,cG,c,ON)
! Compute overlap between old and new MO coefficients
implicit none
! Input variables
integer,intent(in) :: nBas,nO
double precision,intent(in) :: S(nBas,nBas),cG(nBas,nBas),c(nBas,nBas)
! Local variables
integer :: i,j,ploc
double precision,allocatable :: Ov(:,:),pOv(:)
! Output variables
double precision,intent(inout):: ON(nBas)
allocate(Ov(nBas,nBas),pOv(nBas))
Ov = matmul(transpose(cG),matmul(S,c))
pOv(:) = 0d0
do i=1,nBas
do j=1,nBas
pOv(j) = pOv(j) + ON(i)*Ov(i,j)**2
enddo
enddo
pOv(:) = sqrt(pOV(:))
! print*,'--- MOM overlap ---'
! call matout(nBas,1,pOv)
ON(:) = 0d0
do i=1,nO
ploc = maxloc(pOv,nBas)
ON(ploc) = 1d0
pOv(ploc) = 0d0
enddo
! print*,'--- Occupation numbers ---'
! call matout(nBas,1,ON)
end subroutine MOM_overlap

27
src/MCQC/MOtoAO_transform.f90
View File

@ -1,27 +0,0 @@
subroutine MOtoAO_transform(nBas,S,c,A)
! Perform MO to AO transformation of a matrix A for a given metric S
! and coefficients c
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: S(nBas,nBas),c(nBas,nBas)
! Local variables
double precision,allocatable :: Sc(:,:)
! Output variables
double precision,intent(inout):: A(nBas,nBas)
! Memory allocation
allocate(Sc(nBas,nBas))
Sc = matmul(S,c)
A = matmul(Sc,matmul(A,transpose(Sc)))
end subroutine MOtoAO_transform

71
src/MCQC/MP2.f90
View File

@ -1,71 +0,0 @@
subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
! Perform third-order Moller-Plesset calculation
implicit none
! Input variables
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: ENuc,EHF
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),e(nBas)
! Local variables
integer :: i,j,a,b
double precision :: eps,E2a,E2b
! Output variables
double precision,intent(out) :: EcMP2(3)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Moller-Plesset second-order calculation |'
write(*,*)'************************************************'
write(*,*)
! Compute MP2 energy
E2a = 0d0
E2b = 0d0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = e(i) + e(j) - e(a) - e(b)
! Second-order ring diagram
E2a = E2a + ERI(i,j,a,b)*ERI(i,j,a,b)/eps
! Second-order exchange diagram
E2b = E2b + ERI(i,j,a,b)*ERI(i,j,b,a)/eps
enddo
enddo
enddo
enddo
EcMP2(2) = 2d0*E2a
EcMP2(3) = -E2b
EcMP2(1) = EcMP2(2) + EcMP2(3)
write(*,*)
write(*,'(A32)') '-----------------------'
write(*,'(A32)') ' MP2 calculation '
write(*,'(A32)') '-----------------------'
write(*,'(A32,1X,F16.10)') ' MP2 correlation energy',EcMP2(1)
write(*,'(A32,1X,F16.10)') ' Direct part ',EcMP2(2)
write(*,'(A32,1X,F16.10)') ' Exchange part ',EcMP2(3)
write(*,'(A32)') '-----------------------'
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy',EHF + EcMP2(1)
write(*,'(A32,1X,F16.10)') ' MP2 total energy',ENuc + EHF + EcMP2(1)
write(*,'(A32)') '-----------------------'
write(*,*)
end subroutine MP2

167
src/MCQC/MP2F12.f90
View File

@ -1,167 +0,0 @@
subroutine MP2F12(nBas,nC,nO,nV,ERI,F12,Yuk,FC,EHF,e,c)
! Perform MP2-F12 calculation
implicit none
! Input variables
integer,intent(in) :: nBas,nC,nO,nV
double precision,intent(in) :: EHF
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: F12(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Yuk(nBas,nBas,nBas,nBas)
double precision,intent(in) :: FC(nBas,nBas,nBas,nBas,nBas,nBas)
! Local variables
double precision,allocatable :: ooCoo(:,:,:,:)
double precision,allocatable :: ooFoo(:,:,:,:)
double precision,allocatable :: ooYoo(:,:,:,:)
double precision,allocatable :: ooCvv(:,:,:,:)
double precision,allocatable :: ooFvv(:,:,:,:)
double precision,allocatable :: oooFCooo(:,:,:,:,:,:)
double precision,allocatable :: eO(:),eV(:)
double precision,allocatable :: cO(:,:),cV(:,:)
double precision :: E2a,E2b,E3a,E3b,E4a,E4b,E4c,E4d
integer :: i,j,k,l,a,b
double precision :: EcMP2F12(4)
double precision :: EcMP2a
double precision :: EcMP2b
double precision :: EcMP2
double precision :: eps
! Split MOs into occupied and virtual sets
allocate(eO(nO),eV(nV))
eO(1:nO) = e(nC+1:nC+nO)
eV(1:nV) = e(nC+nO+1:nBas)
allocate(cO(nBas,nO),cV(nBas,nV))
cO(1:nBas,1:nO) = c(1:nBas,nC+1:nC+nO)
cV(1:nBas,1:nV) = c(1:nBas,nC+nO+1:nBas)
! Compute conventional MP2 energy
allocate(ooCvv(nO,nO,nV,nV))
call AOtoMO_oovv(nBas,nO,nV,cO,cV,ERI,ooCvv)
EcMP2a = 0d0
EcMP2b = 0d0
do i=1,nO
do j=1,nO
do a=1,nV
do b=1,nV
eps = eO(i) + eO(j) - eV(a) - eV(b)
EcMP2a = EcMP2a + ooCoo(i,j,a,b)/eps
EcMP2b = EcMP2b + ooCoo(i,j,b,a)/eps
enddo
enddo
enddo
enddo
EcMP2 = EcMP2a + EcMP2b
! Compute the two-electron part of the MP2-F12 energy
allocate(ooYoo(nO,nO,nO,nO))
call AOtoMO_oooo(nBas,nO,cO,Yuk,ooYoo)
E2a = 0d0
E2b = 0d0
do i=1,nO
do j=1,nO
E2a = E2a + ooYoo(i,j,i,j)
E2b = E2b + ooYoo(i,j,j,i)
enddo
enddo
deallocate(ooYoo)
! Compute the three-electron part of the MP2-F12 energy
allocate(oooFCooo(nO,nO,nO,nO,nO,nO))
call AOtoMO_oooooo(nBas,nO,cO,FC,oooFCooo)
E3a = 0d0
E3b = 0d0
do i=1,nO
do j=1,nO
do k=1,nO
E3a = E3a + oooFCooo(i,j,k,k,j,i)
E3b = E3b + oooFCooo(i,j,k,k,i,j)
enddo
enddo
enddo
deallocate(oooFCooo)
! Compute the four-electron part of the MP2-F12 energy
allocate(ooCoo(nO,nO,nO,nO),ooFoo(nO,nO,nO,nO))
call AOtoMO_oooo(nBas,nO,cO,ERI,ooCoo)
call AOtoMO_oooo(nBas,nO,cO,F12,ooFoo)
E4a = 0d0
E4b = 0d0
do i=1,nO
do j=1,nO
do k=1,nO
do l=1,nO
E4a = E4a + ooCoo(i,j,k,l)*ooFoo(i,j,k,l)
E4b = E4b + ooCoo(i,j,k,l)*ooFoo(j,i,k,l)
enddo
enddo
enddo
enddo
deallocate(ooCoo,ooFoo)
allocate(ooCvv(nO,nO,nV,nV),ooFvv(nO,nO,nV,nV))
call AOtoMO_oovv(nBas,nO,nV,cO,cV,ERI,ooCvv)
call AOtoMO_oovv(nBas,nO,nV,cO,cV,F12,ooFvv)
E4c = 0d0
E4d = 0d0
do i=1,nO
do j=1,nO
do a=1,nV
do b=1,nV
E4c = E4c + ooCvv(i,j,a,b)*ooFvv(i,j,a,b)
E4d = E4d + ooCvv(i,j,a,b)*ooFvv(j,i,a,b)
enddo
enddo
enddo
enddo
deallocate(ooCvv,ooFvv)
! Final scaling of the various components
EcMP2F12(1) = +0.625d0*E2a - 0.125d0*E2b
EcMP2F12(2) = -1.250d0*E3a + 0.250d0*E3b
EcMP2F12(3) = +0.625d0*E4a - 0.125d0*E4b - 0.625d0*E4c + 0.125d0*E4d
write(*,*)
write(*,'(A32)') '-----------------------'
write(*,'(A32)') ' MP2-F12 calculation '
write(*,'(A32)') '-----------------------'
write(*,'(A32,1X,F16.10)') ' MP2 ',+EcMP2
write(*,'(A32,1X,F16.10)') ' MP2-F12 E(2) ',-EcMP2F12(1)
write(*,'(A32,1X,F16.10)') ' MP2-F12 E(3) ',-EcMP2F12(2)
write(*,'(A32,1X,F16.10)') ' MP2-F12 E(4) ',-EcMP2F12(3)
write(*,'(A32)') '-----------------------'
write(*,'(A32,1X,F16.10)') ' Total ',EcMP2-EcMP2F12(1)-EcMP2F12(2)-EcMP2F12(3)
write(*,'(A32)') '-----------------------'
write(*,*)
deallocate(cO,cV)
end subroutine MP2F12

187
src/MCQC/MP3.f90
View File

@ -1,187 +0,0 @@
subroutine MP3(nBas,nEl,ERI,e,ENuc,EHF)
! Perform third-order Moller-Plesset calculation
implicit none
! Input variables
integer,intent(in) :: nBas,nEl
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 :: nBas2,nO,nV
integer :: i,j,k,l,a,b,c,d
double precision,allocatable :: se(:)
double precision,allocatable :: eO(:)
double precision,allocatable :: eV(:)
double precision,allocatable :: sERI(:,:,:,:)
double precision,allocatable :: dbERI(:,:,:,:)
double precision,allocatable :: OOOO(:,:,:,:)
double precision,allocatable :: OOVV(:,:,:,:)
double precision,allocatable :: OVVO(:,:,:,:)
double precision,allocatable :: VVOO(:,:,:,:)
double precision,allocatable :: VVVV(:,:,:,:)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Moller-Plesset third-order calculation |'
write(*,*)'************************************************'
write(*,*)
! Spatial to spin orbitals
nBas2 = 2*nBas
allocate(se(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
call spatial_to_spin_MO_energy(nBas,e,nBas2,se)
call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI)
! Antysymmetrize ERIs
allocate(dbERI(nBas2,nBas2,nBas2,nBas2))
call antisymmetrize_ERI(2,nBas2,sERI,dbERI)
deallocate(sERI)
! Define occupied and virtual spaces
nO = nEl
nV = nBas2 - nO
! Form energy denominator
allocate(eO(nO),eV(nV))
eO(:) = se(1:nO)
eV(:) = se(nO+1:nBas2)
deallocate(se)
! 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:nBas2,nO+1:nBas2)
OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas2,nO+1:nBas2, 1:nO )
VVOO(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2, 1:nO , 1:nO )
VVVV(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2,nO+1:nBas2,nO+1:nBas2)
deallocate(dbERI)
! Compute MP2 energy
E2 = 0d0
do i=1,nO
do j=1,nO
do a=1,nV
do b=1,nV
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=1,nO
do j=1,nO
do k=1,nO
do l=1,nO
do a=1,nV
do b=1,nV
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=1,nO
do j=1,nO
do a=1,nV
do b=1,nV
do c=1,nV
do d=1,nV
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=1,nO
do j=1,nO
do k=1,nO
do a=1,nV
do b=1,nV
do c=1,nV
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 MP3

36
src/MCQC/Makefile
View File

@ -1,36 +0,0 @@
IDIR =../../include
BDIR =../../bin
ODIR = obj
OODIR = ../IntPak/obj
SDIR =.
FC = gfortran -I$(IDIR)
ifeq ($(DEBUG),1)
FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
else
FFLAGS = -Wall -Wno-unused -Wno-unused-dummy-argument -O2
endif
LIBS = ~/Dropbox/quack/lib/*.a
#LIBS = -lblas -llapack
SRCF90 = $(wildcard *.f90)
SRC = $(wildcard *.f)
OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.f,$(ODIR)/%.o,$(SRC))
$(ODIR)/%.o: %.f90
$(FC) -c -o $@ $< $(FFLAGS)
$(ODIR)/%.o: %.f
$(FC) -c -o $@ $< $(FFLAGS)
$(BDIR)/MCQC: $(OBJ)
$(FC) -o $@ $^ $(FFLAGS) $(LIBS)
debug:
DEBUG=1 make $(BDIR)/MCQC
#DEBUG=1 make clean $(BDIR)/MCQC
clean:
rm -f $(ODIR)/*.o $(BDIR)/MCQC $(BDIR)/debug

121
src/MCQC/MinMCMP2.f90
View File

@ -1,121 +0,0 @@
subroutine MinMCMP2(nBas,nEl,nC,nO,nV,c,e,EcMP2, &
nMC,nEq,nWalk,dt,nPrint, &
nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
TrialType,Norm,cTrial,gradient,hessian)
! Minimize the variance of MC-MP2
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas,nEl,nC,nO,nV,nMC,nEq,nWalk,nPrint
double precision,intent(in) :: EcMP2(3),dt
double precision,intent(in) :: c(nBas,nBas),e(nBas)
integer,intent(in) :: nShell
integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell)
double precision,intent(in) :: CenterShell(maxShell,3),DShell(maxShell,maxK),ExpShell(maxShell,maxK)
! Local variables
logical :: debug,varmin,mincvg
double precision :: thresh
double precision,allocatable :: max_gradient(:),energy_MCMP2(:),variance_MCMP2(:),error_MCMP2(:),NormTr(:)
double precision :: EcMCMP2(3),Err_EcMCMP2(3),Var_EcMCMP2(3)
integer :: it,nIt,i
! Output variables
integer,intent(in) :: TrialType
double precision,intent(inout):: Norm,cTrial(nBas),gradient(nBas),hessian(nBas,nBas)
! Debuging mode
! debug = .true.
debug = .false.
! Minimization parameters
varmin = .true.
mincvg = .false.
nIt = 10
thresh = 1d-5
allocate(max_gradient(nIt),energy_MCMP2(nIt),variance_MCMP2(nIt),error_MCMP2(nIt),normTr(nIt))
if(TrialType == 1) then
! Use HOMO as guess
cTrial(1:nBas) = c(1:nBas,nEl/2)
! Normalization factor will be computed later
endif
!------------------------------------------------------------------------
! Start MC-MP2 variance minimization
!------------------------------------------------------------------------
it = 0
do while (it < nIt .and. .not.(mincvg))
it = it + 1
write(*,*) '**********************************************************************'
write(*,*) ' Variance minimization of MC-MP2 energy '
write(*,*) '**********************************************************************'
write(*,*) ' Iteration n.',it
write(*,*) '**********************************************************************'
write(*,*)
write(*,*) ' Trial wave function coefficients at iteration n.',it
call matout(nBas,1,cTrial)
write(*,*)
call MCMP2(varmin,nBas,nEl,nC,nO,nV,c,e,EcMP2, &
nMC,nEq,nWalk,dt,nPrint, &
nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
TrialType,Norm,cTrial,gradient,hessian, &
EcMCMP2,Err_EcMCMP2,Var_EcMCMP2)
! Newton update of the coefficients
call Newton(nBas,gradient,hessian,cTrial)
! Check for convergence
max_gradient(it) = maxval(abs(gradient))
energy_MCMP2(it) = EcMCMP2(1)
variance_MCMP2(it) = Var_EcMCMP2(1)
error_MCMP2(it) = Err_EcMCMP2(1)
NormTr(it) = Norm
write(*,*)
write(*,*) 'Maximum gradient at iteration n.',it,':',max_gradient(it)
write(*,*)
if(max_gradient(it) < thresh) then
write(*,*) ' Miracle! Variance minimization of MC-MP2 has converged!'
mincvg = .true.
endif
enddo
write(*,*)
write(*,*) '********************************'
write(*,*) 'Summary of variance minimization'
write(*,*) '********************************'
write(*,*)
write(*,'(A3,A20,A20,A20,A20,A20,A20)') &
'It.','Gradient','Ec(MC-MPC2)','Variance','Error','Ec(MC-MP2)-Ec(MP2)','Norm'
write(*,'(I3,4X,F16.10,4X,F16.10,4X,F16.10,4X,F16.10,4X,F16.10,4X,F16.10)') &
(i,max_gradient(i),energy_MCMP2(i),variance_MCMP2(i),error_MCMP2(i),energy_MCMP2(i)-EcMP2(1),NormTr(i),i=1,it)
write(*,*)
!------------------------------------------------------------------------
! End MC-MP2 variance minimization
!------------------------------------------------------------------------
end subroutine MinMCMP2

67
src/MCQC/NDrift.f90
View File

@ -1,67 +0,0 @@
subroutine NDrift(nBas,nShell,nWalk,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,P,r,w,F)
! Compute quantum force numerically
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nWalk,nBas
double precision,intent(in) :: P(nBas,nBas),r(nWalk,2,3),w(nWalk)
integer,intent(in) :: nShell
integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell)
double precision,intent(in) :: CenterShell(maxShell,3),DShell(maxShell,maxK),ExpShell(maxShell,maxK)
! Local variables
integer :: iW,iEl,ixyz
double precision :: delta
double precision :: wp,wm
double precision,allocatable :: rp(:,:,:),rm(:,:,:),r12p(:),r12m(:)
double precision,allocatable :: gAOp(:,:,:),dgAOp(:,:,:,:),gAOm(:,:,:),dgAOm(:,:,:,:)
double precision,allocatable :: gp(:,:),dgp(:,:,:),gm(:,:),dgm(:,:,:)
! Output variables
double precision,intent(out) :: F(nWalk,2,3)
allocate(rp(nWalk,2,3),rm(nWalk,2,3),r12p(nWalk),r12m(nWalk), &
gAOp(nWalk,2,nBas),dgAOp(nWalk,2,3,nBas),gAOm(nWalk,2,nBas),dgAOm(nWalk,2,3,nBas), &
gp(nWalk,2),dgp(nWalk,2,3),gm(nWalk,2),dgm(nWalk,2,3))
do iW=1,nWalk
do iEl=1,2
do ixyz=1,3
delta = 1d-6
rp = r
rm = r
rp(iW,iEl,ixyz) = r(iW,iEl,ixyz) + delta
rm(iW,iEl,ixyz) = r(iW,iEl,ixyz) - delta
call AO_values(.false.,nBas,nShell,nWalk,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,rp,gAOp,dgAOp)
call AO_values(.false.,nBas,nShell,nWalk,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,rm,gAOm,dgAOm)
call Density(.false.,nBas,nWalk,P,gAOp,dgAOp,gp,dgp)
call Density(.false.,nBas,nWalk,P,gAOm,dgAOm,gm,dgm)
call rij(nWalk,rp,r12p)
call rij(nWalk,rm,r12m)
wp = gp(iW,1)*gp(iW,2)/r12p(iW)
wm = gm(iW,1)*gm(iW,2)/r12m(iW)
F(iW,iEl,ixyz) = (wp - wm)/(2d0*delta*w(iw))
enddo
enddo
enddo
! print*,'NF',F
end subroutine NDrift

67
src/MCQC/Newton.f90
View File

@ -1,67 +0,0 @@
subroutine Newton(nWSq,gradient,hessian,cWeight)
! Calculate the Green functions
implicit none
! Input variables
integer,intent(in) :: nWSq
double precision,intent(in) :: gradient(nWSq),hessian(nWSq,nWSq)
! Local variables
integer :: info
integer,allocatable :: ipiv(:)
double precision,allocatable :: scr(:),eigval(:),eigvec(:,:)
! Output variables
double precision,intent(inout):: cWeight(nWSq)
! Memory allocation
allocate(ipiv(nWSq),scr(3*nWsq),eigval(nWSq),eigvec(nWSq,nWSq))
! Compute eigenvectors and eigenvalues
eigvec = hessian
call dsyev('V','U',nWSq,eigvec,nWSq,eigval,scr,3*nWSq,info)
if(info /= 0)then
write(*,*) ' Problem with dsyev!'
stop
endif
write(*,*)
write(*,*) 'Eigenvalues of hessian'
call matout(nWSq,1,eigval)
write(*,*)
! write(*,*) 'Eigenvectors of hessian'
! call matout(nWSq,1,eigval)
! write(*,*)
! Compute inverse of the hessian
call dgetrf(nWSq,nWSq,hessian,nWSq,ipiv,info)
if(info /= 0) then
write(*,*) ' Problem in dgetrf!'
stop
endif
call dgetri(nWSq,hessian,nWSq,ipiv,scr,nWSq,info)
if(info /= 0) then
write(*,*) ' Problem in dgetri!'
stop
endif
print*,'inverse hessian'
call matout(nWSq,nWSq,hessian)
! Compute new coefficients
cWeight = cWeight - matmul(hessian,gradient)
end subroutine Newton

29
src/MCQC/NormCoeff.f90
View File

@ -1,29 +0,0 @@
function NormCoeff(alpha,a)
implicit none
! Input variables
double precision,intent(in) :: alpha
integer,intent(in) :: a(3)
! local variable
double precision :: pi,dfa(3),dfac
integer :: atot
! Output variable
double precision NormCoeff
pi = 4d0*atan(1d0)
atot = a(1) + a(2) + a(3)
dfa(1) = dfac(2*a(1))/(2d0**a(1)*dfac(a(1)))
dfa(2) = dfac(2*a(2))/(2d0**a(2)*dfac(a(2)))
dfa(3) = dfac(2*a(3))/(2d0**a(3)*dfac(a(3)))
NormCoeff = (2d0*alpha/pi)**(3d0/2d0)*(4d0*alpha)**atot
NormCoeff = NormCoeff/(dfa(1)*dfa(2)*dfa(3))
NormCoeff = sqrt(NormCoeff)
end function NormCoeff

171
src/MCQC/RHF.f90
View File

@ -1,171 +0,0 @@
subroutine RHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,c,e,P)
! Perform restricted Hartree-Fock calculation
implicit none
! Input variables
integer,intent(in) :: maxSCF,max_diis,guess_type
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 :: nSCF,nBasSq,n_diis
double precision :: ET,EV,EJ,EK,Conv,Gap
double precision,external :: trace_matrix
double precision,allocatable :: error(:,:),error_diis(:,:),F_diis(:,:)
double precision,allocatable :: J(:,:),K(:,:),cp(:,:),F(:,:),Fp(:,:)
! Output variables
double precision,intent(out) :: ERHF,c(nBas,nBas),e(nBas),P(nBas,nBas)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Restricted Hartree-Fock calculation |'
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), &
error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Guess coefficients and eigenvalues
if(guess_type == 1) then
Fp = matmul(transpose(X),matmul(Hc,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c = matmul(X,cp)
elseif(guess_type == 2) then
call random_number(c)
endif
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! Initialization
n_diis = 0
F_diis(:,:) = 0d0
error_diis(:,:) = 0d0
Conv = 1d0
nSCF = 0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
write(*,*)
write(*,*)'----------------------------------------------------'
write(*,*)'| RHF 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(:,:) + 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(nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
! Diagonalize Fock matrix
Fp = matmul(transpose(X),matmul(F,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c = matmul(X,cp)
! Density matrix
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! 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
! 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 RHF

171
src/MCQC/RHF.f90.x
View File

@ -1,171 +0,0 @@
subroutine RHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,c,e,P)
! Perform restricted Hartree-Fock calculation
implicit none
! Input variables
integer,intent(in) :: maxSCF,max_diis,guess_type
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 :: nSCF,nBasSq,n_diis
double precision :: ET,EV,EJ,EK,Conv,Gap
double precision,external :: trace_matrix
double precision,allocatable :: error(:,:),error_diis(:,:),F_diis(:,:)
double precision,allocatable :: J(:,:),K(:,:),cp(:,:),F(:,:),Fp(:,:)
! Output variables
double precision,intent(out) :: ERHF,c(nBas,nBas),e(nBas),P(nBas,nBas)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Restricted Hartree-Fock calculation |'
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), &
error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Guess coefficients and eigenvalues
if(guess_type == 1) then
Fp = matmul(transpose(X),matmul(Hc,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c = matmul(X,cp)
elseif(guess_type == 2) then
call random_number(c)
endif
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! Initialization
n_diis = 0
F_diis(:,:) = 0d0
error_diis(:,:) = 0d0
Conv = 1d0
nSCF = 0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
write(*,*)
write(*,*)'----------------------------------------------------'
write(*,*)'| RHF 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(:,:) + 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(nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
! Diagonalize Fock matrix
Fp = matmul(transpose(X),matmul(F,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c = matmul(X,cp)
! Density matrix
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! 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
! 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 RHF

170
src/MCQC/SPHF.f90
View File

@ -1,170 +0,0 @@
subroutine SPHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,c,e,P)
! Perform restricted Hartree-Fock calculation
implicit none
! Input variables
integer,intent(in) :: maxSCF,max_diis,guess_type
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 :: nSCF,nBasSq,n_diis
double precision :: ET,EV,EJ,EK,Conv,Gap
double precision,external :: trace_matrix
double precision,allocatable :: error(:,:),error_diis(:,:),F_diis(:,:)
double precision,allocatable :: J(:,:),K(:,:),cp(:,:),F(:,:),Fp(:,:)
! Output variables
double precision,intent(out) :: ERHF,c(nBas,nBas),e(nBas),P(nBas,nBas)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Restricted Hartree-Fock calculation |'
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), &
error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Guess coefficients and eigenvalues
if(guess_type == 1) then
Fp = matmul(transpose(X),matmul(Hc,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c = matmul(X,cp)
elseif(guess_type == 2) then
call random_number(c)
endif
P(:,:) = matmul(c(:,1:nO),transpose(c(:,1:nO)))
! Initialization
n_diis = 0
F_diis(:,:) = 0d0
error_diis(:,:) = 0d0
Conv = 1d0
nSCF = 0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
write(*,*)
write(*,*)'----------------------------------------------------'
write(*,*)'| SPHF 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(:,:) + 2d0*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(nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
! Diagonalize Fock matrix
Fp = matmul(transpose(X),matmul(F,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c = matmul(X,cp)
! Density matrix
P(:,:) = matmul(c(:,1:nO),transpose(c(:,1:nO)))
! Compute HF energy
ERHF = trace_matrix(nBas,matmul(P,Hc)) &
+ 0.5d0*trace_matrix(nBas,matmul(P,J)) &
+ 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
! 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 = 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 SPHF

71
src/MCQC/SPMP2.f90
View File

@ -1,71 +0,0 @@
subroutine SPMP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
! Perform third-order Moller-Plesset calculation
implicit none
! Input variables
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: ENuc,EHF
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),e(nBas)
! Local variables
integer :: i,j,a,b
double precision :: eps,E2a,E2b
! Output variables
double precision,intent(out) :: EcMP2(3)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Moller-Plesset second-order calculation |'
write(*,*)'************************************************'
write(*,*)
! Compute MP2 energy
E2a = 0d0
E2b = 0d0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = e(i) + e(j) - e(a) - e(b)
! Secon-order ring diagram
E2a = E2a + ERI(i,j,a,b)*ERI(i,j,a,b)/eps
! Second-order exchange
E2b = E2b + ERI(i,j,a,b)*ERI(i,j,b,a)/eps
enddo
enddo
enddo
enddo
EcMP2(2) = E2a
EcMP2(3) = -E2b
EcMP2(1) = EcMP2(2) + EcMP2(3)
write(*,*)
write(*,'(A32)') '-----------------------'
write(*,'(A32)') ' MP2 calculation '
write(*,'(A32)') '-----------------------'
write(*,'(A32,1X,F16.10)') ' MP2 correlation energy',EcMP2(1)
write(*,'(A32,1X,F16.10)') ' Direct part ',EcMP2(2)
write(*,'(A32,1X,F16.10)') ' Exchange part ',EcMP2(3)
write(*,'(A32)') '-----------------------'
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy',EHF + EcMP2(1)
write(*,'(A32,1X,F16.10)') ' MP2 total energy',ENuc + EHF + EcMP2(1)
write(*,'(A32)') '-----------------------'
write(*,*)
end subroutine SPMP2

77
src/MCQC/SPTDHF.f90
View File

@ -1,77 +0,0 @@
subroutine SPTDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ERI,e)
! Perform random phase approximation calculation
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: singlet_manifold,triplet_manifold
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),e(nBas)
! Local variables
logical :: dRPA,TDA,BSE
integer :: ispin
double precision,allocatable :: Omega(:,:),XpY(:,:,:)
double precision :: rho
double precision :: EcRPA
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Time-dependent Hartree-Fock calculation |'
write(*,*)'************************************************'
write(*,*)
! Switch on exchange for TDHF
dRPA = .false.
! Switch off Tamm-Dancoff approximation for TDHF
TDA = .false.
! Switch off Bethe-Salpeter equation for TDHF
BSE = .false.
! Memory allocation
allocate(Omega(nS,nspin),XpY(nS,nS,nspin))
! Singlet manifold
if(singlet_manifold) then
ispin = 1
call SP_linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI, &
rho,EcRPA,Omega(:,ispin),XpY(:,:,ispin))
call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
endif
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A27,F15.6)') 'RPA correlation energy =',EcRPA
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Triplet manifold
if(triplet_manifold) then
ispin = 2
call SP_linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI, &
rho,EcRPA,Omega(:,ispin),XpY(:,:,ispin))
call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
endif
end subroutine SPTDHF

81
src/MCQC/SP_linear_response.f90
View File

@ -1,81 +0,0 @@
subroutine SP_linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,EcRPA,Omega,XpY)
! Compute linear response
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dRPA,TDA,BSE
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas),rho(nBas,nBas,nS)
! Local variables
double precision :: trace_matrix
double precision,allocatable :: A(:,:),B(:,:),ApB(:,:),AmB(:,:),AmBSq(:,:),Z(:,:)
! Output variables
double precision,intent(out) :: EcRPA
double precision,intent(out) :: Omega(nS),XpY(nS,nS)
! Memory allocation
allocate(A(nS,nS),B(nS,nS),ApB(nS,nS),AmB(nS,nS),AmBSq(nS,nS),Z(nS,nS))
! Build A and B matrices
call SP_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,e,ERI,A)
if(BSE) call Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,A)
! Tamm-Dancoff approximation
B = 0d0
if(.not. TDA) then
call SP_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,B)
if(BSE) call Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,B)
endif
! Build A + B and A - B matrices
AmB = A - B
ApB = A + B
! print*,'A+B'
! call matout(nS,nS,ApB)
! print*,'A-B'
! call matout(nS,nS,AmB)
! Diagonalize TD-HF matrix
call diagonalize_matrix(nS,AmB,Omega)
if(minval(Omega) < 0d0) &
call print_warning('You may have instabilities in linear response!!')
call ADAt(nS,AmB,sqrt(Omega),AmBSq)
Z = matmul(AmBSq,matmul(ApB,AmBSq))
call diagonalize_matrix(nS,Z,Omega)
if(minval(Omega) < 0d0) &
call print_warning('You may have instabilities in linear response!!')
Omega = sqrt(Omega)
XpY = matmul(transpose(Z),AmBSq)
call DA(nS,1d0/sqrt(Omega),XpY)
! print*,'RPA excitations'
! call matout(nS,1,Omega)
! Compute the RPA correlation energy
EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nS,A))
end subroutine SP_linear_response

56
src/MCQC/SP_linear_response_A_matrix.f90
View File

@ -1,56 +0,0 @@
subroutine SP_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,e,ERI,A_lr)
! Compute linear response
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dRPA
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas)
! Local variables
double precision :: delta_spin,delta_dRPA
double precision :: Kronecker_delta
integer :: i,j,a,b,ia,jb
! Output variables
double precision,intent(out) :: A_lr(nS,nS)
! Singlet or triplet manifold?
delta_spin = 0d0
if(ispin == 1) delta_spin = +1d0
if(ispin == 2) delta_spin = -1d0
! Direct RPA
delta_dRPA = 0d0
if(dRPA) delta_dRPA = 1d0
! Build A matrix
ia = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
ia = ia + 1
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+ 0.5d0*(1d0 + delta_spin)*ERI(i,b,a,j) &
- (1d0 - delta_dRPA)*ERI(i,b,j,a)
enddo
enddo
enddo
enddo
end subroutine SP_linear_response_A_matrix

54
src/MCQC/SP_linear_response_B_matrix.f90
View File

@ -1,54 +0,0 @@
subroutine SP_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,B_lr)
! Compute linear response
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dRPA
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
double precision :: delta_spin,delta_dRPA
integer :: i,j,a,b,ia,jb
! Output variables
double precision,intent(out) :: B_lr(nS,nS)
! Singlet or triplet manifold?
delta_spin = 0d0
if(ispin == 1) delta_spin = +1d0
if(ispin == 2) delta_spin = -1d0
! Direct RPA
delta_dRPA = 0d0
if(dRPA) delta_dRPA = 1d0
! Build A matrix
ia = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
ia = ia + 1
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
B_lr(ia,jb) = 0.5d0*(1d0 + delta_spin)*ERI(i,j,a,b) &
- (1d0 - delta_dRPA)*ERI(i,j,b,a)
enddo
enddo
enddo
enddo
end subroutine SP_linear_response_B_matrix

77
src/MCQC/TDHF.f90
View File

@ -1,77 +0,0 @@
subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ERI,e)
! Perform random phase approximation calculation
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: singlet_manifold,triplet_manifold
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),e(nBas)
! Local variables
logical :: dRPA,TDA,BSE
integer :: ispin
double precision,allocatable :: Omega(:,:),XpY(:,:,:)
double precision :: rho
double precision :: EcRPA
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Time-dependent Hartree-Fock calculation |'
write(*,*)'************************************************'
write(*,*)
! Switch on exchange for TDHF
dRPA = .false.
! Switch off Tamm-Dancoff approximation for TDHF
TDA = .false.
! Switch off Bethe-Salpeter equation for TDHF
BSE = .false.
! Memory allocation
allocate(Omega(nS,nspin),XpY(nS,nS,nspin))
! Singlet manifold
if(singlet_manifold) then
ispin = 1
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI, &
rho,EcRPA,Omega(:,ispin),XpY(:,:,ispin))
call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
endif
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A27,F15.6)') 'RPA correlation energy =',EcRPA
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Triplet manifold
if(triplet_manifold) then
ispin = 2
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI, &
rho,EcRPA,Omega(:,ispin),XpY(:,:,ispin))
call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
endif
end subroutine TDHF

237
src/MCQC/UHF.f90
View File

@ -1,237 +0,0 @@
subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUHF)
! Perform unrestricted Hartree-Fock calculation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
double precision,intent(in) :: thresh
integer,intent(in) :: nBas
integer,intent(in) :: nO(nspin)
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) :: ENuc
! Local variables
integer :: nSCF
integer :: nBasSq
integer :: n_diis
double precision :: conv
double precision :: rcond(nspin)
double precision :: ET(nspin)
double precision :: EV(nspin)
double precision :: EJ(nsp)
double precision :: Ex(nspin)
double precision :: Ec(nsp)
double precision :: EUHF
double precision,allocatable :: eps(:,:)
double precision,allocatable :: c(:,:,:)
double precision,allocatable :: cp(:,:,:)
double precision,allocatable :: J(:,:,:)
double precision,allocatable :: F(:,:,:)
double precision,allocatable :: Fp(:,:,:)
double precision,allocatable :: Fx(:,:,:)
double precision,allocatable :: err(:,:,:)
double precision,allocatable :: err_diis(:,:,:)
double precision,allocatable :: F_diis(:,:,:)
double precision,external :: trace_matrix
double precision,allocatable :: P(:,:,:)
integer :: ispin
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'* Unrestricted Hartree-Fock calculation *'
write(*,*)'************************************************'
write(*,*)
! Useful stuff
nBasSq = nBas*nBas
! Memory allocation
allocate(eps(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin), &
J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas,nspin), &
Fx(nBas,nBas,nspin),err(nBas,nBas,nspin),P(nBas,nBas,nspin), &
err_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin))
! Guess coefficients and eigenvalues
if(guess_type == 1) then
do ispin=1,nspin
F(:,:,ispin) = Hc(:,:)
end do
else if(guess_type == 2) then
do ispin=1,nspin
call random_number(F(:,:,ispin))
end do
end if
! Initialization
nSCF = 0
conv = 1d0
n_diis = 0
F_diis(:,:,:) = 0d0
err_diis(:,:,:) = 0d0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
write(*,*)
write(*,*)'------------------------------------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') &
'|','#','|','E(KS)','|','Ex(KS)','|','Ec(KS)','|','Conv','|'
write(*,*)'------------------------------------------------------------------------------------------'
do while(conv > thresh .and. nSCF < maxSCF)
! Increment
nSCF = nSCF + 1
! Transform Fock matrix in orthogonal basis
do ispin=1,nspin
Fp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(F(:,:,ispin),X(:,:)))
end do
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
cp(:,:,:) = Fp(:,:,:)
do ispin=1,nspin
call diagonalize_matrix(nBas,cp(:,:,ispin),eps(:,ispin))
end do
! Back-transform eigenvectors in non-orthogonal basis
do ispin=1,nspin
c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
end do
! Compute density matrix
do ispin=1,nspin
P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
end do
! Build Coulomb repulsion
do ispin=1,nspin
call Coulomb_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin))
end do
! Compute exchange potential
do ispin=1,nspin
call exchange_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),Fx(:,:,ispin))
end do
! Build Fock operator
do ispin=1,nspin
F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + Fx(:,:,ispin)
end do
! Check convergence
do ispin=1,nspin
err(:,:,ispin) = matmul(F(:,:,ispin),matmul(P(:,:,ispin),S(:,:))) - matmul(matmul(S(:,:),P(:,:,ispin)),F(:,:,ispin))
end do
conv = maxval(abs(err(:,:,:)))
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
do ispin=1,nspin
call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin))
end do
! Reset DIIS if required
if(minval(rcond(:)) < 1d-15) n_diis = 0
!------------------------------------------------------------------------
! Compute UHF energy
!------------------------------------------------------------------------
! Kinetic energy
do ispin=1,nspin
ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:)))
end do
! Potential energy
do ispin=1,nspin
EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:)))
end do
! Coulomb energy
EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1)))
EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2)))
EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2)))
! Exchange energy
do ispin=1,nspin
Ex(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),Fx(:,:,ispin)))
end do
! Total energy
EUHF = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) + sum(Ec(:))
! Dump results
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') &
'|',nSCF,'|',EUHF + ENuc,'|',sum(Ex(:)),'|',sum(Ec(:)),'|',conv,'|'
end do
write(*,*)'------------------------------------------------------------------------------------------'
!------------------------------------------------------------------------
! End of SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
stop
end if
! Compute final UHF energy
call print_UHF(nBas,nO(:),eps(:,:),c(:,:,:),ENuc,ET(:),EV(:),EJ(:),Ex(:),Ec(:),EUHF)
end subroutine UHF

46
src/MCQC/antisymmetrize_ERI.f90
View File

@ -1,46 +0,0 @@
subroutine antisymmetrize_ERI(ispin,nBas,ERI,db_ERI)
! Antisymmetrize ERIs
implicit none
! Input variables
integer,intent(in) :: ispin,nBas
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: i,j,k,l
! Output variables
double precision,intent(out) :: db_ERI(nBas,nBas,nBas,nBas)
if(ispin == 1) then
do i=1,nBas
do j=1,nBas
do k=1,nBas
do l=1,nBas
db_ERI(i,j,k,l) = 2d0*ERI(i,j,k,l) - ERI(i,j,l,k)
enddo
enddo
enddo
enddo
elseif(ispin == 2) then
do i=1,nBas
do j=1,nBas
do k=1,nBas
do l=1,nBas
db_ERI(i,j,k,l) = ERI(i,j,k,l) - ERI(i,j,l,k)
enddo
enddo
enddo
enddo
endif
end subroutine antisymmetrize_ERI

33
src/MCQC/chem_to_phys_ERI.f90
View File

@ -1,33 +0,0 @@
subroutine chem_to_phys_ERI(nBas,ERI)
! Antisymmetrize ERIs
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(inout):: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: p,q,r,s
double precision,allocatable :: pERI(:,:,:,:)
allocate(pERI(nBas,nBas,nBas,nBas))
do p=1,nBas
do q=1,nBas
do r=1,nBas
do s=1,nBas
pERI(p,q,r,s) = ERI(p,r,q,s)
enddo
enddo
enddo
enddo
ERI(:,:,:,:) = pERI(:,:,:,:)
end subroutine chem_to_phys_ERI

84
src/MCQC/dcgw.f90
View File

@ -1,84 +0,0 @@
function SigC_dcgw(x,y) result(SigC)
! Degeneracy-corrected GW
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: x,y
! Local variables
double precision,parameter :: eta = 0.1d0
double precision :: r
! Output variables
double precision :: SigC
! Compute the divergence-free term
r = y/x
! Bare style
SigC = y*r
! DCPT style
! SigC = -0.5d0*x*(1d0-sqrt(1d0+4d0*r*r))
! QDPT style
! SigC = y*r/sqrt(1d0+4d0*r*r)
! Infinitesimal
! SigC = y*y*x/(x*x+eta*eta)
end function SigC_dcgw
function Z_dcgw(x,y) result(Z)
! Derivative of the degeneracy-corrected GW
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: x,y
! Local variables
double precision,parameter :: eta = 1d0
double precision :: r
! Output variables
double precision :: Z
! Compute the divergence-free term
r = y/x
! Bare style
Z = r*r
! DCPT style
! Z = 0.5d0*(1d0-1d0/sqrt(1d0+4d0*r*r))
! QDPT style
! Z = r/sqrt(1d0+4d0*r*r)/(1d0+4d0*r*r)
! Infinitesimal
! Z = y*y*(x*x-eta*eta)/(x*x+eta*eta)**2
end function Z_dcgw

51
src/MCQC/density.f90
View File

@ -1,51 +0,0 @@
subroutine density(doDrift,nBas,nWalk,P,gAO,dgAO,g,dg)
! Calculate the Green functions
implicit none
! Input variables
logical,intent(in) :: doDrift
integer,intent(in) :: nBas,nWalk
double precision,intent(in) :: P(nBas,nBas),gAO(nWalk,2,nBas),dgAO(nWalk,2,3,nBas)
! Local variables
integer :: iW,iEl,ixyz,mu,nu
! Output variables
double precision,intent(out) :: g(nWalk,2),dg(nWalk,2,3)
g = 0d0
do iW=1,nWalk
do iEl=1,2
do mu=1,nBas
do nu=1,nBas
g(iW,iEl) = g(iW,iEl) + gAO(iW,iEl,mu)*P(mu,nu)*gAO(iW,iEl,nu)
enddo
enddo
enddo
enddo
if(doDrift) then
dg = 0d0
do iW=1,nWalk
do iEl=1,2
do ixyz=1,3
do mu=1,nBas
do nu=1,nBas
dg(iW,iEl,ixyz) = dg(iW,iEl,ixyz) &
+ P(mu,nu)*(dgAO(iW,iEl,ixyz,mu)*gAO(iW,iEl,nu) &
+ gAO(iW,iEl,mu)*dgAO(iW,iEl,ixyz,nu))
enddo
enddo
enddo
enddo
enddo
endif
end subroutine density

30
src/MCQC/density_matrix.f90
View File

@ -1,30 +0,0 @@
subroutine density_matrix(nBas,ON,c,P)
! Compute density matrix based on the occupation numbers
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: ON(nBas),c(nBas,nBas)
! Local variables
integer :: mu,nu,i
! Output variables
double precision,intent(out) :: P(nBas,nBas)
P(:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do i=1,nBas
P(mu,nu) = P(mu,nu) + 2d0*ON(i)*c(mu,i)*c(nu,i)
enddo
enddo
enddo
end subroutine density_matrix

50
src/MCQC/drift.f90
View File

@ -1,50 +0,0 @@
subroutine drift(nWalk,r,r12,g,dg,F)
! Compute quantum force
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nWalk
double precision,intent(in) :: r(nWalk,2,3),r12(nWalk),g(nWalk,2),dg(nWalk,2,3)
! Local variables
logical :: smoothDrift
double precision :: rij,rijSq,w,wSq,damp
integer :: iW
! Output variables
double precision,intent(out) :: F(nWalk,2,3)
! Compute
smoothDrift = .false.
w = 0.1d0
wSq = w*w
do iW=1,nWalk
rij = r12(iW)
rijSq = rij*rij
F(iW,1,1:3) = dg(iW,1,1:3)/g(iW,1)
F(iW,2,1:3) = dg(iW,2,1:3)/g(iW,2)
if(smoothDrift) then
damp = 1d0 + 2d0*w/sqrt(pi)*rij*exp(-wSq*rijSq)/erfc(w*rij)
else
damp = 1d0
endif
F(iW,1,1:3) = F(iW,1,1:3) - damp*(r(iW,2,1:3) - r(iW,1,1:3))/rijSq
F(iW,2,1:3) = F(iW,2,1:3) - damp*(r(iW,2,1:3) - r(iW,1,1:3))/rijSq
enddo
! print*,' F',F
end subroutine drift

42
src/MCQC/eNcusp.f90.x
View File

@ -1,42 +0,0 @@
subroutine eNcusp(nEl,nBas,S,T,V,G,X,ENuc,EHF,c,e,P,F)
! Perform restricted Hartree-Fock calculation
implicit none
! Input variables
integer,intent(in) :: nEl,nBas
double precision,intent(in) :: ENuc,EHF
double precision,intent(in) :: S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),G(nBas,nBas,nBas,nBas),X(nBas,nBas)
double precision,intent(out) :: c(nBas,nBas),e(nBas),P(nBas,nBas),F(nBas,nBas)
! Local variables
integer,parameter :: maxSCF = 128
double precision,parameter :: thresh = 1d-6
integer :: nO,nSCF,lwork,info
double precision :: ET,EV,Conv,Gap
double precision,allocatable :: Hc(:,:),cp(:,:),cO(:,:),Fp(:,:),work(:)
integer :: mu,nu,lambda,sigma,i
! Output variables
! Number of occupied orbitals
if(mod(nEl,2) /= 0) then
write(*,*) 'closed-shell system required!'
stop
endif
nO = nEl/2
! Memory allocation
allocate(Hc(nBas,nBas),cp(nBas,nBas),cO(nBas,nO),Fp(nBas,nBas))
lwork = 3*nBas
allocate(work(lwork))
! Core Hamiltonian
Hc = T + V
end subroutine eNcusp

170
src/MCQC/elements.f90
View File

@ -1,170 +0,0 @@
function element_number(element_name)
implicit none
integer,parameter :: nelement_max = 103
character(len=2),intent(in) :: element_name
integer :: element_number
character(len=2),parameter :: element_list(nelement_max) = &
(/' H', 'He', & ! 2
'Li','Be', ' B',' C',' N',' O',' F','Ne', & ! 10
'Na','Mg', 'Al','Si',' P',' S','Cl','Ar', & ! 18
' K','Ca','Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr', & ! 36
'Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te',' I','Xe', & ! 54
'Cs','Ba', & ! 56
'La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb', & ! 70
'Lu','Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn', & ! 86
'Fr','Ra', & ! 88
'Ac','Th','Pa',' U','Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md','No', & ! 102
'Lr' & ! 103
/)
!=====
integer :: ielement
!=====
ielement=1
do while( ADJUSTL(element_name) /= ADJUSTL(element_list(ielement)) )
if( ielement == nelement_max ) then
write(*,'(a,a)') ' Input symbol ',element_name
write(*,'(a,i3,a)') ' Element symbol is not one of first ',nelement_max,' elements'
write(*,*) '!!! element symbol not understood !!!'
stop
endif
ielement = ielement + 1
enddo
element_number = ielement
end function element_number
function element_core(zval,zatom)
implicit none
double precision,intent(in) :: zval
double precision,intent(in) :: zatom
integer :: element_core
!=====
!
! If zval /= zatom, this is certainly an effective core potential
! and no core states should be frozen.
if( ABS(zval - zatom) > 1d0-3 ) then
element_core = 0
else
if( zval <= 4.00001d0 ) then ! up to Be
element_core = 0
else if( zval <= 12.00001d0 ) then ! up to Mg
element_core = 1
else if( zval <= 30.00001d0 ) then ! up to Ca
element_core = 5
else if( zval <= 48.00001d0 ) then ! up to Sr
element_core = 9
else
write(*,*) '!!! not imlemented in element_core !!!'
stop
endif
endif
end function element_core
function element_covalent_radius(zatom)
! Return covalent radius of an atom
implicit none
include 'parameters.h'
integer,intent(in) :: zatom
double precision :: element_covalent_radius
!
! Data from Cambridge Structural Database
! http://en.wikipedia.org/wiki/Covalent_radius
!
! Values are first given in picometer
! They will be converted in bohr later on
select case(zatom)
case( 1)
element_covalent_radius = 31.
case( 2)
element_covalent_radius = 28.
case( 3)
element_covalent_radius = 128.
case( 4)
element_covalent_radius = 96.
case( 5)
element_covalent_radius = 84.
case( 6)
element_covalent_radius = 73.
case( 7)
element_covalent_radius = 71.
case( 8)
element_covalent_radius = 66.
case( 9)
element_covalent_radius = 57.
case(10) ! Ne.
element_covalent_radius = 58.
case(11)
element_covalent_radius = 166.
case(12)
element_covalent_radius = 141.
case(13)
element_covalent_radius = 121.
case(14)
element_covalent_radius = 111.
case(15)
element_covalent_radius = 107.
case(16)
element_covalent_radius = 105.
case(17)
element_covalent_radius = 102.
case(18) ! Ar.
element_covalent_radius = 106.
case(19)
element_covalent_radius = 203.
case(20)
element_covalent_radius = 176.
case(21)
element_covalent_radius = 170.
case(22)
element_covalent_radius = 160.
case(23)
element_covalent_radius = 153.
case(24)
element_covalent_radius = 139.
case(25)
element_covalent_radius = 145.
case(26)
element_covalent_radius = 145.
case(27)
element_covalent_radius = 140.
case(28)
element_covalent_radius = 124.
case(29)
element_covalent_radius = 132.
case(30)
element_covalent_radius = 122.
case(31)
element_covalent_radius = 120.
case(32)
element_covalent_radius = 119.
case(34)
element_covalent_radius = 120.
case(35)
element_covalent_radius = 120.
case(36) ! Kr.
element_covalent_radius = 116.
case default
write(*,*) '!!! covalent radius not available !!!'
stop
end select
! pm to bohr conversion
element_covalent_radius = element_covalent_radius*pmtoau
end function element_covalent_radius

207
src/MCQC/evGW.f90
View File

@ -1,207 +0,0 @@
subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,linearize, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF,eG0W0)
! Perform self-consistent eigenvalue-only GW calculation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: maxSCF,max_diis
double precision,intent(in) :: thresh,ENuc,ERHF
logical,intent(in) :: COHSEX,SOSEX,BSE,TDA,G0W,GW0
logical,intent(in) :: singlet_manifold,triplet_manifold,linearize
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: cHF(nBas,nBas),eHF(nBas),eG0W0(nBas),Hc(nBas,nBas),PHF(nBas,nBas)
double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas),ERI_MO_basis(nBas,nBas,nBas,nBas)
! Local variables
logical :: dRPA,linear_mixing
integer :: ispin,nSCF,n_diis
double precision :: Conv,EcRPA,lambda
double precision,allocatable :: error_diis(:,:),e_diis(:,:)
double precision,allocatable :: eGW(:),eOld(:),Z(:)
double precision,allocatable :: H(:,:),SigmaC(:)
double precision,allocatable :: Omega(:,:),XpY(:,:,:),rho(:,:,:,:),rhox(:,:,:,:)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Self-consistent evGW calculation |'
write(*,*)'************************************************'
write(*,*)
! SOSEX correction
if(SOSEX) write(*,*) 'SOSEX correction activated!'
write(*,*)
! Switch off exchange for G0W0
dRPA = .true.
! Linear mixing
linear_mixing = .false.
lambda = 0.2d0
! Memory allocation
allocate(eGW(nBas),eOld(nBas),Z(nBas), &
H(nBas,nBas),SigmaC(nBas), &
Omega(nS,nspin),XpY(nS,nS,nspin), &
rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), &
error_diis(nBas,max_diis),e_diis(nBas,max_diis))
! Initialization
nSCF = 0
ispin = 1
n_diis = 0
Conv = 1d0
e_diis(:,:) = 0d0
error_diis(:,:) = 0d0
eGW(:) = eG0W0(:)
eOld(:) = eGW(:)
Z(:) = 1d0
! Compute Hartree Hamiltonian in the MO basis
call Hartree_matrix_MO_basis(nBas,cHF,PHF,Hc,ERI_AO_basis,H)
!------------------------------------------------------------------------
! Main loop
!------------------------------------------------------------------------
do while(Conv > thresh .and. nSCF <= maxSCF)
! Compute linear response
if(.not. GW0 .or. nSCF == 0) then
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin))
endif
! Compute correlation part of the self-energy
call excitation_density(nBas,nC,nO,nR,nS,cHF,ERI_AO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,cHF,ERI_AO_basis,XpY(:,:,ispin),rhox(:,:,:,ispin))
! Correlation self-energy
if(G0W) then
call self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,eHF, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),SigmaC)
call renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,eHF,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z)
else
call self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,eGW, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),SigmaC)
call renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z)
endif
! Solve the quasi-particle equation (linearized or not)
if(linearize) then
eGW(:) = eHF(:) + Z(:)*SigmaC(:)
else
eGW(:) = eHF(:) + SigmaC(:)
endif
! Convergence criteria
Conv = maxval(abs(eGW - eOld))
! Print results
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigmaC,Z,eGW,EcRPA)
! Linear mixing or DIIS extrapolation
if(linear_mixing) then
eGW(:) = lambda*eGW(:) + (1d0 - lambda)*eOld(:)
else
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(nBas,nBas,n_diis,error_diis,e_diis,eGW-eOld,eGW)
endif
! Save quasiparticles energy for next cycle
eOld(:) = eGW(:)
! Increment
nSCF = nSCF + 1
enddo
!------------------------------------------------------------------------
! End main loop
!------------------------------------------------------------------------
! Plot stuff
call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin))
! Did it actually converge?
if(nSCF == maxSCF+1) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
if(BSE) stop
endif
! Perform BSE calculation
if(BSE) then
! Singlet manifold
if(singlet_manifold) then
ispin = 1
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
endif
! Triplet manifold
if(triplet_manifold) then
ispin = 2
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,cHF,ERI_AO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
endif
endif
end subroutine evGW

35
src/MCQC/exchange_matrix_AO_basis.f90
View File

@ -1,35 +0,0 @@
subroutine exchange_matrix_AO_basis(nBas,P,G,K)
! Compute exchange matrix in the AO basis
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: G(nBas,nBas,nBas,nBas)
! Local variables
integer :: mu,nu,la,si
! Output variables
double precision,intent(out) :: K(nBas,nBas)
K = 0d0
do nu=1,nBas
do si=1,nBas
do la=1,nBas
do mu=1,nBas
K(mu,nu) = K(mu,nu) + P(la,si)*G(mu,la,si,nu)
enddo
enddo
enddo
enddo
K = -0.5d0*K
end subroutine exchange_matrix_AO_basis

26
src/MCQC/exchange_matrix_MO_basis.f90
View File

@ -1,26 +0,0 @@
subroutine exchange_matrix_MO_basis(nBas,c,P,G,K)
! Compute exchange matrix in the MO basis
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: c(nBas,nBas),P(nBas,nBas)
double precision,intent(in) :: G(nBas,nBas,nBas,nBas)
! Output variables
double precision,intent(out) :: K(nBas,nBas)
! Compute Hartree Hamiltonian in the AO basis
call exchange_matrix_AO_basis(nBas,P,G,K)
! Transform Coulomb matrix in the MO basis
K = matmul(transpose(c),matmul(K,c))
end subroutine exchange_matrix_MO_basis

65
src/MCQC/excitation_density.f90
View File

@ -1,65 +0,0 @@
subroutine excitation_density(nBas,nC,nO,nR,nS,c,G,XpY,rho)
! Compute excitation densities
implicit none
! Input variables
integer,intent(in) :: nBas,nC,nO,nR,nS
double precision,intent(in) :: c(nBas,nBas),G(nBas,nBas,nBas,nBas),XpY(nS,nS)
! Local variables
double precision,allocatable :: scr(:,:,:)
integer :: mu,nu,la,si,ia,jb,x,y,j,b
! Output variables
double precision,intent(out) :: rho(nBas,nBas,nS)
! Memory allocation
allocate(scr(nBas,nBas,nS))
rho(:,:,:) = 0d0
do nu=1,nBas
do si=1,nBas
do ia=1,nS
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
rho(nu,si,ia) = rho(nu,si,ia) + c(nu,j)*XpY(ia,jb)*c(si,b)
enddo
enddo
enddo
enddo
enddo
scr(:,:,:) = 0d0
do mu=1,nBas
do la=1,nBas
do ia=1,nS
do nu=1,nBas
do si=1,nBas
scr(mu,la,ia) = scr(mu,la,ia) + G(mu,nu,la,si)*rho(nu,si,ia)
enddo
enddo
enddo
enddo
enddo
rho(:,:,:) = 0d0
do ia=1,nS
do x=nC+1,nBas-nR
do y=nC+1,nBas-nR
do mu=1,nBas
do la=1,nBas
rho(x,y,ia) = rho(x,y,ia) + c(mu,x)*scr(mu,la,ia)*c(la,y)
enddo
enddo
enddo
enddo
enddo
end subroutine excitation_density

65
src/MCQC/excitation_density_SOSEX.f90
View File

@ -1,65 +0,0 @@
subroutine excitation_density_SOSEX(nBas,nC,nO,nR,nS,c,G,XpY,rho)
! Compute excitation densities
implicit none
! Input variables
integer,intent(in) :: nBas,nC,nO,nR,nS
double precision,intent(in) :: c(nBas,nBas),G(nBas,nBas,nBas,nBas),XpY(nS,nS)
! Local variables
double precision,allocatable :: scr(:,:,:)
integer :: mu,nu,la,si,ia,jb,x,y,j,b
! Output variables
double precision,intent(out) :: rho(nBas,nBas,nS)
! Memory allocation
allocate(scr(nBas,nBas,nS))
rho(:,:,:) = 0d0
do nu=1,nBas
do si=1,nBas
do ia=1,nS
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
rho(nu,si,ia) = rho(nu,si,ia) + c(nu,j)*XpY(ia,jb)*c(si,b)
enddo
enddo
enddo
enddo
enddo
scr(:,:,:) = 0d0
do mu=1,nBas
do la=1,nBas
do ia=1,nS
do nu=1,nBas
do si=1,nBas
scr(mu,la,ia) = scr(mu,la,ia) + G(mu,nu,la,si)*rho(nu,si,ia)
enddo
enddo
enddo
enddo
enddo
rho(:,:,:) = 0d0
do ia=1,nS
do x=nC+1,nBas-nR
do y=nC+1,nBas-nR
do mu=1,nBas
do la=1,nBas
rho(x,y,ia) = rho(x,y,ia) + c(mu,x)*scr(mu,la,ia)*c(la,y)
enddo
enddo
enddo
enddo
enddo
end subroutine excitation_density_SOSEX

35
src/MCQC/excitation_density_SOSEX_from_MO.f90
View File

@ -1,35 +0,0 @@
subroutine excitation_density_SOSEX_from_MO(nBas,nC,nO,nR,nS,G,XpY,rho)
! Compute excitation densities
implicit none
! Input variables
integer,intent(in) :: nBas,nC,nO,nR,nS
double precision,intent(in) :: G(nBas,nBas,nBas,nBas),XpY(nS,nS)
! Local variables
integer :: ia,jb,x,y,j,b
! Output variables
double precision,intent(out) :: rho(nBas,nBas,nS)
rho(:,:,:) = 0d0
do ia=1,nS
do x=nC+1,nBas-nR
do y=nC+1,nBas-nR
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
rho(x,y,ia) = rho(x,y,ia) + G(x,y,b,j)*XpY(ia,jb)
enddo
enddo
enddo
enddo
enddo
end subroutine excitation_density_SOSEX_from_MO

35
src/MCQC/excitation_density_from_MO.f90
View File

@ -1,35 +0,0 @@
subroutine excitation_density_from_MO(nBas,nC,nO,nR,nS,G,XpY,rho)
! Compute excitation densities
implicit none
! Input variables
integer,intent(in) :: nBas,nC,nO,nR,nS
double precision,intent(in) :: G(nBas,nBas,nBas,nBas),XpY(nS,nS)
! Local variables
integer :: ia,jb,x,y,j,b
! Output variables
double precision,intent(out) :: rho(nBas,nBas,nS)
rho(:,:,:) = 0d0
do ia=1,nS
do x=nC+1,nBas-nR
do y=nC+1,nBas-nR
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
rho(x,y,ia) = rho(x,y,ia) + G(x,j,y,b)*XpY(ia,jb)
enddo
enddo
enddo
enddo
enddo
end subroutine excitation_density_from_MO

60
src/MCQC/form_CABS.f90.x
View File

@ -1,60 +0,0 @@
subroutine form_CABS(nBas_OBS,nBas_ABS,c_OBS,c_ABS,S_ABS)
! Perform configuration interaction single calculation`
implicit none
! Input variables
integer,intent(in) :: nBas_OBS,nBas_ABS
double precision,intent(in) :: S_ABS(nBas,nBas),c_OBS(nBas_OBS,nBas_OBS)
! Local variables
integer ::
double precision :: thresh = 1d-07
integer :: i,j,a,b
! Output variables
double precision,intent(out) :: c_ABS(nBas_ABS,nBas_ABS)
allocate(c(nBas_ABS,nBAs_OBS))
c = 0d0
c(1:nBas_OBS,1:nBas_OBS) = c_OBS(1:nBas_OBS,1:nBAs_OBS)
c_ABS = 0d0
do i=1,nBas_ABS
c_ABS(i,i) = 1d0
enddo
v_ABS = S_ABS
call DiagMat(nBas_ABS,v_ABS,e_ABS)
nLD = 0
do i=1,nBas_ABS
if(abs(e_ABS(i)) < thresh) nLD = nLD +1
enddo
write(*,*) 'Number of linear dependencies in ABS',nLD
call DoSVD(nBas_ABS,S_ABS,u,v,w)
! do a SVD of S_ABS to get u, v and w
X_ABS = 0d0
do i=1,nBas_ABS
do j=1,nBas_ABS
do k=1,nBas_ABS
X_ABS(i,k) = X_ABS(i,k) + v_ABS(i,j)*e_ABS(j)*v_ABS(k,j)
enddo
enddo
enddo
cp_ABS = matmul(X_ABS,c_ABS)
S12 = matmul(transpose(c),matmul(S_ABS,cp_ABS))
end subroutine form_CABS

46
src/MCQC/form_T.f90
View File

@ -1,46 +0,0 @@
subroutine form_T(nO,nV,delta_OOOVVV,ub,ubb,EcCCT)
! Compute (T) correction
implicit none
! Input variables
integer,intent(in) :: nO,nV
double precision,intent(in) :: delta_OOOVVV(nO,nO,nO,nV,nV,nV)
double precision,intent(in) :: ub(nO,nO,nO,nV,nV,nV)
double precision,intent(in) :: ubb(nO,nO,nO,nV,nV,nV)
! Local variables
integer :: i,j,k,l
integer :: a,b,c,d
! Output variables
double precision,intent(out) :: EcCCT
EcCCT = 0d0
do i=1,nO
do j=1,nO
do k=1,nO
do a=1,nV
do b=1,nV
do c=1,nV
EcCCT = EcCCT &
+ (ub(i,j,k,a,b,c) + ubb(i,j,k,a,b,c)) &
* ubb(i,j,k,a,b,c)/delta_OOOVVV(i,j,k,a,b,c)
end do
end do
end do
end do
end do
end do
EcCCT = - EcCCT/36d0
end subroutine form_T

92
src/MCQC/form_X.f90
View File

@ -1,92 +0,0 @@
subroutine form_X(nO,nV,OOVV,t2,X1,X2,X3,X4)
! Form intermediate arrays X's in CCD
implicit none
! Input variables
integer,intent(in) :: nO,nV
double precision,intent(in) :: t2(nO,nO,nV,nV)
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
! Local variables
integer :: i,j,k,l
integer :: a,b,c,d
! Output variables
double precision,intent(out) :: X1(nO,nO,nO,nO)
double precision,intent(out) :: X2(nV,nV)
double precision,intent(out) :: X3(nO,nO)
double precision,intent(out) :: X4(nO,nO,nV,nV)
! Initialization
X1(:,:,:,:) = 0d0
X2(:,:) = 0d0
X3(:,:) = 0d0
X4(:,:,:,:) = 0d0
! Build X1
do k=1,nO
do l=1,nO
do i=1,nO
do j=1,nO
do c=1,nV
do d=1,nV
X1(k,l,i,j) = X1(k,l,i,j) + OOVV(k,l,c,d)*t2(i,j,c,d)
enddo
enddo
enddo
enddo
enddo
enddo
! Build X2
do b=1,nV
do c=1,nV
do k=1,nO
do l=1,nO
do d=1,nV
X2(b,c) = X2(b,c) + OOVV(k,l,c,d)*t2(k,l,b,d)
enddo
enddo
enddo
enddo
enddo
! Build X3
do k=1,nO
do j=1,nO
do l=1,nO
do c=1,nV
do d=1,nV
X3(k,j) = X3(k,j) + OOVV(k,l,c,d)*t2(j,l,c,d)
enddo
enddo
enddo
enddo
enddo
! Build X4
do i=1,nO
do l=1,nO
do a=1,nV
do d=1,nV
do k=1,nO
do c=1,nV
X4(i,l,a,d) = X4(i,l,a,d) + OOVV(k,l,c,d)*t2(i,k,a,c)
enddo
enddo
enddo
enddo
enddo
enddo
end subroutine form_X

105
src/MCQC/form_abh.f90
View File

@ -1,105 +0,0 @@
subroutine form_abh(nO,nV,OOOO,OVOO,OOVV,VVVV,VOVV,OVVO,OVVV,t1,tau,aoooo,bvvvv,hovvo)
! Scuseria Eqs. (11),(12) and (13)
implicit none
! Input variables
integer,intent(in) :: nO,nV
double precision,intent(in) :: OOOO(nO,nO,nO,nO)
double precision,intent(in) :: OVOO(nO,nV,nO,nO)
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
double precision,intent(in) :: OVVO(nO,nV,nV,nO)
double precision,intent(in) :: OVVV(nO,nV,nV,nV)
double precision,intent(in) :: VOVV(nV,nO,nV,nV)
double precision,intent(in) :: VVVV(nV,nV,nV,nV)
double precision,intent(in) :: t1(nO,nV)
double precision,intent(in) :: tau(nO,nO,nV,nV)
! Local variables
integer :: i,j,k,l
integer :: a,b,c,d
! Output variables
double precision,intent(out) :: aoooo(nO,nO,nO,nO)
double precision,intent(out) :: bvvvv(nV,nV,nV,nV)
double precision,intent(out) :: hovvo(nO,nV,nV,nO)
aoooo(:,:,:,:) = OOOO(:,:,:,:)
do i=1,nO
do j=1,nO
do k=1,nO
do l=1,nO
do c=1,nV
aoooo(i,j,k,l) = aoooo(i,j,k,l) + OVOO(i,c,k,l)*t1(j,c)
end do
do c=1,nV
aoooo(i,j,k,l) = aoooo(i,j,k,l) - OVOO(j,c,k,l)*t1(i,c)
end do
do c=1,nV
do d=1,nV
aoooo(i,j,k,l) = aoooo(i,j,k,l) + OOVV(k,l,c,d)*tau(i,j,c,d)
end do
end do
end do
end do
end do
end do
bvvvv(:,:,:,:) = VVVV(:,:,:,:)
do c=1,nV
do d=1,nV
do a=1,nV
do b=1,nV
do k=1,nO
bvvvv(c,d,a,b) = bvvvv(c,d,a,b) - VOVV(a,k,c,d)*t1(k,b)
end do
do k=1,nO
bvvvv(c,d,a,b) = bvvvv(c,d,a,b) + VOVV(b,k,c,d)*t1(k,a)
end do
end do
end do
end do
end do
hovvo(:,:,:,:) = OVVO(:,:,:,:)
do i=1,nO
do c=1,nV
do a=1,nV
do k=1,nO
do l=1,nO
hovvo(i,c,a,k) = hovvo(i,c,a,k) - OVOO(i,c,l,k)*t1(l,a)
end do
do d=1,nV
hovvo(i,c,a,k) = hovvo(i,c,a,k) + OVVV(k,a,c,d)*t1(i,d)
end do
do l=1,nO
do d=1,nV
hovvo(i,c,a,k) = hovvo(i,c,a,k) - OOVV(k,l,c,d)*tau(i,l,d,a)
end do
end do
end do
end do
end do
end do
end subroutine form_abh

37
src/MCQC/form_delta_OOOVVV.f90
View File

@ -1,37 +0,0 @@
subroutine form_delta_OOOVVV(nO,nV,eO,eV,delta)
! Form energy denominator for CC
implicit none
! Input variables
integer,intent(in) :: nO,nV
double precision,intent(in) :: eO(nO)
double precision,intent(in) :: eV(nV)
! Local variables
integer :: i,j,k,a,b,c
! Output variables
double precision,intent(out) :: delta(nO,nO,nO,nV,nV,nV)
do i=1,nO
do j=1,nO
do k=1,nO
do a=1,nV
do b=1,nV
do c=1,nV
delta(i,j,k,a,b,c) = eV(a) + eV(b) + eV(c) - eO(i) - eO(j) - eO(k)
enddo
enddo
enddo
enddo
enddo
enddo
end subroutine form_delta_OOOVVV

33
src/MCQC/form_delta_OOVV.f90
View File

@ -1,33 +0,0 @@
subroutine form_delta_OOVV(nO,nV,eO,eV,delta)
! Form energy denominator for CC
implicit none
! Input variables
integer,intent(in) :: nO,nV
double precision,intent(in) :: eO(nO)
double precision,intent(in) :: eV(nV)
! Local variables
integer :: i,j,a,b
! Output variables
double precision,intent(out) :: delta(nO,nO,nV,nV)
do i=1,nO
do j=1,nO
do a=1,nV
do b=1,nV
delta(i,j,a,b) = eV(a) + eV(b) - eO(i) - eO(j)
enddo
enddo
enddo
enddo
end subroutine form_delta_OOVV

27
src/MCQC/form_delta_OV.f90
View File

@ -1,27 +0,0 @@
subroutine form_delta_OV(nO,nV,eO,eV,delta)
! Form energy denominator for CC
implicit none
! Input variables
integer,intent(in) :: nO,nV
double precision,intent(in) :: eO(nO)
double precision,intent(in) :: eV(nV)
! Local variables
integer :: i,a
! Output variables
double precision,intent(out) :: delta(nO,nV)
do i=1,nO
do a=1,nV
delta(i,a) = eV(a) - eO(i)
enddo
enddo
end subroutine form_delta_OV

53
src/MCQC/form_g.f90
View File

@ -1,53 +0,0 @@
subroutine form_g(nO,nV,hvv,hoo,VOVV,OOOV,t1,gvv,goo)
! Scuseria Eqs. (9), (10)
implicit none
! Input variables
integer,intent(in) :: nO,nV
double precision,intent(in) :: hvv(nV,nV)
double precision,intent(in) :: hoo(nO,nO)
double precision,intent(in) :: VOVV(nV,nO,nV,nV)
double precision,intent(in) :: OOOV(nO,nO,nO,nV)
double precision,intent(in) :: t1(nO,nV)
! Local variables
integer :: i,j,k,l
integer :: a,b,c,d
! Output variables
double precision,intent(out) :: gvv(nV,nV)
double precision,intent(out) :: goo(nO,nO)
gvv(:,:) = hvv(:,:)
do c=1,nV
do a=1,nV
do k=1,nO
do d=1,nV
gvv(c,a) = gvv(c,a) + VOVV(a,k,c,d)*t1(k,d)
end do
end do
end do
end do
goo(:,:) = hoo(:,:)
do i=1,nO
do k=1,nO
do l=1,nO
do c=1,nV
goo(i,k) = goo(i,k) + OOOV(k,l,i,c)*t1(l,c)
end do
end do
end do
end do
end subroutine form_g

79
src/MCQC/form_h.f90
View File

@ -1,79 +0,0 @@
subroutine form_h(nO,nV,eO,eV,OOVV,t1,tau,hvv,hoo,hvo)
! Scuseria Eqs. (5), (6) and (7)
implicit none
! Input variables
integer,intent(in) :: nO,nV
double precision,intent(in) :: eO(nO)
double precision,intent(in) :: eV(nV)
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
double precision,intent(in) :: t1(nO,nV)
double precision,intent(in) :: tau(nO,nO,nV,nV)
! Local variables
integer :: i,j,k,l
integer :: a,b,c,d
! Output variables
double precision,intent(out) :: hvv(nV,nV)
double precision,intent(out) :: hoo(nO,nO)
double precision,intent(out) :: hvo(nV,nO)
hvv(:,:) = 0d0
do b=1,nV
hvv(b,b) = eV(b)
do a=1,nV
do j=1,nO
do k=1,nO
do c=1,nV
hvv(b,a) = hvv(b,a) - OOVV(j,k,b,c)*tau(j,k,a,c)
end do
end do
end do
end do
end do
hoo(:,:) = 0d0
do i=1,nO
hoo(i,i) = eO(i)
do j=1,nO
do k=1,nO
do b=1,nV
do c=1,nV
hoo(i,j) = hoo(i,j) + OOVV(j,k,b,c)*tau(i,k,b,c)
end do
end do
end do
end do
end do
hvo(:,:) = 0d0
do b=1,nV
do j=1,nO
do k=1,nO
do c=1,nV
hvo(b,j) = hvo(b,j) + OOVV(j,k,b,c)*t1(k,c)
end do
end do
end do
end do
! print*,'hvv',hvv
end subroutine form_h

77
src/MCQC/form_r1.f90
View File

@ -1,77 +0,0 @@
subroutine form_r1(nO,nV,OVVO,OVVV,OOOV,hvv,hoo,hvo,t1,t2,tau,r1)
! Form tau in CCSD
implicit none
! Input variables
integer,intent(in) :: nO,nV
double precision,intent(in) :: OVVO(nO,nV,nV,nO)
double precision,intent(in) :: OVVV(nO,nV,nV,nV)
double precision,intent(in) :: OOOV(nO,nO,nO,nV)
double precision,intent(in) :: hvv(nV,nV)
double precision,intent(in) :: hoo(nO,nO)
double precision,intent(in) :: hvo(nV,nO)
double precision,intent(in) :: t1(nO,nV)
double precision,intent(in) :: t2(nO,nO,nV,nV)
double precision,intent(in) :: tau(nO,nO,nV,nV)
! Local variables
integer :: i,j,k,l
integer :: a,b,c,d
! Output variables
double precision,intent(out) :: r1(nO,nV)
r1(:,:) = 0d0
do i=1,nO
do a=1,nV
do b=1,nV
r1(i,a) = r1(i,a) + hvv(b,a)*t1(i,b)
end do
do j=1,nO
r1(i,a) = r1(i,a) - hoo(i,j)*t1(j,a)
end do
do j=1,nO
do b=1,nV
r1(i,a) = r1(i,a) + hvo(b,j)*(t2(i,j,a,b) + t1(i,b)*t1(j,a))
end do
end do
do j=1,nO
do b=1,nV
r1(i,a) = r1(i,a) + OVVO(i,b,a,j)*t1(j,b)
end do
end do
do j=1,nO
do b=1,nV
do c=1,nV
r1(i,a) = r1(i,a) - OVVV(j,a,b,c)*tau(i,j,b,c)
end do
end do
end do
do j=1,nO
do k=1,nO
do b=1,nV
r1(i,a) = r1(i,a) - OOOV(j,k,i,b)*tau(j,k,a,b)
end do
end do
end do
end do
end do
end subroutine form_r1

139
src/MCQC/form_r2.f90
View File

@ -1,139 +0,0 @@
subroutine form_r2(nO,nV,OOVV,OVOO,OVVV,OVVO,gvv,goo,aoooo,bvvvv,hovvo,t1,t2,tau,r2)
! Form tau in CCSD
implicit none
! Input variables
integer,intent(in) :: nO,nV
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
double precision,intent(in) :: OVOO(nO,nV,nO,nO)
double precision,intent(in) :: OVVV(nO,nV,nV,nV)
double precision,intent(in) :: OVVO(nO,nV,nV,nO)
double precision,intent(in) :: gvv(nV,nV)
double precision,intent(in) :: goo(nO,nO)
double precision,intent(in) :: aoooo(nO,nO,nO,nO)
double precision,intent(in) :: bvvvv(nV,nV,nV,nV)
double precision,intent(in) :: hovvo(nO,nV,nV,nO)
double precision,intent(in) :: t1(nO,nV)
double precision,intent(in) :: t2(nO,nO,nV,nV)
double precision,intent(in) :: tau(nO,nO,nV,nV)
! Local variables
integer :: i,j,k,l
integer :: a,b,c,d
! Output variables
double precision,intent(out) :: r2(nO,nO,nV,nV)
r2(:,:,:,:) = OOVV(:,:,:,:)
do i=1,nO
do j=1,nO
do a=1,nV
do b=1,nV
do k=1,nO
do l=1,nO
r2(i,j,a,b) = r2(i,j,a,b) + aoooo(i,j,k,l)*tau(k,l,a,b)
end do
end do
do c=1,nV
do d=1,nV
r2(i,j,a,b) = r2(i,j,a,b) + bvvvv(c,d,a,b)*tau(i,j,c,d)
end do
end do
do c=1,nV
r2(i,j,a,b) = r2(i,j,a,b) + gvv(c,a)*t2(i,j,c,b)
end do
do k=1,nO
r2(i,j,a,b) = r2(i,j,a,b) + OVOO(k,a,i,j)*t1(k,b)
end do
do c=1,nV
r2(i,j,a,b) = r2(i,j,a,b) - gvv(c,b)*t2(i,j,c,a)
end do
do k=1,nO
r2(i,j,a,b) = r2(i,j,a,b) - OVOO(k,b,i,j)*t1(k,a)
end do
do k=1,nO
r2(i,j,a,b) = r2(i,j,a,b) - goo(i,k)*t2(k,j,a,b)
end do
do c=1,nV
r2(i,j,a,b) = r2(i,j,a,b) + OVVV(j,c,b,a)*t1(i,c)
end do
do k=1,nO
r2(i,j,a,b) = r2(i,j,a,b) + goo(j,k)*t2(k,i,a,b)
end do
do c=1,nV
r2(i,j,a,b) = r2(i,j,a,b) - OVVV(i,c,b,a)*t1(j,c)
end do
do k=1,nO
do c=1,nV
r2(i,j,a,b) = r2(i,j,a,b) + hovvo(i,c,a,k)*t2(j,k,b,c)
end do
end do
do k=1,nO
do c=1,nV
r2(i,j,a,b) = r2(i,j,a,b) - OVVO(i,c,a,k)*t1(j,c)*t1(k,b)
end do
end do
do k=1,nO
do c=1,nV
r2(i,j,a,b) = r2(i,j,a,b) - hovvo(j,c,a,k)*t2(i,k,b,c)
end do
end do
do k=1,nO
do c=1,nV
r2(i,j,a,b) = r2(i,j,a,b) + OVVO(j,c,a,k)*t1(i,c)*t1(k,b)
end do
end do
do k=1,nO
do c=1,nV
r2(i,j,a,b) = r2(i,j,a,b) - hovvo(i,c,b,k)*t2(j,k,a,c)
end do
end do
do k=1,nO
do c=1,nV
r2(i,j,a,b) = r2(i,j,a,b) + OVVO(i,c,b,k)*t1(j,c)*t1(k,a)
end do
end do
do k=1,nO
do c=1,nV
r2(i,j,a,b) = r2(i,j,a,b) + hovvo(j,c,b,k)*t2(i,k,a,c)
end do
end do
do k=1,nO
do c=1,nV
r2(i,j,a,b) = r2(i,j,a,b) - OVVO(j,c,b,k)*t1(i,c)*t1(k,a)
end do
end do
end do
end do
end do
end do
end subroutine form_r2

34
src/MCQC/form_tau.f90
View File

@ -1,34 +0,0 @@
subroutine form_tau(nO,nV,t1,t2,tau)
! Form tau in CCSD
implicit none
! Input variables
integer,intent(in) :: nO,nV
double precision,intent(in) :: t1(nO,nV)
double precision,intent(in) :: t2(nO,nO,nV,nV)
! Local variables
integer :: i,j,k,l
integer :: a,b,c,d
! Output variables
double precision,intent(out) :: tau(nO,nO,nV,nV)
do i=1,nO
do j=1,nO
do a=1,nV
do b=1,nV
tau(i,j,a,b) = 0.5d0*t2(i,j,a,b) + t1(i,a)*t1(j,b)
enddo
enddo
enddo
enddo
end subroutine form_tau

71
src/MCQC/form_u.f90
View File

@ -1,71 +0,0 @@
subroutine form_u(nO,nV,OOOO,VVVV,OVOV,t2,u)
! Form linear array in CCD
implicit none
! Input variables
integer,intent(in) :: nO,nV
double precision,intent(in) :: t2(nO,nO,nV,nV)
double precision,intent(in) :: OOOO(nO,nO,nO,nO)
double precision,intent(in) :: VVVV(nV,nV,nV,nV)
double precision,intent(in) :: OVOV(nO,nV,nO,nV)
! Local variables
integer :: i,j,k,l
integer :: a,b,c,d
! Output variables
double precision,intent(out) :: u(nO,nO,nV,nV)
u(:,:,:,:) = 0d0
do i=1,nO
do j=1,nO
do a=1,nV
do b=1,nV
do c=1,nV
do d=1,nV
u(i,j,a,b) = u(i,j,a,b) + 0.5d0*VVVV(a,b,c,d)*t2(i,j,c,d)
enddo
enddo
enddo
enddo
enddo
enddo
do i=1,nO
do j=1,nO
do k=1,nO
do l=1,nO
do a=1,nV
do b=1,nV
u(i,j,a,b) = u(i,j,a,b) + 0.5d0*OOOO(k,l,i,j)*t2(k,l,a,b)
enddo
enddo
enddo
enddo
enddo
enddo
do i=1,nO
do j=1,nO
do k=1,nO
do a=1,nV
do b=1,nV
do c=1,nV
u(i,j,a,b) = u(i,j,a,b) - OVOV(k,b,j,c)*t2(i,k,a,c) &
+ OVOV(k,a,j,c)*t2(i,k,b,c) &
- OVOV(k,a,i,c)*t2(j,k,b,c) &
+ OVOV(k,b,i,c)*t2(j,k,a,c)
enddo
enddo
enddo
enddo
enddo
enddo
end subroutine form_u

48
src/MCQC/form_ub.f90
View File

@ -1,48 +0,0 @@
subroutine form_ub(nO,nV,OOVV,t1,ub)
! Form 1st term in (T) correction
implicit none
! Input variables
integer,intent(in) :: nO,nV
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
double precision,intent(in) :: t1(nO,nV)
! Local variables
integer :: i,j,k,l
integer :: a,b,c,d
! Output variables
double precision,intent(out) :: ub(nO,nO,nO,nV,nV,nV)
do i=1,nO
do j=1,nO
do k=1,nO
do a=1,nV
do b=1,nV
do c=1,nV
ub(i,j,k,a,b,c) = t1(i,a)*OOVV(j,k,b,c) &
+ t1(i,b)*OOVV(j,k,c,a) &
+ t1(i,c)*OOVV(j,k,a,b) &
+ t1(j,a)*OOVV(k,i,b,c) &
+ t1(j,b)*OOVV(k,i,c,a) &
+ t1(j,c)*OOVV(k,i,a,b) &
+ t1(k,a)*OOVV(i,j,b,c) &
+ t1(k,b)*OOVV(i,j,c,a) &
+ t1(k,c)*OOVV(i,j,a,b)
end do
end do
end do
end do
end do
end do
end subroutine form_ub

67
src/MCQC/form_ubb.f90
View File

@ -1,67 +0,0 @@
subroutine form_ubb(nO,nV,VVVO,VOOO,t2,ubb)
! Form 2nd term in (T) correction
implicit none
! Input variables
integer,intent(in) :: nO,nV
double precision,intent(in) :: VVVO(nV,nV,nV,nO)
double precision,intent(in) :: VOOO(nV,nO,nO,nO)
double precision,intent(in) :: t2(nO,nO,nV,nV)
! Local variables
integer :: i,j,k,l,m
integer :: a,b,c,d,e
! Output variables
double precision,intent(out) :: ubb(nO,nO,nO,nV,nV,nV)
ubb(:,:,:,:,:,:) = 0d0
do i=1,nO
do j=1,nO
do k=1,nO
do a=1,nV
do b=1,nV
do c=1,nV
do e=1,nV
ubb(i,j,k,a,b,c) = ubb(i,j,k,a,b,c) &
+ t2(i,j,a,e)*VVVO(b,c,e,k) &
+ t2(i,j,b,e)*VVVO(c,a,e,k) &
+ t2(i,j,c,e)*VVVO(a,b,e,k) &
+ t2(k,i,a,e)*VVVO(b,c,e,j) &
+ t2(k,i,b,e)*VVVO(c,a,e,j) &
+ t2(k,i,c,e)*VVVO(a,b,e,j) &
+ t2(j,k,a,e)*VVVO(b,c,e,i) &
+ t2(j,k,b,e)*VVVO(c,a,e,i) &
+ t2(j,k,c,e)*VVVO(a,b,e,i)
end do
do m=1,nO
ubb(i,j,k,a,b,c) = ubb(i,j,k,a,b,c) &
+ t2(i,m,a,b)*VOOO(c,m,j,k) &
+ t2(i,m,b,c)*VOOO(a,m,j,k) &
+ t2(i,m,c,a)*VOOO(b,m,j,k) &
+ t2(j,m,a,b)*VOOO(c,m,k,i) &
+ t2(j,m,b,c)*VOOO(a,m,k,i) &
+ t2(j,m,c,a)*VOOO(b,m,k,i) &
+ t2(k,m,a,b)*VOOO(c,m,i,j) &
+ t2(k,m,b,c)*VOOO(a,m,i,j) &
+ t2(k,m,c,a)*VOOO(b,m,i,j)
end do
end do
end do
end do
end do
end do
end do
end subroutine form_ubb

79
src/MCQC/form_v.f90
View File

@ -1,79 +0,0 @@
subroutine form_v(nO,nV,X1,X2,X3,X4,t2,v)
! Form quadratic array in CCD
implicit none
! Input variables
integer,intent(in) :: nO,nV
double precision,intent(in) :: t2(nO,nO,nV,nV)
double precision,intent(in) :: X1(nO,nO,nO,nO)
double precision,intent(in) :: X2(nV,nV)
double precision,intent(in) :: X3(nO,nO)
double precision,intent(in) :: X4(nO,nO,nV,nV)
! Local variables
integer :: i,j,k,l
integer :: a,b,c,d
! Output variables
double precision,intent(out) :: v(nO,nO,nV,nV)
v(:,:,:,:) = 0d0
do i=1,nO
do j=1,nO
do a=1,nV
do b=1,nV
do k=1,nO
do l=1,nO
v(i,j,a,b) = v(i,j,a,b) + 0.25d0*X1(k,l,i,j)*t2(k,l,a,b)
enddo
enddo
enddo
enddo
enddo
enddo
do i=1,nO
do j=1,nO
do a=1,nV
do b=1,nV
do c=1,nV
v(i,j,a,b) = v(i,j,a,b) - 0.5d0*(X2(b,c)*t2(i,j,a,c) + X2(a,c)*t2(i,j,c,b))
enddo
enddo
enddo
enddo
enddo
do i=1,nO
do j=1,nO
do a=1,nV
do b=1,nV
do k=1,nO
v(i,j,a,b) = v(i,j,a,b) - 0.5d0*(X3(k,j)*t2(i,k,a,b) + X3(k,i)*t2(k,j,a,b))
enddo
enddo
enddo
enddo
enddo
do i=1,nO
do j=1,nO
do a=1,nV
do b=1,nV
do k=1,nO
do c=1,nV
v(i,j,a,b) = v(i,j,a,b) + (X4(i,k,a,c)*t2(j,k,b,c) + X4(i,k,b,c)*t2(k,j,a,c))
enddo
enddo
enddo
enddo
enddo
enddo
end subroutine form_v

30
src/MCQC/generate_shell.f90
View File

@ -1,30 +0,0 @@
subroutine generate_shell(atot,nShellFunction,ShellFunction)
implicit none
! Input variables
integer,intent(in) :: atot,nShellFunction
! Local variables
integer :: ax,ay,az,ia
! Output variables
integer,intent(out) :: ShellFunction(nShellFunction,3)
ia = 0
do ax=atot,0,-1
do az=0,atot
ay = atot - ax - az
if(ay >= 0) then
ia = ia + 1
ShellFunction(ia,1) = ax
ShellFunction(ia,2) = ay
ShellFunction(ia,3) = az
endif
enddo
enddo
end subroutine generate_shell

25
src/MCQC/initialize_random_generator.f90
View File

@ -1,25 +0,0 @@
subroutine initialize_random_generator(iSeed)
! Initialize random number generator
implicit none
! Input variables
integer,intent(in) :: iSeed
! Local variables
integer,allocatable :: Seed(:)
integer :: nSeed
call random_seed(size = nSeed)
allocate(Seed(nSeed))
call random_seed(get=Seed)
if(iSeed /= 0) then
Seed = 0
Seed(1) = iSeed
endif
call random_seed(put=Seed)
end subroutine initialize_random_generator

81
src/MCQC/linear_response.f90
View File

@ -1,81 +0,0 @@
subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,EcRPA,Omega,XpY)
! Compute linear response
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dRPA,TDA,BSE
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas),rho(nBas,nBas,nS)
! Local variables
double precision :: trace_matrix
double precision,allocatable :: A(:,:),B(:,:),ApB(:,:),AmB(:,:),AmBSq(:,:),Z(:,:)
! Output variables
double precision,intent(out) :: EcRPA
double precision,intent(out) :: Omega(nS),XpY(nS,nS)
! Memory allocation
allocate(A(nS,nS),B(nS,nS),ApB(nS,nS),AmB(nS,nS),AmBSq(nS,nS),Z(nS,nS))
! Build A and B matrices
call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,e,ERI,A)
if(BSE) call Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,A)
! Tamm-Dancoff approximation
B = 0d0
if(.not. TDA) then
call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,B)
if(BSE) call Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,B)
endif
! Build A + B and A - B matrices
AmB = A - B
ApB = A + B
! print*,'A+B'
! call matout(nS,nS,ApB)
! print*,'A-B'
! call matout(nS,nS,AmB)
! Diagonalize TD-HF matrix
call diagonalize_matrix(nS,AmB,Omega)
if(minval(Omega) < 0d0) &
call print_warning('You may have instabilities in linear response!!')
call ADAt(nS,AmB,sqrt(Omega),AmBSq)
Z = matmul(AmBSq,matmul(ApB,AmBSq))
call diagonalize_matrix(nS,Z,Omega)
if(minval(Omega) < 0d0) &
call print_warning('You may have instabilities in linear response!!')
Omega = sqrt(Omega)
XpY = matmul(transpose(Z),AmBSq)
call DA(nS,1d0/sqrt(Omega),XpY)
! print*,'RPA excitations'
! call matout(nS,1,Omega)
! Compute the RPA correlation energy
EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nS,A))
end subroutine linear_response

56
src/MCQC/linear_response_A_matrix.f90
View File

@ -1,56 +0,0 @@
subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,e,ERI,A_lr)
! Compute linear response
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dRPA
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas)
! Local variables
double precision :: delta_spin,delta_dRPA
double precision :: Kronecker_delta
integer :: i,j,a,b,ia,jb
! Output variables
double precision,intent(out) :: A_lr(nS,nS)
! Singlet or triplet manifold?
delta_spin = 0d0
if(ispin == 1) delta_spin = +1d0
if(ispin == 2) delta_spin = -1d0
! Direct RPA
delta_dRPA = 0d0
if(dRPA) delta_dRPA = 1d0
! Build A matrix
ia = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
ia = ia + 1
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+ (1d0 + delta_spin)*ERI(i,b,a,j) &
- (1d0 - delta_dRPA)*ERI(i,b,j,a)
enddo
enddo
enddo
enddo
end subroutine linear_response_A_matrix

54
src/MCQC/linear_response_B_matrix.f90
View File

@ -1,54 +0,0 @@
subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,B_lr)
! Compute linear response
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dRPA
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
double precision :: delta_spin,delta_dRPA
integer :: i,j,a,b,ia,jb
! Output variables
double precision,intent(out) :: B_lr(nS,nS)
! Singlet or triplet manifold?
delta_spin = 0d0
if(ispin == 1) delta_spin = +1d0
if(ispin == 2) delta_spin = -1d0
! Direct RPA
delta_dRPA = 0d0
if(dRPA) delta_dRPA = 1d0
! Build A matrix
ia = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
ia = ia + 1
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
B_lr(ia,jb) = (1d0 + delta_spin)*ERI(i,j,a,b) &
- (1d0 - delta_dRPA)*ERI(i,j,b,a)
enddo
enddo
enddo
enddo
end subroutine linear_response_B_matrix

57
src/MCQC/natural_orbital.f90
View File

@ -1,57 +0,0 @@
subroutine natural_orbital(nBas,nO,cHF,c)
! Compute natural orbitals and natural occupancies
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas,nO
double precision,intent(in) :: cHF(nBas,nBas),c(nBas,nBas)
! Local variables
integer :: i,j,k
double precision,allocatable :: eNO(:),cNO(:,:),P(:,:)
! Allocate
allocate(eNO(nBas),cNO(nBas,nBas),P(nBas,nBas))
! Compute density matrix
P = matmul(transpose(cHF),cHF)
call matout(nBas,nBas,P)
cNO = 0d0
do i=1,nBas
do j=1,nBas
do k=1,1
cNO(i,j) = cNO(i,j) + 2d0*P(i,k)*P(j,k)
enddo
enddo
enddo
! cNO(:,:) = matmul(c(:,1:nO),transpose(c(:,1:nO)))
! cNO = matmul(transpose(cHF),matmul(cNO,cHF))
call diagonalize_matrix(nBas,cNO,eNO)
! Print results
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' Natural orbitals '
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,cNO)
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' Natural occupancies'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,1,eNO)
write(*,*)
end subroutine natural_orbital

53
src/MCQC/norm_trial.f90
View File

@ -1,53 +0,0 @@
subroutine norm_trial(nBas,nO,c,P,Norm,NormSq)
! Initialize weight function
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas,nO
double precision,intent(inout):: c(nBas,nO),P(nBas,nBas)
! Local variables
double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),G(:,:,:,:)
integer :: mu,nu,la,si
! Output variables
double precision,intent(inout):: Norm,NormSq
! Memory allocation for one- and two-electron integrals
allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),G(nBas,nBas,nBas,nBas))
! Read integrals
call read_integrals(nBas,S,T,V,Hc,G)
! Compute normalization factor
P = 2d0*matmul(c,transpose(c))
Norm = 0d0
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do si=1,nBas
Norm = Norm + P(mu,nu)*P(la,si)*G(mu,la,nu,si)
enddo
enddo
enddo
enddo
Norm = Norm*Norm
NormSq = Norm*Norm
write(*,*)
write(*,*) 'Normalization of trial wave function: ',Norm
write(*,*)
end subroutine norm_trial

28
src/MCQC/optimize_timestep.f90
View File

@ -1,28 +0,0 @@
subroutine optimize_timestep(nWalk,iMC,Acc,dt)
! Optimize dt to get 50% of accepted moves
implicit none
! Input variables
integer,intent(in) :: nWalk,iMC
double precision,intent(inout):: Acc,dt
! Local variables
double precision :: TotAcc,Current_Acc,Target_Acc,delta
TotAcc = Acc/dble(nWalk)
Current_Acc = 100d0*TotAcc/dble(iMC)
Target_Acc = 50.0d0
delta = dt*abs(Target_Acc - Current_Acc)/100.d0
if(Current_Acc > Target_Acc + 0.5d0)then
dt = dt + delta
elseif(Current_Acc < Target_Acc - 0.5d0)then
dt = dt - delta
endif
end subroutine optimize_timestep

120
src/MCQC/orthogonalization_matrix.f90
View File

@ -1,120 +0,0 @@
subroutine orthogonalization_matrix(ortho_type,nBas,S,X)
! Compute the orthogonalization matrix X
implicit none
! Input variables
integer,intent(in) :: nBas,ortho_type
double precision,intent(in) :: S(nBas,nBas)
! Local variables
logical :: debug
double precision,allocatable :: UVec(:,:),Uval(:)
double precision,parameter :: thresh = 1d-6
integer :: i
! Output variables
double precision,intent(out) :: X(nBas,nBas)
debug = .false.
! Type of orthogonalization ortho_type
!
! 1 = Lowdin
! 2 = Canonical
! 3 = SVD
!
allocate(Uvec(nBas,nBas),Uval(nBas))
if(ortho_type == 1) then
write(*,*)
write(*,*) ' Lowdin orthogonalization'
write(*,*)
Uvec = S
call diagonalize_matrix(nBas,Uvec,Uval)
do i=1,nBas
if(Uval(i) > thresh) then
Uval(i) = 1d0/sqrt(Uval(i))
else
write(*,*) 'Eigenvalue',i,'too small for Lowdin orthogonalization'
endif
enddo
call ADAt(nBas,Uvec,Uval,X)
elseif(ortho_type == 2) then
write(*,*)
write(*,*) 'Canonical orthogonalization'
write(*,*)
Uvec = S
call diagonalize_matrix(nBas,Uvec,Uval)
do i=1,nBas
if(Uval(i) > thresh) then
Uval(i) = 1d0/sqrt(Uval(i))
else
write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization'
endif
enddo
call AD(nBas,Uvec,Uval)
X = Uvec
elseif(ortho_type == 3) then
write(*,*)
write(*,*) ' SVD-based orthogonalization NYI'
write(*,*)
! Uvec = S
! call diagonalize_matrix(nBas,Uvec,Uval)
! do i=1,nBas
! if(Uval(i) > thresh) then
! Uval(i) = 1d0/sqrt(Uval(i))
! else
! write(*,*) 'Eigenvalue',i,'too small for canonical orthogonalization'
! endif
! enddo
!
! call AD(nBas,Uvec,Uval)
! X = Uvec
endif
! Print results
if(debug) then
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Orthogonalization matrix'
write(*,'(A28)') '----------------------'
call matout(nBas,nBas,X)
write(*,*)
endif
end subroutine orthogonalization_matrix

40
src/MCQC/overlap.f90
View File

@ -1,40 +0,0 @@
subroutine overlap(nBas,bra,ket)
! Compute the overlap between two sets of coefficients
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: bra(nBas,nBas),ket(nBas,nBas)
! Local variables
double precision,allocatable :: s(:),Ov(:,:)
! Allocate
allocate(s(nBas),Ov(nBas,nBas))
! Compute overlap
Ov = matmul(transpose(bra),ket)
call diagonalize_matrix(nBas,Ov,s)
! Print results
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' Overlap '
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,Ov)
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' Eigenvalues of overlap matrix'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,1,s)
write(*,*)
end subroutine overlap

113
src/MCQC/plot_GW.f90
View File

@ -1,113 +0,0 @@
subroutine plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega,rho,rhox)
! Dump several GW quantities for external plotting
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eHF(nBas),eGW(nBas),Omega(nS),rho(nBas,nBas,nS),rhox(nBas,nBas,nS)
! Local variables
integer :: i,j,a,b,x,jb,g
integer :: nGrid
double precision :: eps,eta,wmin,wmax,dw
double precision,allocatable :: w(:),SigC(:,:),Z(:,:),S(:,:)
! Infinitesimal
eta = 1d-3
! Construct grid
nGrid = 1000
allocate(w(nGrid),SigC(nBas,nGrid),Z(nBas,nGrid),S(nBas,nGrid))
! Initialize
SigC(:,:) = 0d0
Z(:,:) = 0d0
! Minimum and maximum frequency values
wmin = -5d0
wmax = +5d0
dw = (wmax - wmin)/dble(ngrid)
do g=1,nGrid
w(g) = wmin + dble(g)*dw
enddo
! Occupied part of the self-energy and renormalization factor
do g=1,nGrid
do x=nC+1,nBas-nR
do i=nC+1,nO
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
eps = w(g) - eHF(i) + Omega(jb)
SigC(x,g) = SigC(x,g) + 2d0*rho(x,i,jb)**2*eps/(eps**2 + eta**2)
Z(x,g) = Z(x,g) + 2d0*rho(x,i,jb)**2/eps**2
enddo
enddo
enddo
enddo
enddo
! Virtual part of the self-energy and renormalization factor
do g=1,nGrid
do x=nC+1,nBas-nR
do a=nO+1,nBas-nR
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
eps = w(g) - eHF(a) - Omega(jb)
SigC(x,g) = SigC(x,g) + 2d0*rho(x,a,jb)**2*eps/(eps**2 + eta**2)
Z(x,g) = Z(x,g) + 2d0*rho(x,a,jb)**2/eps**2
enddo
enddo
enddo
enddo
enddo
Z(:,:) = 1d0/(1d0 + Z(:,:))
! Compute spectral function
do g=1,nGrid
do x=nC+1,nBas-nR
S(x,g) = eta/((w(g) - eHF(x) - SigC(x,g))**2 + eta**2)
enddo
enddo
S(:,:) = S(:,:)/pi
! Dump quantities in files as a function of w
open(unit=8 ,file='plot/grid.dat')
open(unit=9 ,file='plot/SigC.dat')
open(unit=10 ,file='plot/Z.dat')
open(unit=11 ,file='plot/A.dat')
do g=1,nGrid
write(8 ,*) w(g)*HaToeV,(SigC(x,g)*HaToeV,x=1,nBas)
write(9 ,*) w(g)*HaToeV,((w(g)-eHF(x))*HaToeV,x=1,nBas)
write(10,*) w(g)*HaToeV,(Z(x,g),x=1,nBas)
write(11,*) w(g)*HaToeV,(S(x,g),x=1,nBas)
enddo
! Closing files
close(unit=8)
close(unit=9)
close(unit=10)
close(unit=11)
end subroutine plot_GW

49
src/MCQC/print_G0W0.f90
View File

@ -1,49 +0,0 @@
subroutine print_G0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM)
! Print one-electron energies and other stuff for G0W0
implicit none
include 'parameters.h'
integer,intent(in) :: nBas,nO
double precision,intent(in) :: ENuc,EHF,EcRPA,EcGM
double precision,intent(in) :: e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas)
integer :: x,HOMO,LUMO
double precision :: Gap
! HOMO and LUMO
HOMO = nO
LUMO = HOMO + 1
Gap = eGW(LUMO)-eGW(HOMO)
! Dump results
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)' One-shot G0W0 calculation'
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
'|','#','|','e_HF (eV)','|','Sigma_c (eV)','|','Z','|','e_QP (eV)','|'
write(*,*)'-------------------------------------------------------------------------------'
do x=1,nBas
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',x,'|',e(x)*HaToeV,'|',SigmaC(x)*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|'
enddo
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A27,F15.6)') 'G0W0 HOMO energy (eV):',eGW(HOMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'G0W0 LUMO energy (eV):',eGW(LUMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A27,F15.6)') 'G0W0 RPA total energy =',ENuc + EHF + EcRPA
write(*,'(2X,A27,F15.6)') 'RPA correlation energy =',EcRPA
write(*,'(2X,A27,F15.6)') 'G0W0 GM total energy =',ENuc + EHF + EcGM
write(*,'(2X,A27,F15.6)') 'GM correlation energy =',EcGM
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end subroutine print_G0W0

44
src/MCQC/print_GF2.f90
View File

@ -1,44 +0,0 @@
subroutine print_GF2(nBas,nO,nSCF,Conv,eHF,eGF2)
! Print one-electron energies and other stuff for GF2
implicit none
include 'parameters.h'
integer,intent(in) :: nBas,nO,nSCF
double precision,intent(in) :: Conv,eHF(nBas),eGF2(nBas)
integer :: x,HOMO,LUMO
double precision :: Gap
! HOMO and LUMO
HOMO = nO
LUMO = HOMO + 1
Gap = eGF2(LUMO)-eGF2(HOMO)
! Dump results
write(*,*)'-------------------------------------------'
write(*,*)' Frequency-dependent diagonal GF2 calculation'
write(*,*)'-------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
'|','#','|','e_HF (eV)','|','e_GF2 (eV)','|'
write(*,*)'-------------------------------------------'
do x=1,nBas
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',x,'|',eHF(x)*HaToeV,'|',eGF2(x)*HaToeV,'|'
enddo
write(*,*)'-------------------------------------------'
write(*,'(2X,A10,I3)') 'Iteration ',nSCF
write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
write(*,*)'-------------------------------------------'
write(*,'(2X,A27,F15.6)') 'GF2 HOMO energy (eV):',eGF2(HOMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'GF2 LUMO energy (eV):',eGF2(LUMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'GF2 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------'
write(*,*)
end subroutine print_GF2

44
src/MCQC/print_GF3.f90
View File

@ -1,44 +0,0 @@
subroutine print_GF3(nBas,nO,nSCF,Conv,eHF,Z,eGF3)
! Print one-electron energies and other stuff for GF3
implicit none
include 'parameters.h'
integer,intent(in) :: nBas,nO,nSCF
double precision,intent(in) :: Conv,eHF(nBas),eGF3(nBas),Z(nBas)
integer :: x,HOMO,LUMO
double precision :: Gap
! HOMO and LUMO
HOMO = nO
LUMO = HOMO + 1
Gap = eGF3(LUMO)-eGF3(HOMO)
! Dump results
write(*,*)'-------------------------------------------------------------'
write(*,*)' Frequency-dependent diagonal GF3 calculation'
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
'|','#','|','e_HF (eV)','|','Z','|','e_GF3 (eV)','|'
write(*,*)'-------------------------------------------------------------'
do x=1,nBas
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',x,'|',eHF(x)*HaToeV,'|',Z(x),'|',eGF3(x)*HaToeV,'|'
enddo
write(*,*)'-------------------------------------------------------------'
write(*,'(2X,A10,I3)') 'Iteration ',nSCF
write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
write(*,*)'-------------------------------------------------------------'
write(*,'(2X,A27,F15.6)') 'GF3 HOMO energy (eV):',eGF3(HOMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'GF3 LUMO energy (eV):',eGF3(LUMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'GF3 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------------------------'
write(*,*)
end subroutine print_GF3

60
src/MCQC/print_RHF.f90
View File

@ -1,60 +0,0 @@
subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF)
! Print one-electron energies and other stuff for G0W0
implicit none
include 'parameters.h'
integer,intent(in) :: nBas,nO
double precision,intent(in) :: eHF(nBas),cHF(nBas,nBas),ENuc,ET,EV,EJ,EK,ERHF
integer :: HOMO,LUMO
double precision :: Gap
! HOMO and LUMO
HOMO = nO
LUMO = HOMO + 1
Gap = eHF(LUMO)-eHF(HOMO)
! Dump results
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' Summary '
write(*,'(A50)') '---------------------------------------'
write(*,'(A32,1X,F16.10)') ' One-electron energy ',ET + EV
write(*,'(A32,1X,F16.10)') ' Kinetic energy ',ET
write(*,'(A32,1X,F16.10)') ' Potential energy ',EV
write(*,'(A50)') '---------------------------------------'
write(*,'(A32,1X,F16.10)') ' Two-electron energy ',EJ + EK
write(*,'(A32,1X,F16.10)') ' Coulomb energy ',EJ
write(*,'(A32,1X,F16.10)') ' Exchange energy ',EK
write(*,'(A50)') '---------------------------------------'
write(*,'(A32,1X,F16.10)') ' Electronic energy ',ERHF
write(*,'(A32,1X,F16.10)') ' Nuclear repulsion ',ENuc
write(*,'(A32,1X,F16.10)') ' Hartree-Fock energy ',ERHF + ENuc
write(*,'(A50)') '---------------------------------------'
write(*,'(A36,F13.6)') ' HF HOMO energy (eV):',eHF(HOMO)*HaToeV
write(*,'(A36,F13.6)') ' HF LUMO energy (eV):',eHF(LUMO)*HaToeV
write(*,'(A36,F13.6)') ' HF HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,'(A50)') '---------------------------------------'
write(*,*)
! Print results
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') 'MO coefficients'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,cHF)
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') 'MO energies'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,1,eHF)
write(*,*)
end subroutine print_RHF

102
src/MCQC/print_UHF.f90
View File

@ -1,102 +0,0 @@
subroutine print_UHF(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew)
! Print one- and two-electron energies and other stuff for UHF calculation
implicit none
include 'parameters.h'
integer,intent(in) :: nBas
integer,intent(in) :: nO(nspin)
double precision,intent(in) :: eps(nBas,nspin)
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: ENuc
double precision,intent(in) :: ET(nspin)
double precision,intent(in) :: EV(nspin)
double precision,intent(in) :: EJ(nsp)
double precision,intent(in) :: Ex(nspin)
double precision,intent(in) :: Ec(nsp)
double precision,intent(in) :: Ew
integer :: HOMO(nspin)
integer :: LUMO(nspin)
double precision :: Gap(nspin)
! HOMO and LUMO
HOMO(:) = nO(:)
LUMO(:) = HOMO(:) + 1
Gap(1) = eps(LUMO(1),1) - eps(HOMO(1),1)
Gap(2) = eps(LUMO(2),2) - eps(HOMO(2),2)
! Dump results
write(*,*)
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40)') ' Summary '
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,1X,F16.10,A3)') ' One-electron energy: ',sum(ET(:)) + sum(EV(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' One-electron a energy: ',ET(1) + EV(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' One-electron b energy: ',ET(2) + EV(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Kinetic energy: ',sum(ET(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Kinetic a energy: ',ET(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Kinetic b energy: ',ET(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Potential energy: ',sum(EV(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Potential a energy: ',EV(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Potential b energy: ',EV(2),' au'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,1X,F16.10,A3)') ' Two-electron a energy: ',sum(EJ(:)) + sum(Ex(:)) + sum(Ec(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Two-electron aa energy: ',EJ(1) + Ex(1) + Ec(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Two-electron ab energy: ',EJ(2) + Ec(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Two-electron bb energy: ',EJ(3) + Ex(2) + Ec(3),' au'
write(*,'(A40,1X,F16.10,A3)') ' Coulomb energy: ',sum(EJ(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Coulomb aa energy: ',EJ(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Coulomb ab energy: ',EJ(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Coulomb bb energy: ',EJ(3),' au'
write(*,'(A40,1X,F16.10,A3)') ' Exchange energy: ',sum(Ex(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Exchange a energy: ',Ex(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Exchange b energy: ',Ex(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Correlation energy: ',sum(Ec(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Correlation aa energy: ',Ec(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Correlation ab energy: ',Ec(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Correlation bb energy: ',Ec(3),' au'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,1X,F16.10,A3)') ' Electronic energy: ',Ew,' au'
write(*,'(A40,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au'
write(*,'(A40,1X,F16.10,A3)') ' UHF energy: ',Ew + ENuc,' au'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,F13.6,A3)') ' UHF HOMO a energy:',eps(HOMO(1),1)*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' UHF LUMO a energy:',eps(LUMO(1),1)*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' UHF HOMOa-LUMOa gap:',Gap(1)*HatoeV,' eV'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,F13.6,A3)') ' UHF HOMO b energy:',eps(HOMO(2),2)*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' UHF LUMO b energy:',eps(LUMO(2),2)*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' UHF HOMOb-LUMOb gap :',Gap(2)*HatoeV,' eV'
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
! Print results
write(*,'(A50)') '-----------------------------------------'
write(*,'(A50)') 'UHF spin-up orbital coefficients '
write(*,'(A50)') '-----------------------------------------'
call matout(nBas,nBas,c(:,:,1))
write(*,'(A50)') '-----------------------------------------'
write(*,'(A50)') 'UHF spin-down orbital coefficients '
write(*,'(A50)') '-----------------------------------------'
call matout(nBas,nBas,c(:,:,2))
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' UHF spin-up orbital energies '
write(*,'(A50)') '---------------------------------------'
call matout(nBas,1,eps(:,1))
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' UHF spin-down orbital energies '
write(*,'(A50)') '---------------------------------------'
call matout(nBas,1,eps(:,2))
write(*,*)
end subroutine print_UHF

54
src/MCQC/print_evGW.f90
View File

@ -1,54 +0,0 @@
subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM)
! Print one-electron energies and other stuff for evGW
implicit none
include 'parameters.h'
integer,intent(in) :: nBas,nO,nSCF
double precision,intent(in) :: ENuc,EHF,EcRPA,EcGM
double precision,intent(in) :: Conv,e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas)
integer :: x,HOMO,LUMO
double precision :: Gap
! HOMO and LUMO
HOMO = nO
LUMO = HOMO + 1
Gap = eGW(LUMO)-eGW(HOMO)
! Dump results
write(*,*)'-------------------------------------------------------------------------------'
if(nSCF < 10) then
write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent evG',nSCF,'W',nSCF,' calculation'
else
write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent evG',nSCF,'W',nSCF,' 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)') &
'|','#','|','e_HF (eV)','|','Sigma_c (eV)','|','Z','|','e_QP (eV)','|'
write(*,*)'-------------------------------------------------------------------------------'
do x=1,nBas
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',x,'|',e(x)*HaToeV,'|',SigmaC(x)*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|'
enddo
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A10,I3)') 'Iteration ',nSCF
write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A27,F15.6)') 'evGW HOMO energy (eV):',eGW(HOMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'evGW LUMO energy (eV):',eGW(LUMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'evGW HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A27,F15.6)') 'evGW RPA total energy =',ENuc + EHF + EcRPA
write(*,'(2X,A27,F15.6)') 'RPA correlation energy =',EcRPA
write(*,'(2X,A27,F15.6)') 'evGW GM total energy =',ENuc + EHF + EcGM
write(*,'(2X,A27,F15.6)') 'GM correlation energy =',EcGM
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end subroutine print_evGW

36
src/MCQC/print_excitation.f90
View File

@ -1,36 +0,0 @@
subroutine print_excitation(method,ispin,nS,Omega)
! Print excitation energies for a given spin manifold
implicit none
include 'parameters.h'
character*5,intent(in) :: method
integer,intent(in) :: ispin,nS
double precision,intent(in) :: Omega(nS)
character*7 :: spin_manifold
integer :: ia
if(ispin == 1) spin_manifold = 'singlet'
if(ispin == 2) spin_manifold = 'triplet'
write(*,*)
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A1,1X,A4,A14,A7,A9,A25)')'|',method,' calculation: ',spin_manifold,' manifold',' |'
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
'|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'
write(*,*)'-------------------------------------------------------------'
do ia=1,nS
write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') &
'|',ia,'|',Omega(ia),'|',Omega(ia)*HaToeV,'|'
enddo
write(*,*)'-------------------------------------------------------------'
write(*,*)
end subroutine print_excitation

Some files were not shown because too many files have changed in this diff Show More