diff --git a/examples/molecule.He b/examples/molecule.He index 1ac370f..c78e87e 100644 --- a/examples/molecule.He +++ b/examples/molecule.He @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 2 0 0 0 + 1 1 1 0 0 # Znuc x y z He 0.0 0.0 0.0 diff --git a/input/methods b/input/methods index cad0336..2c30f9f 100644 --- a/input/methods +++ b/input/methods @@ -1,7 +1,7 @@ # RHF UHF MOM F T F # MP2* MP3 MP2-F12 - F F F + T F F # CCD CCSD CCSD(T) F F F # drCCD rCCD lCCD pCCD @@ -9,11 +9,11 @@ # CIS* CIS(D) CID CISD F F F F # RPA* RPAx* ppRPA - T F F + F F F # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0* evGW* qsGW - F F F + T F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index 5bd9d2e..8abb1c0 100644 --- a/input/options +++ b/input/options @@ -5,7 +5,7 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - T T T T F + F T T T T # 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 @@ -13,6 +13,6 @@ # ACFDT: AC Kx XBS T F T # BSE: BSE dBSE dTDA evDyn - T T T F + T F T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/MBPT/G0W0.f90 b/src/MBPT/G0W0.f90 index a69b76b..f629708 100644 --- a/src/MBPT/G0W0.f90 +++ b/src/MBPT/G0W0.f90 @@ -202,9 +202,9 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & if(doACFDT) then - write(*,*) '------------------------------------------------------' - write(*,*) 'Adiabatic connection version of BSE correlation energy' - write(*,*) '------------------------------------------------------' + write(*,*) '--------------------------------------------------------------' + write(*,*) ' Adiabatic connection version of BSE@UG0W0 correlation energy ' + write(*,*) '--------------------------------------------------------------' write(*,*) if(doXBS) then @@ -218,10 +218,10 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & 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 + ERHF + EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy (singlet) =',EcAC(1) + write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy (triplet) =',EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy =',EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/MBPT/UG0W0.f90 b/src/MBPT/UG0W0.f90 index a4badad..f252496 100644 --- a/src/MBPT/UG0W0.f90 +++ b/src/MBPT/UG0W0.f90 @@ -192,10 +192,10 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev 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(*,'(2X,A50,F20.10)') 'Tr@BSE@UG0W0 correlation energy (spin-conserved) =',EcBSE(1) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@UG0W0 correlation energy (spin-flip) =',EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@UG0W0 correlation energy =',EcBSE(1) + EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@UG0W0 total energy =',ENuc + EUHF + EcBSE(1) + EcBSE(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) @@ -215,22 +215,15 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev end if - 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 + call unrestricted_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) 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(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy (spin-conserved) =',EcAC(1) + write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy (spin-flip) =',EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy =',EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 total energy =',ENuc + EUHF + EcAC(1) + EcAC(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/MBPT/evUGW.f90 b/src/MBPT/evUGW.f90 index f322982..afae05e 100644 --- a/src/MBPT/evUGW.f90 +++ b/src/MBPT/evUGW.f90 @@ -1,6 +1,6 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc, & - ERHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eG0W0) + EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -13,7 +13,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE integer,intent(in) :: max_diis double precision,intent(in) :: thresh double precision,intent(in) :: ENuc - double precision,intent(in) :: ERHF + double precision,intent(in) :: EUHF logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel logical,intent(in) :: doXBS @@ -59,7 +59,6 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE double precision :: EcRPA double precision :: EcBSE(nspin) double precision :: EcAC(nspin) - double precision :: EcGM double precision :: alpha double precision,allocatable :: error_diis(:,:,:) double precision,allocatable :: e_diis(:,:,:) @@ -192,7 +191,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE ! Print results - call print_evUGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA) + call print_evUGW(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA) ! Linear mixing or DIIS extrapolation @@ -256,50 +255,51 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE 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,eGW,eGW,EcBSE) -! if(exchange_kernel) then + if(exchange_kernel) then -! EcBSE(1) = 0.5d0*EcBSE(1) -! EcBSE(2) = 1.5d0*EcBSE(2) + EcBSE(1) = 0.5d0*EcBSE(1) + EcBSE(2) = 1.5d0*EcBSE(2) -! end if + end if -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy (singlet) =',EcBSE(1) -! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy (triplet) =',EcBSE(2) -! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2) -! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@BSE@evUGW correlation energy (spin-conserved) =',EcBSE(1) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@evUGW correlation energy (spin-flip) =',EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@evUGW correlation energy =',EcBSE(1) + EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@evUGW 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@evUGW 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,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcAC) + call unrestricted_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,eGW,eGW,EcAC) -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (singlet) =',EcAC(1) -! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (triplet) =',EcAC(2) -! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy =',EcAC(1) + EcAC(2) -! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'AC@BSE@evUGW correlation energy (spin-conserved) =',EcAC(1) + write(*,'(2X,A50,F20.10)') 'AC@BSE@evUGW correlation energy (spin-flip) =',EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@evUGW correlation energy =',EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@evUGW total energy =',ENuc + EUHF + EcAC(1) + EcAC(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) -! end if + end if endif diff --git a/src/MBPT/unrestricted_Bethe_Salpeter.f90 b/src/MBPT/unrestricted_Bethe_Salpeter.f90 index 22aeb6a..5b0436b 100644 --- a/src/MBPT/unrestricted_Bethe_Salpeter.f90 +++ b/src/MBPT/unrestricted_Bethe_Salpeter.f90 @@ -64,6 +64,10 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, nS_aa = nS(1) nS_bb = nS(2) nS_sc = nS_aa + nS_bb + + 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(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin)) @@ -125,10 +129,6 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, EcBSE(ispin) = 0d0 ! 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(OmBSE_sf(nS_sf),XpY_BSE_sf(nS_sf,nS_sf),XmY_BSE_sf(nS_sf,nS_sf)) diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index f2504bf..a7a976d 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -652,12 +652,12 @@ program QuAcK call cpu_time(start_RPA) if(unrestricted) then - call UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & + call URPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb,eHF,cHF,S) else - call dRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) + call RPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) end if call cpu_time(end_RPA) diff --git a/src/RPA/ACFDT.f90 b/src/RPA/ACFDT.f90 index 9853403..b3f52bc 100644 --- a/src/RPA/ACFDT.f90 +++ b/src/RPA/ACFDT.f90 @@ -126,7 +126,6 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB if(triplet) then ispin = 2 - isp_W = 1 write(*,*) '--------------' write(*,*) 'Triplet states' diff --git a/src/RPA/dRPA.f90 b/src/RPA/RPA.f90 similarity index 97% rename from src/RPA/dRPA.f90 rename to src/RPA/RPA.f90 index 02ce6d3..3d6cb63 100644 --- a/src/RPA/dRPA.f90 +++ b/src/RPA/RPA.f90 @@ -1,4 +1,4 @@ -subroutine dRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) +subroutine RPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform a direct random phase approximation calculation @@ -135,4 +135,4 @@ subroutine dRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR end if -end subroutine dRPA +end subroutine RPA diff --git a/src/RPA/UdRPA.f90 b/src/RPA/URPA.f90 similarity index 85% rename from src/RPA/UdRPA.f90 rename to src/RPA/URPA.f90 index b95e685..1a9f222 100644 --- a/src/RPA/UdRPA.f90 +++ b/src/RPA/URPA.f90 @@ -1,4 +1,4 @@ -subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & +subroutine URPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,e,c,S) ! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism @@ -68,7 +68,7 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n ! Initialization EcRPA(:) = 0d0 - EcAC(:) = 0d0 + EcAC(:) = 0d0 ! Spin-conserved transitions @@ -90,6 +90,7 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & c,S,Omega_sc,XpY_sc,XmY_sc) + deallocate(Omega_sc,XpY_sc,XmY_sc) endif @@ -113,6 +114,8 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & c,S,Omega_sf,XpY_sf,XmY_sf) + deallocate(Omega_sf,XpY_sf,XmY_sf) + endif if(exchange_kernel) then @@ -135,30 +138,23 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n if(doACFDT) then - write(*,*) '-------------------------------------------------------' - write(*,*) 'Adiabatic connection version of RPA correlation energy' - write(*,*) '-------------------------------------------------------' + write(*,*) '---------------------------------------------------------' + write(*,*) ' Adiabatic connection version of URPA correlation energy ' + write(*,*) '---------------------------------------------------------' write(*,*) - call unrestricted_ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,spin_conserved,spin_flip,eta, & + call unrestricted_ACFDT(exchange_kernel,.false.,.true.,.false.,TDA,.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 - - EcAC(1) = 0.5d0*EcAC(1) - EcAC(2) = 1.5d0*EcAC(2) - - end if - 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(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy (spin-conserved) =',EcAC(1) + write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy (spin-flip) =',EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy =',EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@URPA total energy =',ENuc + EUHF + EcAC(1) + EcAC(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) end if -end subroutine UdRPA +end subroutine URPA diff --git a/src/RPA/URPAx.f90 b/src/RPA/URPAx.f90 index b3313cc..90f329c 100644 --- a/src/RPA/URPAx.f90 +++ b/src/RPA/URPAx.f90 @@ -139,27 +139,20 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n if(doACFDT) then - write(*,*) '--------------------------------------------------------' - write(*,*) 'Adiabatic connection version of URPAx correlation energy' - write(*,*) '--------------------------------------------------------' + write(*,*) '----------------------------------------------------------' + write(*,*) ' Adiabatic connection version of URPAx correlation energy ' + write(*,*) '----------------------------------------------------------' write(*,*) - call unrestricted_ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,spin_conserved,spin_flip,eta, & + call unrestricted_ACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.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 - - EcAC(1) = 0.5d0*EcAC(1) - EcAC(2) = 1.5d0*EcAC(2) - - end if - 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(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy (spin-conserved) =',EcAC(1) + write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy (spin-flip) =',EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy =',EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@URPAx total energy =',ENuc + EUHF + EcAC(1) + EcAC(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/RPA/unrestricted_ACFDT.f90 b/src/RPA/unrestricted_ACFDT.f90 index 8735cb5..ded29de 100644 --- a/src/RPA/unrestricted_ACFDT.f90 +++ b/src/RPA/unrestricted_ACFDT.f90 @@ -43,7 +43,7 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons double precision,allocatable :: OmRPA(:) double precision,allocatable :: XpY_RPA(:,:) double precision,allocatable :: XmY_RPA(:,:) - double precision,allocatable :: rho_RPA(:,:,:) + double precision,allocatable :: rho_RPA(:,:,:,:) integer :: nS_aa,nS_bb,nS_sc double precision,allocatable :: Omega_sc(:) @@ -62,7 +62,6 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons ! 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 @@ -88,6 +87,12 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons nS_bb = nS(2) nS_sc = nS_aa + nS_bb + 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(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin)) + 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) @@ -122,10 +127,10 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons 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) + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),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)) + 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_sc,XmY_sc,Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) @@ -149,14 +154,9 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons 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(*,*) '--------------------' @@ -180,10 +180,10 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons 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_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,lambda,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),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, & + 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_sf,XmY_sf,Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) diff --git a/src/RPA/unrestricted_ACFDT_correlation_energy.f90 b/src/RPA/unrestricted_ACFDT_correlation_energy.f90 index 2984ac2..3b1e8d9 100644 --- a/src/RPA/unrestricted_ACFDT_correlation_energy.f90 +++ b/src/RPA/unrestricted_ACFDT_correlation_energy.f90 @@ -28,7 +28,7 @@ subroutine unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,n ! Local variables integer :: i,j,a,b - integer :: ia,jb,kc + integer :: ia,jb double precision :: delta_Kx double precision,allocatable :: Ap(:,:) double precision,allocatable :: Bp(:,:) @@ -183,7 +183,9 @@ subroutine unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,n end do end do - ia = 0 + ! abba block + + ia = 0 do i=nC(1)+1,nO(1) do a=nO(2)+1,nBas-nR(2) ia = ia + 1 @@ -224,8 +226,8 @@ subroutine unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,n 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) + EcAC = trace_matrix(nSt,matmul(X,matmul(Bp,transpose(Y))) + matmul(Y,matmul(Bp,transpose(X)))) & + + trace_matrix(nSt,matmul(X,matmul(Ap,transpose(X))) + matmul(Y,matmul(Ap,transpose(Y)))) & + - trace_matrix(nSt,Ap) end subroutine unrestricted_ACFDT_correlation_energy