From eee1a1c6b496c7b4badea56b2795aaf0dfeb7eaa Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sat, 21 Mar 2020 22:50:43 +0100 Subject: [PATCH] CCD variants --- input/basis | 79 +--------- input/methods | 8 +- input/molecule | 5 +- input/molecule.xyz | 5 +- input/weight | 79 +--------- src/QuAcK/QuAcK.f90 | 62 ++++++-- src/QuAcK/drCCD.f90 | 187 ++++++++++++++++++++++++ src/QuAcK/{ladder_CCD.f90 => lCCD.f90} | 14 +- src/QuAcK/pCCD.f90 | 195 +++++++++++++++++++++++++ src/QuAcK/ppRPA.f90 | 8 +- src/QuAcK/{ring_CCD.f90 => rCCD.f90} | 14 +- src/QuAcK/read_methods.f90 | 42 +++--- 12 files changed, 492 insertions(+), 206 deletions(-) create mode 100644 src/QuAcK/drCCD.f90 rename src/QuAcK/{ladder_CCD.f90 => lCCD.f90} (92%) create mode 100644 src/QuAcK/pCCD.f90 rename src/QuAcK/{ring_CCD.f90 => rCCD.f90} (92%) diff --git a/input/basis b/input/basis index 80e2a23..6796e3b 100644 --- a/input/basis +++ b/input/basis @@ -1,74 +1,9 @@ -1 10 -S 8 - 1 19500.0000000 0.0005070 - 2 2923.0000000 0.0039230 - 3 664.5000000 0.0202000 - 4 187.5000000 0.0790100 - 5 60.6200000 0.2304390 - 6 21.4200000 0.4328720 - 7 7.9500000 0.3499640 - 8 0.8815000 -0.0078920 -S 8 - 1 19500.0000000 -0.0001170 - 2 2923.0000000 -0.0009120 - 3 664.5000000 -0.0047170 - 4 187.5000000 -0.0190860 - 5 60.6200000 -0.0596550 - 6 21.4200000 -0.1400100 - 7 7.9500000 -0.1767820 - 8 0.8815000 0.6050430 +1 3 +S 3 + 1 38.3600000 0.0238090 + 2 5.7700000 0.1548910 + 3 1.2400000 0.4699870 S 1 - 1 2.2570000 1.0000000 -S 1 - 1 0.3041000 1.0000000 -P 3 - 1 43.8800000 0.0166650 - 2 9.9260000 0.1044720 - 3 2.9300000 0.3172600 + 1 0.2976000 1.0000000 P 1 - 1 0.9132000 1.0000000 -P 1 - 1 0.2672000 1.0000000 -D 1 - 1 3.1070000 1.0000000 -D 1 - 1 0.8550000 1.0000000 -F 1 - 1 1.9170000 1.0000000 -2 10 -S 8 - 1 19500.0000000 0.0005070 - 2 2923.0000000 0.0039230 - 3 664.5000000 0.0202000 - 4 187.5000000 0.0790100 - 5 60.6200000 0.2304390 - 6 21.4200000 0.4328720 - 7 7.9500000 0.3499640 - 8 0.8815000 -0.0078920 -S 8 - 1 19500.0000000 -0.0001170 - 2 2923.0000000 -0.0009120 - 3 664.5000000 -0.0047170 - 4 187.5000000 -0.0190860 - 5 60.6200000 -0.0596550 - 6 21.4200000 -0.1400100 - 7 7.9500000 -0.1767820 - 8 0.8815000 0.6050430 -S 1 - 1 2.2570000 1.0000000 -S 1 - 1 0.3041000 1.0000000 -P 3 - 1 43.8800000 0.0166650 - 2 9.9260000 0.1044720 - 3 2.9300000 0.3172600 -P 1 - 1 0.9132000 1.0000000 -P 1 - 1 0.2672000 1.0000000 -D 1 - 1 3.1070000 1.0000000 -D 1 - 1 0.8550000 1.0000000 -F 1 - 1 1.9170000 1.0000000 + 1 1.2750000 1.0000000 diff --git a/input/methods b/input/methods index 2077c1d..5e3dcaa 100644 --- a/input/methods +++ b/input/methods @@ -2,14 +2,14 @@ T F F # MP2 MP3 MP2-F12 T F F -# CCD CCSD CCSD(T) ringCCD ladderCCD - F F F F F +# CCD CCSD CCSD(T) drCCD rCCD lCCD pCCD + F F F T T T T # CIS RPA RPAx ppRPA ADC - T T T F F + F T T T F # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0 evGW qsGW - T F F + F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/molecule b/input/molecule index 9c99be9..c78e87e 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,4 @@ # nAt nEla nElb nCore nRyd - 2 9 9 0 0 + 1 1 1 0 0 # Znuc x y z - F 0. 0. 0. - F 0. 0. 2 + He 0.0 0.0 0.0 diff --git a/input/molecule.xyz b/input/molecule.xyz index f5dba60..797b5fc 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,3 @@ - 2 + 1 - F 0.0000000000 0.0000000000 0.0000000000 - F 0.0000000000 0.0000000000 1.0583544980 + He 0.0000000000 0.0000000000 0.0000000000 diff --git a/input/weight b/input/weight index 80e2a23..6796e3b 100644 --- a/input/weight +++ b/input/weight @@ -1,74 +1,9 @@ -1 10 -S 8 - 1 19500.0000000 0.0005070 - 2 2923.0000000 0.0039230 - 3 664.5000000 0.0202000 - 4 187.5000000 0.0790100 - 5 60.6200000 0.2304390 - 6 21.4200000 0.4328720 - 7 7.9500000 0.3499640 - 8 0.8815000 -0.0078920 -S 8 - 1 19500.0000000 -0.0001170 - 2 2923.0000000 -0.0009120 - 3 664.5000000 -0.0047170 - 4 187.5000000 -0.0190860 - 5 60.6200000 -0.0596550 - 6 21.4200000 -0.1400100 - 7 7.9500000 -0.1767820 - 8 0.8815000 0.6050430 +1 3 +S 3 + 1 38.3600000 0.0238090 + 2 5.7700000 0.1548910 + 3 1.2400000 0.4699870 S 1 - 1 2.2570000 1.0000000 -S 1 - 1 0.3041000 1.0000000 -P 3 - 1 43.8800000 0.0166650 - 2 9.9260000 0.1044720 - 3 2.9300000 0.3172600 + 1 0.2976000 1.0000000 P 1 - 1 0.9132000 1.0000000 -P 1 - 1 0.2672000 1.0000000 -D 1 - 1 3.1070000 1.0000000 -D 1 - 1 0.8550000 1.0000000 -F 1 - 1 1.9170000 1.0000000 -2 10 -S 8 - 1 19500.0000000 0.0005070 - 2 2923.0000000 0.0039230 - 3 664.5000000 0.0202000 - 4 187.5000000 0.0790100 - 5 60.6200000 0.2304390 - 6 21.4200000 0.4328720 - 7 7.9500000 0.3499640 - 8 0.8815000 -0.0078920 -S 8 - 1 19500.0000000 -0.0001170 - 2 2923.0000000 -0.0009120 - 3 664.5000000 -0.0047170 - 4 187.5000000 -0.0190860 - 5 60.6200000 -0.0596550 - 6 21.4200000 -0.1400100 - 7 7.9500000 -0.1767820 - 8 0.8815000 0.6050430 -S 1 - 1 2.2570000 1.0000000 -S 1 - 1 0.3041000 1.0000000 -P 3 - 1 43.8800000 0.0166650 - 2 9.9260000 0.1044720 - 3 2.9300000 0.3172600 -P 1 - 1 0.9132000 1.0000000 -P 1 - 1 0.2672000 1.0000000 -D 1 - 1 3.1070000 1.0000000 -D 1 - 1 0.8550000 1.0000000 -F 1 - 1 1.9170000 1.0000000 + 1 1.2750000 1.0000000 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 6f0ae8f..9a5e06e 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -7,7 +7,7 @@ program QuAcK logical :: doRHF,doUHF,doMOM logical :: doMP2,doMP3,doMP2F12 logical :: doCCD,doCCSD,doCCSDT - logical :: do_ring_CCD,do_ladder_CCD + logical :: do_drCCD,do_rCCD,do_lCCD,do_pCCD logical :: doCIS,doRPA,doRPAx logical :: doppRPA,doADC logical :: doG0F2,doevGF2,doG0F3,doevGF3 @@ -114,15 +114,15 @@ program QuAcK ! Which calculations do you want to do? - call read_methods(doRHF,doUHF,doMOM, & - doMP2,doMP3,doMP2F12, & - doCCD,doCCSD,doCCSDT, & - do_ring_CCD,do_ladder_CCD, & - doCIS,doRPA,doRPAx, & - doppRPA,doADC, & - doG0F2,doevGF2,doG0F3,doevGF3, & - doG0W0,doevGW,doqsGW, & - doG0T0,doevGT,doqsGT, & + call read_methods(doRHF,doUHF,doMOM, & + doMP2,doMP3,doMP2F12, & + doCCD,doCCSD,doCCSDT, & + do_drCCD,do_rCCD,do_lCCD,do_pCCD, & + doCIS,doRPA,doRPAx, & + doppRPA,doADC, & + doG0F2,doevGF2,doG0F3,doevGF3, & + doG0W0,doevGW,doqsGW, & + doG0T0,doevGT,doqsGT, & doMCMP2) ! Read options for methods @@ -392,14 +392,30 @@ program QuAcK end if +!------------------------------------------------------------------------ +! Perform direct ring CCD calculation +!------------------------------------------------------------------------ + + 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 cpu_time(end_CCD) + + t_CCD = end_CCD - start_CCD + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for direct ring CCD = ',t_CCD,' seconds' + write(*,*) + + end if + !------------------------------------------------------------------------ ! Perform ring CCD calculation !------------------------------------------------------------------------ - if(do_ring_CCD) then + if(do_rCCD) then call cpu_time(start_CCD) - call ring_CCD(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,nEl,ERI_MO_basis,ENuc,ERHF,eHF) call cpu_time(end_CCD) t_CCD = end_CCD - start_CCD @@ -412,14 +428,30 @@ program QuAcK ! Perform ladder CCD calculation !------------------------------------------------------------------------ - if(do_ladder_CCD) then + if(do_lCCD) then call cpu_time(start_CCD) - call ladder_CCD(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,nEl,ERI_MO_basis,ENuc,ERHF,eHF) call cpu_time(end_CCD) t_CCD = end_CCD - start_CCD - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CCD,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ladder CCD = ',t_CCD,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform pair CCD calculation +!------------------------------------------------------------------------ + + if(do_pCCD) then + + call cpu_time(start_CCD) + call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF) + call cpu_time(end_CCD) + + t_CCD = end_CCD - start_CCD + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pair CCD = ',t_CCD,' seconds' write(*,*) end if diff --git a/src/QuAcK/drCCD.f90 b/src/QuAcK/drCCD.f90 new file mode 100644 index 0000000..df63a2d --- /dev/null +++ b/src/QuAcK/drCCD.f90 @@ -0,0 +1,187 @@ +subroutine drCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) + +! Direct ring CCD module + + implicit none + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + + integer,intent(in) :: nBas,nEl + double precision,intent(in) :: ENuc,ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: nBas2 + integer :: nO + integer :: nV + integer :: nSCF + double precision :: Conv + double precision :: EcMP2,EcMP3,EcMP4 + double precision :: ECCD,EcCCD + double precision,allocatable :: seHF(:) + double precision,allocatable :: sERI(:,:,:,:) + + double precision,allocatable :: eO(:) + 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(:,:,:,:) + +! Hello world + + write(*,*) + write(*,*)'**************************************' + write(*,*)'| direct ring CCD calculation |' + write(*,*)'**************************************' + write(*,*) + +! Spatial to spin orbitals + + nBas2 = 2*nBas + + allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2)) + + call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF) + call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI) + +! Antysymmetrize ERIs + +! 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) + + call form_delta_OOVV(nO,nV,eO,eV,delta_OOVV) + + deallocate(seHF) + +! Create integral batches + + allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VVVV(nV,nV,nV,nV)) + + 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) + + deallocate(sERI) + +! MP2 guess amplitudes + + allocate(t2(nO,nO,nV,nV)) + + t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) + + EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.)) + EcMP4 = 0d0 + +! Initialization + + allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(nO,nO,nV,nV)) + allocate(X1(nO,nO,nO,nO),X2(nV,nV),X3(nO,nO),X4(nO,nO,nV,nV)) + + Conv = 1d0 + nSCF = 0 + +!------------------------------------------------------------------------ +! Main SCF loop +!------------------------------------------------------------------------ + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)'| direct ring CCD calculation |' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(CCD)','|','Ec(CCD)','|','Conv','|' + write(*,*)'----------------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + +! Increment + + nSCF = nSCF + 1 + +! Compute residual + + call form_ring_r(nO,nV,OVVO,OOVV,t2,r2) + + r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:) + +! Check convergence + + Conv = maxval(abs(r2(:,:,:,:))) + +! 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.)) + +! Dump results + + ECCD = ERHF + EcCCD + + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & + '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|' + + enddo + write(*,*)'----------------------------------------------------' +!------------------------------------------------------------------------ +! End of SCF loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + endif + + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)' direct ring CCD energy ' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A30,1X,F15.10)')' E(drCCD) = ',ECCD + write(*,'(1X,A30,1X,F15.10)')' Ec(drCCD) = ',EcCCD + write(*,*)'----------------------------------------------------' + write(*,*) + +end subroutine drCCD diff --git a/src/QuAcK/ladder_CCD.f90 b/src/QuAcK/lCCD.f90 similarity index 92% rename from src/QuAcK/ladder_CCD.f90 rename to src/QuAcK/lCCD.f90 index 43c3586..7ca2c20 100644 --- a/src/QuAcK/ladder_CCD.f90 +++ b/src/QuAcK/lCCD.f90 @@ -1,4 +1,4 @@ -subroutine ladder_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) +subroutine lCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) ! Ladder CCD module @@ -52,7 +52,7 @@ subroutine ladder_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) write(*,*) write(*,*)'**************************************' - write(*,*)'| ladder-CCD calculation |' + write(*,*)'| ladder CCD calculation |' write(*,*)'**************************************' write(*,*) @@ -123,7 +123,7 @@ subroutine ladder_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) !------------------------------------------------------------------------ write(*,*) write(*,*)'----------------------------------------------------' - write(*,*)'| ladder-CCD calculation |' + write(*,*)'| ladder CCD calculation |' write(*,*)'----------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & '|','#','|','E(CCD)','|','Ec(CCD)','|','Conv','|' @@ -184,11 +184,11 @@ subroutine ladder_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) write(*,*) write(*,*)'----------------------------------------------------' - write(*,*)' ladder-CCD energy ' + write(*,*)' ladder CCD energy ' write(*,*)'----------------------------------------------------' - write(*,'(1X,A30,1X,F15.10)')' E(ladder-CCD) = ',ECCD - write(*,'(1X,A30,1X,F15.10)')' Ec(ladder-CCSD) = ',EcCCD + write(*,'(1X,A30,1X,F15.10)')' E(lCCD) = ',ECCD + write(*,'(1X,A30,1X,F15.10)')' Ec(lCCD) = ',EcCCD write(*,*)'----------------------------------------------------' write(*,*) -end subroutine ladder_CCD +end subroutine lCCD diff --git a/src/QuAcK/pCCD.f90 b/src/QuAcK/pCCD.f90 new file mode 100644 index 0000000..f29aaa4 --- /dev/null +++ b/src/QuAcK/pCCD.f90 @@ -0,0 +1,195 @@ +subroutine pCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) + +! pair CCD module + + implicit none + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + + integer,intent(in) :: nBas,nEl + double precision,intent(in) :: ENuc,ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: nBas2 + integer :: nO + integer :: nV + integer :: nSCF + double precision :: Conv + double precision :: EcMP2,EcMP3,EcMP4 + double precision :: ECCD,EcCCD + double precision,allocatable :: seHF(:) + double precision,allocatable :: sERI(:,:,:,:) + double precision,allocatable :: dbERI(:,:,:,:) + + double precision,allocatable :: eO(:) + double precision,allocatable :: eV(:) + double precision,allocatable :: delta_OOVV(:,:,:,:) + + double precision,allocatable :: OOOO(:,:,:,:) + double precision,allocatable :: OOVV(:,:,:,:) + double precision,allocatable :: OVOV(:,:,:,:) + double precision,allocatable :: VVVV(:,:,:,:) + + double precision,allocatable :: X1(:,:,:,:) + double precision,allocatable :: X2(:,:) + double precision,allocatable :: X3(:,:) + double precision,allocatable :: X4(:,:,:,:) + + double precision,allocatable :: u(:,:,:,:) + double precision,allocatable :: v(:,:,:,:) + + double precision,allocatable :: r2(:,:,:,:) + double precision,allocatable :: t2(:,:,:,:) + +! Hello world + + write(*,*) + write(*,*)'**************************************' + write(*,*)'| CCD calculation |' + write(*,*)'**************************************' + write(*,*) + +! Spatial to spin orbitals + + nBas2 = 2*nBas + + allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2)) + + call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF) + call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI) + +! Antysymmetrize ERIs + + allocate(dbERI(nBas2,nBas2,nBas2,nBas2)) + + call antisymmetrize_ERI(2,nBas2,sERI,dbERI) + + deallocate(sERI) + +! Define occupied and virtual spaces + + nO = 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) + + call form_delta_OOVV(nO,nV,eO,eV,delta_OOVV) + + deallocate(seHF) + +! Create integral batches + + allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVOV(nO,nV,nO,nV),VVVV(nV,nV,nV,nV)) + + OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO ) + OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas2,nO+1:nBas2) + OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas2, 1:nO ,nO+1:nBas2) + VVVV(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2,nO+1:nBas2,nO+1:nBas2) + + deallocate(dbERI) + +! MP2 guess amplitudes + + allocate(t2(nO,nO,nV,nV)) + + t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) + + EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.)) + EcMP4 = 0d0 + +! Initialization + + allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(nO,nO,nV,nV)) + allocate(X1(nO,nO,nO,nO),X2(nV,nV),X3(nO,nO),X4(nO,nO,nV,nV)) + + Conv = 1d0 + nSCF = 0 + +!------------------------------------------------------------------------ +! Main SCF loop +!------------------------------------------------------------------------ + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)'| pair CCD calculation |' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(pCCD)','|','Ec(pCCD)','|','Conv','|' + write(*,*)'----------------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + +! Increment + + nSCF = nSCF + 1 + +! Form linear array + + call form_u(nO,nV,OOOO,VVVV,OVOV,t2,u) + +! Form interemediate arrays + + call form_X(nO,nV,OOVV,t2,X1,X2,X3,X4) + +! Form quadratic array + + call form_v(nO,nV,X1,X2,X3,X4,t2,v) + +! Compute residual + + r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + u(:,:,:,:) + v(:,:,:,:) + +! Check convergence + + Conv = maxval(abs(r2(:,:,:,:))) + +! Update amplitudes + + t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:) + +! Compute correlation energy + + EcCCD = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.)) + + if(nSCF == 1) EcMP3 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2 + v/delta_OOVV,.true.)) + +! Dump results + + ECCD = ERHF + EcCCD + + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & + '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|' + + enddo + write(*,*)'----------------------------------------------------' +!------------------------------------------------------------------------ +! End of SCF loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + endif + +end subroutine pCCD diff --git a/src/QuAcK/ppRPA.f90 b/src/QuAcK/ppRPA.f90 index 642fc23..3c45cc5 100644 --- a/src/QuAcK/ppRPA.f90 +++ b/src/QuAcK/ppRPA.f90 @@ -104,10 +104,10 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A40,F15.6)') 'Tr@ppRPA correlation energy (singlet) =',Ec_ppRPA(1) - write(*,'(2X,A40,F15.6)') 'Tr@ppRPA correlation energy (triplet) =',3d0*Ec_ppRPA(2) - write(*,'(2X,A40,F15.6)') 'Tr@ppRPA correlation energy =',Ec_ppRPA(1) + 3d0*Ec_ppRPA(2) - write(*,'(2X,A40,F15.6)') 'Tr@ppRPA total energy =',ENuc + ERHF + Ec_ppRPA(1) + 3d0*Ec_ppRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (singlet) =',Ec_ppRPA(1) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (triplet) =',3d0*Ec_ppRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy =',Ec_ppRPA(1) + 3d0*Ec_ppRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy =',ENuc + ERHF + Ec_ppRPA(1) + 3d0*Ec_ppRPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/QuAcK/ring_CCD.f90 b/src/QuAcK/rCCD.f90 similarity index 92% rename from src/QuAcK/ring_CCD.f90 rename to src/QuAcK/rCCD.f90 index 717ceda..e53dd52 100644 --- a/src/QuAcK/ring_CCD.f90 +++ b/src/QuAcK/rCCD.f90 @@ -1,4 +1,4 @@ -subroutine ring_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) +subroutine rCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) ! Ring CCD module @@ -52,7 +52,7 @@ subroutine ring_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) write(*,*) write(*,*)'**************************************' - write(*,*)'| ring-CCD calculation |' + write(*,*)'| ring CCD calculation |' write(*,*)'**************************************' write(*,*) @@ -123,7 +123,7 @@ subroutine ring_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) !------------------------------------------------------------------------ write(*,*) write(*,*)'----------------------------------------------------' - write(*,*)'| ring-CCD calculation |' + write(*,*)'| ring CCD calculation |' write(*,*)'----------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & '|','#','|','E(CCD)','|','Ec(CCD)','|','Conv','|' @@ -184,11 +184,11 @@ subroutine ring_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) write(*,*) write(*,*)'----------------------------------------------------' - write(*,*)' ring-CCD energy ' + write(*,*)' ring CCD energy ' write(*,*)'----------------------------------------------------' - write(*,'(1X,A30,1X,F15.10)')' E(ring-CCD) = ',ECCD - write(*,'(1X,A30,1X,F15.10)')' Ec(ring-CCD) = ',EcCCD + write(*,'(1X,A30,1X,F15.10)')' E(rCCD) = ',ECCD + write(*,'(1X,A30,1X,F15.10)')' Ec(rCCD) = ',EcCCD write(*,*)'----------------------------------------------------' write(*,*) -end subroutine ring_CCD +end subroutine rCCD diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index e8e89fa..99d14ef 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -1,12 +1,12 @@ -subroutine read_methods(doRHF,doUHF,doMOM, & - doMP2,doMP3,doMP2F12, & - doCCD,doCCSD,doCCSDT, & - do_ring_CCD,do_ladder_CCD, & - doCIS,doRPA,doRPAx, & - doppRPA,doADC, & - doG0F2,doevGF2,doG0F3,doevGF3, & - doG0W0,doevGW,doqsGW, & - doG0T0,doevGT,doqsGT, & +subroutine read_methods(doRHF,doUHF,doMOM, & + doMP2,doMP3,doMP2F12, & + doCCD,doCCSD,doCCSDT, & + do_drCCD,do_rCCD,do_lCCD,do_pCCD, & + doCIS,doRPA,doRPAx, & + doppRPA,doADC, & + doG0F2,doevGF2,doG0F3,doevGF3, & + doG0W0,doevGW,doqsGW, & + doG0T0,doevGT,doqsGT, & doMCMP2) ! Read desired methods @@ -18,7 +18,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, & logical,intent(out) :: doRHF,doUHF,doMOM logical,intent(out) :: doMP2,doMP3,doMP2F12 logical,intent(out) :: doCCD,doCCSD,doCCSDT - logical,intent(out) :: do_ring_CCD,do_ladder_CCD + logical,intent(out) :: do_drCCD,do_rCCD,do_lCCD,do_pCCD logical,intent(out) :: doCIS,doRPA,doRPAx,doppRPA,doADC logical,intent(out) :: doG0F2,doevGF2,doG0F3,doevGF3 logical,intent(out) :: doG0W0,doevGW,doqsGW @@ -27,7 +27,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, & ! Local variables - character(len=1) :: answer1,answer2,answer3,answer4,answer5 + character(len=1) :: answer1,answer2,answer3,answer4,answer5,answer6,answer7 ! Open file with method specification @@ -47,8 +47,10 @@ subroutine read_methods(doRHF,doUHF,doMOM, & doCCSD = .false. doCCSDT = .false. - do_ring_CCD = .false. - do_ladder_CCD = .false. + do_drCCD = .false. + do_rCCD = .false. + do_lCCD = .false. + do_pCCD = .false. doCIS = .false. doRPA = .false. @@ -90,12 +92,14 @@ subroutine read_methods(doRHF,doUHF,doMOM, & ! Read CC methods read(1,*) - read(1,*) answer1,answer2,answer3,answer4,answer5 - if(answer1 == 'T') doCCD = .true. - if(answer2 == 'T') doCCSD = .true. - if(answer3 == 'T') doCCSDT = .true. - if(answer4 == 'T') do_ring_CCD = .true. - if(answer5 == 'T') do_ladder_CCD = .true. + read(1,*) answer1,answer2,answer3,answer4,answer5,answer6,answer7 + if(answer1 == 'T') doCCD = .true. + if(answer2 == 'T') doCCSD = .true. + if(answer3 == 'T') doCCSDT = .true. + if(answer4 == 'T') do_drCCD = .true. + if(answer5 == 'T') do_rCCD = .true. + if(answer6 == 'T') do_lCCD = .true. + if(answer7 == 'T') do_pCCD = .true. ! Read excited state methods