mirror of
https://github.com/pfloos/quack
synced 2024-10-31 19:23:52 +01:00
core in CC
This commit is contained in:
parent
fb012773cb
commit
4983c8c16b
35
input/basis
35
input/basis
@ -1,9 +1,30 @@
|
||||
1 3
|
||||
S 3
|
||||
1 38.3600000 0.0238090
|
||||
2 5.7700000 0.1548910
|
||||
3 1.2400000 0.4699870
|
||||
1 6
|
||||
S 8
|
||||
1 17880.0000000 0.0007380
|
||||
2 2683.0000000 0.0056770
|
||||
3 611.5000000 0.0288830
|
||||
4 173.5000000 0.1085400
|
||||
5 56.6400000 0.2909070
|
||||
6 20.4200000 0.4483240
|
||||
7 7.8100000 0.2580260
|
||||
8 1.6530000 0.0150630
|
||||
S 8
|
||||
1 17880.0000000 -0.0001720
|
||||
2 2683.0000000 -0.0013570
|
||||
3 611.5000000 -0.0067370
|
||||
4 173.5000000 -0.0276630
|
||||
5 56.6400000 -0.0762080
|
||||
6 20.4200000 -0.1752270
|
||||
7 7.8100000 -0.1070380
|
||||
8 1.6530000 0.5670500
|
||||
S 1
|
||||
1 0.2976000 1.0000000
|
||||
1 0.4869000 1.0000000
|
||||
P 3
|
||||
1 28.3900000 0.0460870
|
||||
2 6.2700000 0.2401810
|
||||
3 1.6950000 0.5087440
|
||||
P 1
|
||||
1 1.2750000 1.0000000
|
||||
1 0.4317000 1.0000000
|
||||
D 1
|
||||
1 2.2020000 1.0000000
|
||||
|
||||
|
@ -3,11 +3,11 @@
|
||||
# Chemist notation for two-electron integral?
|
||||
F
|
||||
# Exposant of the Slater geminal
|
||||
1.0
|
||||
0.5
|
||||
# One-electron integrals: Ov Kin Nuc
|
||||
T T T
|
||||
# Two-electron integrals: ERI F12 Yuk Erf
|
||||
T F F F
|
||||
T F F T
|
||||
# Three-electron integrals: Type1 Type2 Type3
|
||||
F F F
|
||||
# Four-electron integrals: Type1 Type2 Type3
|
||||
|
@ -13,6 +13,6 @@
|
||||
# G0W0 evGW qsGW
|
||||
F F F
|
||||
# G0T0 evGT qsGT
|
||||
T F F
|
||||
F F F
|
||||
# MCMP2
|
||||
F
|
||||
|
@ -1,4 +1,4 @@
|
||||
# nAt nEla nElb nCore nRyd
|
||||
1 1 1 0 0
|
||||
1 5 5 0 0
|
||||
# Znuc x y z
|
||||
He 0.0 0.0 0.0
|
||||
Ne 0.0 0.0 0.0
|
||||
|
@ -1,3 +1,3 @@
|
||||
1
|
||||
|
||||
He 0.0000000000 0.0000000000 0.0000000000
|
||||
Ne 0.0000000000 0.0000000000 0.0000000000
|
||||
|
@ -11,6 +11,6 @@
|
||||
# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta
|
||||
256 0.00001 T 5 F F F F F F F 0.000
|
||||
# ACFDT: AC Kx XBS
|
||||
T F T
|
||||
T T T
|
||||
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
|
||||
1000000 100000 10 0.3 10000 1234 T
|
||||
|
35
input/weight
35
input/weight
@ -1,9 +1,30 @@
|
||||
1 3
|
||||
S 3
|
||||
1 38.3600000 0.0238090
|
||||
2 5.7700000 0.1548910
|
||||
3 1.2400000 0.4699870
|
||||
1 6
|
||||
S 8
|
||||
1 17880.0000000 0.0007380
|
||||
2 2683.0000000 0.0056770
|
||||
3 611.5000000 0.0288830
|
||||
4 173.5000000 0.1085400
|
||||
5 56.6400000 0.2909070
|
||||
6 20.4200000 0.4483240
|
||||
7 7.8100000 0.2580260
|
||||
8 1.6530000 0.0150630
|
||||
S 8
|
||||
1 17880.0000000 -0.0001720
|
||||
2 2683.0000000 -0.0013570
|
||||
3 611.5000000 -0.0067370
|
||||
4 173.5000000 -0.0276630
|
||||
5 56.6400000 -0.0762080
|
||||
6 20.4200000 -0.1752270
|
||||
7 7.8100000 -0.1070380
|
||||
8 1.6530000 0.5670500
|
||||
S 1
|
||||
1 0.2976000 1.0000000
|
||||
1 0.4869000 1.0000000
|
||||
P 3
|
||||
1 28.3900000 0.0460870
|
||||
2 6.2700000 0.2401810
|
||||
3 1.6950000 0.5087440
|
||||
P 1
|
||||
1 1.2750000 1.0000000
|
||||
1 0.4317000 1.0000000
|
||||
D 1
|
||||
1 2.2020000 1.0000000
|
||||
|
||||
|
@ -1,6 +1,6 @@
|
||||
subroutine CalcOmErf(maxm,ExpY,fG,NormYSq,Om)
|
||||
|
||||
! Compute the 0^m for the long-range Coulomb operator: (00|erf(r)/r|00)^m
|
||||
! Compute the 0^m for the long-range Coulomb operator: (00|erf(mu r)/r|00)^m
|
||||
|
||||
implicit none
|
||||
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
subroutine CCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! CCD module
|
||||
|
||||
@ -10,16 +10,22 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
|
||||
integer,intent(in) :: nBas,nEl
|
||||
integer,intent(in) :: nBasin
|
||||
integer,intent(in) :: nCin
|
||||
integer,intent(in) :: nOin
|
||||
integer,intent(in) :: nVin
|
||||
integer,intent(in) :: nRin
|
||||
double precision,intent(in) :: ENuc,ERHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: eHF(nBasin)
|
||||
double precision,intent(in) :: ERI(nBasin,nBasin,nBasin,nBasin)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: nBas2
|
||||
integer :: nBas
|
||||
integer :: nC
|
||||
integer :: nO
|
||||
integer :: nV
|
||||
integer :: nR
|
||||
integer :: nSCF
|
||||
double precision :: Conv
|
||||
double precision :: EcMP2,EcMP3,EcMP4
|
||||
@ -63,35 +69,34 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Spatial to spin orbitals
|
||||
|
||||
nBas2 = 2*nBas
|
||||
nBas = 2*nBasin
|
||||
nC = 2*nCin
|
||||
nO = 2*nOin
|
||||
nV = 2*nVin
|
||||
nR = 2*nRin
|
||||
|
||||
allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
|
||||
allocate(seHF(nBas),sERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
|
||||
call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI)
|
||||
call spatial_to_spin_MO_energy(nBasin,eHF,nBas,seHF)
|
||||
call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI)
|
||||
|
||||
! Antysymmetrize ERIs
|
||||
|
||||
allocate(dbERI(nBas2,nBas2,nBas2,nBas2))
|
||||
allocate(dbERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call antisymmetrize_ERI(2,nBas2,sERI,dbERI)
|
||||
call antisymmetrize_ERI(2,nBas,sERI,dbERI)
|
||||
|
||||
deallocate(sERI)
|
||||
|
||||
! Define occupied and virtual spaces
|
||||
|
||||
nO = 2*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)
|
||||
eV(:) = seHF(nO+1:nBas)
|
||||
|
||||
call form_delta_OOVV(nO,nV,eO,eV,delta_OOVV)
|
||||
call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV)
|
||||
|
||||
deallocate(seHF)
|
||||
|
||||
@ -99,10 +104,10 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
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)
|
||||
OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO )
|
||||
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
|
||||
OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas, 1:nO ,nO+1:nBas)
|
||||
VVVV(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas,nO+1:nBas,nO+1:nBas)
|
||||
|
||||
deallocate(dbERI)
|
||||
|
||||
@ -112,7 +117,7 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.))
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcMP2)
|
||||
EcMP4 = 0d0
|
||||
|
||||
! Memory allocation for DIIS
|
||||
@ -150,15 +155,15 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Form linear array
|
||||
|
||||
call form_u(nO,nV,OOOO,VVVV,OVOV,t2,u)
|
||||
call form_u(nC,nO,nV,nR,OOOO,VVVV,OVOV,t2,u)
|
||||
|
||||
! Form interemediate arrays
|
||||
|
||||
call form_X(nO,nV,OOVV,t2,X1,X2,X3,X4)
|
||||
call form_X(nC,nO,nV,nR,OOVV,t2,X1,X2,X3,X4)
|
||||
|
||||
! Form quadratic array
|
||||
|
||||
call form_v(nO,nV,X1,X2,X3,X4,t2,v)
|
||||
call form_v(nC,nO,nV,nR,X1,X2,X3,X4,t2,v)
|
||||
|
||||
! Compute residual
|
||||
|
||||
@ -166,7 +171,7 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Check convergence
|
||||
|
||||
Conv = maxval(abs(r2(:,:,:,:)))
|
||||
Conv = maxval(abs(r2(nC+1:nO,nC+1:nO,1:nV-nR,1:nV-nR)))
|
||||
|
||||
! Update amplitudes
|
||||
|
||||
@ -174,9 +179,9 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Compute correlation energy
|
||||
|
||||
EcCCD = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.))
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCCD)
|
||||
|
||||
if(nSCF == 1) EcMP3 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2 + v/delta_OOVV,.true.))
|
||||
if(nSCF == 1) call MP3_correlation_energy(nC,nO,nV,nR,OOVV,t2,v,delta_OOVV,EcMP3)
|
||||
|
||||
! Dump results
|
||||
|
||||
|
36
src/QuAcK/CCD_correlation_energy.f90
Normal file
36
src/QuAcK/CCD_correlation_energy.f90
Normal file
@ -0,0 +1,36 @@
|
||||
subroutine CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCCD)
|
||||
|
||||
! Compute the CCD energy
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
|
||||
double precision,intent(in) :: t2(nO,nO,nV,nV)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,j,a,b
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: EcCCD
|
||||
|
||||
EcCCD = 0d0
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
|
||||
EcCCD = EcCCD + OOVV(i,j,a,b)*t2(i,j,a,b)
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
EcCCD = 0.25d0*EcCCD
|
||||
|
||||
end subroutine CCD_correlation_energy
|
@ -1,4 +1,4 @@
|
||||
subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! CCSD module
|
||||
|
||||
@ -11,17 +11,24 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
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)
|
||||
integer,intent(in) :: nBasin
|
||||
integer,intent(in) :: nCin
|
||||
integer,intent(in) :: nOin
|
||||
integer,intent(in) :: nVin
|
||||
integer,intent(in) :: nRin
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: ERHF
|
||||
double precision,intent(in) :: eHF(nBasin)
|
||||
double precision,intent(in) :: ERI(nBasin,nBasin,nBasin,nBasin)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: start_CCSDT,end_CCSDT,t_CCSDT
|
||||
integer :: nBas2
|
||||
integer :: nBas
|
||||
integer :: nC
|
||||
integer :: nO
|
||||
integer :: nV
|
||||
integer :: nR
|
||||
integer :: nSCF
|
||||
double precision :: Conv
|
||||
double precision :: EcMP2
|
||||
@ -81,36 +88,35 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Spatial to spin orbitals
|
||||
|
||||
nBas2 = 2*nBas
|
||||
nBas = 2*nBasin
|
||||
nC = 2*nCin
|
||||
nO = 2*nOin
|
||||
nV = 2*nVin
|
||||
nR = 2*nRin
|
||||
|
||||
allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
|
||||
allocate(seHF(nBas),sERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
|
||||
call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI)
|
||||
call spatial_to_spin_MO_energy(nBasin,eHF,nBas,seHF)
|
||||
call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI)
|
||||
|
||||
! Antysymmetrize ERIs
|
||||
|
||||
allocate(dbERI(nBas2,nBas2,nBas2,nBas2))
|
||||
allocate(dbERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call antisymmetrize_ERI(2,nBas2,sERI,dbERI)
|
||||
call antisymmetrize_ERI(2,nBas,sERI,dbERI)
|
||||
|
||||
deallocate(sERI)
|
||||
|
||||
! Define occupied and virtual spaces
|
||||
|
||||
nO = 2*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)
|
||||
eV(:) = seHF(nO+1:nBas)
|
||||
|
||||
call form_delta_OV(nO,nV,eO,eV,delta_OV)
|
||||
call form_delta_OOVV(nO,nV,eO,eV,delta_OOVV)
|
||||
call form_delta_OV(nC,nO,nV,nR,eO,eV,delta_OV)
|
||||
call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV)
|
||||
|
||||
deallocate(seHF)
|
||||
|
||||
@ -122,16 +128,16 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
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)
|
||||
OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO )
|
||||
OOOV(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO ,nO+1:nBas)
|
||||
OVOO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas, 1:nO , 1:nO )
|
||||
VOOO(:,:,:,:) = dbERI(nO+1:nBas, 1:nO , 1:nO , 1:nO )
|
||||
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
|
||||
OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO )
|
||||
OVVV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas,nO+1:nBas)
|
||||
VOVV(:,:,:,:) = dbERI(nO+1:nBas, 1:nO ,nO+1:nBas,nO+1:nBas)
|
||||
VVVO(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas,nO+1:nBas, 1:nO )
|
||||
VVVV(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas,nO+1:nBas,nO+1:nBas)
|
||||
|
||||
deallocate(dbERI)
|
||||
|
||||
@ -141,7 +147,7 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
t1(:,:) = 0d0
|
||||
t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
call form_tau(nO,nV,t1,t2,tau)
|
||||
call form_tau(nC,nO,nV,nR,t1,t2,tau)
|
||||
|
||||
EcMP2 = 0.5d0*dot_product(pack(OOVV,.true.),pack(tau,.true.))
|
||||
write(*,'(1X,A20,1X,F10.6)') 'Ec(MP2) = ',EcMP2
|
||||
@ -186,34 +192,34 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Scuseria Eqs. (5), (6) and (7)
|
||||
|
||||
call form_h(nO,nV,eO,eV,OOVV,t1,tau,hvv,hoo,hvo)
|
||||
call form_h(nC,nO,nV,nR,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_g(nC,nO,nV,nR,hvv,hoo,VOVV,OOOV,t1,gvv,goo)
|
||||
|
||||
call form_abh(nO,nV,OOOO,OVOO,OOVV,VVVV,VOVV,OVVO,OVVV,t1,tau,aoooo,bvvvv,hovvo)
|
||||
call form_abh(nC,nO,nV,nR,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_r1(nC,nO,nV,nR,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)
|
||||
call form_r2(nC,nO,nV,nR,OOVV,OVOO,OVVV,OVVO,gvv,goo,aoooo,bvvvv,hovvo,t1,t2,tau,r2)
|
||||
|
||||
! Check convergence
|
||||
|
||||
Conv = max(maxval(abs(r1(:,:))),maxval(abs(r2(:,:,:,:))))
|
||||
Conv = max(maxval(abs(r1(nC+1:nO,1:nV-nR))),maxval(abs(r2(nC+1:nO,nC+1:nO,1:nV-nR,1:nV-nR))))
|
||||
|
||||
! Update
|
||||
|
||||
t1(:,:) = t1(:,:) - r1(:,:) /delta_OV (:,:)
|
||||
t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
call form_tau(nO,nV,t1,t2,tau)
|
||||
call form_tau(nC,nO,nV,nR,t1,t2,tau)
|
||||
|
||||
! Compute correlation energy
|
||||
|
||||
EcCCSD = 0.5d0*dot_product(pack(OOVV,.true.),pack(tau,.true.))
|
||||
call CCSD_correlation_energy(nC,nO,nV,nR,OOVV,tau,EcCCSD)
|
||||
|
||||
! Dump results
|
||||
|
||||
@ -222,12 +228,12 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
! DIIS extrapolation
|
||||
|
||||
n_diis = min(n_diis+1,max_diis)
|
||||
call DIIS_extrapolation(rcond1,nO*nV ,nO*nV ,n_diis,err1_diis,t1_diis,-r1/delta_OV ,t1)
|
||||
call DIIS_extrapolation(rcond2,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,err2_diis,t2_diis,-r2/delta_OOVV,t2)
|
||||
! call DIIS_extrapolation(rcond1,nO*nV ,nO*nV ,n_diis,err1_diis,t1_diis,-r1/delta_OV ,t1)
|
||||
! call DIIS_extrapolation(rcond2,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,err2_diis,t2_diis,-r2/delta_OOVV,t2)
|
||||
|
||||
! Reset DIIS if required
|
||||
|
||||
if(min(abs(rcond1),abs(rcond2)) < 1d-15) n_diis = 0
|
||||
! if(min(abs(rcond1),abs(rcond2)) < 1d-15) n_diis = 0
|
||||
|
||||
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,'|'
|
||||
@ -276,7 +282,7 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
if(doCCSDT) then
|
||||
|
||||
call cpu_time(start_CCSDT)
|
||||
call CCSDT(nO,nV,eO,eV,OOVV,VVVO,VOOO,t1,t2,EcCCT)
|
||||
call CCSDT(nC,nO,nV,nR,eO,eV,OOVV,VVVO,VOOO,t1,t2,EcCCT)
|
||||
call cpu_time(end_CCSDT)
|
||||
|
||||
t_CCSDT = end_CCSDT - start_CCSDT
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine CCSDT(nO,nV,eO,eV,OOVV,VVVO,VOOO,t1,t2,EcCCT)
|
||||
subroutine CCSDT(nC,nO,nV,nR,eO,eV,OOVV,VVVO,VOOO,t1,t2,EcCCT)
|
||||
|
||||
! Compute the (T) correction of the CCSD(T) energy
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine CCSDT(nO,nV,eO,eV,OOVV,VVVO,VOOO,t1,t2,EcCCT)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
|
||||
double precision,intent(in) :: eO(nO)
|
||||
double precision,intent(in) :: eV(nV)
|
||||
@ -34,12 +34,12 @@ subroutine CCSDT(nO,nV,eO,eV,OOVV,VVVO,VOOO,t1,t2,EcCCT)
|
||||
|
||||
! Form CCSD(T) quantities
|
||||
|
||||
call form_delta_OOOVVV(nO,nV,eO,eV,delta_OOOVVV)
|
||||
call form_delta_OOOVVV(nC,nO,nV,nR,eO,eV,delta_OOOVVV)
|
||||
|
||||
call form_ub(nO,nV,OOVV,t1,ub)
|
||||
call form_ub(nC,nO,nV,nR,OOVV,t1,ub)
|
||||
|
||||
call form_ubb(nO,nV,VVVO,VOOO,t2,ubb)
|
||||
call form_ubb(nC,nO,nV,nR,VVVO,VOOO,t2,ubb)
|
||||
|
||||
call form_T(nO,nV,delta_OOOVVV,ub,ubb,EcCCT)
|
||||
call form_T(nC,nO,nV,nR,delta_OOOVVV,ub,ubb,EcCCT)
|
||||
|
||||
end subroutine CCSDT
|
||||
|
36
src/QuAcK/CCSD_correlation_energy.f90
Normal file
36
src/QuAcK/CCSD_correlation_energy.f90
Normal file
@ -0,0 +1,36 @@
|
||||
subroutine CCSD_correlation_energy(nC,nO,nV,nR,OOVV,tau,EcCCSD)
|
||||
|
||||
! Compute the CCSD energy
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
|
||||
double precision,intent(in) :: tau(nO,nO,nV,nV)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,j,a,b
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: EcCCSD
|
||||
|
||||
EcCCSD = 0d0
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
|
||||
EcCCSD = EcCCSD + OOVV(i,j,a,b)*tau(i,j,a,b)
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
EcCCSD = 0.5d0*EcCCSD
|
||||
|
||||
end subroutine CCSD_correlation_energy
|
@ -131,6 +131,9 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m
|
||||
! Compute T-matrix version of the self-energy
|
||||
!----------------------------------------------
|
||||
|
||||
rho2s(:,:,:) = 0d0
|
||||
rho2t(:,:,:) = 0d0
|
||||
|
||||
call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF(:), &
|
||||
Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:), &
|
||||
Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:), &
|
||||
|
38
src/QuAcK/MP3_correlation_energy.f90
Normal file
38
src/QuAcK/MP3_correlation_energy.f90
Normal file
@ -0,0 +1,38 @@
|
||||
subroutine MP3_correlation_energy(nC,nO,nV,nR,OOVV,t2,v,delta_OOVV,EcMP3)
|
||||
|
||||
! Compute the MP3 correlation energy
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
|
||||
double precision,intent(in) :: t2(nO,nO,nV,nV)
|
||||
double precision,intent(in) :: v(nO,nO,nV,nV)
|
||||
double precision,intent(in) :: delta_OOVV(nO,nO,nV,nV)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,j,a,b
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: EcMP3
|
||||
|
||||
EcMP3 = 0d0
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
|
||||
EcMP3 = EcMP3 + OOVV(i,j,a,b)*(t2(i,j,a,b) + v(i,j,a,b)/delta_OOVV(i,j,a,b))
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
EcMP3 = 0.25d0*EcMP3
|
||||
|
||||
end subroutine MP3_correlation_energy
|
@ -375,7 +375,8 @@ program QuAcK
|
||||
if(doCCD) then
|
||||
|
||||
call cpu_time(start_CCD)
|
||||
call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), &
|
||||
ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
call cpu_time(end_CCD)
|
||||
|
||||
t_CCD = end_CCD - start_CCD
|
||||
@ -393,7 +394,8 @@ program QuAcK
|
||||
if(doCCSD) then
|
||||
|
||||
call cpu_time(start_CCSD)
|
||||
call CCSD(maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
call CCSD(maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC(1),nO(1),nV(1),nR(1), &
|
||||
ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
call cpu_time(end_CCSD)
|
||||
|
||||
t_CCSD = end_CCSD - start_CCSD
|
||||
@ -409,7 +411,8 @@ program QuAcK
|
||||
if(do_drCCD) then
|
||||
|
||||
call cpu_time(start_CCD)
|
||||
call drCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
call drCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), &
|
||||
ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
call cpu_time(end_CCD)
|
||||
|
||||
t_CCD = end_CCD - start_CCD
|
||||
@ -425,7 +428,8 @@ program QuAcK
|
||||
if(do_rCCD) then
|
||||
|
||||
call cpu_time(start_CCD)
|
||||
call rCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
call rCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), &
|
||||
ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
call cpu_time(end_CCD)
|
||||
|
||||
t_CCD = end_CCD - start_CCD
|
||||
@ -441,7 +445,8 @@ program QuAcK
|
||||
if(do_lCCD) then
|
||||
|
||||
call cpu_time(start_CCD)
|
||||
call lCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
call lCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), &
|
||||
ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
call cpu_time(end_CCD)
|
||||
|
||||
t_CCD = end_CCD - start_CCD
|
||||
@ -457,7 +462,8 @@ program QuAcK
|
||||
if(do_pCCD) then
|
||||
|
||||
call cpu_time(start_CCD)
|
||||
call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), &
|
||||
ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
call cpu_time(end_CCD)
|
||||
|
||||
t_CCD = end_CCD - start_CCD
|
||||
|
@ -61,7 +61,7 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
|
||||
|
||||
call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
|
||||
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
|
||||
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
|
||||
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
|
||||
|
||||
endif
|
||||
|
||||
@ -73,7 +73,7 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
|
||||
|
||||
call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
|
||||
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
|
||||
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
|
||||
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
|
||||
|
||||
endif
|
||||
|
||||
|
@ -61,7 +61,7 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
|
||||
|
||||
call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
|
||||
EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
|
||||
call print_excitation('RPAx ',ispin,nS,Omega(:,ispin))
|
||||
call print_excitation('RPAx ',ispin,nS,Omega(:,ispin))
|
||||
|
||||
endif
|
||||
|
||||
@ -73,7 +73,7 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
|
||||
|
||||
call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
|
||||
EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
|
||||
call print_excitation('RPAx ',ispin,nS,Omega(:,ispin))
|
||||
call print_excitation('RPAx ',ispin,nS,Omega(:,ispin))
|
||||
|
||||
endif
|
||||
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine drCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Direct ring CCD module
|
||||
|
||||
@ -10,19 +10,25 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
|
||||
integer,intent(in) :: nBas,nEl
|
||||
integer,intent(in) :: nBasin
|
||||
integer,intent(in) :: nCin
|
||||
integer,intent(in) :: nOin
|
||||
integer,intent(in) :: nVin
|
||||
integer,intent(in) :: nRin
|
||||
double precision,intent(in) :: ENuc,ERHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: eHF(nBasin)
|
||||
double precision,intent(in) :: ERI(nBasin,nBasin,nBasin,nBasin)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: nBas2
|
||||
integer :: nBas
|
||||
integer :: nC
|
||||
integer :: nO
|
||||
integer :: nV
|
||||
integer :: nR
|
||||
integer :: nSCF
|
||||
double precision :: Conv
|
||||
double precision :: EcMP2,EcMP3,EcMP4
|
||||
double precision :: EcMP2
|
||||
double precision :: ECCD,EcCCD
|
||||
double precision,allocatable :: seHF(:)
|
||||
double precision,allocatable :: sERI(:,:,:,:)
|
||||
@ -31,18 +37,8 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
double precision,allocatable :: eV(:)
|
||||
double precision,allocatable :: delta_OOVV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: OOOO(:,:,:,:)
|
||||
double precision,allocatable :: OOVV(:,:,:,:)
|
||||
double precision,allocatable :: OVVO(:,:,:,:)
|
||||
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(:,:,:,:)
|
||||
@ -63,19 +59,16 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Spatial to spin orbitals
|
||||
|
||||
nBas2 = 2*nBas
|
||||
nBas = 2*nBasin
|
||||
nC = 2*nCin
|
||||
nO = 2*nOin
|
||||
nV = 2*nVin
|
||||
nR = 2*nRin
|
||||
|
||||
allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
|
||||
allocate(seHF(nBas),sERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
|
||||
call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI)
|
||||
|
||||
! Antysymmetrize ERIs
|
||||
|
||||
! Define occupied and virtual spaces
|
||||
|
||||
nO = 2*nEl
|
||||
nV = nBas2 - nO
|
||||
call spatial_to_spin_MO_energy(nBasin,eHF,nBas,seHF)
|
||||
call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI)
|
||||
|
||||
! Form energy denominator
|
||||
|
||||
@ -83,20 +76,18 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
allocate(delta_OOVV(nO,nO,nV,nV))
|
||||
|
||||
eO(:) = seHF(1:nO)
|
||||
eV(:) = seHF(nO+1:nBas2)
|
||||
eV(:) = seHF(nO+1:nBas)
|
||||
|
||||
call form_delta_OOVV(nO,nV,eO,eV,delta_OOVV)
|
||||
call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV)
|
||||
|
||||
deallocate(seHF)
|
||||
|
||||
! Create integral batches
|
||||
|
||||
allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VVVV(nV,nV,nV,nV))
|
||||
allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO))
|
||||
|
||||
OOOO(:,:,:,:) = sERI( 1:nO , 1:nO , 1:nO , 1:nO )
|
||||
OOVV(:,:,:,:) = sERI( 1:nO , 1:nO ,nO+1:nBas2,nO+1:nBas2)
|
||||
OVVO(:,:,:,:) = sERI( 1:nO ,nO+1:nBas2,nO+1:nBas2, 1:nO )
|
||||
VVVV(:,:,:,:) = sERI(nO+1:nBas2,nO+1:nBas2,nO+1:nBas2,nO+1:nBas2)
|
||||
OOVV(:,:,:,:) = sERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
|
||||
OVVO(:,:,:,:) = sERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO )
|
||||
|
||||
deallocate(sERI)
|
||||
|
||||
@ -106,8 +97,7 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.))
|
||||
EcMP4 = 0d0
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcMP2)
|
||||
|
||||
! Memory allocation for DIIS
|
||||
|
||||
@ -115,8 +105,7 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! 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))
|
||||
allocate(r2(nO,nO,nV,nV))
|
||||
|
||||
Conv = 1d0
|
||||
nSCF = 0
|
||||
@ -144,23 +133,22 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Compute residual
|
||||
|
||||
call form_ring_r(nO,nV,OVVO,OOVV,t2,r2)
|
||||
call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
|
||||
|
||||
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)
|
||||
|
||||
! Check convergence
|
||||
|
||||
Conv = maxval(abs(r2(:,:,:,:)))
|
||||
|
||||
Conv = maxval(abs(r2(nC+1:nO,nC+1:nO,1:nV-nR,1:nV-nR)))
|
||||
|
||||
! Update amplitudes
|
||||
|
||||
t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
! Compute correlation energy
|
||||
|
||||
EcCCD = 0.5d0*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.))
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCCD)
|
||||
EcCCD = 2d0*EcCCD
|
||||
|
||||
! Dump results
|
||||
|
||||
|
@ -45,7 +45,10 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
|
||||
do c=nO+1,nBas-nR
|
||||
do d=c,nBas-nR
|
||||
cd = cd + 1
|
||||
rho1(p,i,ab) = rho1(p,i,ab) + ERI(p,i,c,d)*X1(cd,ab)
|
||||
rho1(p,i,ab) = rho1(p,i,ab) + ERI(p,i,c,d)*X1(cd,ab)
|
||||
! rho1(p,i,ab) = rho1(p,i,ab) &
|
||||
! + (ERI(p,i,c,d) + ERI(p,i,d,c))*X1(cd,ab) &
|
||||
! /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d)))
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -54,6 +57,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
|
||||
do l=k,nO
|
||||
kl = kl + 1
|
||||
rho1(p,i,ab) = rho1(p,i,ab) + ERI(p,i,k,l)*Y1(kl,ab)
|
||||
! rho1(p,i,ab) = rho1(p,i,ab) &
|
||||
! + (ERI(p,i,k,l) + ERI(p,i,l,k))*Y1(kl,ab) &
|
||||
! /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l)))
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -68,6 +74,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
|
||||
do d=c,nBas-nR
|
||||
cd = cd + 1
|
||||
rho2(p,a,ij) = rho2(p,a,ij) + ERI(p,nO+a,c,d)*X2(cd,ij)
|
||||
! rho2(p,a,ij) = rho2(p,a,ij) &
|
||||
! + (ERI(p,nO+a,c,d) + ERI(p,nO+a,d,c))*X2(cd,ij) &
|
||||
! /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(c,d)))
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -76,6 +85,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
|
||||
do l=k,nO
|
||||
kl = kl + 1
|
||||
rho2(p,a,ij) = rho2(p,a,ij) + ERI(p,nO+a,k,l)*Y2(kl,ij)
|
||||
! rho2(p,a,ij) = rho2(p,a,ij) &
|
||||
! + (ERI(p,nO+a,k,l) + ERI(p,nO+a,l,k))*Y2(kl,ij) &
|
||||
! /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(k,l)))
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -101,7 +113,8 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
|
||||
do c=nO+1,nBas-nR
|
||||
do d=c+1,nBas-nR
|
||||
cd = cd + 1
|
||||
rho1(p,i,ab) = rho1(p,i,ab) + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab)
|
||||
rho1(p,i,ab) = rho1(p,i,ab) + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab)
|
||||
! *3d0/2d0*sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d)))
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -110,6 +123,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
|
||||
do l=k+1,nO
|
||||
kl = kl + 1
|
||||
rho1(p,i,ab) = rho1(p,i,ab) + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab)
|
||||
! *3d0/2d0/sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l)))
|
||||
end do
|
||||
end do
|
||||
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_T(nO,nV,delta_OOOVVV,ub,ubb,EcCCT)
|
||||
subroutine form_T(nC,nO,nV,nR,delta_OOOVVV,ub,ubb,EcCCT)
|
||||
|
||||
! Compute (T) correction
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_T(nO,nV,delta_OOOVVV,ub,ubb,EcCCT)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
|
||||
double precision,intent(in) :: delta_OOOVVV(nO,nO,nO,nV,nV,nV)
|
||||
double precision,intent(in) :: ub(nO,nO,nO,nV,nV,nV)
|
||||
@ -23,12 +23,12 @@ subroutine form_T(nO,nV,delta_OOOVVV,ub,ubb,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
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
do c=1,nV-nR
|
||||
|
||||
EcCCT = EcCCT &
|
||||
+ (ub(i,j,k,a,b,c) + ubb(i,j,k,a,b,c)) &
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_X(nO,nV,OOVV,t2,X1,X2,X3,X4)
|
||||
subroutine form_X(nC,nO,nV,nR,OOVV,t2,X1,X2,X3,X4)
|
||||
|
||||
! Form intermediate arrays X's in CCD
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_X(nO,nV,OOVV,t2,X1,X2,X3,X4)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
double precision,intent(in) :: t2(nO,nO,nV,nV)
|
||||
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
|
||||
|
||||
@ -31,12 +31,12 @@ subroutine form_X(nO,nV,OOVV,t2,X1,X2,X3,X4)
|
||||
|
||||
! 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
|
||||
do k=nC+1,nO
|
||||
do l=nC+1,nO
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
do d=1,nV-nR
|
||||
X1(k,l,i,j) = X1(k,l,i,j) + OOVV(k,l,c,d)*t2(i,j,c,d)
|
||||
enddo
|
||||
enddo
|
||||
@ -47,11 +47,11 @@ subroutine form_X(nO,nV,OOVV,t2,X1,X2,X3,X4)
|
||||
|
||||
! Build X2
|
||||
|
||||
do b=1,nV
|
||||
do c=1,nV
|
||||
do k=1,nO
|
||||
do l=1,nO
|
||||
do d=1,nV
|
||||
do b=1,nV-nR
|
||||
do c=1,nV-nR
|
||||
do k=nC+1,nO
|
||||
do l=nC+1,nO
|
||||
do d=1,nV-nR
|
||||
X2(b,c) = X2(b,c) + OOVV(k,l,c,d)*t2(k,l,b,d)
|
||||
enddo
|
||||
enddo
|
||||
@ -61,11 +61,11 @@ subroutine form_X(nO,nV,OOVV,t2,X1,X2,X3,X4)
|
||||
|
||||
! Build X3
|
||||
|
||||
do k=1,nO
|
||||
do j=1,nO
|
||||
do l=1,nO
|
||||
do c=1,nV
|
||||
do d=1,nV
|
||||
do k=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do l=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
do d=1,nV-nR
|
||||
X3(k,j) = X3(k,j) + OOVV(k,l,c,d)*t2(j,l,c,d)
|
||||
enddo
|
||||
enddo
|
||||
@ -75,12 +75,12 @@ subroutine form_X(nO,nV,OOVV,t2,X1,X2,X3,X4)
|
||||
|
||||
! 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
|
||||
do i=nC+1,nO
|
||||
do l=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do d=1,nV-nR
|
||||
do k=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
X4(i,l,a,d) = X4(i,l,a,d) + OOVV(k,l,c,d)*t2(i,k,a,c)
|
||||
enddo
|
||||
enddo
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_abh(nO,nV,OOOO,OVOO,OOVV,VVVV,VOVV,OVVO,OVVV,t1,tau,aoooo,bvvvv,hovvo)
|
||||
subroutine form_abh(nC,nO,nV,nR,OOOO,OVOO,OOVV,VVVV,VOVV,OVVO,OVVV,t1,tau,aoooo,bvvvv,hovvo)
|
||||
|
||||
! Scuseria Eqs. (11),(12) and (13)
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_abh(nO,nV,OOOO,OVOO,OOVV,VVVV,VOVV,OVVO,OVVV,t1,tau,aoooo,bvvvv,
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
|
||||
double precision,intent(in) :: OOOO(nO,nO,nO,nO)
|
||||
double precision,intent(in) :: OVOO(nO,nV,nO,nO)
|
||||
@ -32,21 +32,21 @@ subroutine form_abh(nO,nV,OOOO,OVOO,OOVV,VVVV,VOVV,OVVO,OVVV,t1,tau,aoooo,bvvvv,
|
||||
|
||||
aoooo(:,:,:,:) = OOOO(:,:,:,:)
|
||||
|
||||
do i=1,nO
|
||||
do j=1,nO
|
||||
do k=1,nO
|
||||
do l=1,nO
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do l=nC+1,nO
|
||||
|
||||
do c=1,nV
|
||||
do c=1,nV-nR
|
||||
aoooo(i,j,k,l) = aoooo(i,j,k,l) + OVOO(i,c,k,l)*t1(j,c)
|
||||
end do
|
||||
|
||||
do c=1,nV
|
||||
do c=1,nV-nR
|
||||
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
|
||||
do c=1,nV-nR
|
||||
do d=1,nV-nR
|
||||
aoooo(i,j,k,l) = aoooo(i,j,k,l) + OOVV(k,l,c,d)*tau(i,j,c,d)
|
||||
end do
|
||||
end do
|
||||
@ -58,16 +58,16 @@ subroutine form_abh(nO,nV,OOOO,OVOO,OOVV,VVVV,VOVV,OVVO,OVVV,t1,tau,aoooo,bvvvv,
|
||||
|
||||
bvvvv(:,:,:,:) = VVVV(:,:,:,:)
|
||||
|
||||
do c=1,nV
|
||||
do d=1,nV
|
||||
do a=1,nV
|
||||
do b=1,nV
|
||||
do c=1,nV-nR
|
||||
do d=1,nV-nR
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
|
||||
do k=1,nO
|
||||
do k=nC+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
|
||||
do k=nC+1,nO
|
||||
bvvvv(c,d,a,b) = bvvvv(c,d,a,b) + VOVV(b,k,c,d)*t1(k,a)
|
||||
end do
|
||||
|
||||
@ -78,21 +78,21 @@ subroutine form_abh(nO,nV,OOOO,OVOO,OOVV,VVVV,VOVV,OVVO,OVVV,t1,tau,aoooo,bvvvv,
|
||||
|
||||
hovvo(:,:,:,:) = OVVO(:,:,:,:)
|
||||
|
||||
do i=1,nO
|
||||
do c=1,nV
|
||||
do a=1,nV
|
||||
do k=1,nO
|
||||
do i=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
do a=1,nV-nR
|
||||
do k=nC+1,nO
|
||||
|
||||
do l=1,nO
|
||||
do l=nC+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
|
||||
do d=1,nV-nR
|
||||
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
|
||||
do l=nC+1,nO
|
||||
do d=1,nV-nR
|
||||
hovvo(i,c,a,k) = hovvo(i,c,a,k) - OOVV(k,l,c,d)*tau(i,l,d,a)
|
||||
end do
|
||||
end do
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_delta_OOOVVV(nO,nV,eO,eV,delta)
|
||||
subroutine form_delta_OOOVVV(nC,nO,nV,nR,eO,eV,delta)
|
||||
|
||||
! Form energy denominator for CC
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_delta_OOOVVV(nO,nV,eO,eV,delta)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
double precision,intent(in) :: eO(nO)
|
||||
double precision,intent(in) :: eV(nV)
|
||||
|
||||
@ -18,12 +18,12 @@ subroutine form_delta_OOOVVV(nO,nV,eO,eV,delta)
|
||||
|
||||
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
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
do c=1,nV-nR
|
||||
|
||||
delta(i,j,k,a,b,c) = eV(a) + eV(b) + eV(c) - eO(i) - eO(j) - eO(k)
|
||||
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_delta_OOVV(nO,nV,eO,eV,delta)
|
||||
subroutine form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta)
|
||||
|
||||
! Form energy denominator for CC
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_delta_OOVV(nO,nV,eO,eV,delta)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
double precision,intent(in) :: eO(nO)
|
||||
double precision,intent(in) :: eV(nV)
|
||||
|
||||
@ -18,10 +18,10 @@ subroutine form_delta_OOVV(nO,nV,eO,eV,delta)
|
||||
|
||||
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
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
|
||||
delta(i,j,a,b) = eV(a) + eV(b) - eO(i) - eO(j)
|
||||
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_delta_OV(nO,nV,eO,eV,delta)
|
||||
subroutine form_delta_OV(nC,nO,nV,nR,eO,eV,delta)
|
||||
|
||||
! Form energy denominator for CC
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_delta_OV(nO,nV,eO,eV,delta)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
double precision,intent(in) :: eO(nO)
|
||||
double precision,intent(in) :: eV(nV)
|
||||
|
||||
@ -18,8 +18,8 @@ subroutine form_delta_OV(nO,nV,eO,eV,delta)
|
||||
|
||||
double precision,intent(out) :: delta(nO,nV)
|
||||
|
||||
do i=1,nO
|
||||
do a=1,nV
|
||||
do i=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
delta(i,a) = eV(a) - eO(i)
|
||||
enddo
|
||||
enddo
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_g(nO,nV,hvv,hoo,VOVV,OOOV,t1,gvv,goo)
|
||||
subroutine form_g(nC,nO,nV,nR,hvv,hoo,VOVV,OOOV,t1,gvv,goo)
|
||||
|
||||
! Scuseria Eqs. (9), (10)
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_g(nO,nV,hvv,hoo,VOVV,OOOV,t1,gvv,goo)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
|
||||
double precision,intent(in) :: hvv(nV,nV)
|
||||
double precision,intent(in) :: hoo(nO,nO)
|
||||
@ -28,10 +28,10 @@ subroutine form_g(nO,nV,hvv,hoo,VOVV,OOOV,t1,gvv,goo)
|
||||
|
||||
gvv(:,:) = hvv(:,:)
|
||||
|
||||
do c=1,nV
|
||||
do a=1,nV
|
||||
do k=1,nO
|
||||
do d=1,nV
|
||||
do c=1,nV-nR
|
||||
do a=1,nV-nR
|
||||
do k=nC+1,nO
|
||||
do d=1,nV-nR
|
||||
gvv(c,a) = gvv(c,a) + VOVV(a,k,c,d)*t1(k,d)
|
||||
end do
|
||||
end do
|
||||
@ -40,10 +40,10 @@ subroutine form_g(nO,nV,hvv,hoo,VOVV,OOOV,t1,gvv,goo)
|
||||
|
||||
goo(:,:) = hoo(:,:)
|
||||
|
||||
do i=1,nO
|
||||
do k=1,nO
|
||||
do l=1,nO
|
||||
do c=1,nV
|
||||
do i=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do l=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
goo(i,k) = goo(i,k) + OOOV(k,l,i,c)*t1(l,c)
|
||||
end do
|
||||
end do
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_h(nO,nV,eO,eV,OOVV,t1,tau,hvv,hoo,hvo)
|
||||
subroutine form_h(nC,nO,nV,nR,eO,eV,OOVV,t1,tau,hvv,hoo,hvo)
|
||||
|
||||
! Scuseria Eqs. (5), (6) and (7)
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_h(nO,nV,eO,eV,OOVV,t1,tau,hvv,hoo,hvo)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
|
||||
double precision,intent(in) :: eO(nO)
|
||||
double precision,intent(in) :: eV(nV)
|
||||
@ -28,12 +28,12 @@ subroutine form_h(nO,nV,eO,eV,OOVV,t1,tau,hvv,hoo,hvo)
|
||||
|
||||
hvv(:,:) = 0d0
|
||||
|
||||
do b=1,nV
|
||||
do b=1,nV-nR
|
||||
hvv(b,b) = eV(b)
|
||||
do a=1,nV
|
||||
do j=1,nO
|
||||
do k=1,nO
|
||||
do c=1,nV
|
||||
do a=1,nV-nR
|
||||
do j=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
|
||||
hvv(b,a) = hvv(b,a) - OOVV(j,k,b,c)*tau(j,k,a,c)
|
||||
|
||||
@ -45,12 +45,12 @@ subroutine form_h(nO,nV,eO,eV,OOVV,t1,tau,hvv,hoo,hvo)
|
||||
|
||||
hoo(:,:) = 0d0
|
||||
|
||||
do i=1,nO
|
||||
do i=nC+1,nO
|
||||
hoo(i,i) = eO(i)
|
||||
do j=1,nO
|
||||
do k=1,nO
|
||||
do b=1,nV
|
||||
do c=1,nV
|
||||
do j=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do b=1,nV-nR
|
||||
do c=1,nV-nR
|
||||
|
||||
hoo(i,j) = hoo(i,j) + OOVV(j,k,b,c)*tau(i,k,b,c)
|
||||
|
||||
@ -62,10 +62,10 @@ subroutine form_h(nO,nV,eO,eV,OOVV,t1,tau,hvv,hoo,hvo)
|
||||
|
||||
hvo(:,:) = 0d0
|
||||
|
||||
do b=1,nV
|
||||
do j=1,nO
|
||||
do k=1,nO
|
||||
do c=1,nV
|
||||
do b=1,nV-nR
|
||||
do j=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
|
||||
hvo(b,j) = hvo(b,j) + OOVV(j,k,b,c)*t1(k,c)
|
||||
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_ladder_r(nO,nV,OOOO,OOVV,VVVV,t2,r2)
|
||||
subroutine form_ladder_r(nC,nO,nV,nR,OOOO,OOVV,VVVV,t2,r2)
|
||||
|
||||
! Form residuals for ladder CCD
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_ladder_r(nO,nV,OOOO,OOVV,VVVV,t2,r2)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
double precision,intent(in) :: t2(nO,nO,nV,nV)
|
||||
double precision,intent(in) :: OOOO(nO,nO,nO,nO)
|
||||
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
|
||||
@ -16,6 +16,7 @@ subroutine form_ladder_r(nO,nV,OOOO,OOVV,VVVV,t2,r2)
|
||||
|
||||
integer :: i,j,k,l
|
||||
integer :: a,b,c,d
|
||||
double precision,allocatable :: y(:,:,:,:)
|
||||
|
||||
! Output variables
|
||||
|
||||
@ -23,31 +24,44 @@ subroutine form_ladder_r(nO,nV,OOOO,OOVV,VVVV,t2,r2)
|
||||
|
||||
r2(:,:,:,:) = 0d0
|
||||
|
||||
do i=1,nO
|
||||
do j=1,nO
|
||||
do a=1,nV
|
||||
do b=1,nV
|
||||
allocate(y(nO,nO,nO,nO))
|
||||
|
||||
do k=1,nO
|
||||
do l=1,nO
|
||||
y(:,:,:,:) = 0d0
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do l=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
do d=1,nV-nR
|
||||
y(i,j,k,l) = y(i,j,k,l) + t2(i,j,c,d)*OOVV(k,l,c,d)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
|
||||
do k=nC+1,nO
|
||||
do l=nC+1,nO
|
||||
r2(i,j,a,b) = r2(i,j,a,b) + 0.5d0*t2(k,l,a,b)*OOOO(i,j,k,l)
|
||||
end do
|
||||
end do
|
||||
|
||||
do c=1,nV
|
||||
do d=1,nV
|
||||
do c=1,nV-nR
|
||||
do d=1,nV-nR
|
||||
r2(i,j,a,b) = r2(i,j,a,b) + 0.5d0*VVVV(c,d,a,b)*t2(i,j,c,d)
|
||||
end do
|
||||
end do
|
||||
|
||||
do k=1,nO
|
||||
do l=1,nO
|
||||
do c=1,nV
|
||||
do d=1,nV
|
||||
r2(i,j,a,b) = r2(i,j,a,b) + 0.25d0*t2(i,j,c,d)*OOVV(k,l,c,d)*t2(k,l,a,b)
|
||||
|
||||
end do
|
||||
end do
|
||||
do k=nC+1,nO
|
||||
do l=nC+1,nO
|
||||
r2(i,j,a,b) = r2(i,j,a,b) + 0.25d0*y(i,j,k,l)*t2(k,l,a,b)
|
||||
end do
|
||||
end do
|
||||
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_r1(nO,nV,OVVO,OVVV,OOOV,hvv,hoo,hvo,t1,t2,tau,r1)
|
||||
subroutine form_r1(nC,nO,nV,nR,OVVO,OVVV,OOOV,hvv,hoo,hvo,t1,t2,tau,r1)
|
||||
|
||||
! Form tau in CCSD
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_r1(nO,nV,OVVO,OVVV,OOOV,hvv,hoo,hvo,t1,t2,tau,r1)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
|
||||
double precision,intent(in) :: OVVO(nO,nV,nV,nO)
|
||||
double precision,intent(in) :: OVVV(nO,nV,nV,nV)
|
||||
@ -32,40 +32,40 @@ subroutine form_r1(nO,nV,OVVO,OVVV,OOOV,hvv,hoo,hvo,t1,t2,tau,r1)
|
||||
|
||||
r1(:,:) = 0d0
|
||||
|
||||
do i=1,nO
|
||||
do a=1,nV
|
||||
do i=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
|
||||
do b=1,nV
|
||||
do b=1,nV-nR
|
||||
r1(i,a) = r1(i,a) + hvv(b,a)*t1(i,b)
|
||||
end do
|
||||
|
||||
do j=1,nO
|
||||
do j=nC+1,nO
|
||||
r1(i,a) = r1(i,a) - hoo(i,j)*t1(j,a)
|
||||
end do
|
||||
|
||||
do j=1,nO
|
||||
do b=1,nV
|
||||
do j=nC+1,nO
|
||||
do b=1,nV-nR
|
||||
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
|
||||
do j=nC+1,nO
|
||||
do b=1,nV-nR
|
||||
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
|
||||
do j=nC+1,nO
|
||||
do b=1,nV-nR
|
||||
do c=1,nV-nR
|
||||
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
|
||||
do j=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do b=1,nV-nR
|
||||
r1(i,a) = r1(i,a) - OOOV(j,k,i,b)*tau(j,k,a,b)
|
||||
end do
|
||||
end do
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_r2(nO,nV,OOVV,OVOO,OVVV,OVVO,gvv,goo,aoooo,bvvvv,hovvo,t1,t2,tau,r2)
|
||||
subroutine form_r2(nC,nO,nV,nR,OOVV,OVOO,OVVV,OVVO,gvv,goo,aoooo,bvvvv,hovvo,t1,t2,tau,r2)
|
||||
|
||||
! Form tau in CCSD
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_r2(nO,nV,OOVV,OVOO,OVVV,OVVO,gvv,goo,aoooo,bvvvv,hovvo,t1,t2,tau
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
|
||||
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
|
||||
double precision,intent(in) :: OVOO(nO,nV,nO,nO)
|
||||
@ -34,99 +34,99 @@ subroutine form_r2(nO,nV,OOVV,OVOO,OVVV,OVVO,gvv,goo,aoooo,bvvvv,hovvo,t1,t2,tau
|
||||
|
||||
r2(:,:,:,:) = OOVV(:,:,:,:)
|
||||
|
||||
do i=1,nO
|
||||
do j=1,nO
|
||||
do a=1,nV
|
||||
do b=1,nV
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
|
||||
do k=1,nO
|
||||
do l=1,nO
|
||||
do k=nC+1,nO
|
||||
do l=nC+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
|
||||
do c=1,nV-nR
|
||||
do d=1,nV-NR
|
||||
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
|
||||
do c=1,nV-nR
|
||||
r2(i,j,a,b) = r2(i,j,a,b) + gvv(c,a)*t2(i,j,c,b)
|
||||
end do
|
||||
|
||||
do k=1,nO
|
||||
do k=nC+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
|
||||
do c=1,nV-NR
|
||||
r2(i,j,a,b) = r2(i,j,a,b) - gvv(c,b)*t2(i,j,c,a)
|
||||
end do
|
||||
|
||||
do k=1,nO
|
||||
do k=nC+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
|
||||
do k=nC+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
|
||||
do c=1,nV-nR
|
||||
r2(i,j,a,b) = r2(i,j,a,b) + OVVV(j,c,b,a)*t1(i,c)
|
||||
end do
|
||||
|
||||
do k=1,nO
|
||||
do k=nC+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
|
||||
do c=1,nV-nR
|
||||
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
|
||||
do k=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
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
|
||||
do k=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
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
|
||||
do k=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
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
|
||||
do k=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
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
|
||||
do k=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
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
|
||||
do k=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
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
|
||||
do k=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
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
|
||||
do k=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
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
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_ring_r(nO,nV,OVVO,OOVV,t2,r2)
|
||||
subroutine form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
|
||||
|
||||
! Form residuals for ring CCD
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_ring_r(nO,nV,OVVO,OOVV,t2,r2)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
double precision,intent(in) :: t2(nO,nO,nV,nV)
|
||||
double precision,intent(in) :: OVVO(nO,nV,nV,nO)
|
||||
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
|
||||
@ -15,6 +15,7 @@ subroutine form_ring_r(nO,nV,OVVO,OOVV,t2,r2)
|
||||
|
||||
integer :: i,j,k,l
|
||||
integer :: a,b,c,d
|
||||
double precision,allocatable :: y(:,:,:,:)
|
||||
|
||||
! Output variables
|
||||
|
||||
@ -22,28 +23,42 @@ subroutine form_ring_r(nO,nV,OVVO,OOVV,t2,r2)
|
||||
|
||||
r2(:,:,:,:) = 0d0
|
||||
|
||||
do i=1,nO
|
||||
do j=1,nO
|
||||
do a=1,nV
|
||||
do b=1,nV
|
||||
allocate(y(nO,nV,nO,nV))
|
||||
|
||||
do k=1,nO
|
||||
do c=1,nV
|
||||
y(:,:,:,:) = 0d0
|
||||
|
||||
do i=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do l=nC+1,nO
|
||||
do d=1,nV-nR
|
||||
do k=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
y(i,a,l,d) = y(i,a,l,d) + t2(i,k,a,c)*OOVV(k,l,c,d)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
|
||||
do k=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
|
||||
r2(i,j,a,b) = r2(i,j,a,b) + OVVO(i,c,a,k)*t2(k,j,c,b) + OVVO(j,c,b,k)*t2(i,k,a,c)
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
do k=1,nO
|
||||
do l=1,nO
|
||||
do c=1,nV
|
||||
do d=1,nV
|
||||
do l=nC+1,nO
|
||||
do d=1,nV-nR
|
||||
|
||||
r2(i,j,a,b) = r2(i,j,a,b) + t2(i,k,a,c)*OOVV(k,l,c,d)*t2(l,j,d,b)
|
||||
r2(i,j,a,b) = r2(i,j,a,b) + y(i,a,l,d)*t2(l,j,d,b)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_tau(nO,nV,t1,t2,tau)
|
||||
subroutine form_tau(nC,nO,nV,nR,t1,t2,tau)
|
||||
|
||||
! Form tau in CCSD
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_tau(nO,nV,t1,t2,tau)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
double precision,intent(in) :: t1(nO,nV)
|
||||
double precision,intent(in) :: t2(nO,nO,nV,nV)
|
||||
|
||||
@ -19,10 +19,10 @@ subroutine form_tau(nO,nV,t1,t2,tau)
|
||||
|
||||
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
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
|
||||
tau(i,j,a,b) = 0.5d0*t2(i,j,a,b) + t1(i,a)*t1(j,b)
|
||||
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_u(nO,nV,OOOO,VVVV,OVOV,t2,u)
|
||||
subroutine form_u(nC,nO,nV,nR,OOOO,VVVV,OVOV,t2,u)
|
||||
|
||||
! Form linear array in CCD
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_u(nO,nV,OOOO,VVVV,OVOV,t2,u)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
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)
|
||||
@ -23,12 +23,12 @@ subroutine form_u(nO,nV,OOOO,VVVV,OVOV,t2,u)
|
||||
|
||||
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
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
do c=1,nV-nR
|
||||
do d=1,nV-nR
|
||||
u(i,j,a,b) = u(i,j,a,b) + 0.5d0*VVVV(a,b,c,d)*t2(i,j,c,d)
|
||||
enddo
|
||||
enddo
|
||||
@ -37,12 +37,12 @@ subroutine form_u(nO,nV,OOOO,VVVV,OVOV,t2,u)
|
||||
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
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do l=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
u(i,j,a,b) = u(i,j,a,b) + 0.5d0*OOOO(k,l,i,j)*t2(k,l,a,b)
|
||||
enddo
|
||||
enddo
|
||||
@ -51,12 +51,12 @@ subroutine form_u(nO,nV,OOOO,VVVV,OVOV,t2,u)
|
||||
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
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
do c=1,nV-nR
|
||||
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) &
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_ub(nO,nV,OOVV,t1,ub)
|
||||
subroutine form_ub(nC,nO,nV,nR,OOVV,t1,ub)
|
||||
|
||||
! Form 1st term in (T) correction
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_ub(nO,nV,OOVV,t1,ub)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
|
||||
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
|
||||
|
||||
@ -21,12 +21,12 @@ subroutine form_ub(nO,nV,OOVV,t1,ub)
|
||||
|
||||
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
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
do c=1,nV-nR
|
||||
|
||||
ub(i,j,k,a,b,c) = t1(i,a)*OOVV(j,k,b,c) &
|
||||
+ t1(i,b)*OOVV(j,k,c,a) &
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_ubb(nO,nV,VVVO,VOOO,t2,ubb)
|
||||
subroutine form_ubb(nC,nO,nV,nR,VVVO,VOOO,t2,ubb)
|
||||
|
||||
! Form 2nd term in (T) correction
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_ubb(nO,nV,VVVO,VOOO,t2,ubb)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
|
||||
double precision,intent(in) :: VVVO(nV,nV,nV,nO)
|
||||
double precision,intent(in) :: VOOO(nV,nO,nO,nO)
|
||||
@ -24,14 +24,14 @@ subroutine form_ubb(nO,nV,VVVO,VOOO,t2,ubb)
|
||||
|
||||
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 i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
do c=1,nV-nR
|
||||
|
||||
do e=1,nV
|
||||
do e=1,nV-nR
|
||||
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) &
|
||||
@ -44,7 +44,7 @@ subroutine form_ubb(nO,nV,VVVO,VOOO,t2,ubb)
|
||||
+ t2(j,k,c,e)*VVVO(a,b,e,i)
|
||||
end do
|
||||
|
||||
do m=1,nO
|
||||
do m=nC+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) &
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine form_v(nO,nV,X1,X2,X3,X4,t2,v)
|
||||
subroutine form_v(nC,nO,nV,nR,X1,X2,X3,X4,t2,v)
|
||||
|
||||
! Form quadratic array in CCD
|
||||
|
||||
@ -6,7 +6,7 @@ subroutine form_v(nO,nV,X1,X2,X3,X4,t2,v)
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
integer,intent(in) :: nC,nO,nV,nR
|
||||
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)
|
||||
@ -24,12 +24,12 @@ subroutine form_v(nO,nV,X1,X2,X3,X4,t2,v)
|
||||
|
||||
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
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
do k=nC+1,nO
|
||||
do l=nC+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
|
||||
@ -38,11 +38,11 @@ subroutine form_v(nO,nV,X1,X2,X3,X4,t2,v)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do i=1,nO
|
||||
do j=1,nO
|
||||
do a=1,nV
|
||||
do b=1,nV
|
||||
do c=1,nV
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
do c=1,nV-nR
|
||||
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
|
||||
@ -50,11 +50,11 @@ subroutine form_v(nO,nV,X1,X2,X3,X4,t2,v)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do i=1,nO
|
||||
do j=1,nO
|
||||
do a=1,nV
|
||||
do b=1,nV
|
||||
do k=1,nO
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
do k=nC+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
|
||||
@ -62,12 +62,12 @@ subroutine form_v(nO,nV,X1,X2,X3,X4,t2,v)
|
||||
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
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
do k=nC+1,nO
|
||||
do c=1,nV-nR
|
||||
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
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine lCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
subroutine lCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Ladder CCD module
|
||||
|
||||
@ -10,19 +10,25 @@ subroutine lCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
|
||||
integer,intent(in) :: nBas,nEl
|
||||
integer,intent(in) :: nBasin
|
||||
integer,intent(in) :: nCin
|
||||
integer,intent(in) :: nOin
|
||||
integer,intent(in) :: nVin
|
||||
integer,intent(in) :: nRin
|
||||
double precision,intent(in) :: ENuc,ERHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: eHF(nBasin)
|
||||
double precision,intent(in) :: ERI(nBasin,nBasin,nBasin,nBasin)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: nBas2
|
||||
integer :: nBas
|
||||
integer :: nC
|
||||
integer :: nO
|
||||
integer :: nV
|
||||
integer :: nR
|
||||
integer :: nSCF
|
||||
double precision :: Conv
|
||||
double precision :: EcMP2,EcMP3,EcMP4
|
||||
double precision :: EcMP2
|
||||
double precision :: ECCD,EcCCD
|
||||
double precision,allocatable :: seHF(:)
|
||||
double precision,allocatable :: sERI(:,:,:,:)
|
||||
@ -63,35 +69,34 @@ subroutine lCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Spatial to spin orbitals
|
||||
|
||||
nBas2 = 2*nBas
|
||||
nBas = 2*nBasin
|
||||
nC = 2*nCin
|
||||
nO = 2*nOin
|
||||
nV = 2*nVin
|
||||
nR = 2*nRin
|
||||
|
||||
allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
|
||||
allocate(seHF(nBas),sERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
|
||||
call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI)
|
||||
call spatial_to_spin_MO_energy(nBasin,eHF,nBas,seHF)
|
||||
call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI)
|
||||
|
||||
! Antysymmetrize ERIs
|
||||
|
||||
allocate(dbERI(nBas2,nBas2,nBas2,nBas2))
|
||||
allocate(dbERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call antisymmetrize_ERI(2,nBas2,sERI,dbERI)
|
||||
call antisymmetrize_ERI(2,nBas,sERI,dbERI)
|
||||
|
||||
deallocate(sERI)
|
||||
|
||||
! Define occupied and virtual spaces
|
||||
|
||||
nO = 2*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)
|
||||
eV(:) = seHF(nO+1:nBas)
|
||||
|
||||
call form_delta_OOVV(nO,nV,eO,eV,delta_OOVV)
|
||||
call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV)
|
||||
|
||||
deallocate(seHF)
|
||||
|
||||
@ -99,10 +104,10 @@ subroutine lCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
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)
|
||||
OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO )
|
||||
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
|
||||
OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas, 1:nO ,nO+1:nBas)
|
||||
VVVV(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas,nO+1:nBas,nO+1:nBas)
|
||||
|
||||
deallocate(dbERI)
|
||||
|
||||
@ -112,8 +117,7 @@ subroutine lCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.))
|
||||
EcMP4 = 0d0
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcMP2)
|
||||
|
||||
! Memory allocation for DIIS
|
||||
|
||||
@ -150,13 +154,13 @@ subroutine lCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Compute residual
|
||||
|
||||
call form_ladder_r(nO,nV,OOOO,OOVV,VVVV,t2,r2)
|
||||
call form_ladder_r(nC,nO,nV,nR,OOOO,OOVV,VVVV,t2,r2)
|
||||
|
||||
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)
|
||||
|
||||
! Check convergence
|
||||
|
||||
Conv = maxval(abs(r2(:,:,:,:)))
|
||||
Conv = maxval(abs(r2(nC+1:nO,nC+1:nO,1:nV-nR,1:nV-nR)))
|
||||
|
||||
! Update amplitudes
|
||||
|
||||
@ -164,9 +168,7 @@ subroutine lCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! 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.))
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCCD)
|
||||
|
||||
! Dump results
|
||||
|
||||
|
@ -1,4 +1,4 @@
|
||||
subroutine rCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Ring CCD module
|
||||
|
||||
@ -10,19 +10,25 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
|
||||
integer,intent(in) :: nBas,nEl
|
||||
integer,intent(in) :: nBasin
|
||||
integer,intent(in) :: nCin
|
||||
integer,intent(in) :: nOin
|
||||
integer,intent(in) :: nVin
|
||||
integer,intent(in) :: nRin
|
||||
double precision,intent(in) :: ENuc,ERHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: eHF(nBasin)
|
||||
double precision,intent(in) :: ERI(nBasin,nBasin,nBasin,nBasin)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: nBas2
|
||||
integer :: nBas
|
||||
integer :: nC
|
||||
integer :: nO
|
||||
integer :: nV
|
||||
integer :: nR
|
||||
integer :: nSCF
|
||||
double precision :: Conv
|
||||
double precision :: EcMP2,EcMP3,EcMP4
|
||||
double precision :: EcMP2
|
||||
double precision :: ECCD,EcCCD
|
||||
double precision,allocatable :: seHF(:)
|
||||
double precision,allocatable :: sERI(:,:,:,:)
|
||||
@ -32,18 +38,8 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
double precision,allocatable :: eV(:)
|
||||
double precision,allocatable :: delta_OOVV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: OOOO(:,:,:,:)
|
||||
double precision,allocatable :: OOVV(:,:,:,:)
|
||||
double precision,allocatable :: OVVO(:,:,:,:)
|
||||
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(:,:,:,:)
|
||||
@ -63,46 +59,43 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Spatial to spin orbitals
|
||||
|
||||
nBas2 = 2*nBas
|
||||
nBas = 2*nBasin
|
||||
nC = 2*nCin
|
||||
nO = 2*nOin
|
||||
nV = 2*nVin
|
||||
nR = 2*nRin
|
||||
|
||||
allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
|
||||
allocate(seHF(nBas),sERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
|
||||
call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI)
|
||||
call spatial_to_spin_MO_energy(nBasin,eHF,nBas,seHF)
|
||||
call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI)
|
||||
|
||||
! Antysymmetrize ERIs
|
||||
|
||||
allocate(dbERI(nBas2,nBas2,nBas2,nBas2))
|
||||
allocate(dbERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call antisymmetrize_ERI(2,nBas2,sERI,dbERI)
|
||||
call antisymmetrize_ERI(2,nBas,sERI,dbERI)
|
||||
|
||||
deallocate(sERI)
|
||||
|
||||
! Define occupied and virtual spaces
|
||||
|
||||
nO = 2*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)
|
||||
eV(:) = seHF(nO+1:nBas)
|
||||
|
||||
call form_delta_OOVV(nO,nV,eO,eV,delta_OOVV)
|
||||
call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV)
|
||||
|
||||
deallocate(seHF)
|
||||
|
||||
! Create integral batches
|
||||
|
||||
allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VVVV(nV,nV,nV,nV))
|
||||
allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO))
|
||||
|
||||
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 )
|
||||
VVVV(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2,nO+1:nBas2,nO+1:nBas2)
|
||||
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
|
||||
OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO )
|
||||
|
||||
deallocate(dbERI)
|
||||
|
||||
@ -112,9 +105,7 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.))
|
||||
EcMP4 = 0d0
|
||||
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcMP2)
|
||||
|
||||
! Memory allocation for DIIS
|
||||
|
||||
@ -122,8 +113,7 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! 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))
|
||||
allocate(r2(nO,nO,nV,nV))
|
||||
|
||||
Conv = 1d0
|
||||
nSCF = 0
|
||||
@ -151,13 +141,13 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Compute residual
|
||||
|
||||
call form_ring_r(nO,nV,OVVO,OOVV,t2,r2)
|
||||
call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
|
||||
|
||||
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)
|
||||
|
||||
! Check convergence
|
||||
|
||||
Conv = maxval(abs(r2(:,:,:,:)))
|
||||
Conv = maxval(abs(r2(nC+1:nO,nC+1:nO,1:nV-nR,1:nV-nR)))
|
||||
|
||||
! Update amplitudes
|
||||
|
||||
@ -165,9 +155,7 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! 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.))
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCCD)
|
||||
|
||||
! Dump results
|
||||
|
||||
|
@ -100,10 +100,10 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
|
||||
! Read weird CC methods
|
||||
read(1,*)
|
||||
read(1,*) answer1,answer2,answer3,answer4
|
||||
if(answer4 == 'T') do_drCCD = .true.
|
||||
if(answer5 == 'T') do_rCCD = .true.
|
||||
if(answer6 == 'T') do_lCCD = .true.
|
||||
if(answer7 == 'T') do_pCCD = .true.
|
||||
if(answer1 == 'T') do_drCCD = .true.
|
||||
if(answer2 == 'T') do_rCCD = .true.
|
||||
if(answer3 == 'T') do_lCCD = .true.
|
||||
if(answer4 == 'T') do_pCCD = .true.
|
||||
|
||||
! Read excited state methods
|
||||
|
||||
|
@ -95,6 +95,8 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
|
||||
! Compute T-matrix version of the self-energy
|
||||
!----------------------------------------------
|
||||
|
||||
rho2(:,:,:) = 0d0
|
||||
|
||||
call self_energy_Tmatrix_diag_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), &
|
||||
Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), &
|
||||
SigT(:))
|
||||
|
@ -127,75 +127,75 @@ subroutine sort_ppRPA(ortho_eigvec,nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
|
||||
|
||||
if(ortho_eigvec) then
|
||||
|
||||
! Find degenerate eigenvalues
|
||||
! ! Find degenerate eigenvalues
|
||||
|
||||
deg1 = 1
|
||||
ab_start = 1
|
||||
! deg1 = 1
|
||||
! ab_start = 1
|
||||
|
||||
do ab=1,nVV
|
||||
! do ab=1,nVV
|
||||
|
||||
if(ab < nVV .and. abs(Omega1(ab) - Omega1(ab+1)) < 1d-6) then
|
||||
! if(ab < nVV .and. abs(Omega1(ab) - Omega1(ab+1)) < 1d-6) then
|
||||
|
||||
if(deg1 == 1) ab_start = ab
|
||||
deg1 = deg1 + 1
|
||||
! if(deg1 == 1) ab_start = ab
|
||||
! deg1 = deg1 + 1
|
||||
|
||||
else
|
||||
! else
|
||||
|
||||
ab_end = ab
|
||||
! ab_end = ab
|
||||
! print*,'deg = ',deg1,ab_start,ab_end
|
||||
|
||||
allocate(S1(deg1,deg1),O1(deg1,deg1))
|
||||
! allocate(S1(deg1,deg1),O1(deg1,deg1))
|
||||
|
||||
S1 = matmul(transpose(Z1(:,ab_start:ab_end)),matmul(M,Z1(:,ab_start:ab_end)))
|
||||
call orthogonalization_matrix(1,deg1,S1,O1)
|
||||
Z1(:,ab_start:ab_end) = matmul(Z1(:,ab_start:ab_end),O1)
|
||||
! S1 = matmul(transpose(Z1(:,ab_start:ab_end)),matmul(M,Z1(:,ab_start:ab_end)))
|
||||
! call orthogonalization_matrix(1,deg1,S1,O1)
|
||||
! Z1(:,ab_start:ab_end) = matmul(Z1(:,ab_start:ab_end),O1)
|
||||
|
||||
deallocate(S1,O1)
|
||||
! deallocate(S1,O1)
|
||||
|
||||
deg1 = 1
|
||||
ab_start = ab + 1
|
||||
! deg1 = 1
|
||||
! ab_start = ab + 1
|
||||
|
||||
end if
|
||||
end do
|
||||
! end if
|
||||
! end do
|
||||
|
||||
deg2 = 1
|
||||
ij_start = 1
|
||||
! deg2 = 1
|
||||
! ij_start = 1
|
||||
|
||||
do ij=1,nOO
|
||||
! do ij=1,nOO
|
||||
|
||||
if(ij < nOO .and. abs(Omega2(ij) - Omega2(ij+1)) < 1d-6) then
|
||||
! if(ij < nOO .and. abs(Omega2(ij) - Omega2(ij+1)) < 1d-6) then
|
||||
|
||||
if(deg2 == 1) ij_start = ij
|
||||
deg2 = deg2 + 1
|
||||
! if(deg2 == 1) ij_start = ij
|
||||
! deg2 = deg2 + 1
|
||||
|
||||
else
|
||||
! else
|
||||
|
||||
ij_end = ij
|
||||
! ij_end = ij
|
||||
! print*,'deg = ',deg2,ij_start,ij_end
|
||||
|
||||
allocate(S2(deg2,deg2),O2(deg2,deg2))
|
||||
! allocate(S2(deg2,deg2),O2(deg2,deg2))
|
||||
|
||||
S2 = - matmul(transpose(Z2(:,ij_start:ij_end)),matmul(M,Z2(:,ij_start:ij_end)))
|
||||
call orthogonalization_matrix(1,deg2,S2,O2)
|
||||
Z2(:,ij_start:ij_end) = matmul(Z2(:,ij_start:ij_end),O2)
|
||||
! S2 = - matmul(transpose(Z2(:,ij_start:ij_end)),matmul(M,Z2(:,ij_start:ij_end)))
|
||||
! call orthogonalization_matrix(1,deg2,S2,O2)
|
||||
! Z2(:,ij_start:ij_end) = matmul(Z2(:,ij_start:ij_end),O2)
|
||||
|
||||
deallocate(S2,O2)
|
||||
! deallocate(S2,O2)
|
||||
|
||||
deg2 = 1
|
||||
ij_start = ij + 1
|
||||
! deg2 = 1
|
||||
! ij_start = ij + 1
|
||||
|
||||
end if
|
||||
end do
|
||||
! end if
|
||||
! end do
|
||||
|
||||
! allocate(S1(nVV,nVV),S2(nOO,nOO),O1(nVV,nVV),O2(nOO,nOO))
|
||||
! S1 = + matmul(transpose(Z1),matmul(M,Z1))
|
||||
! S2 = - matmul(transpose(Z2),matmul(M,Z2))
|
||||
allocate(S1(nVV,nVV),S2(nOO,nOO),O1(nVV,nVV),O2(nOO,nOO))
|
||||
S1 = + matmul(transpose(Z1),matmul(M,Z1))
|
||||
S2 = - matmul(transpose(Z2),matmul(M,Z2))
|
||||
|
||||
! if(nVV > 0) call orthogonalization_matrix(1,nVV,S1,O1)
|
||||
! if(nOO > 0) call orthogonalization_matrix(1,nOO,S2,O2)
|
||||
if(nVV > 0) call orthogonalization_matrix(1,nVV,S1,O1)
|
||||
if(nOO > 0) call orthogonalization_matrix(1,nOO,S2,O2)
|
||||
|
||||
! Z1 = matmul(Z1,O1)
|
||||
! Z2 = matmul(Z2,O2)
|
||||
Z1 = matmul(Z1,O1)
|
||||
Z2 = matmul(Z2,O2)
|
||||
|
||||
end if
|
||||
|
||||
|
Loading…
Reference in New Issue
Block a user