From cec0fe92056dbe5c5fd7386b3fbe631b10fe5bb2 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 10 Nov 2021 09:42:30 +0100 Subject: [PATCH] working on RPA third channel --- input/methods | 10 +- input/options | 2 +- mol/h2.xyz | 2 +- src/CC/CCD.f90 | 9 + src/CC/crCCD.f90 | 237 ++++++++++++++++++ src/CC/drCCD.f90 | 6 +- src/CC/form_ring_r.f90 | 5 +- src/CC/rCCD.f90 | 10 +- src/MBPT/Bethe_Salpeter_A_matrix_dynamic.f90 | 19 +- ..._Salpeter_Tmatrix_dynamic_perturbation.f90 | 6 +- .../Bethe_Salpeter_dynamic_perturbation.f90 | 6 +- src/MBPT/dynamic_Tmatrix_A.f90 | 26 +- src/MBPT/dynamic_Tmatrix_ZA.f90 | 73 ------ src/MBPT/static_Tmatrix_TA.f90 | 13 +- src/MBPT/static_Tmatrix_TB.f90 | 13 +- src/QuAcK/QuAcK.f90 | 23 +- src/RPA/ehRPA.f90 | 130 ++++++++++ 17 files changed, 468 insertions(+), 122 deletions(-) create mode 100644 src/CC/crCCD.f90 delete mode 100644 src/MBPT/dynamic_Tmatrix_ZA.f90 create mode 100644 src/RPA/ehRPA.f90 diff --git a/input/methods b/input/methods index efad4bb..c22322a 100644 --- a/input/methods +++ b/input/methods @@ -5,17 +5,17 @@ # CCD DCD CCSD CCSD(T) F F F F # drCCD rCCD lCCD pCCD - F F F F + T T F F # CIS* CIS(D) CID CISD FCI F F F F F # RPA* RPAx* ppRPA - F F F + T T F # G0F2* evGF2* qsGF2* G0F3 evGF3 - T F F F F + F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW - T F F F F + 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 88c4acd..a4093ba 100644 --- a/input/options +++ b/input/options @@ -5,7 +5,7 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000000001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - F T F T T + F T T T T # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.00367493 3 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 diff --git a/mol/h2.xyz b/mol/h2.xyz index 73b19fa..a4e936a 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0. 0. 0. -H 0. -0.740848 0.740848 +H 0. 0. 2.000000 diff --git a/src/CC/CCD.f90 b/src/CC/CCD.f90 index fad21b4..f5d06b1 100644 --- a/src/CC/CCD.f90 +++ b/src/CC/CCD.f90 @@ -220,6 +220,15 @@ subroutine CCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,e endif + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)' CCD energy ' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A30,1X,F15.10)')' E(CCD) = ',ECCD + write(*,'(1X,A30,1X,F15.10)')' Ec(CCD) = ',EcCCD + write(*,*)'----------------------------------------------------' + write(*,*) + ! Moller-Plesset energies write(*,*) diff --git a/src/CC/crCCD.f90 b/src/CC/crCCD.f90 new file mode 100644 index 0000000..36cd10d --- /dev/null +++ b/src/CC/crCCD.f90 @@ -0,0 +1,237 @@ +subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF) + +! Crossed-ring CCD module + + implicit none + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + + integer,intent(in) :: nBasin + integer,intent(in) :: nCin + integer,intent(in) :: nOin + integer,intent(in) :: nVin + integer,intent(in) :: nRin + double precision,intent(in) :: ENuc,ERHF + double precision,intent(in) :: eHF(nBasin) + double precision,intent(in) :: ERI(nBasin,nBasin,nBasin,nBasin) + +! Local variables + + integer :: nBas + integer :: nC + integer :: nO + integer :: nV + integer :: nR + integer :: nSCF + double precision :: Conv + double precision :: EcMP2 + 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 :: 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 :: r2r(:,:,:,:) + double precision,allocatable :: r2l(:,:,:,:) + + 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(*,*)'| crossed-ring CCD calculation |' + write(*,*)'**************************************' + write(*,*) + +! Spatial to spin orbitals + + nBas = 2*nBasin + nC = 2*nCin + nO = 2*nOin + nV = 2*nVin + nR = 2*nRin + + allocate(seHF(nBas),sERI(nBas,nBas,nBas,nBas)) + + call spatial_to_spin_MO_energy(nBasin,eHF,nBas,seHF) + call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI) + +! Antysymmetrize ERIs + + allocate(dbERI(nBas,nBas,nBas,nBas)) + + call antisymmetrize_ERI(2,nBas,sERI,dbERI) + + deallocate(sERI) + +! Form energy denominator + + allocate(eO(nO),eV(nV)) + allocate(delta_OOVV(nO,nO,nV,nV)) + + eO(:) = seHF(1:nO) + eV(:) = seHF(nO+1:nBas) + + call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV) + + deallocate(seHF) + +! Create integral batches + + allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVOV(nO,nV,nO,nV),VVVV(nV,nV,nV,nV),OVVO(nO,nV,nV,nO)) + + OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO ) + OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas) + OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas, 1:nO ,nO+1:nBas) + OVVO(:,:,:,:) = dbERI( 1:nO ,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(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(r2r(nO,nO,nV,nV),r2l(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(*,*)'| 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 + +! Form linear array + + call form_u(nC,nO,nV,nR,OOOO,VVVV,OVOV,t2,u) + +! Form interemediate arrays + + call form_X(nC,nO,nV,nR,OOVV,t2,X1,X2,X3,X4) + +! Form quadratic array + + call form_v(nC,nO,nV,nR,X1,X2,X3,X4,t2,v) + + call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2r) + + call form_ladder_r(nC,nO,nV,nR,OOOO,OOVV,VVVV,t2,r2l) + + r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + u(:,:,:,:) + v(:,:,:,:) - r2r(:,:,:,:) - r2l(:,:,:,:) + +! 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,EcCCD) + +! Dump results + + ECCD = ERHF + EcCCD + + ! 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,'|',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(*,*)' crossed-ring CCD energy ' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A30,1X,F15.10)')' E(crCCD) = ',ECCD + write(*,'(1X,A30,1X,F15.10)')' Ec(crCCD) = ',EcCCD + write(*,*)'----------------------------------------------------' + write(*,*) + +end subroutine crCCD diff --git a/src/CC/drCCD.f90 b/src/CC/drCCD.f90 index 063c4b0..0d21165 100644 --- a/src/CC/drCCD.f90 +++ b/src/CC/drCCD.f90 @@ -39,6 +39,7 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF double precision,allocatable :: OOVV(:,:,:,:) double precision,allocatable :: OVVO(:,:,:,:) + double precision,allocatable :: VOOV(:,:,:,:) double precision,allocatable :: r2(:,:,:,:) double precision,allocatable :: t2(:,:,:,:) @@ -84,10 +85,11 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF ! Create integral batches - allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO)) + allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VOOV(nV,nO,nO,nV)) OOVV(:,:,:,:) = sERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas) OVVO(:,:,:,:) = sERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO ) + VOOV(:,:,:,:) = sERI(nO+1:nBas, 1:nO , 1:nO ,nO+1:nBas) deallocate(sERI) @@ -133,7 +135,7 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF ! Compute residual - call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2) + call form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2) r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:) diff --git a/src/CC/form_ring_r.f90 b/src/CC/form_ring_r.f90 index 162cff5..2757650 100644 --- a/src/CC/form_ring_r.f90 +++ b/src/CC/form_ring_r.f90 @@ -1,4 +1,4 @@ -subroutine form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2) +subroutine form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2) ! Form residuals for ring CCD @@ -9,6 +9,7 @@ subroutine form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2) integer,intent(in) :: nC,nO,nV,nR double precision,intent(in) :: t2(nO,nO,nV,nV) double precision,intent(in) :: OVVO(nO,nV,nV,nO) + double precision,intent(in) :: VOOV(nV,nO,nO,nV) double precision,intent(in) :: OOVV(nO,nO,nV,nV) ! Local variables @@ -49,7 +50,7 @@ subroutine form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2) do k=nC+1,nO do c=1,nV-nR - r2(i,j,a,b) = r2(i,j,a,b) + OVVO(i,c,a,k)*t2(k,j,c,b) + OVVO(j,c,b,k)*t2(i,k,a,c) + r2(i,j,a,b) = r2(i,j,a,b) + VOOV(a,k,i,c)*t2(k,j,c,b) + OVVO(k,b,c,j)*t2(i,k,a,c) end do end do diff --git a/src/CC/rCCD.f90 b/src/CC/rCCD.f90 index bc8c10d..346f715 100644 --- a/src/CC/rCCD.f90 +++ b/src/CC/rCCD.f90 @@ -40,6 +40,7 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF, double precision,allocatable :: OOVV(:,:,:,:) double precision,allocatable :: OVVO(:,:,:,:) + double precision,allocatable :: VOOV(:,:,:,:) double precision,allocatable :: r2(:,:,:,:) double precision,allocatable :: t2(:,:,:,:) @@ -92,10 +93,11 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF, ! Create integral batches - allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO)) + allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VOOV(nV,nO,nO,nV)) - OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas) - OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas, 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 ) + VOOV(:,:,:,:) = dbERI(nO+1:nBas, 1:nO , 1:nO ,nO+1:nBas) deallocate(dbERI) @@ -141,7 +143,7 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF, ! Compute residual - call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2) + call form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2) r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:) diff --git a/src/MBPT/Bethe_Salpeter_A_matrix_dynamic.f90 b/src/MBPT/Bethe_Salpeter_A_matrix_dynamic.f90 index 634323f..9e0a2bf 100644 --- a/src/MBPT/Bethe_Salpeter_A_matrix_dynamic.f90 +++ b/src/MBPT/Bethe_Salpeter_A_matrix_dynamic.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,A_dyn) +subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,A_dyn,ZA_dyn) ! Compute the dynamic part of the Bethe-Salpeter equation matrices @@ -25,10 +25,12 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om ! Output variables double precision,intent(out) :: A_dyn(nS,nS) + double precision,intent(out) :: ZA_dyn(nS,nS) ! Initialization - A_dyn(:,:) = 0d0 + A_dyn(:,:) = 0d0 + ZA_dyn(:,:) = 0d0 ! Number of poles taken into account @@ -67,6 +69,19 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi + chi = 0d0 + do kc=1,maxS + + eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) + chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) + chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + enddo + + ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + 2d0*lambda*chi + enddo enddo enddo diff --git a/src/MBPT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90 b/src/MBPT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90 index 05b9fc2..75eae8b 100644 --- a/src/MBPT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90 +++ b/src/MBPT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90 @@ -88,11 +88,7 @@ subroutine Bethe_Salpeter_Tmatrix_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR ! Resonant part of the BSE correction for dynamical TDA - call dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGT,Omega1,Omega2,rho1,rho2,OmBSE(ia),Ap_dyn) - - ! Renormalization factor of the resonant parts for dynamical TDA - - call dynamic_Tmatrix_ZA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGT,Omega1,Omega2,rho1,rho2,OmBSE(ia),ZAp_dyn) + call dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGT,Omega1,Omega2,rho1,rho2,OmBSE(ia),Ap_dyn,Zap_dyn) ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) OmDyn(ia) = dot_product(X,matmul( Ap_dyn,X)) diff --git a/src/MBPT/Bethe_Salpeter_dynamic_perturbation.f90 b/src/MBPT/Bethe_Salpeter_dynamic_perturbation.f90 index 88dcc08..0af6209 100644 --- a/src/MBPT/Bethe_Salpeter_dynamic_perturbation.f90 +++ b/src/MBPT/Bethe_Salpeter_dynamic_perturbation.f90 @@ -82,11 +82,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,e ! Resonant part of the BSE correction for dynamical TDA - call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),Ap_dyn) - - ! Renormalization factor of the resonant parts for dynamical TDA - - call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),ZAp_dyn) + call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),Ap_dyn,ZAp_dyn) ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) OmDyn(ia) = dot_product(X,matmul( Ap_dyn,X)) diff --git a/src/MBPT/dynamic_Tmatrix_A.f90 b/src/MBPT/dynamic_Tmatrix_A.f90 index 7ab7ae5..31e75e2 100644 --- a/src/MBPT/dynamic_Tmatrix_A.f90 +++ b/src/MBPT/dynamic_Tmatrix_A.f90 @@ -1,4 +1,4 @@ -subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,Omega2,rho1,rho2,OmBSE,A_dyn) +subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,Omega2,rho1,rho2,OmBSE,A_dyn,ZA_dyn) ! Compute the dynamic part of the Bethe-Salpeter equation matrices for GT @@ -36,11 +36,13 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O ! Output variables - double precision,intent(out) :: A_dyn(nS,nS) + double precision,intent(out) :: A_dyn(nS,nS) + double precision,intent(out) :: ZA_dyn(nS,nS) ! Initialization - A_dyn(:,:) = 0d0 + A_dyn(:,:) = 0d0 + ZA_dyn(:,:) = 0d0 ! Build dynamic A matrix @@ -63,7 +65,7 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O chi = chi - rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2) end do - A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi + A_dyn(ia,jb) = A_dyn(ia,jb) - 1d0*lambda*chi chi = 0d0 @@ -77,7 +79,21 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) end do - A_dyn(ia,jb) = A_dyn(ia,jb) + 2d0*lambda*chi + A_dyn(ia,jb) = A_dyn(ia,jb) + 1d0*lambda*chi + + chi = 0d0 + + do cd=1,nVV + eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j)) + chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + + do kl=1,nOO + eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) + chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + + ZA_dyn(ia,jb) = ZA_dyn(ia,jb) - 1d0*lambda*chi end do end do diff --git a/src/MBPT/dynamic_Tmatrix_ZA.f90 b/src/MBPT/dynamic_Tmatrix_ZA.f90 deleted file mode 100644 index d3792a5..0000000 --- a/src/MBPT/dynamic_Tmatrix_ZA.f90 +++ /dev/null @@ -1,73 +0,0 @@ -subroutine dynamic_Tmatrix_ZA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,Omega2,rho1,rho2,OmBSE,ZA_dyn) - -! Compute the dynamic part of the Bethe-Salpeter equation matrices - - implicit none - include 'parameters.h' - -! Input variables - - double precision,intent(in) :: eta - integer,intent(in) :: nBas - integer,intent(in) :: nC - integer,intent(in) :: nO - integer,intent(in) :: nV - integer,intent(in) :: nR - integer,intent(in) :: nS - - integer,intent(in) :: nOO - integer,intent(in) :: nVV - - double precision,intent(in) :: lambda - double precision,intent(in) :: eGT(nBas) - double precision,intent(in) :: OmBSE - - double precision,intent(in) :: Omega1(nVV) - double precision,intent(in) :: Omega2(nOO) - double precision,intent(in) :: rho1(nBas,nBas,nVV) - double precision,intent(in) :: rho2(nBas,nBas,nOO) - -! Local variables - - double precision :: chi - double precision :: eps - integer :: i,j,a,b,ia,jb,cd,kl - -! Output variables - - double precision,intent(out) :: ZA_dyn(nS,nS) - -! Initialization - - ZA_dyn(:,:) = 0d0 - -! Build dynamic A matrix - - ia = 0 - do i=nC+1,nO - do a=nO+1,nBas-nR - ia = ia + 1 - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - - chi = 0d0 - do cd=1,nVV - eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j)) - chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - end do - - do kl=1,nOO - eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) - chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - end do - - ZA_dyn(ia,jb) = ZA_dyn(ia,jb) - 2d0*lambda*chi - - end do - end do - end do - end do - -end subroutine dynamic_Tmatrix_ZA diff --git a/src/MBPT/static_Tmatrix_TA.f90 b/src/MBPT/static_Tmatrix_TA.f90 index d6e48aa..ac4ed6e 100644 --- a/src/MBPT/static_Tmatrix_TA.f90 +++ b/src/MBPT/static_Tmatrix_TA.f90 @@ -26,7 +26,6 @@ subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r ! Local variables double precision :: chi - double precision :: eps integer :: i,j,a,b,ia,jb,kl,cd ! Output variables @@ -45,18 +44,16 @@ subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r chi = 0d0 do cd=1,nVV - eps = Omega1(cd)**2 + eta**2 -! chi = chi + lambda*rho1(i,j,cd)*rho1(a,b,cd)*Omega1(cd)/eps - chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*Omega1(cd)/eps +! chi = chi + lambda*rho1(i,j,cd)*rho1(a,b,cd)*Omega1(cd)/(Omega1(cd)**2 + eta**2) + chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*Omega1(cd)/(Omega1(cd)**2 + eta**2) enddo do kl=1,nOO - eps = Omega2(kl)**2 + eta**2 -! chi = chi - lambda*rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/eps - chi = chi - rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/eps +! chi = chi - lambda*rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2) + chi = chi - rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2) enddo - TA(ia,jb) = TA(ia,jb) + 2d0*lambda*chi + TA(ia,jb) = TA(ia,jb) + 1d0*lambda*chi enddo enddo diff --git a/src/MBPT/static_Tmatrix_TB.f90 b/src/MBPT/static_Tmatrix_TB.f90 index 9c00e8a..3613980 100644 --- a/src/MBPT/static_Tmatrix_TB.f90 +++ b/src/MBPT/static_Tmatrix_TB.f90 @@ -26,7 +26,6 @@ subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r ! Local variables double precision :: chi - double precision :: eps integer :: i,j,a,b,ia,jb,kl,cd ! Output variables @@ -45,18 +44,16 @@ subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r chi = 0d0 do cd=1,nVV - eps = Omega1(cd)**2 + eta**2 -! chi = chi + lambda*rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/eps - chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/eps +! chi = chi + lambda*rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/Omega1(cd)**2 + eta**2 + chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/Omega1(cd)**2 + eta**2 enddo do kl=1,nOO - eps = Omega2(kl)**2 + eta**2 -! chi = chi - lambda*rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/eps - chi = chi - rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/eps +! chi = chi - lambda*rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2 + chi = chi - rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2 enddo - TB(ia,jb) = TB(ia,jb) + 2d0*lambda*chi + TB(ia,jb) = TB(ia,jb) + 1d0*lambda*chi enddo enddo diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 8fd4f99..984550f 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -10,7 +10,7 @@ program QuAcK logical :: doKS logical :: doMP2,doMP3,doMP2F12 logical :: doCCD,doDCD,doCCSD,doCCSDT - logical :: do_drCCD,do_rCCD,do_lCCD,do_pCCD + logical :: do_drCCD,do_rCCD,do_lCCD,do_crCCD,do_pCCD logical :: doCIS,doCIS_D,doCID,doCISD,doFCI logical :: doRPA,doRPAx,doppRPA logical :: doADC @@ -662,6 +662,27 @@ program QuAcK end if +!------------------------------------------------------------------------ +! Perform crossed-ring CCD calculation +!------------------------------------------------------------------------ + + do_crCCD = .false. + + if(do_crCCD) then + + call cpu_time(start_CCD) + call ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,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) + + t_CCD = end_CCD - start_CCD + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CCD,' seconds' + write(*,*) + + end if + !------------------------------------------------------------------------ ! Perform pair CCD calculation !------------------------------------------------------------------------ diff --git a/src/RPA/ehRPA.f90 b/src/RPA/ehRPA.f90 new file mode 100644 index 0000000..103f5c3 --- /dev/null +++ b/src/RPA/ehRPA.f90 @@ -0,0 +1,130 @@ +subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & + ERI,dipole_int,eHF) + +! Perform random phase approximation calculation with exchange (aka TDHF) + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + logical,intent(in) :: TDA + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: singlet + double precision,intent(in) :: eta + logical,intent(in) :: triplet + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + 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) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + integer :: ispin + double precision,allocatable :: Omega(:,:) + double precision,allocatable :: XpY(:,:,:) + double precision,allocatable :: XmY(:,:,:) + + double precision :: rho + double precision :: EcRPAx(nspin) + double precision :: EcAC(nspin) + +! Hello world + + write(*,*) + write(*,*)'***********************************************************' + write(*,*)'| Random phase approximation calculation: eh channel |' + write(*,*)'***********************************************************' + write(*,*) + +! TDA + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if + +! Initialization + + EcRPAx(:) = 0d0 + EcAC(:) = 0d0 + +! Memory allocation + + allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) + +! Singlet manifold + + if(singlet) then + + ispin = 1 + + call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Omega(:,ispin),rho, & + EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_excitation('ehRPA@HF ',ispin,nS,Omega(:,ispin)) + call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + + endif + +! Triplet manifold + + if(triplet) then + + ispin = 2 + + call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Omega(:,ispin),rho, & + EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_excitation('ehRPA@HF ',ispin,nS,Omega(:,ispin)) + call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + + endif + + if(exchange_kernel) then + + EcRPAx(1) = 0.5d0*EcRPAx(1) + EcRPAx(2) = 1.5d0*EcRPAx(2) + + end if + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@ehRPA correlation energy (singlet) =',EcRPAx(1) + write(*,'(2X,A50,F20.10)') 'Tr@ehRPA correlation energy (triplet) =',EcRPAx(2) + write(*,'(2X,A50,F20.10)') 'Tr@ehRPA correlation energy =',EcRPAx(1) + EcRPAx(2) + write(*,'(2X,A50,F20.10)') 'Tr@ehRPA total energy =',ENuc + ERHF + EcRPAx(1) + EcRPAx(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! Compute the correlation energy via the adiabatic connection + +! if(doACFDT) then + +! write(*,*) '-------------------------------------------------------' +! write(*,*) 'Adiabatic connection version of ehRPA correlation energy' +! write(*,*) '-------------------------------------------------------' +! write(*,*) + +! call ACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, & +! nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (singlet) =',EcAC(1) +! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (triplet) =',EcAC(2) +! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy =',EcAC(1) + EcAC(2) +! write(*,'(2X,A50,F20.10)') 'AC@RPAx total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + +end subroutine ehRPA