From 1f2f6634bb6e1f58a8ec890aa5e4aa601040dfde Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 27 Sep 2022 14:07:31 +0200 Subject: [PATCH] DIP- and DEA-EOM-CCD --- input/methods | 4 +- input/options | 4 +- src/CC/CCD.f90 | 60 ++++++++++++---- src/CC/DEA_EOM_CCD_2p.f90 | 116 ++++++++++++++++++++++++++++++ src/CC/DIP_EOM_CCD_2h.f90 | 116 ++++++++++++++++++++++++++++++ src/CC/rCCD.f90 | 13 +++- src/GT/static_Tmatrix_C_pp.f90 | 8 ++- src/GW/static_screening_WC_pp.f90 | 2 +- src/LR/print_excitation.f90 | 2 +- src/QuAcK/QuAcK.f90 | 50 +++++++------ src/RPA/ppRPA.f90 | 12 ++-- 11 files changed, 333 insertions(+), 54 deletions(-) create mode 100644 src/CC/DEA_EOM_CCD_2p.f90 create mode 100644 src/CC/DIP_EOM_CCD_2h.f90 diff --git a/input/methods b/input/methods index 511a496..5ef1b53 100644 --- a/input/methods +++ b/input/methods @@ -3,7 +3,7 @@ # MP2* MP3 MP2-F12 F F F # CCD pCCD DCD CCSD CCSD(T) - F F F F F + T F F F F # drCCD rCCD crCCD lCCD F F F F # CIS* CIS(D) CID CISD FCI @@ -15,7 +15,7 @@ # G0W0* evGW* qsGW* ufG0W0 ufGW F F F F F # G0T0 evGT qsGT - T F F + F F F # MCMP2 F # * unrestricted version available diff --git a/input/options b/input/options index b445287..5990b1c 100644 --- a/input/options +++ b/input/options @@ -5,7 +5,7 @@ # CC: maxSCF thresh DIIS n_diis 64 0.00001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - T T F T T + T T T T T # GF: maxSCF thresh DIIS n_diis lin eta renorm reg 256 0.00001 T 5 T 0.0 3 F # GW: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 reg @@ -15,6 +15,6 @@ # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn ppBSE - F F T F T + T F T F F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/CC/CCD.f90 b/src/CC/CCD.f90 index 7b69337..3653e7f 100644 --- a/src/CC/CCD.f90 +++ b/src/CC/CCD.f90 @@ -52,14 +52,20 @@ subroutine CCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ER double precision,allocatable :: u(:,:,:,:) double precision,allocatable :: v(:,:,:,:) - double precision,allocatable :: r2(:,:,:,:) - double precision,allocatable :: t2(:,:,:,:) + double precision,allocatable :: r(:,:,:,:) + double precision,allocatable :: t(:,:,:,:) integer :: n_diis,i,j,a,b 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 = .true. + logical :: do_DIP_EOM_CC_2h = .true. + ! Hello world write(*,*) @@ -123,11 +129,11 @@ subroutine CCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ER ! MP2 guess amplitudes - allocate(t2(nO-nC,nO-nC,nV-nR,nV-nR)) + allocate(t(nO-nC,nO-nC,nV-nR,nV-nR)) - t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) + t(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) - call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcMP2) + call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t,EcMP2) EcMP4 = 0d0 ! Memory allocation for DIIS @@ -136,7 +142,7 @@ subroutine CCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ER ! Initialization - allocate(r2(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(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 @@ -165,33 +171,33 @@ subroutine CCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ER ! Form linear array - call form_u(nC,nO,nV,nR,OOOO,VVVV,OVOV,t2,u) + call form_u(nC,nO,nV,nR,OOOO,VVVV,OVOV,t,u) ! Form interemediate arrays - call form_X(nC,nO,nV,nR,OOVV,t2,X1,X2,X3,X4) + 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,t2,v) + call form_v(nC,nO,nV,nR,X1,X2,X3,X4,t,v) ! Compute residual - r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + u(:,:,:,:) + v(:,:,:,:) + r(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t(:,:,:,:) + u(:,:,:,:) + v(:,:,:,:) ! Check convergence - Conv = maxval(abs(r2(:,:,:,:))) + Conv = maxval(abs(r(:,:,:,:))) ! Update amplitudes - t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:) + t(:,:,:,:) = t(:,:,:,:) - r(:,:,:,:)/delta_OOVV(:,:,:,:) ! Compute correlation energy - call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCCD) + call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t,EcCCD) - if(nSCF == 1) call MP3_correlation_energy(nC,nO,nV,nR,OOVV,t2,v,delta_OOVV,EcMP3) + if(nSCF == 1) call MP3_correlation_energy(nC,nO,nV,nR,OOVV,t,v,delta_OOVV,EcMP3) ! Dump results @@ -200,7 +206,7 @@ subroutine CCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ER ! 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,-r2/delta_OOVV,t2) + 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 @@ -246,4 +252,28 @@ subroutine CCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ER write(*,'(1X,A15,1X,F10.6)') 'Ec(MP4-SDQ) = ',EcMP4 write(*,*) +!------------------------------------------------------------------------ +! EOM section +!------------------------------------------------------------------------ + +! EE-EOM-CCD (1h1p) + +! if(do_EE-EOM-CC_1h1p) call EE-EOM-CCD_1h1p() + +! 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) + end subroutine CCD diff --git a/src/CC/DEA_EOM_CCD_2p.f90 b/src/CC/DEA_EOM_CCD_2p.f90 new file mode 100644 index 0000000..4760530 --- /dev/null +++ b/src/CC/DEA_EOM_CCD_2p.f90 @@ -0,0 +1,116 @@ +subroutine DEA_EOM_CCD_2p(nC,nO,nV,nR,eV,OOVV,VVVV,t) + +! DEA-EOM-CCD calculation up to 2p + + implicit none + +! Input variables + + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: eV(nV) + double precision,intent(in) :: OOVV(nO,nO,nV,nV) + double precision,intent(in) :: VVVV(nV,nV,nV,nV) + double precision,intent(in) :: t(nO,nO,nV,nV) + +! Local variables + + integer :: a,b,c,d,ab,cd + integer :: i,j + integer :: nVV + double precision,external :: Kronecker_delta + double precision,allocatable :: F(:,:) + double precision,allocatable :: W(:,:,:,:) + double precision,allocatable :: H(:,:) + double precision,allocatable :: Om(:) + + +! Hello world + + write(*,*) + write(*,*)'********************' + write(*,*)'| DEA-EOM-CCD (2p) |' + write(*,*)'********************' + write(*,*) + +! Size of the EOM Hamiltonian + + nVV = nV*(nV-1)/2 + +! Memory allocation + + allocate(F(nV,nV),W(nV,nV,nV,nV),H(nVV,nVV),Om(nVV)) + +! Form one-body terms + + do a=1,nV-nR + do b=1,nV-nR + + F(a,b) = eV(a)*Kronecker_delta(a,b) + + do i=1,nO-nR + do j=1,nO-nR + do c=1,nV-nC + + F(a,b) = F(a,b) - 0.5d0*OOVV(i,j,b,c)*t(i,j,a,c) + + end do + end do + end do + + end do + end do + +! Form two-body terms + + do a=1,nV-nR + do b=1,nV-nR + do c=1,nV-nR + do d=1,nV-nR + + W(a,b,c,d) = VVVV(a,b,c,d) + + do i=1,nO-nC + do j=i+1,nO-nC + + W(a,b,c,d) = W(a,b,c,d) + OOVV(i,j,a,b)*t(i,j,a,b) + + end do + end do + + end do + end do + end do + end do + +! Form EOM Hamiltonian + + ab = 0 + do a=1,nV-nR + do b=a+1,nV-nR + ab = ab + 1 + + cd = 0 + do c=1,nV-nR + do d=c+1,nV-nR + cd = cd + 1 + + H(ab,cd) = F(a,c)*Kronecker_delta(b,d) + Kronecker_delta(a,c)*F(b,d) + W(a,b,c,d) + + end do + end do + + end do + end do + +! Diagonalize EOM Hamiltonian + + if(nVV > 0) call diagonalize_matrix(nVV,H,Om) + +! Dump results + + call print_excitation('DEA-EOM-CCD ',3,nVV,Om) + +end subroutine DEA_EOM_CCD_2p diff --git a/src/CC/DIP_EOM_CCD_2h.f90 b/src/CC/DIP_EOM_CCD_2h.f90 new file mode 100644 index 0000000..c48c5f7 --- /dev/null +++ b/src/CC/DIP_EOM_CCD_2h.f90 @@ -0,0 +1,116 @@ +subroutine DIP_EOM_CCD_2h(nC,nO,nV,nR,eO,OOVV,OOOO,t) + +! DIP-EOM-CCD calculation up to 2h + + implicit none + +! Input variables + + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: eO(nO) + double precision,intent(in) :: OOVV(nO,nO,nV,nV) + double precision,intent(in) :: OOOO(nO,nO,nO,nO) + double precision,intent(in) :: t(nO,nO,nV,nV) + +! Local variables + + integer :: i,j,k,l,ij,kl + integer :: a,b + integer :: nOO + double precision,external :: Kronecker_delta + double precision,allocatable :: F(:,:) + double precision,allocatable :: W(:,:,:,:) + double precision,allocatable :: H(:,:) + double precision,allocatable :: Om(:) + + +! Hello world + + write(*,*) + write(*,*)'********************' + write(*,*)'| DIP-EOM-CCD (2h) |' + write(*,*)'********************' + write(*,*) + +! Size of the EOM Hamiltonian + + nOO = nO*(nO-1)/2 + +! Memory allocation + + allocate(F(nO,nO),W(nO,nO,nO,nO),H(nOO,nOO),Om(nOO)) + +! Form one-body terms + + do i=1,nO-nR + do j=1,nO-nR + + F(i,j) = eO(i)*Kronecker_delta(i,j) + + do k=1,nO-nR + do a=1,nV-nC + do b=1,nV-nC + + F(i,j) = F(i,j) + 0.5d0*OOVV(i,k,a,b)*t(j,k,a,b) + + end do + end do + end do + + end do + end do + +! Form two-body terms + + do i=1,nO-nC + do j=1,nO-nC + do k=1,nO-nC + do l=1,nO-nC + + W(i,j,k,l) = OOOO(i,j,k,l) + + do a=1,nV-nR + do b=a+1,nV-nR + + W(i,j,k,l) = W(i,j,k,l) + OOVV(i,j,a,b)*t(i,j,a,b) + + end do + end do + + end do + end do + end do + end do + +! Form EOM Hamiltonian + + ij = 0 + do i=1,nO-nC + do j=i+1,nO-nC + ij = ij + 1 + + kl = 0 + do k=1,nO-nC + do l=k+1,nO-nC + kl = kl + 1 + + H(ij,kl) = - F(i,k)*Kronecker_delta(j,l) - Kronecker_delta(i,k)*F(j,l) + W(i,j,k,l) + + end do + end do + + end do + end do + +! Diagonalize EOM Hamiltonian + + if(nOO > 0) call diagonalize_matrix(nOO,H,Om) + +! Dump results + + call print_excitation('DIP-EOM-CCD ',3,nOO,Om) + +end subroutine DIP_EOM_CCD_2h diff --git a/src/CC/rCCD.f90 b/src/CC/rCCD.f90 index 14dcb0c..7b4318a 100644 --- a/src/CC/rCCD.f90 +++ b/src/CC/rCCD.f90 @@ -1,4 +1,4 @@ -subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF) +subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF) ! Ring CCD module @@ -6,6 +6,7 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF, ! Input variables + logical,intent(in) :: BSE integer,intent(in) :: maxSCF integer,intent(in) :: max_diis double precision,intent(in) :: thresh @@ -74,7 +75,15 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF, allocate(dbERI(nBas,nBas,nBas,nBas)) - call antisymmetrize_ERI(2,nBas,sERI,dbERI) + if(BSE) then + + call static_screening(nBas,nC,nO,nV,nR,seHF,sERI,dbERI) + + else + + call antisymmetrize_ERI(2,nBas,sERI,dbERI) + + end if deallocate(sERI) diff --git a/src/GT/static_Tmatrix_C_pp.f90 b/src/GT/static_Tmatrix_C_pp.f90 index 1904bda..0151b9b 100644 --- a/src/GT/static_Tmatrix_C_pp.f90 +++ b/src/GT/static_Tmatrix_C_pp.f90 @@ -55,15 +55,17 @@ subroutine static_Tmatrix_C_pp(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,nOOx,nVVx,lamb do ef=1,nVV eps = + Om1(ef) - chi = chi + rho1(a,b,ef)*rho1(c,d,ef)*eps/(eps**2 + eta**2) + chi = chi + rho1(a,b,ef)*rho1(c,d,ef)*eps/(eps**2 + eta**2) & + + rho1(a,b,ef)*rho1(d,c,ef)*eps/(eps**2 + eta**2) end do do mn=1,nOO eps = - Om2(mn) - chi = chi + rho2(a,b,mn)*rho2(c,d,mn)*eps/(eps**2 + eta**2) + chi = chi + rho2(a,b,mn)*rho2(c,d,mn)*eps/(eps**2 + eta**2) & + + rho2(a,b,mn)*rho2(d,c,mn)*eps/(eps**2 + eta**2) end do - TC(ab,cd) = lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + TC(ab,cd) = 0.5d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) end do end do diff --git a/src/GW/static_screening_WC_pp.f90 b/src/GW/static_screening_WC_pp.f90 index dd5205f..f9b4ccc 100644 --- a/src/GW/static_screening_WC_pp.f90 +++ b/src/GW/static_screening_WC_pp.f90 @@ -56,7 +56,7 @@ subroutine static_screening_WC_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E do m=1,nS eps = Omega(m)**2 + eta**2 chi = chi + rho(a,c,m)*rho(b,d,m)*Omega(m)/eps & - - rho(a,d,m)*rho(b,c,m)*Omega(m)/eps + + rho(a,d,m)*rho(b,c,m)*Omega(m)/eps enddo WC(ab,cd) = + 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) diff --git a/src/LR/print_excitation.f90 b/src/LR/print_excitation.f90 index 66f72e2..898f9f3 100644 --- a/src/LR/print_excitation.f90 +++ b/src/LR/print_excitation.f90 @@ -33,7 +33,7 @@ subroutine print_excitation(method,ispin,nS,Omega) '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' write(*,*)'-------------------------------------------------------------' - do ia=1,maxS + do ia=1,min(maxS,nS) write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & '|',ia,'|',Omega(ia),'|',Omega(ia)*HaToeV,'|' enddo diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index c1fafd6..14b88cf 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -582,15 +582,15 @@ program QuAcK write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CCD,' seconds' write(*,*) - call cpu_time(start_CCD) - call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, & - linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) - call CCD(.true.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eG0W0) - call cpu_time(end_CCD) +! call cpu_time(start_CCD) +! call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, & +! linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) +! call CCD(.true.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eG0W0) +! call cpu_time(end_CCD) - t_CCD = end_CCD - start_CCD - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD@BSE = ',t_CCD,' seconds' - write(*,*) +! t_CCD = end_CCD - start_CCD +! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD@BSE = ',t_CCD,' seconds' +! write(*,*) end if @@ -627,16 +627,16 @@ program QuAcK write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CCSD,' seconds' write(*,*) - call cpu_time(start_CCSD) - call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, & - linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) +! call cpu_time(start_CCSD) +! call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, & +! linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) - call CCSD(.true.,maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eG0W0) - call cpu_time(end_CCSD) +! call CCSD(.true.,maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eG0W0) +! call cpu_time(end_CCSD) - t_CCSD = end_CCSD - start_CCSD - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CCSD,' seconds' - write(*,*) +! t_CCSD = end_CCSD - start_CCSD +! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CCSD,' seconds' +! write(*,*) end if @@ -647,8 +647,7 @@ program QuAcK if(do_drCCD) then call cpu_time(start_CCD) - call drCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & - ERI_MO,ENuc,ERHF,eHF) + call drCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call cpu_time(end_CCD) t_CCD = end_CCD - start_CCD @@ -664,8 +663,16 @@ program QuAcK if(do_rCCD) then call cpu_time(start_CCD) - call rCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & - ERI_MO,ENuc,ERHF,eHF) + call rCCD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) + call cpu_time(end_CCD) + + t_CCD = end_CCD - start_CCD + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ring CCD = ',t_CCD,' seconds' + write(*,*) + call cpu_time(start_CCD) + call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, & + linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) + call rCCD(.true.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eG0W0) call cpu_time(end_CCD) t_CCD = end_CCD - start_CCD @@ -681,8 +688,7 @@ program QuAcK if(do_crCCD) then call cpu_time(start_CCD) - call crCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & - ERI_MO,ENuc,ERHF,eHF) + call crCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call cpu_time(end_CCD) diff --git a/src/RPA/ppRPA.f90 b/src/RPA/ppRPA.f90 index 4b968c2..e04a09c 100644 --- a/src/RPA/ppRPA.f90 +++ b/src/RPA/ppRPA.f90 @@ -74,10 +74,10 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI, & Omega1,X1,Y1,Omega2,X2,Y2,Ec_ppRPA(ispin)) - call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2) +! call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2) -! call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1) -! call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2) + call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1) + call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2) deallocate(Omega1,X1,Y1,Omega2,X2,Y2) @@ -102,10 +102,10 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI, & Omega1,X1,Y1,Omega2,X2,Y2,Ec_ppRPA(ispin)) - call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2) +! call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2) -! call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1) -! call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2) + call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1) + call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2) deallocate(Omega1,X1,Y1,Omega2,X2,Y2)