From 0dfa3b00717d5dfdf0a0badc522b6796c605d251 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 31 Jul 2023 16:01:12 +0200 Subject: [PATCH] clean up in RPA --- input/methods | 2 +- input/options | 2 +- src/CC/CC.f90 | 32 ++++++++++++++++---------------- src/CC/CCSD.f90 | 4 ++-- src/CI/CI.f90 | 16 ++++++++-------- src/GF/GF.f90 | 20 ++++++++++---------- src/GT/GT.f90 | 4 ++-- src/MP/MP.f90 | 8 ++++---- src/RPA/RPA.f90 | 16 ++++++++-------- src/RPA/crRPA.f90 | 33 ++++++++++++++++----------------- src/RPA/phRPA.f90 | 32 +++++++++++++++----------------- src/RPA/phRPAx.f90 | 32 +++++++++++++++----------------- src/RPA/ppRPA.f90 | 12 +++++------- src/RPA/sort_ppRPA.f90 | 30 +++++++++++++++--------------- 14 files changed, 118 insertions(+), 125 deletions(-) diff --git a/input/methods b/input/methods index 91124e8..6c28410 100644 --- a/input/methods +++ b/input/methods @@ -15,5 +15,5 @@ # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW F F F F F F # G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh - T F F F F F + F F F F F T # * unrestricted version available diff --git a/input/options b/input/options index 0b1a49b..316358c 100644 --- a/input/options +++ b/input/options @@ -11,7 +11,7 @@ # GW: maxSCF thresh DIIS n_diis lin eta TDA_W reg 256 0.00001 T 5 T 0.0 F F # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg - 256 0.00001 T 5 T 0.0 F F + 256 0.00001 T 5 T 0.1 T F # ACFDT: AC Kx XBS F F T # BSE: phBSE phBSE2 ppBSE dBSE dTDA diff --git a/src/CC/CC.f90 b/src/CC/CC.f90 index 70460c0..2a90da6 100644 --- a/src/CC/CC.f90 +++ b/src/CC/CC.f90 @@ -42,9 +42,9 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC if(doCCD) then - call cpu_time(start_CC) + call wall_time(start_CC) call CCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) - call cpu_time(end_CC) + call wall_time(end_CC) t_CC = end_CC - start_CC write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CC,' seconds' @@ -58,10 +58,10 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC if(doDCD) then - call cpu_time(start_CC) + call wall_time(start_CC) call DCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR, & ERI,ENuc,EHF,epsHF) - call cpu_time(end_CC) + call wall_time(end_CC) t_CC = end_CC - start_CC write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for DCD = ',t_CC,' seconds' @@ -77,9 +77,9 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC if(doCCSD) then - call cpu_time(start_CC) + call wall_time(start_CC) call CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) - call cpu_time(end_CC) + call wall_time(end_CC) t_CC = end_CC - start_CC write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CC,' seconds' @@ -93,9 +93,9 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC if(do_drCCD) then - call cpu_time(start_CC) + call wall_time(start_CC) call drCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) - call cpu_time(end_CC) + call wall_time(end_CC) t_CC = end_CC - start_CC write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for direct ring CCD = ',t_CC,' seconds' @@ -109,9 +109,9 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC if(do_rCCD) then - call cpu_time(start_CC) + call wall_time(start_CC) call rCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF,epsHF) - call cpu_time(end_CC) + call wall_time(end_CC) t_CC = end_CC - start_CC write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for rCCD = ',t_CC,' seconds' @@ -125,9 +125,9 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC if(do_crCCD) then - call cpu_time(start_CC) + call wall_time(start_CC) call crCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) - call cpu_time(end_CC) + call wall_time(end_CC) t_CC = end_CC - start_CC write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CC,' seconds' @@ -141,9 +141,9 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC if(do_lCCD) then - call cpu_time(start_CC) + call wall_time(start_CC) call lCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) - call cpu_time(end_CC) + call wall_time(end_CC) t_CC = end_CC - start_CC write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ladder CCD = ',t_CC,' seconds' @@ -157,9 +157,9 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC if(dopCCD) then - call cpu_time(start_CC) + call wall_time(start_CC) call pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) - call cpu_time(end_CC) + call wall_time(end_CC) t_CC = end_CC - start_CC write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pair CCD = ',t_CC,' seconds' diff --git a/src/CC/CCSD.f90 b/src/CC/CCSD.f90 index 4ffb335..42d0a90 100644 --- a/src/CC/CCSD.f90 +++ b/src/CC/CCSD.f90 @@ -283,9 +283,9 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBasin,nCin,nOin,nVin,nRin,ERI,EN !------------------------------------------------------------------------ if(doCCSDT) then - call cpu_time(start_CCSDT) + call wall_time(start_CCSDT) call CCSDT(nC,nO,nV,nR,eO,eV,OOVV,VVVO,VOOO,t1,t2,EcCCT) - call cpu_time(end_CCSDT) + call wall_time(end_CCSDT) t_CCSDT = end_CCSDT - start_CCSDT write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for (T) = ',t_CCSDT,' seconds' diff --git a/src/CI/CI.f90 b/src/CI/CI.f90 index 54a0c8e..0dfa91a 100644 --- a/src/CI/CI.f90 +++ b/src/CI/CI.f90 @@ -49,14 +49,14 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin if(doCIS) then - call cpu_time(start_CI) + call wall_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 call CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,epsHF) end if - call cpu_time(end_CI) + call wall_time(end_CI) t_CI = end_CI - start_CI write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CIS = ',t_CI,' seconds' @@ -70,9 +70,9 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin if(doCID) then - call cpu_time(start_CI) + call wall_time(start_CI) call CID(singlet,triplet,nBas,nC,nO,nV,nR,ERI,F,EHF) - call cpu_time(end_CI) + call wall_time(end_CI) t_CI = end_CI - start_CI write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CID = ',t_CI,' seconds' @@ -86,9 +86,9 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin if(doCISD) then - call cpu_time(start_CI) + call wall_time(start_CI) call CISD(singlet,triplet,nBas,nC,nO,nV,nR,ERI,F,EHF) - call cpu_time(end_CI) + call wall_time(end_CI) t_CI = end_CI - start_CI write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CISD = ',t_CI,' seconds' @@ -102,10 +102,10 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin if(doFCI) then - call cpu_time(start_CI) + call wall_time(start_CI) write(*,*) ' FCI is not yet implemented! Sorry.' ! call FCI(nBas,nC,nO,nV,nR,ERI,epsHF) - call cpu_time(end_CI) + call wall_time(end_CI) t_CI = end_CI - start_CI write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for FCI = ',t_CI,' seconds' diff --git a/src/GF/GF.f90 b/src/GF/GF.f90 index 5282856..90174d2 100644 --- a/src/GF/GF.f90 +++ b/src/GF/GF.f90 @@ -75,7 +75,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t if(doG0F2) then - call cpu_time(start_GF) + call wall_time(start_GF) if(unrestricted) then call UG0F2(dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linearize,eta,regularize, & nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, & @@ -84,7 +84,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t call G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) end if - call cpu_time(end_GF) + call wall_time(end_GF) t_GF = end_GF - start_GF write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds' @@ -98,7 +98,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t if(doevGF2) then - call cpu_time(start_GF) + call wall_time(start_GF) if(unrestricted) then call evUGF2(maxSCF,thresh,max_diis,dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip, & eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, & @@ -108,7 +108,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF, & ERI,dipole_int,epsHF) end if - call cpu_time(end_GF) + call wall_time(end_GF) t_GF = end_GF - start_GF write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds' @@ -122,7 +122,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t if(doqsGF2) then - call cpu_time(start_GF) + call wall_time(start_GF) if(unrestricted) then call qsUGF2(maxSCF,thresh,max_diis,dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,regularize, & nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO, & @@ -131,7 +131,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t call qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc, & nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF) end if - call cpu_time(end_GF) + call wall_time(end_GF) t_GF = end_GF - start_GF write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGF2 = ',t_GF,' seconds' @@ -145,13 +145,13 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t if(doG0F3) then - call cpu_time(start_GF) + call wall_time(start_GF) if(unrestricted) then print*,'Unrestricted version of G0F3 not yet implemented! Sorry.' else call G0F3(renorm,nBas,nC,nO,nV,nR,ERI,epsHF) end if - call cpu_time(end_GF) + call wall_time(end_GF) t_GF = end_GF - start_GF write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF,' seconds' @@ -165,13 +165,13 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t if(doevGF3) then - call cpu_time(start_GF) + call wall_time(start_GF) if(unrestricted) then print*,'Unrestricted version of evGF3 not yet implemented! Sorry.' else call evGF3(maxSCF,thresh,max_diis,renorm,nBas,nC,nO,nV,nR,ERI,epsHF) end if - call cpu_time(end_GF) + call wall_time(end_GF) t_GF = end_GF - start_GF write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF,' seconds' diff --git a/src/GT/GT.f90 b/src/GT/GT.f90 index eb48af6..0aad9b5 100644 --- a/src/GT/GT.f90 +++ b/src/GT/GT.f90 @@ -158,8 +158,8 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted else call G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) - call soG0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) +! call soG0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & +! linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) end if call wall_time(end_GT) diff --git a/src/MP/MP.f90 b/src/MP/MP.f90 index 23ffc09..c08a52b 100644 --- a/src/MP/MP.f90 +++ b/src/MP/MP.f90 @@ -37,13 +37,13 @@ subroutine MP(doMP2,doMP3,unrestricted,regularize,nBas,nC,nO,nV,nR,ERI,ERI_aaaa, if(doMP2) then - call cpu_time(start_MP) + call wall_time(start_MP) if(unrestricted) then call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF) else call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) end if - call cpu_time(end_MP) + call wall_time(end_MP) t_MP = end_MP - start_MP write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds' @@ -57,7 +57,7 @@ subroutine MP(doMP2,doMP3,unrestricted,regularize,nBas,nC,nO,nV,nR,ERI,ERI_aaaa, if(doMP3) then - call cpu_time(start_MP) + call wall_time(start_MP) if(unrestricted) then write(*,*) 'MP3 NYI for UHF reference' @@ -65,7 +65,7 @@ subroutine MP(doMP2,doMP3,unrestricted,regularize,nBas,nC,nO,nV,nR,ERI,ERI_aaaa, else call MP3(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) end if - call cpu_time(end_MP) + call wall_time(end_MP) t_MP = end_MP - start_MP write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds' diff --git a/src/RPA/RPA.f90 b/src/RPA/RPA.f90 index bacf31f..a67550d 100644 --- a/src/RPA/RPA.f90 +++ b/src/RPA/RPA.f90 @@ -52,14 +52,14 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted, if(dophRPA) then - call cpu_time(start_RPA) + call wall_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) + call wall_time(end_RPA) t_RPA = end_RPA - start_RPA write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPA = ',t_RPA,' seconds' @@ -73,14 +73,14 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted, if(dophRPAx) then - call cpu_time(start_RPA) + call wall_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) + call wall_time(end_RPA) t_RPA = end_RPA - start_RPA write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPA,' seconds' @@ -94,13 +94,13 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted, if(docrRPA) then - call cpu_time(start_RPA) + call wall_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) + call wall_time(end_RPA) t_RPA = end_RPA - start_RPA write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds' @@ -114,13 +114,13 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted, if(doppRPA) then - call cpu_time(start_RPA) + call wall_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) + call wall_time(end_RPA) t_RPA = end_RPA - start_RPA write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds' diff --git a/src/RPA/crRPA.f90 b/src/RPA/crRPA.f90 index 277a2dc..b5283d7 100644 --- a/src/RPA/crRPA.f90 +++ b/src/RPA/crRPA.f90 @@ -35,8 +35,7 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS double precision,allocatable :: XpY(:,:) double precision,allocatable :: XmY(:,:) - double precision :: EcTr(nspin) - double precision :: EcAC(nspin) + double precision :: EcRPA(nspin) ! Hello world @@ -57,8 +56,8 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS dRPA = .false. - EcTr(:) = 0d0 - EcAC(:) = 0d0 + EcRPA(:) = 0d0 + EcRPA(:) = 0d0 ! Memory allocation @@ -74,7 +73,7 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,e,ERI,Aph) if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph) - call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY) + call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY) call print_excitation_energies('crRPA@HF',ispin,nS,Om) call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) @@ -89,7 +88,7 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,e,ERI,Aph) if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph) - call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY) + call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY) call print_excitation_energies('crRPA@HF',ispin,nS,Om) call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) @@ -97,17 +96,17 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS if(exchange_kernel) then - EcTr(1) = 0.5d0*EcTr(1) - EcTr(2) = 1.5d0*EcTr(2) + EcRPA(1) = 0.5d0*EcRPA(1) + EcRPA(2) = 1.5d0*EcRPA(2) end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (singlet) =',EcTr(1) - write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (triplet) =',EcTr(2) - write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy =',EcTr(1) + EcTr(2) - write(*,'(2X,A50,F20.10)') 'Tr@crRPA total energy =',ENuc + EHF + EcTr(1) + EcTr(2) + write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (singlet) =',EcRPA(1) + write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (triplet) =',EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy =',EcRPA(1) + EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@crRPA total energy =',ENuc + EHF + EcRPA(1) + EcRPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) @@ -120,14 +119,14 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS write(*,*) '-------------------------------------------------------' write(*,*) - call crACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) + call crACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcRPA) write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (singlet) =',EcAC(1) - write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (triplet) =',EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy =',EcAC(1) + EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@crRPA total energy =',ENuc + EHF + EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (singlet) =',EcRPA(1) + write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (triplet) =',EcRPA(2) + write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy =',EcRPA(1) + EcRPA(2) + write(*,'(2X,A50,F20.10)') 'AC@crRPA total energy =',ENuc + EHF + EcRPA(1) + EcRPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/RPA/phRPA.f90 b/src/RPA/phRPA.f90 index b69b18b..7eeb51f 100644 --- a/src/RPA/phRPA.f90 +++ b/src/RPA/phRPA.f90 @@ -35,8 +35,7 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS double precision,allocatable :: XpY(:,:) double precision,allocatable :: XmY(:,:) - double precision :: EcTr(nspin) - double precision :: EcAC(nspin) + double precision :: EcRPA(nspin) ! Hello world @@ -57,8 +56,7 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS dRPA = .true. - EcTr(:) = 0d0 - EcAC(:) = 0d0 + EcRPA(:) = 0d0 ! Memory allocation @@ -74,7 +72,7 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph) if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) - call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY) + call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY) call print_excitation_energies('phRPA@HF',ispin,nS,Om) call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) @@ -89,7 +87,7 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph) if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) - call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY) + call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY) call print_excitation_energies('phRPA@HF',ispin,nS,Om) call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) @@ -97,17 +95,17 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS if(exchange_kernel) then - EcTr(1) = 0.5d0*EcTr(1) - EcTr(2) = 1.5d0*EcTr(2) + EcRPA(1) = 0.5d0*EcRPA(1) + EcRPA(2) = 1.5d0*EcRPA(2) end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy (singlet) =',EcTr(1) - write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy (triplet) =',EcTr(2) - write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy =',EcTr(1) + EcTr(2) - write(*,'(2X,A50,F20.10)') 'Tr@phRPA total energy =',ENuc + EHF + EcTr(1) + EcTr(2) + write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy (singlet) =',EcRPA(1) + write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy (triplet) =',EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy =',EcRPA(1) + EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@phRPA total energy =',ENuc + EHF + EcRPA(1) + EcRPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) @@ -122,14 +120,14 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS write(*,*) '--------------------------------------------------------' write(*,*) - call phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) + call phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcRPA) write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy (singlet) =',EcAC(1) - write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy (triplet) =',EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy =',EcAC(1) + EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@phRPA total energy =',ENuc + EHF + EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy (singlet) =',EcRPA(1) + write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy (triplet) =',EcRPA(2) + write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy =',EcRPA(1) + EcRPA(2) + write(*,'(2X,A50,F20.10)') 'AC@phRPA total energy =',ENuc + EHF + EcRPA(1) + EcRPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/RPA/phRPAx.f90 b/src/RPA/phRPAx.f90 index 676f486..1ae61b7 100644 --- a/src/RPA/phRPAx.f90 +++ b/src/RPA/phRPAx.f90 @@ -35,8 +35,7 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n double precision,allocatable :: XpY(:,:) double precision,allocatable :: XmY(:,:) - double precision :: EcTr(nspin) - double precision :: EcAC(nspin) + double precision :: EcRPA(nspin) ! Hello world @@ -58,8 +57,7 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n dRPA = .false. - EcTr(:) = 0d0 - EcAC(:) = 0d0 + EcRPA(:) = 0d0 ! Memory allocation @@ -75,7 +73,7 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph) if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) - call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY) + call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY) call print_excitation_energies('phRPAx@HF',ispin,nS,Om) call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) @@ -90,7 +88,7 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph) if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) - call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY) + call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY) call print_excitation_energies('phRPAx@HF',ispin,nS,Om) call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) @@ -98,17 +96,17 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n if(exchange_kernel) then - EcTr(1) = 0.5d0*EcTr(1) - EcTr(2) = 1.5d0*EcTr(2) + EcRPA(1) = 0.5d0*EcRPA(1) + EcRPA(2) = 1.5d0*EcRPA(2) end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy (singlet) =',EcTr(1) - write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy (triplet) =',EcTr(2) - write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy =',EcTr(1) + EcTr(2) - write(*,'(2X,A50,F20.10)') 'Tr@phRPAx total energy =',ENuc + EHF + EcTr(1) + EcTr(2) + write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy (singlet) =',EcRPA(1) + write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy (triplet) =',EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy =',EcRPA(1) + EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@phRPAx total energy =',ENuc + EHF + EcRPA(1) + EcRPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) @@ -125,14 +123,14 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n write(*,*) '---------------------------------------------------------' write(*,*) - call phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) + call phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcRPA) write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy (singlet) =',EcAC(1) - write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy (triplet) =',EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy =',EcAC(1) + EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@phRPAx total energy =',ENuc + EHF + EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy (singlet) =',EcRPA(1) + write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy (triplet) =',EcRPA(2) + write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy =',EcRPA(1) + EcRPA(2) + write(*,'(2X,A50,F20.10)') 'AC@phRPAx total energy =',ENuc + EHF + EcRPA(1) + EcRPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/RPA/ppRPA.f90 b/src/RPA/ppRPA.f90 index a63f3cc..9da9357 100644 --- a/src/RPA/ppRPA.f90 +++ b/src/RPA/ppRPA.f90 @@ -38,7 +38,6 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipol double precision,allocatable :: Y2(:,:) double precision :: EcRPA(nspin) - double precision :: EcAC(nspin) ! Hello world @@ -51,7 +50,6 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipol ! Initialization EcRPA(:) = 0d0 - EcAC(:) = 0d0 ! Singlet manifold @@ -135,14 +133,14 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipol write(*,*) '--------------------------------------------------------' write(*,*) - call ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ERI,e,EcAC) + call ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ERI,e,EcRPA) write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (singlet) =',EcAC(1),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (triplet) =',EcAC(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy =',EcAC(1) + EcAC(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy =',ENuc + EHF + EcAC(1) + EcAC(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (singlet) =',EcRPA(1),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (triplet) =',EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy =',EcRPA(1) + EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy =',ENuc + EHF + EcRPA(1) + EcRPA(2),' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/RPA/sort_ppRPA.f90 b/src/RPA/sort_ppRPA.f90 index c5e2ca6..9b61219 100644 --- a/src/RPA/sort_ppRPA.f90 +++ b/src/RPA/sort_ppRPA.f90 @@ -1,4 +1,4 @@ -subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) +subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) ! Compute the metric matrix for pp-RPA @@ -9,7 +9,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) integer,intent(in) :: nOO integer,intent(in) :: nVV - double precision,intent(in) :: Omega(nOO+nVV) + double precision,intent(in) :: Om(nOO+nVV) double precision,intent(in) :: Z(nOO+nVV,nOO+nVV) ! Local variables @@ -32,10 +32,10 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) ! Output variables - double precision,intent(out) :: Omega1(nVV) + double precision,intent(out) :: Om1(nVV) double precision,intent(out) :: X1(nVV,nVV) double precision,intent(out) :: Y1(nOO,nVV) - double precision,intent(out) :: Omega2(nOO) + double precision,intent(out) :: Om2(nOO) double precision,intent(out) :: X2(nVV,nOO) double precision,intent(out) :: Y2(nOO,nOO) @@ -47,11 +47,11 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) ! Initializatiom - Omega1(:) = 0d0 + Om1(:) = 0d0 X1(:,:) = 0d0 Y1(:,:) = 0d0 - Omega2(:) = 0d0 + Om2(:) = 0d0 X2(:,:) = 0d0 Y2(:,:) = 0d0 @@ -74,24 +74,24 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) do pq=1,nOO+nVV - if(Omega(pq) > 0d0) then + if(Om(pq) > 0d0) then ab = ab + 1 - Omega1(ab) = Omega(pq) + Om1(ab) = Om(pq) Z1(1:nOO+nVV,ab) = Z(1:nOO+nVV,pq) else ij = ij + 1 - Omega2(ij) = Omega(pq) + Om2(ij) = Om(pq) Z2(1:nOO+nVV,ij) = Z(1:nOO+nVV,pq) end if end do - if(minval(Omega1(:)) < 0d0 .or. ab /= nVV) call print_warning('You may have instabilities in pp-RPA!!') - if(maxval(Omega2(:)) > 0d0 .or. ij /= nOO) call print_warning('You may have instabilities in pp-RPA!!') + if(minval(Om1(:)) < 0d0 .or. ab /= nVV) call print_warning('You may have instabilities in pp-RPA!!') + if(maxval(Om2(:)) > 0d0 .or. ij /= nOO) call print_warning('You may have instabilities in pp-RPA!!') if(nVV > 0) then @@ -100,7 +100,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) order1(ab) = ab end do - call quick_sort(Omega1(:),order1(:),nVV) + call quick_sort(Om1(:),order1(:),nVV) call set_order(Z1(:,:),order1(:),nOO+nVV,nVV) end if @@ -111,7 +111,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) order2(ij) = ij end do - call quick_sort(Omega2(:),order2(:),nOO) + call quick_sort(Om2(:),order2(:),nOO) call set_order(Z2(:,:),order2(:),nOO+nVV,nOO) end if @@ -123,7 +123,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) ! do ab=2,nVV -! if(abs(Omega1(ab-1) - Omega1(ab)) < 1d-6) then +! if(abs(Om1(ab-1) - Om1(ab)) < 1d-6) then ! deg1 = deg1 + 1 @@ -164,7 +164,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) ! do ij=2,nOO -! if(abs(Omega2(ij-1) - Omega2(ij)) < 1d-6) then +! if(abs(Om2(ij-1) - Om2(ij)) < 1d-6) then ! deg2 = deg2 + 1