From 4983c8c16b056793944f4a5f95463f7911a22543 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 26 Mar 2020 16:57:46 +0100 Subject: [PATCH] core in CC --- input/basis | 35 +++++++-- input/int | 4 +- input/methods | 2 +- input/molecule | 4 +- input/molecule.xyz | 2 +- input/options | 2 +- input/weight | 35 +++++++-- src/IntPak/CalcOmErf.f90 | 2 +- src/QuAcK/CCD.f90 | 63 ++++++++-------- src/QuAcK/CCD_correlation_energy.f90 | 36 ++++++++++ src/QuAcK/CCSD.f90 | 92 +++++++++++++----------- src/QuAcK/CCSDT.f90 | 12 ++-- src/QuAcK/CCSD_correlation_energy.f90 | 36 ++++++++++ src/QuAcK/G0T0.f90 | 3 + src/QuAcK/MP3_correlation_energy.f90 | 38 ++++++++++ src/QuAcK/QuAcK.f90 | 18 +++-- src/QuAcK/RPA.f90 | 4 +- src/QuAcK/RPAx.f90 | 4 +- src/QuAcK/drCCD.f90 | 76 +++++++++----------- src/QuAcK/excitation_density_Tmatrix.f90 | 18 ++++- src/QuAcK/form_T.f90 | 16 ++--- src/QuAcK/form_X.f90 | 48 ++++++------- src/QuAcK/form_abh.f90 | 48 ++++++------- src/QuAcK/form_delta_OOOVVV.f90 | 16 ++--- src/QuAcK/form_delta_OOVV.f90 | 12 ++-- src/QuAcK/form_delta_OV.f90 | 8 +-- src/QuAcK/form_g.f90 | 20 +++--- src/QuAcK/form_h.f90 | 32 ++++----- src/QuAcK/form_ladder_r.f90 | 50 ++++++++----- src/QuAcK/form_r1.f90 | 32 ++++----- src/QuAcK/form_r2.f90 | 68 +++++++++--------- src/QuAcK/form_ring_r.f90 | 45 ++++++++---- src/QuAcK/form_tau.f90 | 12 ++-- src/QuAcK/form_u.f90 | 40 +++++------ src/QuAcK/form_ub.f90 | 16 ++--- src/QuAcK/form_ubb.f90 | 20 +++--- src/QuAcK/form_v.f90 | 48 ++++++------- src/QuAcK/lCCD.f90 | 62 ++++++++-------- src/QuAcK/rCCD.f90 | 76 +++++++++----------- src/QuAcK/read_methods.f90 | 8 +-- src/QuAcK/soG0T0.f90 | 2 + src/QuAcK/sort_ppRPA.f90 | 84 +++++++++++----------- 42 files changed, 722 insertions(+), 527 deletions(-) create mode 100644 src/QuAcK/CCD_correlation_energy.f90 create mode 100644 src/QuAcK/CCSD_correlation_energy.f90 create mode 100644 src/QuAcK/MP3_correlation_energy.f90 diff --git a/input/basis b/input/basis index 6796e3b..f19a2d0 100644 --- a/input/basis +++ b/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 + diff --git a/input/int b/input/int index ed14d90..9341cab 100644 --- a/input/int +++ b/input/int @@ -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 diff --git a/input/methods b/input/methods index c1d2dda..52e59fd 100644 --- a/input/methods +++ b/input/methods @@ -13,6 +13,6 @@ # G0W0 evGW qsGW F F F # G0T0 evGT qsGT - T F F + F F F # MCMP2 F diff --git a/input/molecule b/input/molecule index c78e87e..edeba31 100644 --- a/input/molecule +++ b/input/molecule @@ -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 diff --git a/input/molecule.xyz b/input/molecule.xyz index 797b5fc..1c70680 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,3 +1,3 @@ 1 - He 0.0000000000 0.0000000000 0.0000000000 + Ne 0.0000000000 0.0000000000 0.0000000000 diff --git a/input/options b/input/options index 2de160f..f20ae20 100644 --- a/input/options +++ b/input/options @@ -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 diff --git a/input/weight b/input/weight index 6796e3b..f19a2d0 100644 --- a/input/weight +++ b/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 + diff --git a/src/IntPak/CalcOmErf.f90 b/src/IntPak/CalcOmErf.f90 index 26adc72..ef1535a 100644 --- a/src/IntPak/CalcOmErf.f90 +++ b/src/IntPak/CalcOmErf.f90 @@ -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 diff --git a/src/QuAcK/CCD.f90 b/src/QuAcK/CCD.f90 index 3929415..b662148 100644 --- a/src/QuAcK/CCD.f90 +++ b/src/QuAcK/CCD.f90 @@ -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 diff --git a/src/QuAcK/CCD_correlation_energy.f90 b/src/QuAcK/CCD_correlation_energy.f90 new file mode 100644 index 0000000..b8bd52b --- /dev/null +++ b/src/QuAcK/CCD_correlation_energy.f90 @@ -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 diff --git a/src/QuAcK/CCSD.f90 b/src/QuAcK/CCSD.f90 index 3240b44..8936721 100644 --- a/src/QuAcK/CCSD.f90 +++ b/src/QuAcK/CCSD.f90 @@ -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 diff --git a/src/QuAcK/CCSDT.f90 b/src/QuAcK/CCSDT.f90 index 2af9b3b..5b58b64 100644 --- a/src/QuAcK/CCSDT.f90 +++ b/src/QuAcK/CCSDT.f90 @@ -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 diff --git a/src/QuAcK/CCSD_correlation_energy.f90 b/src/QuAcK/CCSD_correlation_energy.f90 new file mode 100644 index 0000000..3a678cb --- /dev/null +++ b/src/QuAcK/CCSD_correlation_energy.f90 @@ -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 diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90 index f44ffb4..17eada4 100644 --- a/src/QuAcK/G0T0.f90 +++ b/src/QuAcK/G0T0.f90 @@ -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(:,:,:), & diff --git a/src/QuAcK/MP3_correlation_energy.f90 b/src/QuAcK/MP3_correlation_energy.f90 new file mode 100644 index 0000000..739335f --- /dev/null +++ b/src/QuAcK/MP3_correlation_energy.f90 @@ -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 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index d8ad4dc..c2041c3 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -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 diff --git a/src/QuAcK/RPA.f90 b/src/QuAcK/RPA.f90 index 03547a2..a01a30e 100644 --- a/src/QuAcK/RPA.f90 +++ b/src/QuAcK/RPA.f90 @@ -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 diff --git a/src/QuAcK/RPAx.f90 b/src/QuAcK/RPAx.f90 index 6674ddf..c57d638 100644 --- a/src/QuAcK/RPAx.f90 +++ b/src/QuAcK/RPAx.f90 @@ -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 diff --git a/src/QuAcK/drCCD.f90 b/src/QuAcK/drCCD.f90 index b3cbca0..063c4b0 100644 --- a/src/QuAcK/drCCD.f90 +++ b/src/QuAcK/drCCD.f90 @@ -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 diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/QuAcK/excitation_density_Tmatrix.f90 index 1ac2d24..32e74fa 100644 --- a/src/QuAcK/excitation_density_Tmatrix.f90 +++ b/src/QuAcK/excitation_density_Tmatrix.f90 @@ -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 diff --git a/src/QuAcK/form_T.f90 b/src/QuAcK/form_T.f90 index dc66b63..54a47c6 100644 --- a/src/QuAcK/form_T.f90 +++ b/src/QuAcK/form_T.f90 @@ -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)) & diff --git a/src/QuAcK/form_X.f90 b/src/QuAcK/form_X.f90 index 9f2eeb5..699e418 100644 --- a/src/QuAcK/form_X.f90 +++ b/src/QuAcK/form_X.f90 @@ -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 diff --git a/src/QuAcK/form_abh.f90 b/src/QuAcK/form_abh.f90 index 4b0efcc..5e6300c 100644 --- a/src/QuAcK/form_abh.f90 +++ b/src/QuAcK/form_abh.f90 @@ -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 diff --git a/src/QuAcK/form_delta_OOOVVV.f90 b/src/QuAcK/form_delta_OOOVVV.f90 index 973b579..e29989d 100644 --- a/src/QuAcK/form_delta_OOOVVV.f90 +++ b/src/QuAcK/form_delta_OOOVVV.f90 @@ -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) diff --git a/src/QuAcK/form_delta_OOVV.f90 b/src/QuAcK/form_delta_OOVV.f90 index 15bb6d9..fbcd2d6 100644 --- a/src/QuAcK/form_delta_OOVV.f90 +++ b/src/QuAcK/form_delta_OOVV.f90 @@ -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) diff --git a/src/QuAcK/form_delta_OV.f90 b/src/QuAcK/form_delta_OV.f90 index e785a0c..02b737d 100644 --- a/src/QuAcK/form_delta_OV.f90 +++ b/src/QuAcK/form_delta_OV.f90 @@ -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 diff --git a/src/QuAcK/form_g.f90 b/src/QuAcK/form_g.f90 index ebb8427..b9475d1 100644 --- a/src/QuAcK/form_g.f90 +++ b/src/QuAcK/form_g.f90 @@ -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 diff --git a/src/QuAcK/form_h.f90 b/src/QuAcK/form_h.f90 index 4bc16f1..b1405c7 100644 --- a/src/QuAcK/form_h.f90 +++ b/src/QuAcK/form_h.f90 @@ -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) diff --git a/src/QuAcK/form_ladder_r.f90 b/src/QuAcK/form_ladder_r.f90 index b66e1bd..678e25f 100644 --- a/src/QuAcK/form_ladder_r.f90 +++ b/src/QuAcK/form_ladder_r.f90 @@ -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 diff --git a/src/QuAcK/form_r1.f90 b/src/QuAcK/form_r1.f90 index ab1bd57..1392c9f 100644 --- a/src/QuAcK/form_r1.f90 +++ b/src/QuAcK/form_r1.f90 @@ -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 diff --git a/src/QuAcK/form_r2.f90 b/src/QuAcK/form_r2.f90 index c57b8a8..c0508c3 100644 --- a/src/QuAcK/form_r2.f90 +++ b/src/QuAcK/form_r2.f90 @@ -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 diff --git a/src/QuAcK/form_ring_r.f90 b/src/QuAcK/form_ring_r.f90 index 099af66..162cff5 100644 --- a/src/QuAcK/form_ring_r.f90 +++ b/src/QuAcK/form_ring_r.f90 @@ -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 diff --git a/src/QuAcK/form_tau.f90 b/src/QuAcK/form_tau.f90 index d666383..918c219 100644 --- a/src/QuAcK/form_tau.f90 +++ b/src/QuAcK/form_tau.f90 @@ -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) diff --git a/src/QuAcK/form_u.f90 b/src/QuAcK/form_u.f90 index c8fc76d..739f5bd 100644 --- a/src/QuAcK/form_u.f90 +++ b/src/QuAcK/form_u.f90 @@ -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) & diff --git a/src/QuAcK/form_ub.f90 b/src/QuAcK/form_ub.f90 index 14f8d9c..abbb11d 100644 --- a/src/QuAcK/form_ub.f90 +++ b/src/QuAcK/form_ub.f90 @@ -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) & diff --git a/src/QuAcK/form_ubb.f90 b/src/QuAcK/form_ubb.f90 index 5beaa6e..252b5a6 100644 --- a/src/QuAcK/form_ubb.f90 +++ b/src/QuAcK/form_ubb.f90 @@ -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) & diff --git a/src/QuAcK/form_v.f90 b/src/QuAcK/form_v.f90 index f589d31..12c6150 100644 --- a/src/QuAcK/form_v.f90 +++ b/src/QuAcK/form_v.f90 @@ -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 diff --git a/src/QuAcK/lCCD.f90 b/src/QuAcK/lCCD.f90 index c92c5e7..1b748c7 100644 --- a/src/QuAcK/lCCD.f90 +++ b/src/QuAcK/lCCD.f90 @@ -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 diff --git a/src/QuAcK/rCCD.f90 b/src/QuAcK/rCCD.f90 index 5331be2..bc8c10d 100644 --- a/src/QuAcK/rCCD.f90 +++ b/src/QuAcK/rCCD.f90 @@ -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 diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index 56a5d67..ae8ede1 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -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 diff --git a/src/QuAcK/soG0T0.f90 b/src/QuAcK/soG0T0.f90 index fc218d4..30c9f00 100644 --- a/src/QuAcK/soG0T0.f90 +++ b/src/QuAcK/soG0T0.f90 @@ -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(:)) diff --git a/src/QuAcK/sort_ppRPA.f90 b/src/QuAcK/sort_ppRPA.f90 index cbdb16c..258dd9c 100644 --- a/src/QuAcK/sort_ppRPA.f90 +++ b/src/QuAcK/sort_ppRPA.f90 @@ -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