diff --git a/basis/aug-cc-pvtz b/basis/aug-cc-pvtz index 5a5fd36..c8b16ad 100644 --- a/basis/aug-cc-pvtz +++ b/basis/aug-cc-pvtz @@ -121,7 +121,7 @@ S 1 S 1 1 4.409000E-02 1.000000E+00 S 1 -1 1.470000E-02 1.000000E+00 +1 1.503000E-02 1.000000E+00 P 5 1 7.436000E+00 1.073600E-02 2 1.577000E+00 6.285400E-02 @@ -133,17 +133,17 @@ P 1 P 1 1 4.994000E-02 1.000000E+00 P 1 -1 9.300000E-03 1.000000E+00 +1 7.060000E-03 1.000000E+00 D 1 -1 3.493000E-01 1.000000E+00 +1 3.480000E-01 1.000000E+00 D 1 -1 1.724000E-01 1.000000E+00 +1 1.803000E-01 1.000000E+00 D 1 -1 5.880000E-02 1.000000E+00 +1 6.540000E-02 1.000000E+00 F 1 -1 3.423000E-01 1.0000000 +1 3.250000E-01 1.000000E+00 F 1 -1 1.188000E-01 1.000000E+00 +1 1.533000E-01 1.000000E+00 BORON S 8 diff --git a/input/methods b/input/methods index d16be27..cad0336 100644 --- a/input/methods +++ b/input/methods @@ -9,11 +9,11 @@ # CIS* CIS(D) CID CISD F F F F # RPA* RPAx* ppRPA - F F F + T F F # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0* evGW* qsGW - T F F + F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index 02d5453..5bd9d2e 100644 --- a/input/options +++ b/input/options @@ -1,17 +1,17 @@ # HF: maxSCF thresh DIIS n_diis guess_type ortho_type - 128 0.00001 T 3 1 1 + 128 0.0000001 T 5 1 1 # MP: # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - T T T T T + T T T T F # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.0 3 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 256 0.00001 T 5 T 0.0 F F F F F # ACFDT: AC Kx XBS - F F T + T F T # BSE: BSE dBSE dTDA evDyn T T T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift diff --git a/src/MBPT/UG0W0.f90 b/src/MBPT/UG0W0.f90 index ad33756..a4badad 100644 --- a/src/MBPT/UG0W0.f90 +++ b/src/MBPT/UG0W0.f90 @@ -183,58 +183,58 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,S, & ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eGW,EcBSE) -! if(exchange_kernel) then -! -! EcRPA(1) = 0.5d0*EcRPA(1) -! EcRPA(2) = 1.5d0*EcRPA(1) -! -! end if + if(exchange_kernel) then + + EcBSE(1) = 0.5d0*EcBSE(1) + EcBSE(2) = 1.5d0*EcBSE(1) + + end if -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1) -! write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy (triplet) =',EcBSE(2) -! write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2) -! write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 total energy =',ENuc + EUHF + EcBSE(1) + EcBSE(2) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy (spin-conserved) =',EcBSE(1) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy (spin-flip) =',EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 total energy =',ENuc + EUHF + EcBSE(1) + EcBSE(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) ! Compute the BSE correlation energy via the adiabatic connection -! if(doACFDT) then + if(doACFDT) then -! write(*,*) '------------------------------------------------------' -! write(*,*) 'Adiabatic connection version of BSE correlation energy' -! write(*,*) '------------------------------------------------------' -! write(*,*) + write(*,*) '------------------------------------------------------' + write(*,*) 'Adiabatic connection version of BSE correlation energy' + write(*,*) '------------------------------------------------------' + write(*,*) -! if(doXBS) then + if(doXBS) then -! write(*,*) '*** scaled screening version (XBS) ***' -! write(*,*) + write(*,*) '*** scaled screening version (XBS) ***' + write(*,*) -! end if + end if -! call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, & -! nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcAC) + call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip,eta, & + nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGW,EcAC) -! if(exchange_kernel) then -! -! EcAC(1) = 0.5d0*EcAC(1) -! EcAC(2) = 1.5d0*EcAC(1) -! -! end if + if(exchange_kernel) then + + EcAC(1) = 0.5d0*EcAC(1) + EcAC(2) = 1.5d0*EcAC(1) + + end if -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1) -! write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2) -! write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy =',EcAC(1) + EcAC(2) -! write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 total energy =',ENuc + EUHF + EcAC(1) + EcAC(2) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (spin-conserved) =',EcAC(1) + write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (spin-flip) =',EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy =',EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 total energy =',ENuc + EUHF + EcAC(1) + EcAC(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) -! end if + end if end if diff --git a/src/RPA/URPAx.f90 b/src/RPA/URPAx.f90 index 42cc633..b3313cc 100644 --- a/src/RPA/URPAx.f90 +++ b/src/RPA/URPAx.f90 @@ -119,12 +119,12 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n endif -! if(exchange_kernel) then + if(exchange_kernel) then -! EcRPAx(1) = 0.5d0*EcRPAx(1) -! EcRPAx(2) = 1.5d0*EcRPAx(2) + EcRPAx(1) = 0.5d0*EcRPAx(1) + EcRPAx(2) = 1.5d0*EcRPAx(2) -! end if + end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -137,32 +137,32 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n ! Compute the correlation energy via the adiabatic connection -! if(doACFDT) then + if(doACFDT) then -! write(*,*) '-------------------------------------------------------' -! write(*,*) 'Adiabatic connection version of RPAx correlation energy' -! write(*,*) '-------------------------------------------------------' -! write(*,*) + write(*,*) '--------------------------------------------------------' + write(*,*) 'Adiabatic connection version of URPAx correlation energy' + write(*,*) '--------------------------------------------------------' + write(*,*) -! call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,singlet,triplet,eta, & -! nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC) + call unrestricted_ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,spin_conserved,spin_flip,eta, & + nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,e,e,EcAC) -! if(exchange_kernel) then + if(exchange_kernel) then -! EcAC(1) = 0.5d0*EcAC(1) -! EcAC(2) = 1.5d0*EcAC(2) + EcAC(1) = 0.5d0*EcAC(1) + EcAC(2) = 1.5d0*EcAC(2) -! end if + end if -! 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 + EUHF + EcAC(1) + EcAC(2) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (spin-conserved) =',EcAC(1) + write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (spin-flip) =',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 + EUHF + EcAC(1) + EcAC(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) -! end if + end if end subroutine URPAx diff --git a/src/RPA/UdRPA.f90 b/src/RPA/UdRPA.f90 index 22f29ec..b95e685 100644 --- a/src/RPA/UdRPA.f90 +++ b/src/RPA/UdRPA.f90 @@ -115,12 +115,12 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n endif -! if(exchange_kernel) then + if(exchange_kernel) then -! EcRPA(1) = 0.5d0*EcRPA(1) -! EcRPA(2) = 1.5d0*EcRPA(2) + EcRPA(1) = 0.5d0*EcRPA(1) + EcRPA(2) = 1.5d0*EcRPA(2) -! end if + end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -133,32 +133,32 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n ! Compute the correlation energy via the adiabatic connection -! if(doACFDT) then + if(doACFDT) then -! write(*,*) '-------------------------------------------------------' -! write(*,*) 'Adiabatic connection version of RPA correlation energy' -! write(*,*) '-------------------------------------------------------' -! write(*,*) + write(*,*) '-------------------------------------------------------' + write(*,*) 'Adiabatic connection version of RPA correlation energy' + write(*,*) '-------------------------------------------------------' + write(*,*) -! call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,singlet,triplet,eta, & -! nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC) + call unrestricted_ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,spin_conserved,spin_flip,eta, & + nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,e,e,EcAC) -! if(exchange_kernel) then + if(exchange_kernel) then -! EcAC(1) = 0.5d0*EcAC(1) -! EcAC(2) = 1.5d0*EcAC(2) + EcAC(1) = 0.5d0*EcAC(1) + EcAC(2) = 1.5d0*EcAC(2) -! end if + end if -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy (singlet) =',EcAC(1) -! write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy (triplet) =',EcAC(2) -! write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy =',EcAC(1) + EcAC(2) -! write(*,'(2X,A50,F20.10)') 'AC@RPA total energy =',ENuc + EUHF + EcAC(1) + EcAC(2) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy (spin-conserved) =',EcAC(1) + write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy (spin-flip) =',EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy =',EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@RPA total energy =',ENuc + EUHF + EcAC(1) + EcAC(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) -! end if + end if end subroutine UdRPA diff --git a/src/RPA/unrestricted_ACFDT.f90 b/src/RPA/unrestricted_ACFDT.f90 new file mode 100644 index 0000000..8735cb5 --- /dev/null +++ b/src/RPA/unrestricted_ACFDT.f90 @@ -0,0 +1,206 @@ +subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin_flip,eta, & + nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eW,e,EcAC) + +! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + logical,intent(in) :: doXBS + logical,intent(in) :: dRPA + logical,intent(in) :: TDA_W + logical,intent(in) :: TDA + logical,intent(in) :: BSE + logical,intent(in) :: exchange_kernel + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nS(nspin) + double precision,intent(in) :: eW(nBas,nspin) + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: ispin + integer :: isp_W + integer :: iAC + double precision :: lambda + double precision,allocatable :: Ec(:,:) + + double precision :: EcRPA + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) + + integer :: nS_aa,nS_bb,nS_sc + double precision,allocatable :: Omega_sc(:) + double precision,allocatable :: XpY_sc(:,:) + double precision,allocatable :: XmY_sc(:,:) + + integer :: nS_ab,nS_ba,nS_sf + double precision,allocatable :: Omega_sf(:) + double precision,allocatable :: XpY_sf(:,:) + double precision,allocatable :: XmY_sf(:,:) + +! Output variables + + double precision,intent(out) :: EcAC(nspin) + +! Memory allocation + + allocate(Ec(nAC,nspin)) + allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc)) + +! Antisymmetrized kernel version + + if(exchange_kernel) then + + write(*,*) + write(*,*) '*** Exchange kernel version ***' + write(*,*) + + end if + + EcAC(:) = 0d0 + Ec(:,:) = 0d0 + +! Compute (singlet) RPA screening + + isp_W = 1 + EcRPA = 0d0 + + ! Memory allocation + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,eW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + +! Spin-conserved manifold + + if(spin_conserved) then + + ispin = 1 + + allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) + + write(*,*) '------------------------' + write(*,*) 'Spin-conserved manifold ' + write(*,*) '------------------------' + write(*,*) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)' + write(*,*) '-----------------------------------------------------------------------------------' + + do iAC=1,nAC + + lambda = rAC(iAC) + + if(doXBS) then + + call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + + end if + + call unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,Omega_sc,XpY_sc,XmY_sc) + + call unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf, & + ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sc,XmY_sc,Ec(iAC,ispin)) + + write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) + + end do + + EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin)) + + if(exchange_kernel) EcAC(ispin) = 0.5d0*EcAC(ispin) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) + write(*,*) '-----------------------------------------------------------------------------------' + write(*,*) + + deallocate(Omega_sc,XpY_sc,XmY_sc) + + end if + +! spin-flip manifold + + if(spin_flip) then + + ispin = 2 + isp_W = 1 + + ! Memory allocation + + nS_ab = (nO(1) - nC(1))*(nV(2) - nR(2)) + nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1)) + nS_sf = nS_ab + nS_ba + + allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) + + write(*,*) '--------------------' + write(*,*) ' Spin-flip manifold ' + write(*,*) '--------------------' + write(*,*) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)' + write(*,*) '-----------------------------------------------------------------------------------' + + do iAC=1,nAC + + lambda = rAC(iAC) + + if(doXBS) then + + call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + + end if + + call unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sc,lambda,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,Omega_sf,XpY_sf,XmY_sf) + + call unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc, & + ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sf,XmY_sf,Ec(iAC,ispin)) + + write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) + + end do + + EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin)) + + if(exchange_kernel) EcAC(ispin) = 1.5d0*EcAC(ispin) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) + write(*,*) '-----------------------------------------------------------------------------------' + write(*,*) + + deallocate(Omega_sf,XpY_sf,XmY_sf) + + end if + +end subroutine unrestricted_ACFDT diff --git a/src/RPA/unrestricted_ACFDT_correlation_energy.f90 b/src/RPA/unrestricted_ACFDT_correlation_energy.f90 new file mode 100644 index 0000000..2984ac2 --- /dev/null +++ b/src/RPA/unrestricted_ACFDT_correlation_energy.f90 @@ -0,0 +1,231 @@ +subroutine unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt, & + ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,XmY,EcAC) + +! Compute the correlation energy via the adiabatic connection formula + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + logical,intent(in) :: exchange_kernel + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nS(nspin) + integer,intent(in) :: nSa + integer,intent(in) :: nSb + integer,intent(in) :: nSt + double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: XpY(nSt,nSt) + double precision,intent(in) :: XmY(nSt,nSt) + +! Local variables + + integer :: i,j,a,b + integer :: ia,jb,kc + double precision :: delta_Kx + double precision,allocatable :: Ap(:,:) + double precision,allocatable :: Bp(:,:) + double precision,allocatable :: X(:,:) + double precision,allocatable :: Y(:,:) + double precision,external :: trace_matrix + +! Output variables + + double precision,intent(out) :: EcAC + +! Exchange kernel + + delta_Kx = 0d0 + if(exchange_kernel) delta_Kx = 1d0 + +! Memory allocation + + allocate(Ap(nSt,nSt),Bp(nSt,nSt),X(nSt,nSt),Y(nSt,nSt)) + +! Compute Aiajb = (ia|bj) and Biajb = (ia|jb) + +! Initialization + + Ap(:,:) = 0d0 + Bp(:,:) = 0d0 + +!----------------------------------------------- +! Build kernel for spin-conserving transitions +!----------------------------------------------- + + if(ispin == 1) then + + ! aaaa block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + Ap(ia,jb) = ERI_aaaa(i,b,a,j) - delta_Kx*ERI_aaaa(i,b,j,a) + Bp(ia,jb) = ERI_aaaa(i,j,a,b) - delta_Kx*ERI_aaaa(i,j,b,a) + + end do + end do + end do + end do + + ! aabb block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + Ap(ia,nSa+jb) = ERI_aabb(i,b,a,j) + Bp(ia,nSa+jb) = ERI_aabb(i,j,a,b) + + + end do + end do + end do + end do + + ! bbaa block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + Ap(nSa+ia,jb) = ERI_aabb(b,i,j,a) + Bp(nSa+ia,jb) = ERI_aabb(j,i,b,a) + + end do + end do + end do + end do + + ! bbbb block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + Ap(nSa+ia,nSa+jb) = ERI_bbbb(i,b,a,j) - delta_Kx*ERI_bbbb(i,b,j,a) + Bp(nSa+ia,nSa+jb) = ERI_bbbb(i,j,a,b) - delta_Kx*ERI_bbbb(i,j,b,a) + + + end do + end do + end do + end do + + end if + +!----------------------------------------------- +! Build A matrix for spin-flip transitions +!----------------------------------------------- + + if(ispin == 2) then + + ! abab block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + Ap(ia,jb) = - delta_Kx*ERI_aabb(i,b,j,a) + + end do + end do + end do + end do + + ! baba block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + Ap(nSa+ia,nSa+jb) = - delta_Kx*ERI_aabb(b,i,a,j) + + end do + end do + end do + end do + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + Bp(ia,nSa+jb) = - delta_Kx*ERI_aabb(i,j,b,a) + + end do + end do + end do + end do + + ! baab block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + Bp(nSa+ia,jb) = - delta_Kx*ERI_aabb(j,i,a,b) + + end do + end do + end do + end do + + end if + +! Compute Tr(K x P_lambda) + + X(:,:) = 0.5d0*(XpY(:,:) + XmY(:,:)) + Y(:,:) = 0.5d0*(XpY(:,:) - XmY(:,:)) + + EcAC = trace_matrix(nS,matmul(X,matmul(Bp,transpose(Y))) + matmul(Y,matmul(Bp,transpose(X)))) & + + trace_matrix(nS,matmul(X,matmul(Ap,transpose(X))) + matmul(Y,matmul(Ap,transpose(Y)))) & + - trace_matrix(nS,Ap) + +end subroutine unrestricted_ACFDT_correlation_energy