From ed084506616ae6e6d9fb34fbd1b12388f61e2b1c Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sun, 23 Jul 2023 11:16:42 +0200 Subject: [PATCH] RPA module --- input/methods | 6 +- src/CI/CI.f90 | 8 +- src/GW/UG0W0.f90 | 4 +- src/{RPA/UACFDT.f90 => GW/UGW_phACFDT.f90} | 22 +-- src/GW/evUGW.f90 | 4 +- src/GW/qsUGW.f90 | 4 +- src/LR/phULR.f90 | 7 +- src/QuAcK/QuAcK.f90 | 119 +++------------- src/RPA/RPA.f90 | 131 ++++++++++++++++++ src/{GW/UGW_ACFDT.f90 => RPA/phUACFDT.f90} | 23 ++- ...gy.f90 => phUACFDT_correlation_energy.f90} | 4 +- src/RPA/{URPA.f90 => phURPA.f90} | 13 +- src/RPA/{URPAx.f90 => phURPAx.f90} | 13 +- 13 files changed, 200 insertions(+), 158 deletions(-) rename src/{RPA/UACFDT.f90 => GW/UGW_phACFDT.f90} (84%) create mode 100644 src/RPA/RPA.f90 rename src/{GW/UGW_ACFDT.f90 => RPA/phUACFDT.f90} (84%) rename src/RPA/{UACFDT_correlation_energy.f90 => phUACFDT_correlation_energy.f90} (96%) rename src/RPA/{URPA.f90 => phURPA.f90} (90%) rename src/RPA/{URPAx.f90 => phURPAx.f90} (89%) diff --git a/input/methods b/input/methods index 6878a71..4da1670 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF RMOM UMOM KS - F F F F F + F T F F F # MP2* MP3 F F # CCD pCCD DCD CCSD CCSD(T) @@ -7,9 +7,9 @@ # drCCD rCCD crCCD lCCD F F F F # CIS* CIS(D) CID CISD FCI - T F F F F + F F F F F # phRPA* phRPAx* crRPA ppRPA - F F F F + T T T T # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW diff --git a/src/CI/CI.f90 b/src/CI/CI.f90 index 48b268c..54a0c8e 100644 --- a/src/CI/CI.f90 +++ b/src/CI/CI.f90 @@ -50,18 +50,12 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin if(doCIS) then call cpu_time(start_CI) - if(unrestricted) then - call UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb, & ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S) - - else - + else call CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,epsHF) - end if - call cpu_time(end_CI) t_CI = end_CI - start_CI diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index 20aa140..79600f1 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -217,8 +217,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons end if - call UGW_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) + call UGW_phACFDT(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(*,*)'-------------------------------------------------------------------------------' diff --git a/src/RPA/UACFDT.f90 b/src/GW/UGW_phACFDT.f90 similarity index 84% rename from src/RPA/UACFDT.f90 rename to src/GW/UGW_phACFDT.f90 index dacfed8..63fe4df 100644 --- a/src/RPA/UACFDT.f90 +++ b/src/GW/UGW_phACFDT.f90 @@ -1,5 +1,5 @@ -subroutine UACFDT(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) +subroutine UGW_phACFDT(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 @@ -93,7 +93,7 @@ subroutine UACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin_f allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin)) - call phULR(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,eW, & + call phULR(isp_W,.true.,TDA_W,.false.,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) @@ -120,17 +120,17 @@ subroutine UACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin_f if(doXBS) then - call phULR(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & + call phULR(isp_W,.true.,TDA_W,.false.,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 phULR(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,e, & + call phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Om_sc,XpY_sc,XmY_sc) - call UACFDT_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)) + call phUACFDT_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) @@ -174,17 +174,17 @@ subroutine UACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin_f if(doXBS) then - call phULR(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & + call phULR(isp_W,.true.,TDA_W,.false.,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 phULR(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,lambda,e, & + call phULR(ispin,dRPA,TDA,BSE,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),Om_sf,XpY_sf,XmY_sf) - call UACFDT_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)) + call phUACFDT_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/GW/evUGW.f90 b/src/GW/evUGW.f90 index 19690e7..e36bf0f 100644 --- a/src/GW/evUGW.f90 +++ b/src/GW/evUGW.f90 @@ -262,8 +262,8 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, end if - call UGW_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) + call UGW_phACFDT(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(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90 index 2f0151b..a835f52 100644 --- a/src/GW/qsUGW.f90 +++ b/src/GW/qsUGW.f90 @@ -401,8 +401,8 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, end if - call UGW_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) + call UGW_phACFDT(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(*,*)'-------------------------------------------------------------------------------' diff --git a/src/LR/phULR.f90 b/src/LR/phULR.f90 index d3bcd82..53080b1 100644 --- a/src/LR/phULR.f90 +++ b/src/LR/phULR.f90 @@ -1,4 +1,4 @@ -subroutine phULR(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, & +subroutine phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,Om,XpY,XmY) ! Compute linear response for unrestricted formalism @@ -12,7 +12,6 @@ subroutine phULR(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambd logical,intent(in) :: dRPA logical,intent(in) :: TDA logical,intent(in) :: BSE - double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC(nspin) integer,intent(in) :: nO(nspin) @@ -58,7 +57,7 @@ subroutine phULR(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambd call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) if(BSE) & - call unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, & + call unrestricted_Bethe_Salpeter_A_matrix(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,Aph) ! Tamm-Dancoff approximation @@ -76,7 +75,7 @@ subroutine phULR(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambd call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) if(BSE) & - call unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & + call unrestricted_Bethe_Salpeter_B_matrix(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,Bph) ! Build A + B and A - B matrices diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 146ae4f..f0839d7 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -9,9 +9,9 @@ program QuAcK logical :: doKS logical :: doMP,doMP2,doMP3 logical :: doCC,doCCD,dopCCD,doDCD,doCCSD,doCCSDT - logical :: do_drCCD,do_rCCD,do_crCCD,do_lCCD + logical :: dodrCCD,dorCCD,docrCCD,dolCCD logical :: doCI,doCIS,doCIS_D,doCID,doCISD,doFCI - logical :: dophRPA,dophRPAx,docrRPA,doppRPA + logical :: doRPA,dophRPA,dophRPAx,docrRPA,doppRPA logical :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW logical :: doG0T0pp,doevGTpp,doqsGTpp @@ -122,7 +122,7 @@ program QuAcK call read_methods(doRHF,doUHF,doRMOM,doUMOM,doKS, & doMP2,doMP3, & doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & - do_drCCD,do_rCCD,do_crCCD,do_lCCD, & + dodrCCD,dorCCD,docrCCD,dolCCD, & doCIS,doCIS_D,doCID,doCISD,doFCI, & dophRPA,dophRPAx,docrRPA,doppRPA, & doG0F2,doevGF2,doqsGF2, & @@ -165,9 +165,6 @@ program QuAcK call read_geometry(nNuc,ZNuc,rNuc,ENuc) -! allocate(CenterShell(maxShell,ncart),TotAngMomShell(maxShell),KShell(maxShell),DShell(maxShell,maxK), & -! ExpShell(maxShell,maxK),max_ang_mom(nNuc),min_exponent(nNuc,maxL+1),max_exponent(nNuc)) - !------------------------------------------------------------------------ ! Read basis set information from PySCF !------------------------------------------------------------------------ @@ -181,7 +178,7 @@ program QuAcK ! Memory allocation for one- and two-electron integrals - allocate(cHF(nBas,nBas,nspin),epsHF(nBas,nspin),PHF(nBas,nBas,nspin),S(nBas,nBas),T(nBas,nBas), & + allocate(cHF(nBas,nBas,nspin),epsHF(nBas,nspin),PHF(nBas,nBas,nspin),S(nBas,nBas),T(nBas,nBas), & V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas),dipole_int_AO(nBas,nBas,ncart), & dipole_int_MO(nBas,nBas,ncart),F_AO(nBas,nBas)) @@ -194,10 +191,10 @@ program QuAcK call wall_time(end_int) - t_int = end_int - start_int - write(*,*) - write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading integrals = ',t_int,' seconds' - write(*,*) + t_int = end_int - start_int + write(*,*) + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading integrals = ',t_int,' seconds' + write(*,*) ! Compute orthogonalization matrix @@ -337,17 +334,11 @@ program QuAcK if(dostab) then call cpu_time(start_stab) - if(unrestricted) then - call UHF_stability(nBas,nC,nO,nV,nR,nS,epsHF,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb) - else - call RHF_stability(nBas,nC,nO,nV,nR,nS,epsHF,ERI_MO) - end if - call cpu_time(end_stab) t_stab = end_stab - start_stab @@ -378,12 +369,13 @@ program QuAcK ! Coupled-cluster module !------------------------------------------------------------------------ - doCC = doCCD .or. dopCCD .or. doDCD .or. doCCSD .or. doCCSDT + doCC = doCCD .or. dopCCD .or. doDCD .or. doCCSD .or. doCCSDT .or. & + dodrCCD .or. dorCCD .or. docrCCD .or. dolCCD if(doCC) then call cpu_time(start_CC) - call CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCCD, & + call CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) call cpu_time(end_CC) @@ -402,7 +394,7 @@ program QuAcK if(doCI) then call cpu_time(start_CI) - call CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin_conserved,spin_flip, & + call CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin_conserved,spin_flip, & nBas,nC,nO,nV,nR,nS,ERI_MO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_MO,dipole_int_aa,dipole_int_bb, & epsHF,EHF,cHF,S,F_MO) call cpu_time(end_CI) @@ -414,22 +406,18 @@ program QuAcK end if !------------------------------------------------------------------------ -! Compute (direct) RPA excitations +! Random-phase approximation module !------------------------------------------------------------------------ - if(dophRPA) then + doRPA = dophRPA .or. dophRPAx .or. docrRPA .or. doppRPA + + if(doRPA) then call cpu_time(start_RPA) - if(unrestricted) then - - call URPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EHF, & - ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S) - - else - - call phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF) - - end if + call RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted, & + TDA,doACFDT,exchange_kernel,singlet,triplet,spin_conserved,spin_flip, & + nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & + dipole_int_MO,dipole_int_aa,dipole_int_bb,epsHF,cHF,S) call cpu_time(end_RPA) t_RPA = end_RPA - start_RPA @@ -438,73 +426,6 @@ program QuAcK end if -!------------------------------------------------------------------------ -! Compute RPAx (RPA with exchange) excitations -!------------------------------------------------------------------------ - - if(dophRPAx) then - - call cpu_time(start_RPA) - if(unrestricted) then - - call URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EHF, & - ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S) - - else - - call phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF) - - end if - call cpu_time(end_RPA) - - t_RPA = end_RPA - start_RPA - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPA,' seconds' - write(*,*) - - end if - -!------------------------------------------------------------------------ -! Compute crRPA excitations -!------------------------------------------------------------------------ - - if(docrRPA) then - - call cpu_time(start_RPA) - call crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF) - call cpu_time(end_RPA) - - t_RPA = end_RPA - start_RPA - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds' - write(*,*) - - end if - -!------------------------------------------------------------------------ -! Compute ppRPA excitations -!------------------------------------------------------------------------ - - if(doppRPA) then - - call cpu_time(start_RPA) - - if(unrestricted) then - - call ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EHF,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,epsHF) - - else - - call ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF) - - end if - - call cpu_time(end_RPA) - - t_RPA = end_RPA - start_RPA - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds' - write(*,*) - - end if - !------------------------------------------------------------------------ ! Compute G0F2 electronic binding energies !------------------------------------------------------------------------ diff --git a/src/RPA/RPA.f90 b/src/RPA/RPA.f90 new file mode 100644 index 0000000..bacf31f --- /dev/null +++ b/src/RPA/RPA.f90 @@ -0,0 +1,131 @@ +subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted, & + TDA,doACFDT,exchange_kernel,singlet,triplet,spin_conserved,spin_flip, & + nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb, & + dipole_int,dipole_int_aa,dipole_int_bb,epsHF,cHF,S) + +! Random-phase approximation module + + implicit none + include 'parameters.h' + +! Input variables + + logical :: dophRPA + logical :: dophRPAx + logical :: docrRPA + logical :: doppRPA + logical :: unrestricted + + logical,intent(in) :: TDA + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: singlet + logical,intent(in) :: triplet + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + 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) :: ENuc + double precision,intent(in) :: EHF + double precision,intent(in) :: epsHF(nBas,nspin) + double precision,intent(in) :: cHF(nBas,nBas,nspin) + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + 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) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) + +! Local variables + + double precision :: start_RPA ,end_RPA ,t_RPA + +!------------------------------------------------------------------------ +! Compute (direct) RPA excitations +!------------------------------------------------------------------------ + + if(dophRPA) then + + call cpu_time(start_RPA) + if(unrestricted) then + call phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, & + ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S) + else + call phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) + end if + call cpu_time(end_RPA) + + t_RPA = end_RPA - start_RPA + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPA = ',t_RPA,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Compute RPAx (RPA with exchange) excitations +!------------------------------------------------------------------------ + + if(dophRPAx) then + + call cpu_time(start_RPA) + if(unrestricted) then + call phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, & + ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S) + else + call phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) + end if + call cpu_time(end_RPA) + + t_RPA = end_RPA - start_RPA + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPA,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Compute crRPA excitations +!------------------------------------------------------------------------ + + if(docrRPA) then + + call cpu_time(start_RPA) + if(unrestricted) then + write(*,*) 'Unrestricted version of crRPA not yet implemented! Sorry.' + else + call crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) + end if + call cpu_time(end_RPA) + + t_RPA = end_RPA - start_RPA + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Compute ppRPA excitations +!------------------------------------------------------------------------ + + if(doppRPA) then + + call cpu_time(start_RPA) + if(unrestricted) then + call ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,epsHF) + else + call ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,epsHF) + end if + call cpu_time(end_RPA) + + t_RPA = end_RPA - start_RPA + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds' + write(*,*) + + end if + +end subroutine diff --git a/src/GW/UGW_ACFDT.f90 b/src/RPA/phUACFDT.f90 similarity index 84% rename from src/GW/UGW_ACFDT.f90 rename to src/RPA/phUACFDT.f90 index d7eaf5c..e913d55 100644 --- a/src/GW/UGW_ACFDT.f90 +++ b/src/RPA/phUACFDT.f90 @@ -1,5 +1,5 @@ -subroutine UGW_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) +subroutine phUACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin_flip, & + 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 @@ -18,7 +18,6 @@ subroutine UGW_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spi 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) @@ -93,7 +92,7 @@ subroutine UGW_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spi allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin)) - call phULR(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,eW, & + call phULR(isp_W,.true.,TDA_W,.false.,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) @@ -120,17 +119,17 @@ subroutine UGW_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spi if(doXBS) then - call phULR(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & + call phULR(isp_W,.true.,TDA_W,.false.,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 phULR(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,e, & + call phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Om_sc,XpY_sc,XmY_sc) - call UACFDT_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)) + call phUACFDT_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) @@ -174,17 +173,17 @@ subroutine UGW_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spi if(doXBS) then - call phULR(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & + call phULR(isp_W,.true.,TDA_W,.false.,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 phULR(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,lambda,e, & + call phULR(ispin,dRPA,TDA,BSE,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),Om_sf,XpY_sf,XmY_sf) - call UACFDT_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)) + call phUACFDT_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/UACFDT_correlation_energy.f90 b/src/RPA/phUACFDT_correlation_energy.f90 similarity index 96% rename from src/RPA/UACFDT_correlation_energy.f90 rename to src/RPA/phUACFDT_correlation_energy.f90 index b788960..a1b89b5 100644 --- a/src/RPA/UACFDT_correlation_energy.f90 +++ b/src/RPA/phUACFDT_correlation_energy.f90 @@ -1,5 +1,5 @@ -subroutine UACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt, & - ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,XmY,EcAC) +subroutine phUACFDT_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 diff --git a/src/RPA/URPA.f90 b/src/RPA/phURPA.f90 similarity index 90% rename from src/RPA/URPA.f90 rename to src/RPA/phURPA.f90 index eb8b81d..b457f32 100644 --- a/src/RPA/URPA.f90 +++ b/src/RPA/phURPA.f90 @@ -1,5 +1,5 @@ -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) +subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,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 @@ -20,7 +20,6 @@ subroutine URPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC integer,intent(in) :: nV(nspin) integer,intent(in) :: nR(nspin) integer,intent(in) :: nS(nspin) - double precision,intent(in) :: eta double precision,intent(in) :: ENuc double precision,intent(in) :: EUHF double precision,intent(in) :: e(nBas,nspin) @@ -83,7 +82,7 @@ subroutine URPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) - call phULR(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & + call phULR(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) call print_excitation('URPA ',5,nS_sc,Omega_sc) call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & @@ -107,7 +106,7 @@ subroutine URPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) - call phULR(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, & + call phULR(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sf,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf) call print_excitation('URPA ',6,nS_sf,Omega_sf) call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & @@ -142,8 +141,8 @@ subroutine URPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC write(*,*) '---------------------------------------------------------' write(*,*) - call UACFDT(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,EcRPA) + call phUACFDT(exchange_kernel,.false.,.true.,.false.,TDA,.false.,spin_conserved,spin_flip, & + nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,e,e,EcRPA) if(exchange_kernel) then diff --git a/src/RPA/URPAx.f90 b/src/RPA/phURPAx.f90 similarity index 89% rename from src/RPA/URPAx.f90 rename to src/RPA/phURPAx.f90 index eb1b081..3ed0867 100644 --- a/src/RPA/URPAx.f90 +++ b/src/RPA/phURPAx.f90 @@ -1,5 +1,5 @@ -subroutine URPAx(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) +subroutine phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,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 @@ -20,7 +20,6 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n integer,intent(in) :: nV(nspin) integer,intent(in) :: nR(nspin) integer,intent(in) :: nS(nspin) - double precision,intent(in) :: eta double precision,intent(in) :: ENuc double precision,intent(in) :: EUHF double precision,intent(in) :: e(nBas,nspin) @@ -84,7 +83,7 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) - call phULR(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & + call phULR(ispin,.false.,TDA,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) call print_excitation('URPAx ',5,nS_sc,Omega_sc) call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & @@ -108,7 +107,7 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) - call phULR(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, & + call phULR(ispin,.false.,TDA,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sf,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf) call print_excitation('URPAx ',6,nS_sf,Omega_sf) call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & @@ -147,8 +146,8 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n write(*,*) '----------------------------------------------------------' write(*,*) - call UACFDT(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,EcRPA) + call phUACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,spin_conserved,spin_flip, & + nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,e,e,EcRPA) write(*,*) write(*,*)'-------------------------------------------------------------------------------'