From 3318daf241fd1672baa537e6ea3e33253473b918 Mon Sep 17 00:00:00 2001 From: pfloos Date: Mon, 13 Nov 2023 22:15:00 +0100 Subject: [PATCH] plenty of CC for GHF reference --- input/methods | 18 +-- input/options | 2 +- src/CC/GCC.f90 | 172 +++++++++++++++++++++++++ src/CC/GCCD.f90 | 262 +++++++++++++++++++++++++++++++++++++ src/CC/GCCSD.f90 | 291 ++++++++++++++++++++++++++++++++++++++++++ src/CC/RCC.f90 | 23 ++-- src/CC/crGCCD.f90 | 189 +++++++++++++++++++++++++++ src/CC/drGCCD.f90 | 181 ++++++++++++++++++++++++++ src/CC/lGCCD.f90 | 202 +++++++++++++++++++++++++++++ src/CC/rCCD.f90 | 13 +- src/CC/rGCCD.f90 | 208 ++++++++++++++++++++++++++++++ src/HF/GHF_search.f90 | 2 +- src/HF/RHF_search.f90 | 2 +- src/HF/UHF_search.f90 | 2 +- src/QuAcK/GQuAcK.f90 | 63 ++++++--- src/QuAcK/QuAcK.f90 | 5 +- src/QuAcK/RQuAcK.f90 | 34 ++--- src/QuAcK/UQuAcK.f90 | 32 ++--- 18 files changed, 1615 insertions(+), 86 deletions(-) create mode 100644 src/CC/GCC.f90 create mode 100644 src/CC/GCCD.f90 create mode 100644 src/CC/GCCSD.f90 create mode 100644 src/CC/crGCCD.f90 create mode 100644 src/CC/drGCCD.f90 create mode 100644 src/CC/lGCCD.f90 create mode 100644 src/CC/rGCCD.f90 diff --git a/input/methods b/input/methods index 956bcdf..f99a66d 100644 --- a/input/methods +++ b/input/methods @@ -1,20 +1,20 @@ # RHF UHF GHF ROHF - T F F F + F F T F # MP2 MP3 T T # CCD pCCD DCD CCSD CCSD(T) - T T T T F + F F F F F # drCCD rCCD crCCD lCCD - T T T T + F F F F # CIS CIS(D) CID CISD FCI - T F F F F + F F F F F # phRPA phRPAx crRPA ppRPA - T T T T + F F F F # G0F2 evGF2 qsGF2 G0F3 evGF3 - T T F F F + F F F F F # G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW - T T F F F F + F F F F F F # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh - T T F F F F + F F F F F F # Rtest Utest Gtest - T F F + F F F diff --git a/input/options b/input/options index 92084cd..b124ee7 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ # HF: maxSCF thresh DIIS guess mix shift stab search - 10000 0.0000001 5 1 0.0 0.0 F F + 10000 0.0000001 5 2 0.0 0.0 F T # MP: reg F # CC: maxSCF thresh DIIS diff --git a/src/CC/GCC.f90 b/src/CC/GCC.f90 new file mode 100644 index 0000000..4ce7985 --- /dev/null +++ b/src/CC/GCC.f90 @@ -0,0 +1,172 @@ +subroutine GCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & + maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF) + +! Generalized Coupled-cluster module + + implicit none + include 'parameters.h' + +! Input variables + + logical :: dotest + + logical :: doCCD + logical :: dopCCD + logical :: doDCD + logical :: doCCSD + logical :: doCCSDT + logical :: dodrCCD + logical :: dorCCD + logical :: docrCCD + logical :: dolCCD + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + double precision,intent(in) :: ENuc + double precision,intent(in) :: EGHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: start_CC ,end_CC ,t_CC + +!------------------------------------------------------------------------ +! Perform CCD calculation +!------------------------------------------------------------------------ + + if(doCCD) then + + call wall_time(start_CC) + call GCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF) + call wall_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CC,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform DCD calculation +!------------------------------------------------------------------------ + + if(doDCD) then + + call wall_time(start_CC) +! call DCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR, & +! ERI,ENuc,EGHF,eHF) + call wall_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for DCD = ',t_CC,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform CCSD or CCSD(T) calculation +!------------------------------------------------------------------------ + + if(doCCSDT) doCCSD = .true. + + if(doCCSD) then + + call wall_time(start_CC) + call GCCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF) + call wall_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CC,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform direct ring CCD calculation +!------------------------------------------------------------------------ + + if(dodrCCD) then + + call wall_time(start_CC) + call drGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF) + call wall_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for direct ring CCD = ',t_CC,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform ring CCD calculation +!------------------------------------------------------------------------ + + if(dorCCD) then + + call wall_time(start_CC) + call rGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF) + call wall_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for rCCD = ',t_CC,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform crossed-ring CCD calculation +!------------------------------------------------------------------------ + + if(docrCCD) then + + call wall_time(start_CC) + call crGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF) + call wall_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CC,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform ladder CCD calculation +!------------------------------------------------------------------------ + + if(dolCCD) then + + call wall_time(start_CC) + call lGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF) + call wall_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ladder CCD = ',t_CC,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform pair CCD calculation +!------------------------------------------------------------------------ + + if(dopCCD) then + + call wall_time(start_CC) +! call pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF) + call wall_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pair CCD = ',t_CC,' seconds' + write(*,*) + + end if + +end subroutine diff --git a/src/CC/GCCD.f90 b/src/CC/GCCD.f90 new file mode 100644 index 0000000..2f5869a --- /dev/null +++ b/src/CC/GCCD.f90 @@ -0,0 +1,262 @@ +subroutine GCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF) + +! Generalized CCD module + + implicit none + +! Input variables + + logical,intent(in) :: dotest + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: ENuc + double precision,intent(in) :: EGHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: nSCF + double precision :: Conv + double precision :: EcMP2,EcMP3,EcMP4 + double precision :: ECC,EcCC + 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 :: OVVO(:,:,:,:) + + 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 :: r(:,:,:,:) + double precision,allocatable :: t(:,:,:,:) + + integer :: n_diis + double precision :: rcond + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: t_diis(:,:) + + logical :: do_EE_EOM_CC_1h1p = .false. + logical :: do_EA_EOM_CC_1p = .false. + logical :: do_IP_EOM_CC_1h = .false. + logical :: do_DEA_EOM_CC_2p = .false. + logical :: do_DIP_EOM_CC_2h = .false. + +! Hello world + + write(*,*) + write(*,*)'*******************************' + write(*,*)'* Generalized CCD Calculation *' + write(*,*)'*******************************' + write(*,*) + +! Antysymmetrize ERIs + + allocate(dbERI(nBas,nBas,nBas,nBas)) + + call antisymmetrize_ERI(2,nBas,ERI,dbERI) + +! Form energy denominator + + allocate(eO(nO-nC),eV(nV-nR)) + allocate(delta_OOVV(nO-nC,nO-nC,nV-nR,nV-nR)) + + eO(:) = eHF(nC+1:nO) + eV(:) = eHF(nO+1:nBas-nR) + + call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV) + +! Create integral batches + + allocate(OOOO(nO-nC,nO-nC,nO-nC,nO-nC),OOVV(nO-nC,nO-nC,nV-nR,nV-nR), & + OVOV(nO-nC,nV-nR,nO-nC,nV-nR),VVVV(nV-nR,nV-nR,nV-nR,nV-nR)) + + OOOO(:,:,:,:) = dbERI(nC+1:nO ,nC+1:nO ,nC+1:nO ,nC+1:nO ) + OOVV(:,:,:,:) = dbERI(nC+1:nO ,nC+1:nO ,nO+1:nBas-nR,nO+1:nBas-nR) + OVOV(:,:,:,:) = dbERI(nC+1:nO ,nO+1:nBas-nR,nC+1:nO ,nO+1:nBas-nR) + VVVV(:,:,:,:) = dbERI(nO+1:nBas-nR,nO+1:nBas-nR,nO+1:nBas-nR,nO+1:nBas-nR) + + allocate(OVVO(nO-nC,nV-nR,nV-nR,nO-nC)) + OVVO(:,:,:,:) = dbERI(nC+1:nO,nO+1:nBas-nR,nO+1:nBas-nR,nC+1:nO) + + deallocate(dbERI) + +! MP2 guess amplitudes + + allocate(t(nO-nC,nO-nC,nV-nR,nV-nR)) + + t(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) + + call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t,EcMP2) + EcMP4 = 0d0 + +! Memory allocation for DIIS + + allocate(error_diis((nO-nR)**2*(nV-nR)**2,max_diis),t_diis((nO-nR)**2*(nV-nR)**2,max_diis)) + +! Initialization + + allocate(r(nO-nC,nO-nC,nV-nR,nV-nR),u(nO-nC,nO-nC,nV-nR,nV-nR),v(nO-nC,nO-nC,nV-nR,nV-nR)) + allocate(X1(nO-nC,nO-nC,nO-nC,nO-nC),X2(nV-nR,nV-nR),X3(nO-nC,nO-nC),X4(nO-nC,nO-nC,nV-nR,nV-nR)) + + Conv = 1d0 + nSCF = 0 + + n_diis = 0 + t_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + +!------------------------------------------------------------------------ +! Main SCF loop +!------------------------------------------------------------------------ + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)'| GCCD calculation |' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(GCCD)','|','Ec(GCCD)','|','Conv','|' + write(*,*)'----------------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + +! Increment + + nSCF = nSCF + 1 + +! Form linear array + + call form_u(nC,nO,nV,nR,OOOO,VVVV,OVOV,t,u) + +! Form interemediate arrays + + call form_X(nC,nO,nV,nR,OOVV,t,X1,X2,X3,X4) + +! Form quadratic array + + call form_v(nC,nO,nV,nR,X1,X2,X3,X4,t,v) + +! Compute residual + + r(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t(:,:,:,:) + u(:,:,:,:) + v(:,:,:,:) + +! Check convergence + + Conv = maxval(abs(r(:,:,:,:))) + +! Update amplitudes + + t(:,:,:,:) = t(:,:,:,:) - r(:,:,:,:)/delta_OOVV(:,:,:,:) + +! Compute correlation energy + + call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t,EcCC) + + if(nSCF == 1) call MP3_correlation_energy(nC,nO,nV,nR,OOVV,t,v,delta_OOVV,EcMP3) + +! Dump results + + ECC = EGHF + EcCC + + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,(nO-nC)**2*(nV-nR)**2,(nO-nC)**2*(nV-nR)**2,n_diis,error_diis,t_diis,-r/delta_OOVV,t) + + ! Reset DIIS if required + + if(abs(rcond) < 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,'|',ECC+ENuc,'|',EcCC,'|',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(*,*)' GCCD energy ' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A30,1X,F15.10)')' E(GCCD) = ',ECC + write(*,'(1X,A30,1X,F15.10)')' Ec(GCCD) = ',EcCC + write(*,*)'----------------------------------------------------' + write(*,*) + +! Moller-Plesset energies + + write(*,*) + write(*,'(1X,A15,1X,F10.6)') 'Ec(GMP2) = ',EcMP2 + write(*,'(1X,A15,1X,F10.6)') 'Ec(GMP3) = ',EcMP3 + write(*,'(1X,A15,1X,F10.6)') 'Ec(GMP4-SDQ) = ',EcMP4 + write(*,*) + +!------------------------------------------------------------------------ +! EOM section +!------------------------------------------------------------------------ + +! EE-EOM-CCD (1h1p) + + if(do_EE_EOM_CC_1h1p) call EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t) + +! EA-EOM (1p) + +! if(do_EA-EOM-CC_1p) call EA-EOM-CCD_1p() + +! IP-EOM-CCD(1h) + +! if(do_IP-EOM-CC_1h) call IP-EOM-CCD_1h() + +! DEA-EOM (2p) + + if(do_DEA_EOM_CC_2p) call DEA_EOM_CCD_2p(nC,nO,nV,nR,eV,OOVV,VVVV,t) + +! DIP-EOM-CCD(2h) + + if(do_DIP_EOM_CC_2h) call DIP_EOM_CCD_2h(nC,nO,nV,nR,eO,OOVV,OOOO,t) + +! Testing zone + + if(dotest) then + + call dump_test_value('G','GCCD correlation energy',EcCC) + + end if + +end subroutine diff --git a/src/CC/GCCSD.f90 b/src/CC/GCCSD.f90 new file mode 100644 index 0000000..31bc9ae --- /dev/null +++ b/src/CC/GCCSD.f90 @@ -0,0 +1,291 @@ +subroutine GCCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF) + +! Generalized CCSD module + + implicit none + +! Input variables + + logical,intent(in) :: dotest + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + + logical,intent(in) :: doCCSDT + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: ENuc + double precision,intent(in) :: EGHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: start_CCSDT,end_CCSDT,t_CCSDT + integer :: nSCF + double precision :: Conv + double precision :: EcMP2 + double precision :: ECC + double precision :: EcCC + double precision :: EcCCT + + double precision,allocatable :: dbERI(:,:,:,:) + double precision,allocatable :: delta_OV(:,:) + double precision,allocatable :: delta_OOVV(:,:,:,:) + + double precision,allocatable :: OOOO(:,:,:,:) + double precision,allocatable :: OOOV(:,:,:,:) + double precision,allocatable :: OVOO(:,:,:,:) + double precision,allocatable :: VOOO(:,:,:,:) + double precision,allocatable :: OOVV(:,:,:,:) + double precision,allocatable :: OVVO(:,:,:,:) + double precision,allocatable :: OVVV(:,:,:,:) + double precision,allocatable :: VOVV(:,:,:,:) + double precision,allocatable :: VVVO(:,:,:,:) + double precision,allocatable :: VVVV(:,:,:,:) + + double precision,allocatable :: eO(:) + double precision,allocatable :: eV(:) + double precision,allocatable :: hvv(:,:) + double precision,allocatable :: hoo(:,:) + double precision,allocatable :: hvo(:,:) + double precision,allocatable :: gvv(:,:) + double precision,allocatable :: goo(:,:) + double precision,allocatable :: aoooo(:,:,:,:) + double precision,allocatable :: bvvvv(:,:,:,:) + double precision,allocatable :: hovvo(:,:,:,:) + + double precision,allocatable :: r1(:,:) + double precision,allocatable :: r2(:,:,:,:) + + double precision,allocatable :: t1(:,:) + double precision,allocatable :: t2(:,:,:,:) + double precision,allocatable :: tau(:,:,:,:) + + integer :: n_diis + double precision :: rcond1 + double precision :: rcond2 + double precision,allocatable :: err1_diis(:,:) + double precision,allocatable :: err2_diis(:,:) + double precision,allocatable :: t1_diis(:,:) + double precision,allocatable :: t2_diis(:,:) + +! Hello world + + write(*,*) + write(*,*)'********************************' + write(*,*)'* Generalized CCSD Calculation *' + write(*,*)'********************************' + write(*,*) + +! Antysymmetrize ERIs + + allocate(dbERI(nBas,nBas,nBas,nBas)) + + call antisymmetrize_ERI(2,nBas,ERI,dbERI) + +! Form energy denominator + + allocate(eO(nO),eV(nV)) + allocate(delta_OV(nO,nV),delta_OOVV(nO,nO,nV,nV)) + + eO(:) = eHF(1:nO) + eV(:) = eHF(nO+1:nBas) + + call form_delta_OV(nC,nO,nV,nR,eO,eV,delta_OV) + call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV) + +! Create integral batches + + allocate(OOOO(nO,nO,nO,nO), & + OOOV(nO,nO,nO,nV),OVOO(nO,nV,nO,nO),VOOO(nV,nO,nO,nO), & + OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO), & + OVVV(nO,nV,nV,nV),VOVV(nV,nO,nV,nV),VVVO(nV,nV,nV,nO), & + VVVV(nV,nV,nV,nV)) + + OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO ) + OOOV(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO ,nO+1: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) + +! MP2 guess amplitudes + + allocate(t1(nO,nV),t2(nO,nO,nV,nV),tau(nO,nO,nV,nV)) + + t1(:,:) = 0d0 + t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) + call form_tau(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 + +! Initialization + + allocate(hvv(nV,nV),hoo(nO,nO),hvo(nV,nO), & + gvv(nV,nV),goo(nO,nO), & + aoooo(nO,nO,nO,nO),bvvvv(nV,nV,nV,nV),hovvo(nO,nV,nV,nO), & + r1(nO,nV),r2(nO,nO,nV,nV)) + +! Memory allocation for DIIS + + allocate(err1_diis(nO*nV ,max_diis),t1_diis(nO*nV ,max_diis), & + err2_diis(nO*nO*nV*nV,max_diis),t2_diis(nO*nO*nV*nV,max_diis)) + + Conv = 1d0 + nSCF = 0 + ECC = EGHF + + n_diis = 0 + t1_diis(:,:) = 0d0 + t2_diis(:,:) = 0d0 + err1_diis(:,:) = 0d0 + err2_diis(:,:) = 0d0 + +!------------------------------------------------------------------------ +! Main SCF loop +!------------------------------------------------------------------------ + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)'| GCCSD calculation |' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(GCCSD)','|','Ec(GCCSD)','|','Conv','|' + write(*,*)'----------------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + +! Increment + + nSCF = nSCF + 1 + +! Scuseria Eqs. (5), (6) and (7) + + 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(nC,nO,nV,nR,hvv,hoo,VOVV,OOOV,t1,gvv,goo) + + call form_abh(nC,nO,nV,nR,OOOO,OVOO,OOVV,VVVV,VOVV,OVVO,OVVV,t1,tau,aoooo,bvvvv,hovvo) + +! Compute residuals + + call form_r1(nC,nO,nV,nR,OVVO,OVVV,OOOV,hvv,hoo,hvo,t1,t2,tau,r1) + + 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(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(nC,nO,nV,nR,t1,t2,tau) + +! Compute correlation energy + + call CCSD_correlation_energy(nC,nO,nV,nR,OOVV,tau,EcCC) + +! Dump results + + ECC = EGHF + EcCC + + ! 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) + + ! Reset DIIS if required + +! 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,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|' + + end do + write(*,*)'----------------------------------------------------' +!------------------------------------------------------------------------ +! End of SCF loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + end if + + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)' CCSD energy ' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A20,1X,F15.10)')' E(CCSD) = ',ENuc+ECC + write(*,'(1X,A20,1X,F10.6)')' Ec(CCSD) = ',EcCC + write(*,*)'----------------------------------------------------' + write(*,*) + +! Deallocate memory + + deallocate(hvv,hoo,hvo, & + delta_OV,delta_OOVV, & + gvv,goo, & + aoooo,bvvvv,hovvo, & + tau, & + r1,r2) + +!------------------------------------------------------------------------ +! (T) correction +!------------------------------------------------------------------------ + if(doCCSDT) then + + call wall_time(start_CCSDT) + call CCSDT(nC,nO,nV,nR,eO,eV,OOVV,VVVO,VOOO,t1,t2,EcCCT) + call wall_time(end_CCSDT) + + t_CCSDT = end_CCSDT - start_CCSDT + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for (T) = ',t_CCSDT,' seconds' + write(*,*) + + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)' GCCSD(T) energy ' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A20,1X,F15.10)')' E(GCCSD(T)) = ',ENuc + ECC + EcCCT + write(*,'(1X,A20,1X,F10.6)')' Ec(GCCSD(T)) = ',EcCC + EcCCT + write(*,*)'----------------------------------------------------' + write(*,*) + + end if + +! Testing zone + + if(dotest) then + + call dump_test_value('R','GCCSD correlation energy',EcCC) + + end if + +end subroutine diff --git a/src/CC/RCC.f90 b/src/CC/RCC.f90 index 5f2e531..142d8fc 100644 --- a/src/CC/RCC.f90 +++ b/src/CC/RCC.f90 @@ -1,5 +1,5 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & - maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) + maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) ! Coupled-cluster module @@ -30,8 +30,8 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d integer,intent(in) :: nV(nspin) integer,intent(in) :: nR(nspin) double precision,intent(in) :: ENuc - double precision,intent(in) :: EHF - double precision,intent(in) :: epsHF(nBas) + double precision,intent(in) :: ERHF + double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables @@ -45,7 +45,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(doCCD) then call wall_time(start_CC) - call CCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) + call CCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -62,7 +62,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d call wall_time(start_CC) call DCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR, & - ERI,ENuc,EHF,epsHF) + ERI,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -80,7 +80,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(doCCSD) then call wall_time(start_CC) - call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) + call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -96,7 +96,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(dodrCCD) then call wall_time(start_CC) - call drCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) + call drCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -112,7 +112,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(dorCCD) then call wall_time(start_CC) - call rCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF,epsHF) + call rCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -128,7 +128,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(docrCCD) then call wall_time(start_CC) - call crCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) + call crCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -144,7 +144,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(dolCCD) then call wall_time(start_CC) - call lCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) + call lCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -160,8 +160,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(dopCCD) then call wall_time(start_CC) - call pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) - + call pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC diff --git a/src/CC/crGCCD.f90 b/src/CC/crGCCD.f90 new file mode 100644 index 0000000..b1058aa --- /dev/null +++ b/src/CC/crGCCD.f90 @@ -0,0 +1,189 @@ +subroutine crGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + +! Generalized crossed-ring CCD module + + implicit none + +! Input variables + + logical,intent(in) :: dotest + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: ENuc,ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: nSCF + double precision :: Conv + double precision :: EcMP2 + double precision :: ECC,EcCC + double precision,allocatable :: dbERI(:,:,:,:) + + double precision,allocatable :: eO(:) + double precision,allocatable :: eV(:) + double precision,allocatable :: delta_OOVV(:,:,:,:) + + double precision,allocatable :: OOVV(:,:,:,:) + double precision,allocatable :: OVOV(:,:,:,:) + + double precision,allocatable :: r2(:,:,:,:) + double precision,allocatable :: t2(:,:,:,:) + + integer :: n_diis + double precision :: rcond + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: t_diis(:,:) + +! Hello world + + write(*,*) + write(*,*)'*********************************' + write(*,*)'* Generalized crCCD Calculation *' + write(*,*)'*********************************' + write(*,*) + +! Antysymmetrize ERIs + + allocate(dbERI(nBas,nBas,nBas,nBas)) + + call antisymmetrize_ERI(2,nBas,ERI,dbERI) + +! Form energy denominator + + allocate(eO(nO),eV(nV)) + allocate(delta_OOVV(nO,nO,nV,nV)) + + eO(:) = eHF(1:nO) + eV(:) = eHF(nO+1:nBas) + + call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV) + +! Create integral batches + + allocate(OOVV(nO,nO,nV,nV),OVOV(nO,nV,nO,nV)) + + OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas) + OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas, 1:nO ,nO+1:nBas) + + deallocate(dbERI) + +! MP2 guess amplitudes + + allocate(t2(nO,nO,nV,nV)) + + t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) + + call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcMP2) + +! Memory allocation for DIIS + + allocate(error_diis(nO*nO*nV*nV,max_diis),t_diis(nO*nO*nV*nV,max_diis)) + +! Initialization + + allocate(r2(nO,nO,nV,nV)) + + Conv = 1d0 + nSCF = 0 + + n_diis = 0 + t_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + +!------------------------------------------------------------------------ +! Main SCF loop +!------------------------------------------------------------------------ + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)'| crossed-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_crossed_ring_r(nC,nO,nV,nR,OVOV,OOVV,t2,r2) + + r2(:,:,:,:) = OOVV(:,:,:,:) - delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:) + +! Check convergence + + 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 + + call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCC) + +! Dump results + + ECC = ERHF + EcCC + + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,error_diis,t_diis,-r2/delta_OOVV,t2) + + ! Reset DIIS if required + + if(abs(rcond) < 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,'|',ECC+ENuc,'|',EcCC,'|',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(*,*)' crossed-ring CCD energy ' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A30,1X,F15.10)')' E(crCCD) = ',ECC + write(*,'(1X,A30,1X,F15.10)')' Ec(crCCD) = ',EcCC + write(*,*)'----------------------------------------------------' + write(*,*) + + if(dotest) then + + call dump_test_value('R','crCCD correlation energy',EcCC) + + end if + +end subroutine diff --git a/src/CC/drGCCD.f90 b/src/CC/drGCCD.f90 new file mode 100644 index 0000000..0692d63 --- /dev/null +++ b/src/CC/drGCCD.f90 @@ -0,0 +1,181 @@ +subroutine drGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + +! Generalized Direct ring CCD module + + implicit none + +! Input variables + + logical,intent(in) :: dotest + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: ENuc,ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: nSCF + double precision :: Conv + double precision :: EcMP2 + double precision :: ECC,EcCC + + double precision,allocatable :: eO(:) + double precision,allocatable :: eV(:) + double precision,allocatable :: delta_OOVV(:,:,:,:) + + double precision,allocatable :: OOVV(:,:,:,:) + double precision,allocatable :: OVVO(:,:,:,:) + + double precision,allocatable :: r2(:,:,:,:) + double precision,allocatable :: t2(:,:,:,:) + + + integer :: n_diis + double precision :: rcond + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: t_diis(:,:) + +! Hello world + + write(*,*) + write(*,*)'*********************************' + write(*,*)'* Generalized drCCD Calculation *' + write(*,*)'*********************************' + write(*,*) +! Form energy denominator + + allocate(eO(nO),eV(nV)) + allocate(delta_OOVV(nO,nO,nV,nV)) + + eO(:) = eHF(1:nO) + eV(:) = eHF(nO+1:nBas) + + call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV) + +! Create integral batches + + allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO)) + + OOVV(:,:,:,:) = ERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas) + OVVO(:,:,:,:) = ERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO ) + +! MP2 guess amplitudes + + allocate(t2(nO,nO,nV,nV)) + + t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) + + call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcMP2) + +! Memory allocation for DIIS + + allocate(error_diis(nO*nO*nV*nV,max_diis),t_diis(nO*nO*nV*nV,max_diis)) + +! Initialization + + allocate(r2(nO,nO,nV,nV)) + + Conv = 1d0 + nSCF = 0 + + n_diis = 0 + t_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + +!------------------------------------------------------------------------ +! Main SCF loop +!------------------------------------------------------------------------ + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)'| direct ring GCCD calculation |' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(GCCD)','|','Ec(GCCD)','|','Conv','|' + write(*,*)'----------------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + +! Increment + + nSCF = nSCF + 1 + +! Compute residual + + call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2) + + r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:) + +! Check convergence + + 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 + + call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCC) + EcCC = 2d0*EcCC + +! Dump results + + ECC = ERHF + EcCC + + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,error_diis,t_diis,-r2/delta_OOVV,t2) + + ! Reset DIIS if required + + if(abs(rcond) < 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,'|',ECC+ENuc,'|',EcCC,'|',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 GCCD energy ' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A30,1X,F15.10)')' E(drGCCD) = ',ECC + write(*,'(1X,A30,1X,F15.10)')' Ec(drGCCD) = ',EcCC + write(*,*)'----------------------------------------------------' + write(*,*) + + if(dotest) then + + call dump_test_value('G','drCCD correlation energy',EcCC) + + end if + +end subroutine diff --git a/src/CC/lGCCD.f90 b/src/CC/lGCCD.f90 new file mode 100644 index 0000000..3c4a3d3 --- /dev/null +++ b/src/CC/lGCCD.f90 @@ -0,0 +1,202 @@ +subroutine lGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + +! Generalized Ladder CCD module + + implicit none + +! Input variables + + logical,intent(in) :: dotest + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: ENuc,ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: nSCF + double precision :: Conv + double precision :: EcMP2 + double precision :: ECC,EcCC + 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(:,:,:,:) + + integer :: n_diis + double precision :: rcond + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: t_diis(:,:) + +! Hello world + + write(*,*) + write(*,*)'********************************' + write(*,*)'* Generalized lCCD Calculation *' + write(*,*)'********************************' + write(*,*) + +! Antysymmetrize ERIs + + allocate(dbERI(nBas,nBas,nBas,nBas)) + + call antisymmetrize_ERI(2,nBas,ERI,dbERI) + +! Form energy denominator + + allocate(eO(nO),eV(nV)) + allocate(delta_OOVV(nO,nO,nV,nV)) + + eO(:) = eHF(1:nO) + eV(:) = eHF(nO+1:nBas) + + call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV) + +! 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: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) + +! MP2 guess amplitudes + + allocate(t2(nO,nO,nV,nV)) + + t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) + + call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcMP2) + +! Memory allocation for DIIS + + allocate(error_diis(nO*nO*nV*nV,max_diis),t_diis(nO*nO*nV*nV,max_diis)) + +! 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 + + n_diis = 0 + t_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + +!------------------------------------------------------------------------ +! Main SCF loop +!------------------------------------------------------------------------ + write(*,*) + write(*,*)'----------------------------------------------------' + 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','|' + write(*,*)'----------------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + +! Increment + + nSCF = nSCF + 1 + +! Compute residual + + 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(nC+1:nO,nC+1:nO,1:nV-nR,1:nV-nR))) + +! Update amplitudes + + t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:) + +! Compute correlation energy + + call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCC) + +! Dump results + + ECC = ERHF + EcCC + + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,error_diis,t_diis,-r2/delta_OOVV,t2) + + ! Reset DIIS if required + + if(abs(rcond) < 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,'|',ECC+ENuc,'|',EcCC,'|',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(*,*)' ladder CCD energy ' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A30,1X,F15.10)')' E(lCCD) = ',ECC + write(*,'(1X,A30,1X,F15.10)')' Ec(lCCD) = ',EcCC + write(*,*)'----------------------------------------------------' + write(*,*) + + if(dotest) then + + call dump_test_value('R','lCCD correlation energy',EcCC) + + end if + +end subroutine diff --git a/src/CC/rCCD.f90 b/src/CC/rCCD.f90 index ce61a66..bbd6444 100644 --- a/src/CC/rCCD.f90 +++ b/src/CC/rCCD.f90 @@ -1,4 +1,4 @@ -subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF,eGW) +subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF) ! Ring CCD module @@ -20,7 +20,6 @@ subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENu double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nBasin) - double precision,intent(in) :: eGW(nBasin) double precision,intent(in) :: ERI(nBasin,nBasin,nBasin,nBasin) ! Local variables @@ -35,7 +34,6 @@ subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENu double precision :: EcMP2 double precision :: ECC,EcCC double precision,allocatable :: seHF(:) - double precision,allocatable :: seGW(:) double precision,allocatable :: sERI(:,:,:,:) double precision,allocatable :: dbERI(:,:,:,:) @@ -72,10 +70,9 @@ subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENu nV = 2*nVin nR = 2*nRin - allocate(seHF(nBas),seGW(nBas),sERI(nBas,nBas,nBas,nBas)) + allocate(seHF(nBas),sERI(nBas,nBas,nBas,nBas)) call spatial_to_spin_MO_energy(nBasin,eHF,nBas,seHF) - call spatial_to_spin_MO_energy(nBasin,eGW,nBas,seGW) call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI) ! Antysymmetrize ERIs @@ -91,12 +88,12 @@ subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENu allocate(eO(nO),eV(nV)) allocate(delta_OOVV(nO,nO,nV,nV)) - eO(:) = seGW(1:nO) - eV(:) = seGW(nO+1:nBas) + eO(:) = seHF(1:nO) + eV(:) = seHF(nO+1:nBas) call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV) -! deallocate(seHF,seGW) + deallocate(seHF) ! Create integral batches diff --git a/src/CC/rGCCD.f90 b/src/CC/rGCCD.f90 new file mode 100644 index 0000000..6aea413 --- /dev/null +++ b/src/CC/rGCCD.f90 @@ -0,0 +1,208 @@ +subroutine rGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + +! Generalized ring CCD module + + implicit none + +! Input variables + + logical,intent(in) :: dotest + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: nSCF + double precision :: Conv + double precision :: EcMP2 + double precision :: ECC,EcCC + double precision,allocatable :: dbERI(:,:,:,:) + + double precision,allocatable :: eO(:) + double precision,allocatable :: eV(:) + double precision,allocatable :: delta_OOVV(:,:,:,:) + + double precision,allocatable :: OOVV(:,:,:,:) + double precision,allocatable :: OVVO(:,:,:,:) + + double precision,allocatable :: r(:,:,:,:) + double precision,allocatable :: t(:,:,:,:) + + integer :: n_diis + double precision :: rcond + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: t_diis(:,:) + + logical :: do_EE_EOM_CC_1h1p = .true. + +! Hello world + + write(*,*) + write(*,*)'********************************' + write(*,*)'* Generalized rCCD Calculation *' + write(*,*)'********************************' + write(*,*) + +! Antysymmetrize ERIs + + allocate(dbERI(nBas,nBas,nBas,nBas)) + + call antisymmetrize_ERI(2,nBas,ERI,dbERI) + +! Form energy denominator + + allocate(eO(nO),eV(nV)) + allocate(delta_OOVV(nO,nO,nV,nV)) + + eO(:) = eHF(1:nO) + eV(:) = eHF(nO+1:nBas) + + call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV) + +! Create integral batches + + allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,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 ) + + deallocate(dbERI) + +! MP2 guess amplitudes + + allocate(t(nO,nO,nV,nV)) + + t(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) + + call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t,EcMP2) + +! Memory allocation for DIIS + + allocate(error_diis(nO*nO*nV*nV,max_diis),t_diis(nO*nO*nV*nV,max_diis)) + +! Initialization + + allocate(r(nO,nO,nV,nV)) + + Conv = 1d0 + nSCF = 0 + + n_diis = 0 + t_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + +!------------------------------------------------------------------------ +! Main SCF loop +!------------------------------------------------------------------------ + write(*,*) + write(*,*)'----------------------------------------------------' + 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','|' + write(*,*)'----------------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + +! Increment + + nSCF = nSCF + 1 + +! Compute residual + + call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t,r) + + r(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t(:,:,:,:) + r(:,:,:,:) + +! Check convergence + + Conv = maxval(abs(r(nC+1:nO,nC+1:nO,1:nV-nR,1:nV-nR))) + +! Update amplitudes + + t(:,:,:,:) = t(:,:,:,:) - r(:,:,:,:)/delta_OOVV(:,:,:,:) + +! Compute correlation energy + + call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t,EcCC) + +! Dump results + + ECC = ERHF + EcCC + + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,error_diis,t_diis,-r/delta_OOVV,t) + + ! Reset DIIS if required + + if(abs(rcond) < 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,'|',ECC+ENuc,'|',EcCC,'|',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(*,*)' ring CCD energy ' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A30,1X,F15.10)')' E(rCCD) = ',ECC + write(*,'(1X,A30,1X,F15.10)')' Ec(rCCD) = ',EcCC + write(*,*)'----------------------------------------------------' + write(*,*) + +! write(*,*) +! write(*,*)'----------------------------------------------------' +! write(*,*)' ring CCD amplitudes ' +! write(*,*)'----------------------------------------------------' +! call matout(nO*nO,nV*nV,t) +! write(*,*) + +!------------------------------------------------------------------------ +! EOM section +!------------------------------------------------------------------------ + +! EE-EOM-CCD (1h1p) + + if(do_EE_EOM_CC_1h1p) call EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t) + + + if(dotest) then + + call dump_test_value('R','rCCD correlation energy',EcCC) + + end if + +end subroutine diff --git a/src/HF/GHF_search.f90 b/src/HF/GHF_search.f90 index e8082f7..5d43822 100644 --- a/src/HF/GHF_search.f90 +++ b/src/HF/GHF_search.f90 @@ -98,7 +98,7 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu !---------------------! call wall_time(start_HF) - call GHF(maxSCF,thresh,max_diis,guess,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & + call GHF(.false.,maxSCF,thresh,max_diis,guess,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & nBas,nBas2,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,e,c,P) call wall_time(end_HF) diff --git a/src/HF/RHF_search.f90 b/src/HF/RHF_search.f90 index 337b244..eb12b94 100644 --- a/src/HF/RHF_search.f90 +++ b/src/HF/RHF_search.f90 @@ -91,7 +91,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN !---------------------! call wall_time(start_HF) - call RHF(maxSCF,thresh,max_diis,guess,level_shift,nNuc,ZNuc,rNuc,ENuc, & + call RHF(.false.,maxSCF,thresh,max_diis,guess,level_shift,nNuc,ZNuc,rNuc,ENuc, & nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,e,c,P) call wall_time(end_HF) diff --git a/src/HF/UHF_search.f90 b/src/HF/UHF_search.f90 index 0ea95c2..ae4de45 100644 --- a/src/HF/UHF_search.f90 +++ b/src/HF/UHF_search.f90 @@ -102,7 +102,7 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu !---------------------! call wall_time(start_HF) - call UHF(maxSCF,thresh,max_diis,guess,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & + call UHF(.false.,maxSCF,thresh,max_diis,guess,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,e,c,P) call wall_time(end_HF) diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90 index 4514402..4635b5c 100644 --- a/src/QuAcK/GQuAcK.f90 +++ b/src/QuAcK/GQuAcK.f90 @@ -1,9 +1,11 @@ -subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,doppRPA, & - doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & - nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & - maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, & - TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & - maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & +subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,doppRPA, & + doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & + nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & + maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, & + maxSCF_CC,max_diis_CC,thresh_CC, & + TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & + maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) implicit none @@ -16,6 +18,8 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp logical,intent(in) :: dosearch logical,intent(in) :: doMP2 logical,intent(in) :: doMP3 + logical,intent(in) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT + logical,intent(in) :: dodrCCD,dorCCD,docrCCD,dolCCD logical,intent(in) :: dophRPA,dophRPAx,doppRPA logical,intent(in) :: doG0F2,doevGF2,doqsGF2 logical,intent(in) :: doG0W0,doevGW,doqsGW @@ -43,6 +47,9 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp logical,intent(in) :: reg_MP + integer,intent(in) :: maxSCF_CC,max_diis_CC + double precision,intent(in) :: thresh_CC + logical,intent(in) :: TDA integer,intent(in) :: maxSCF_GF,max_diis_GF @@ -60,18 +67,19 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp ! Local variables - logical :: doMP,doRPA,doGF,doGW + logical :: doMP,doCC,doRPA,doGF,doGW double precision :: start_HF ,end_HF ,t_HF double precision :: start_stab ,end_stab ,t_stab double precision :: start_AOtoMO ,end_AOtoMO ,t_AOtoMO double precision :: start_MP ,end_MP ,t_MP + double precision :: start_CC ,end_CC ,t_CC double precision :: start_RPA ,end_RPA ,t_RPA double precision :: start_GF ,end_GF ,t_GF double precision :: start_GW ,end_GW ,t_GW - double precision,allocatable :: cHF(:,:),epsHF(:),PHF(:,:) - double precision :: EHF + double precision,allocatable :: cHF(:,:),eHF(:),PHF(:,:) + double precision :: EGHF double precision,allocatable :: dipole_int_MO(:,:,:) double precision,allocatable :: ERI_MO(:,:,:,:) double precision,allocatable :: ERI_tmp(:,:,:,:) @@ -93,7 +101,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp nBas2 = 2*nBas - allocate(cHF(nBas2,nBas2),epsHF(nBas2),PHF(nBas2,nBas2), & + allocate(cHF(nBas2,nBas2),eHF(nBas2),PHF(nBas2,nBas2), & dipole_int_MO(nBas2,nBas2,ncart),ERI_MO(nBas2,nBas2,nBas2,nBas2)) !---------------------! @@ -104,7 +112,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp call wall_time(start_HF) call GHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nBas2,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF) + nBas,nBas2,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EGHF,eHF,cHF,PHF) call wall_time(end_HF) t_HF = end_HF - start_HF @@ -165,7 +173,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp if(dostab) then call wall_time(start_stab) - call GHF_stability(nBas2,nC,nO,nV,nR,nS,epsHF,ERI_MO) + call GHF_stability(nBas2,nC,nO,nV,nR,nS,eHF,ERI_MO) call wall_time(end_stab) t_stab = end_stab - start_stab @@ -178,7 +186,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp call wall_time(start_stab) call GHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nBas2,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF) + nBas,nBas2,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EGHF,eHF,cHF,PHF) call wall_time(end_stab) t_stab = end_stab - start_stab @@ -196,7 +204,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp if(doMP) then call wall_time(start_MP) - call GMP(dotest,doMP2,doMP3,reg_MP,nBas2,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) + call GMP(dotest,doMP2,doMP3,reg_MP,nBas2,nC,nO,nV,nR,ERI_MO,ENuc,EGHF,eHF) call wall_time(end_MP) t_MP = end_MP - start_MP @@ -205,6 +213,25 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp end if +!------------------------! +! Coupled-cluster module ! +!------------------------! + + doCC = doCCD .or. doCCSD .or. doCCSDT .or. dodrCCD .or. dorCCD .or. docrCCD .or. dolCCD + + if(doCC) then + + call wall_time(start_CC) + call GCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & + maxSCF_CC,thresh_CC,max_diis_CC,nBas2,nC,nO,nV,nR,ERI_MO,ENuc,EGHF,eHF) + call wall_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CC = ',t_CC,' seconds' + write(*,*) + + end if + !-----------------------------------! ! Random-phase approximation module ! !-----------------------------------! @@ -214,8 +241,8 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp if(doRPA) then call wall_time(start_RPA) - call GRPA(dotest,dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas2,nC,nO,nV,nR,nS,ENuc,EHF, & - ERI_MO,dipole_int_MO,epsHF,cHF,S) + call GRPA(dotest,dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF, & + ERI_MO,dipole_int_MO,eHF,cHF,S) call wall_time(end_RPA) t_RPA = end_RPA - start_RPA @@ -234,7 +261,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp call wall_time(start_GF) call GGF(dotest,doG0F2,doevGF2,doqsGF2,maxSCF_GF,thresh_GF,max_diis_GF,dophBSE,doppBSE,TDA,dBSE,dTDA,lin_GF,eta_GF,reg_GF, & - nNuc,ZNuc,rNuc,ENuc,nBas2,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF) + nNuc,ZNuc,rNuc,ENuc,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) call wall_time(end_GF) t_GF = end_GF - start_GF @@ -254,7 +281,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dopp call wall_time(start_GW) call GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT,exchange_kernel,doXBS, & dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc, & - nBas,nBas2,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF) + nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 8fb8048..2b9fcdb 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -226,11 +226,12 @@ program QuAcK !--------------------------! if(doGQuAcK) & - call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,doppRPA, & + call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,doppRPA, & doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, & - TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & + maxSCF_CC,max_diis_CC,thresh_CC,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index f2d1d97..f82e48c 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -90,8 +90,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d double precision :: start_GW ,end_GW ,t_GW double precision :: start_GT ,end_GT ,t_GT - double precision,allocatable :: cHF(:,:),epsHF(:),PHF(:,:) - double precision :: EHF + double precision,allocatable :: cHF(:,:),eHF(:),PHF(:,:) + double precision :: ERHF double precision,allocatable :: dipole_int_MO(:,:,:) double precision,allocatable :: ERI_MO(:,:,:,:) integer :: ixyz @@ -107,7 +107,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d ! Memory allocation ! !-------------------! - allocate(cHF(nBas,nBas),epsHF(nBas),PHF(nBas,nBas), & + allocate(cHF(nBas,nBas),eHF(nBas),PHF(nBas,nBas), & dipole_int_MO(nBas,nBas,ncart),ERI_MO(nBas,nBas,nBas,nBas)) !---------------------! @@ -118,7 +118,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_HF) call RHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF) + nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF) call wall_time(end_HF) t_HF = end_HF - start_HF @@ -131,7 +131,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_HF) call ROHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF) + nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF) call wall_time(end_HF) t_HF = end_HF - start_HF @@ -175,7 +175,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(dostab) then call wall_time(start_stab) - call RHF_stability(nBas,nC,nO,nV,nR,nS,epsHF,ERI_MO) + call RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_MO) call wall_time(end_stab) t_stab = end_stab - start_stab @@ -188,7 +188,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_stab) call RHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF) + nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF) call wall_time(end_stab) t_stab = end_stab - start_stab @@ -206,7 +206,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doMP) then call wall_time(start_MP) - call RMP(dotest,doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) + call RMP(dotest,doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_MP) t_MP = end_MP - start_MP @@ -226,7 +226,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_CC) call RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & - maxSCF_CC,thresh_CC,max_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) + maxSCF_CC,thresh_CC,max_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -245,7 +245,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_CI) call RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, & - epsHF,EHF,cHF,S) + eHF,ERHF,cHF,S) call wall_time(end_CI) t_CI = end_CI - start_CI @@ -264,7 +264,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_RPA) call RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, & - nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF,cHF,S) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF,cHF,S) call wall_time(end_RPA) t_RPA = end_RPA - start_RPA @@ -284,8 +284,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_GF) call RGF(dotest,doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, & dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF,eta_GF,reg_GF, & - nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO, & - dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF) + nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO, & + dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) call wall_time(end_GF) t_GF = end_GF - start_GF @@ -305,8 +305,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_GW) call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, & exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, & - lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, & - ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF) + lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, & + ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -326,8 +326,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_GT) call RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF_GT,thresh_GT,max_diis_GT,doACFDT, & exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & - lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, & - ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF) + lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, & + ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT diff --git a/src/QuAcK/UQuAcK.f90 b/src/QuAcK/UQuAcK.f90 index 4cd189f..7675dc1 100644 --- a/src/QuAcK/UQuAcK.f90 +++ b/src/QuAcK/UQuAcK.f90 @@ -88,8 +88,8 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do double precision :: start_GW ,end_GW ,t_GW double precision :: start_GT ,end_GT ,t_GT - double precision,allocatable :: cHF(:,:,:),epsHF(:,:),PHF(:,:,:) - double precision :: EHF + double precision,allocatable :: cHF(:,:,:),eHF(:,:),PHF(:,:,:) + double precision :: EUHF double precision,allocatable :: dipole_int_aa(:,:,:),dipole_int_bb(:,:,:) double precision,allocatable :: ERI_aaaa(:,:,:,:),ERI_aabb(:,:,:,:),ERI_bbbb(:,:,:,:) integer :: ixyz @@ -105,7 +105,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do ! Memory allocation ! !-------------------! - allocate(cHF(nBas,nBas,nspin),epsHF(nBas,nspin),PHF(nBas,nBas,nspin), & + allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),PHF(nBas,nBas,nspin), & dipole_int_aa(nBas,nBas,ncart),dipole_int_bb(nBas,nBas,ncart), & ERI_aaaa(nBas,nBas,nBas,nBas),ERI_aabb(nBas,nBas,nBas,nBas),ERI_bbbb(nBas,nBas,nBas,nBas)) @@ -117,7 +117,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do call wall_time(start_HF) call UHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF) + nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,eHF,cHF,PHF) call wall_time(end_HF) t_HF = end_HF - start_HF @@ -138,7 +138,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do ! write(*,*) ! call eDFT(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nC, & ! nO,nV,nR,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, & -! max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int_AO,EHF,epsHF,cHF,PHF,Vxc) +! max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int_AO,EUHF,eHF,cHF,PHF,Vxc) ! call wall_time(end_KS) ! t_KS = end_KS - start_KS @@ -191,7 +191,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do if(dostab) then call wall_time(start_stab) - call UHF_stability(nBas,nC,nO,nV,nR,nS,epsHF,ERI_aaaa,ERI_aabb,ERI_bbbb) + call UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb) call wall_time(end_stab) t_stab = end_stab - start_stab @@ -204,7 +204,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do call wall_time(start_stab) call UHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF) + nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,eHF,cHF,PHF) call wall_time(end_stab) @@ -223,7 +223,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do if(doMP) then call wall_time(start_MP) - call UMP(dotest,doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF) + call UMP(dotest,doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EUHF,eHF) call wall_time(end_MP) t_MP = end_MP - start_MP @@ -260,7 +260,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do call wall_time(start_CI) call UCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS, & - ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,EHF,cHF,S) + ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,eHF,EUHF,cHF,S) call wall_time(end_CI) t_CI = end_CI - start_CI @@ -279,7 +279,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do call wall_time(start_RPA) call URPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip, & - nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S) + nBas,nC,nO,nV,nR,nS,ENuc,EUHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,eHF,cHF,S) call wall_time(end_RPA) t_RPA = end_RPA - start_RPA @@ -299,8 +299,8 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do call wall_time(start_GF) call UGF(dotest,doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, & dophBSE,doppBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,lin_GF,eta_GF,reg_GF, & - nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, & - dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) + nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, & + dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) call wall_time(end_GF) t_GF = end_GF - start_GF @@ -320,8 +320,8 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do call wall_time(start_GW) call UGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, & exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & - lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, & - ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) + lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc, & + ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -341,8 +341,8 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do call wall_time(start_GT) call UGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF_GT,thresh_GT,max_diis_GT,doACFDT, & exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, & - lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, & - ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) + lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc, & + ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT