From 96621bd038666ed6e85c996e4ce4d205574865ba Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 28 Sep 2022 16:09:09 +0200 Subject: [PATCH] adding routines for MBPT vs CC --- input/methods | 6 +- input/options | 8 +- src/CC/CCD.f90 | 2 +- src/CC/EE_EOM_CCD_1h1p.f90 | 25 ++++-- src/CC/rCCD.f90 | 61 ++++++++++----- src/GW/soBSE.f90 | 122 ++++++++++++++++++++++++++++++ src/GW/static_screening_WA_so.f90 | 53 +++++++++++++ src/GW/static_screening_WB_so.f90 | 58 ++++++++++++++ src/QuAcK/QuAcK.f90 | 19 ++--- src/RPA/soRPAx.f90 | 113 +++++++++++++++++++++++++++ src/RPA/sort_ppRPA.f90 | 4 +- 11 files changed, 428 insertions(+), 43 deletions(-) create mode 100644 src/GW/soBSE.f90 create mode 100644 src/GW/static_screening_WA_so.f90 create mode 100644 src/GW/static_screening_WB_so.f90 create mode 100644 src/RPA/soRPAx.f90 diff --git a/input/methods b/input/methods index 55ea358..e248d6e 100644 --- a/input/methods +++ b/input/methods @@ -3,13 +3,13 @@ # MP2* MP3 MP2-F12 F F F # CCD pCCD DCD CCSD CCSD(T) - T F F F F + F F F F F # drCCD rCCD crCCD lCCD - F F F F + F T F F # CIS* CIS(D) CID CISD FCI F F F F F # RPA* RPAx* crRPA ppRPA - F T F F + F F F F # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW diff --git a/input/options b/input/options index 5990b1c..d83b3d9 100644 --- a/input/options +++ b/input/options @@ -3,9 +3,9 @@ # MP: # CC: maxSCF thresh DIIS n_diis - 64 0.00001 T 5 + 64 0.0000001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - T T T T T + F 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 @@ -13,8 +13,8 @@ # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg 10 0.00001 T 5 T 0.0 T F # ACFDT: AC Kx XBS - F F T + F T T # BSE: BSE dBSE dTDA evDyn ppBSE - T F T F F + F 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 81d4ab1..e418dac 100644 --- a/src/CC/CCD.f90 +++ b/src/CC/CCD.f90 @@ -62,7 +62,7 @@ subroutine CCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ER double precision,allocatable :: error_diis(:,:) double precision,allocatable :: t_diis(:,:) - logical :: do_EE_EOM_CC_1h1p = .true. + 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. diff --git a/src/CC/EE_EOM_CCD_1h1p.f90 b/src/CC/EE_EOM_CCD_1h1p.f90 index 1d8fbdc..9ef2200 100644 --- a/src/CC/EE_EOM_CCD_1h1p.f90 +++ b/src/CC/EE_EOM_CCD_1h1p.f90 @@ -28,7 +28,9 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t) double precision,allocatable :: Wovvo(:,:,:,:) double precision,allocatable :: H(:,:) double precision,allocatable :: Om(:) + double precision,allocatable :: Z(:,:) + integer,allocatable :: order(:) ! Hello world @@ -44,7 +46,9 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t) ! Memory allocation - allocate(Foo(nO,nO),Fvv(nV,nV),Wovvo(nO,nV,nV,nO),H(nS,nS),Om(nS)) + allocate(Foo(nO,nO),Fvv(nV,nV),Wovvo(nO,nV,nV,nO),H(nS,nS),Om(nS),Z(nS,nS)) + allocate(order(nS)) + ! Form one-body terms @@ -57,7 +61,7 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t) do j=1,nO-nC do c=1,nV-nR - Fvv(a,b) = Fvv(a,b) - 0.5d0*OOVV(i,j,b,c)*t(i,j,a,c) +! Fvv(a,b) = Fvv(a,b) - 0.5d0*OOVV(i,j,b,c)*t(i,j,a,c) end do end do @@ -75,7 +79,7 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t) do a=1,nV-nR do b=1,nV-nR - Foo(i,j) = Foo(i,j) + 0.5d0*OOVV(i,k,a,b)*t(j,k,a,b) +! Foo(i,j) = Foo(i,j) + 0.5d0*OOVV(i,k,a,b)*t(j,k,a,b) end do end do @@ -128,10 +132,19 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t) ! Diagonalize EOM Hamiltonian - if(nS > 0) call diagonalize_matrix(nS,H,Om) + if(nS > 0) then -! Dump results + call diagonalize_general_matrix(nS,H,Om,Z) - call print_excitation('EE-EOM-CCD ',3,nS,Om) + do ia=1,nS + order(ia) = ia + end do + + call quick_sort(Om,order,nS) + call set_order(Z,order,nS,nS) + + call print_excitation('EE-EOM-CCD ',3,nS,Om) + + end if end subroutine EE_EOM_CCD_1h1p diff --git a/src/CC/rCCD.f90 b/src/CC/rCCD.f90 index 7b4318a..f30ebd8 100644 --- a/src/CC/rCCD.f90 +++ b/src/CC/rCCD.f90 @@ -1,4 +1,4 @@ -subroutine rCCD(BSE,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,eGW) ! Ring CCD module @@ -16,8 +16,10 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E integer,intent(in) :: nOin integer,intent(in) :: nVin integer,intent(in) :: nRin - double precision,intent(in) :: ENuc,ERHF + 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 @@ -32,6 +34,7 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E double precision :: EcMP2 double precision :: ECCD,EcCCD double precision,allocatable :: seHF(:) + double precision,allocatable :: seGW(:) double precision,allocatable :: sERI(:,:,:,:) double precision,allocatable :: dbERI(:,:,:,:) @@ -42,14 +45,16 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E double precision,allocatable :: OOVV(:,:,:,:) double precision,allocatable :: OVVO(:,:,:,:) - double precision,allocatable :: r2(:,:,:,:) - double precision,allocatable :: t2(:,:,:,:) + 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(*,*) @@ -66,11 +71,16 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E nV = 2*nVin nR = 2*nRin - allocate(seHF(nBas),sERI(nBas,nBas,nBas,nBas)) + allocate(seHF(nBas),seGW(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) +! call soRPAx(.false.,nBas,nC,nO,nV,nR,nO*nV,ENuc,ERHF,sERI,seHF) + + call soBSE(.false.,.false.,0.0,nBas,nC,nO,nV,nR,nO*nV,sERI,seHF,seGW) + ! Antysymmetrize ERIs allocate(dbERI(nBas,nBas,nBas,nBas)) @@ -92,12 +102,12 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E allocate(eO(nO),eV(nV)) allocate(delta_OOVV(nO,nO,nV,nV)) - eO(:) = seHF(1:nO) - eV(:) = seHF(nO+1:nBas) + eO(:) = seGW(1:nO) + eV(:) = seGW(nO+1:nBas) call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV) - deallocate(seHF) +! deallocate(seHF,seGW) ! Create integral batches @@ -110,11 +120,11 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E ! MP2 guess amplitudes - allocate(t2(nO,nO,nV,nV)) + allocate(t(nO,nO,nV,nV)) - 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) ! Memory allocation for DIIS @@ -122,7 +132,7 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E ! Initialization - allocate(r2(nO,nO,nV,nV)) + allocate(r(nO,nO,nV,nV)) Conv = 1d0 nSCF = 0 @@ -150,21 +160,21 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E ! Compute residual - call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2) + call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t,r) - r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:) + r(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t(:,:,:,:) + r(:,:,:,:) ! Check convergence - Conv = maxval(abs(r2(nC+1:nO,nC+1:nO,1:nV-nR,1:nV-nR))) + Conv = maxval(abs(r(nC+1:nO,nC+1:nO,1:nV-nR,1:nV-nR))) ! 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) ! Dump results @@ -173,7 +183,7 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E ! 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) + 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 @@ -211,4 +221,19 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E 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) + end subroutine rCCD diff --git a/src/GW/soBSE.f90 b/src/GW/soBSE.f90 new file mode 100644 index 0000000..28266df --- /dev/null +++ b/src/GW/soBSE.f90 @@ -0,0 +1,122 @@ +subroutine soBSE(TDA_W,TDA,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,eGW) + +! Compute the Bethe-Salpeter excitation energies + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: TDA_W + logical,intent(in) :: TDA + + 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 + double precision,intent(in) :: eW(nBas) + double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: ispin + integer :: isp_W + + double precision :: EcRPA + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) + + double precision,allocatable :: OmBSE(:) + double precision,allocatable :: XpY_BSE(:,:) + double precision,allocatable :: XmY_BSE(:,:) + + double precision,allocatable :: WA_sta(:,:) + double precision,allocatable :: WB_sta(:,:) + + double precision,allocatable :: X(:,:) + double precision,allocatable :: Y(:,:) + double precision,allocatable :: Xinv(:,:) + double precision,allocatable :: t(:,:,:,:) + + integer ::i,a,j,b,k,c,ia,jb,kc + + double precision :: EcBSE + +! Memory allocation + + allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & + WA_sta(nS,nS),WB_sta(nS,nS),OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS)) + +!--------------------------------- +! Compute (singlet) RPA screening +!--------------------------------- + + isp_W = 3 + EcRPA = 0d0 + + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + + call static_screening_WA_so(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA_sta) + call static_screening_WB_so(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB_sta) + + + ispin = 3 + EcBSE = 0d0 + + ! Compute BSE excitation energies + + call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,WA_sta,WB_sta, & + EcBSE,OmBSE,XpY_BSE,XmY_BSE) + call print_excitation('soBSE@GW ',ispin,nS,OmBSE) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy =',0.5d0*EcBSE,' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! allocate(X(nS,nS),Y(nS,nS),Xinv(nS,nS),t(nO,nO,nV,nV)) + +! X(:,:) = transpose(0.5d0*(XpY_BSE(:,:) + XmY_BSE(:,:))) +! Y(:,:) = transpose(0.5d0*(XpY_BSE(:,:) - XmY_BSE(:,:))) + +! call matout(nS,nS,matmul(X,transpose(X))-matmul(Y,transpose(Y))) + +! call inverse_matrix(nS,X,Xinv) + +! t = 0d0 +! ia = 0 +! do i=1,nO +! do a=1,nV +! ia = ia + 1 + +! jb = 0 +! do j=1,nO +! do b=1,nV +! jb = jb + 1 + +! kc = 0 +! do k=1,nO +! do c=1,nV +! kc = kc + 1 + +! t(i,j,a,b) = t(i,j,a,b) + Y(ia,kc)*Xinv(kc,jb) + +! end do +! end do +! end do +! end do +! end do +! end do + +! call matout(nO*nO,nV*nV,t) + +end subroutine soBSE diff --git a/src/GW/static_screening_WA_so.f90 b/src/GW/static_screening_WA_so.f90 new file mode 100644 index 0000000..951a013 --- /dev/null +++ b/src/GW/static_screening_WA_so.f90 @@ -0,0 +1,53 @@ +subroutine static_screening_WA_so(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WA) + +! Compute the OOVV block of the static screening W for the resonant block + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + double precision :: chi + double precision :: eps + integer :: i,j,a,b,ia,jb,kc + +! Output variables + + double precision,intent(out) :: WA(nS,nS) + +! Initialize + + WA(:,:) = 0d0 + + 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 kc=1,nS + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,j,kc)*rho(a,b,kc)*Omega(kc)/eps + enddo + + WA(ia,jb) = WA(ia,jb) + lambda*ERI(i,b,j,a) - 2d0*lambda*chi + + enddo + enddo + enddo + enddo + +end subroutine static_screening_WA_so diff --git a/src/GW/static_screening_WB_so.f90 b/src/GW/static_screening_WB_so.f90 new file mode 100644 index 0000000..849f489 --- /dev/null +++ b/src/GW/static_screening_WB_so.f90 @@ -0,0 +1,58 @@ +subroutine static_screening_WB_so(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WB) + +! Compute the static screening W for the coupling block + + implicit none + include 'parameters.h' + +! Input variables + + 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) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + double precision :: chi + double precision :: eps + integer :: i,j,a,b,ia,jb,kc + +! Output variables + + double precision,intent(out) :: WB(nS,nS) + +! Initialize + + WB(:,:) = 0d0 + + 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 kc=1,nS + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,b,kc)*rho(a,j,kc)*Omega(kc)/eps + enddo + + WB(ia,jb) = WB(ia,jb) + lambda*ERI(i,j,b,a) - 2d0*lambda*chi + + enddo + enddo + enddo + enddo + +end subroutine static_screening_WB_so diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 14b88cf..1cc2354 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -662,21 +662,22 @@ program QuAcK if(do_rCCD) then - call cpu_time(start_CCD) - 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) +! call cpu_time(start_CCD) +! call rCCD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,eHF) +! call cpu_time(end_CCD) + +! t_CCD = end_CCD - start_CCD +! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for rCCD = ',t_CCD,' seconds' +! write(*,*) - 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 rCCD(.true.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,eG0W0) 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(*,'(A65,1X,F9.3,A8)') 'Total CPU time for rCCD@BSE = ',t_CCD,' seconds' write(*,*) end if @@ -829,7 +830,7 @@ program QuAcK else - call RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) + call RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) end if call cpu_time(end_RPA) diff --git a/src/RPA/soRPAx.f90 b/src/RPA/soRPAx.f90 new file mode 100644 index 0000000..8b37398 --- /dev/null +++ b/src/RPA/soRPAx.f90 @@ -0,0 +1,113 @@ +subroutine soRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) + +! Perform random phase approximation calculation with exchange (aka TDHF) in the +! spinorbital basis + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + logical,intent(in) :: TDA + 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) + +! Local variables + + integer :: ispin + double precision,allocatable :: Omega(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) + double precision,allocatable :: X(:,:) + double precision,allocatable :: Y(:,:) + double precision,allocatable :: Xinv(:,:) + double precision,allocatable :: t(:,:,:,:) + + double precision :: EcRPAx + + integer ::i,a,j,b,k,c,ia,jb,kc + +! Hello world + + write(*,*) + write(*,*)'***********************************************************' + write(*,*)'| Random phase approximation calculation with exchange |' + write(*,*)'***********************************************************' + write(*,*) + +! TDA + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) ' => RPAx + TDA = CIS ' + write(*,*) + end if + +! Initialization + + EcRPAx = 0d0 + +! Memory allocation + + allocate(Omega(nS),XpY(nS,nS),XmY(nS,nS)) + + ispin = 3 + + call linear_response(ispin,.false.,TDA,0d0,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,EcRPAx,Omega,XpY,XmY) + call print_excitation('soRPAx@HF ',ispin,nS,Omega) + + EcRPAx = 0.5d0*EcRPAx + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@RPAx correlation energy =',EcRPAx + write(*,'(2X,A50,F20.10)') 'Tr@RPAx total energy =',ENuc + ERHF + EcRPAx + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! allocate(X(nS,nS),Y(nS,nS),Xinv(nS,nS),t(nO,nO,nV,nV)) + +! X(:,:) = transpose(0.5d0*(XpY(:,:) + XmY(:,:))) +! Y(:,:) = transpose(0.5d0*(XpY(:,:) - XmY(:,:))) + +! call matout(nS,nS,matmul(transpose(X),X)-matmul(transpose(Y),Y)) + +! call inverse_matrix(nS,X,Xinv) + +! t = 0d0 +! ia = 0 +! do i=1,nO +! do a=1,nV +! ia = ia + 1 + +! jb = 0 +! do j=1,nO +! do b=1,nV +! jb = jb + 1 + +! kc = 0 +! do k=1,nO +! do c=1,nV +! kc = kc + 1 + +! t(i,j,a,b) = t(i,j,a,b) + Y(ia,kc)*Xinv(kc,jb) + +! end do +! end do +! end do +! end do +! end do +! end do + +! call matout(nO*nO,nV*nV,t) + +end subroutine soRPAx diff --git a/src/RPA/sort_ppRPA.f90 b/src/RPA/sort_ppRPA.f90 index ef5ec93..3d2cace 100644 --- a/src/RPA/sort_ppRPA.f90 +++ b/src/RPA/sort_ppRPA.f90 @@ -15,8 +15,8 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) ! Local variables integer :: pq,ab,ij - integer :: deg1,ab_start,ab_end - integer :: deg2,ij_start,ij_end +! integer :: deg1,ab_start,ab_end +! integer :: deg2,ij_start,ij_end double precision,allocatable :: M(:,:) double precision,allocatable :: Z1(:,:) double precision,allocatable :: Z2(:,:)