From c25e934e8b1a010f166e4f06388bc6063bb2ecb8 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Thu, 29 Aug 2024 00:47:12 +0200 Subject: [PATCH] introduce nBas_MOs in RGT --- src/GT/RGT.f90 | 75 ++++++++++---------- src/GT/print_qsRGTeh.f90 | 24 ++++--- src/GT/print_qsRGTpp.f90 | 24 ++++--- src/GT/qsRGTeh.f90 | 110 ++++++++++++++++++----------- src/GT/qsRGTpp.f90 | 147 +++++++++++++++++++++++---------------- src/QuAcK/RQuAcK.f90 | 10 +-- 6 files changed, 229 insertions(+), 161 deletions(-) diff --git a/src/GT/RGT.f90 b/src/GT/RGT.f90 index 4225d36..a4f46b6 100644 --- a/src/GT/RGT.f90 +++ b/src/GT/RGT.f90 @@ -1,7 +1,11 @@ -subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis, & - doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & - linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, & - ERI_AO,ERI_MO,dipole_int_AO,dipole_int,PHF,cHF,eHF) + +! --- + +subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevGTeh, doqsGTeh, & + maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, & + doppBSE, TDA_T, TDA, dBSE, dTDA, singlet, triplet, linearize, eta, regularize, & + nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, S, X, T, & + V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! T-matrix module @@ -44,7 +48,7 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas + integer,intent(in) :: nBas_AOs, nBas_MOs integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -52,18 +56,18 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do integer,intent(in) :: nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: PHF(nBas,nBas) - double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: T(nBas,nBas) - double precision,intent(in) :: V(nBas,nBas) - double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) - double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: eHF(nBas_MOs) + double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs) + double precision,intent(in) :: PHF(nBas_AOs,nBas_AOs) + double precision,intent(in) :: S(nBas_AOs,nBas_AOs) + double precision,intent(in) :: T(nBas_AOs,nBas_AOs) + double precision,intent(in) :: V(nBas_AOs,nBas_AOs) + double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs) + double precision,intent(in) :: X(nBas_AOs,nBas_MOs) + double precision,intent(in) :: ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs) + double precision,intent(in) :: ERI_MO(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs) + double precision,intent(in) :: dipole_int_AO(nBas_AOs,nBas_AOs,ncart) + double precision,intent(in) :: dipole_int_MO(nBas_MOs,nBas_MOs,ncart) ! Local variables @@ -78,11 +82,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do call wall_time(start_GT) call RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) + linearize,eta,regularize,nBas_MOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0pp = ',t_GT,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for G0T0pp = ',t_GT,' seconds' write(*,*) end if @@ -95,11 +99,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do call wall_time(start_GT) call evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & - linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) + linearize,eta,regularize,nBas_MOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGTpp = ',t_GT,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for evGTpp = ',t_GT,' seconds' write(*,*) end if @@ -111,13 +115,13 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do if(doqsGTpp) then call wall_time(start_GT) - call qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & - eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int, & - PHF,cHF,eHF) + call qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, TDA_T, TDA, dBSE, & + dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, & + nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) call wall_time(end_GT) t_GT = end_GT - start_GT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGTpp = ',t_GT,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGTpp = ',t_GT,' seconds' write(*,*) end if @@ -129,11 +133,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do if(doufG0T0pp) then call wall_time(start_GT) - call ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + call ufG0T0pp(dotest,TDA_T,nBas_MOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0T0pp = ',t_GT,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ufG0T0pp = ',t_GT,' seconds' write(*,*) end if @@ -146,11 +150,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do call wall_time(start_GT) call RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) + linearize,eta,regularize,nBas_MOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0eh = ',t_GT,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for G0T0eh = ',t_GT,' seconds' write(*,*) end if @@ -163,11 +167,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do call wall_time(start_GT) call evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & - singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) + singlet,triplet,linearize,eta,regularize,nBas_MOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGTeh = ',t_GT,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for evGTeh = ',t_GT,' seconds' write(*,*) end if @@ -179,13 +183,14 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do if(doqsGTeh) then call wall_time(start_GT) - call qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & - eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int, & - PHF,cHF,eHF) + call qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, & + dophBSE2, TDA_T, TDA, dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, & + ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, S, X, T, V, & + Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) call wall_time(end_GT) t_GT = end_GT - start_GT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGTeh = ',t_GT,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGTeh = ',t_GT,' seconds' write(*,*) end if diff --git a/src/GT/print_qsRGTeh.f90 b/src/GT/print_qsRGTeh.f90 index 9f8f266..72ffe99 100644 --- a/src/GT/print_qsRGTeh.f90 +++ b/src/GT/print_qsRGTeh.f90 @@ -1,4 +1,8 @@ -subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) + +! --- + +subroutine print_qsRGTeh(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c, SigC, & + Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGT, dipole) ! Print one-electron energies and other stuff for qsGTeh @@ -7,7 +11,7 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas_AOs, nBas_MOs integer,intent(in) :: nO integer,intent(in) :: nSCF double precision,intent(in) :: ENuc @@ -19,11 +23,11 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ double precision,intent(in) :: EcRPA(nspin) double precision,intent(in) :: Conv double precision,intent(in) :: thresh - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: eGT(nBas) - double precision,intent(in) :: c(nBas) - double precision,intent(in) :: SigC(nBas,nBas) - double precision,intent(in) :: Z(nBas) + double precision,intent(in) :: eHF(nBas_MOs) + double precision,intent(in) :: eGT(nBas_MOs) + double precision,intent(in) :: c(nBas_AOs,nBas_MOs) + double precision,intent(in) :: SigC(nBas_MOs,nBas_MOs) + double precision,intent(in) :: Z(nBas_MOs) double precision,intent(in) :: EqsGT double precision,intent(in) :: dipole(ncart) @@ -58,7 +62,7 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ '|','#','|','e_HF (eV)','|','Sig_GTeh (eV)','|','Z','|','e_GTeh (eV)','|' write(*,*)'-------------------------------------------------------------------------------' - do p=1,nBas + do p=1,nBas_MOs write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|' end do @@ -109,13 +113,13 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGTeh MO coefficients' write(*,'(A50)') '---------------------------------------' - call matout(nBas,nBas,c) + call matout(nBas_AOs, nBas_MOs, c) write(*,*) end if write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGTeh MO energies' write(*,'(A50)') '---------------------------------------' - call vecout(nBas,eGT) + call vecout(nBas_MOs, eGT) write(*,*) end if diff --git a/src/GT/print_qsRGTpp.f90 b/src/GT/print_qsRGTpp.f90 index c3bebfa..9d1479d 100644 --- a/src/GT/print_qsRGTpp.f90 +++ b/src/GT/print_qsRGTpp.f90 @@ -1,4 +1,8 @@ -subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) + +! --- + +subroutine print_qsRGTpp(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c, SigC, Z, & + ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGT, dipole) ! Print one-electron energies and other stuff for qsGT @@ -7,7 +11,7 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas_AOs, nBas_MOs integer,intent(in) :: nO integer,intent(in) :: nSCF double precision,intent(in) :: ENuc @@ -19,11 +23,11 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ double precision,intent(in) :: EcRPA(nspin) double precision,intent(in) :: Conv double precision,intent(in) :: thresh - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: eGT(nBas) - double precision,intent(in) :: c(nBas) - double precision,intent(in) :: SigC(nBas,nBas) - double precision,intent(in) :: Z(nBas) + double precision,intent(in) :: eHF(nBas_MOs) + double precision,intent(in) :: eGT(nBas_MOs) + double precision,intent(in) :: c(nBas_AOs,nBas_MOs) + double precision,intent(in) :: SigC(nBas_MOs,nBas_MOs) + double precision,intent(in) :: Z(nBas_MOs) double precision,intent(in) :: EqsGT double precision,intent(in) :: dipole(ncart) @@ -58,7 +62,7 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ '|','#','|','e_HF (eV)','|','Sig_GTpp (eV)','|','Z','|','e_GTpp (eV)','|' write(*,*)'-------------------------------------------------------------------------------' - do p=1,nBas + do p=1,nBas_MOs write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|' end do @@ -109,13 +113,13 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGTpp MO coefficients' write(*,'(A50)') '---------------------------------------' - call matout(nBas,nBas,c) + call matout(nBas_AOs, nBas_MOs, c) write(*,*) end if write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGTpp MO energies' write(*,'(A50)') '---------------------------------------' - call vecout(nBas,eGT) + call vecout(nBas_MOs, eGT) write(*,*) end if diff --git a/src/GT/qsRGTeh.f90 b/src/GT/qsRGTeh.f90 index 889e752..e608d9a 100644 --- a/src/GT/qsRGTeh.f90 +++ b/src/GT/qsRGTeh.f90 @@ -1,6 +1,10 @@ -subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA, & - dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & - S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + +! --- + +subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, & + dophBSE2, TDA_T, TDA, dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, & + ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, S, X, T, V, & + Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! Perform a quasiparticle self-consistent GTeh calculation @@ -33,31 +37,31 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas + integer,intent(in) :: nBas_AOs, nBas_MOs integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: PHF(nBas,nBas) - double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: T(nBas,nBas) - double precision,intent(in) :: V(nBas,nBas) - double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) - double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) - double precision,intent(in) :: dipole_int_MO(nBas,nBas,ncart) + double precision,intent(in) :: eHF(nBas_MOs) + double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs) + double precision,intent(in) :: PHF(nBas_AOs,nBas_AOs) + double precision,intent(in) :: S(nBas_AOs,nBas_AOs) + double precision,intent(in) :: T(nBas_AOs,nBas_AOs) + double precision,intent(in) :: V(nBas_AOs,nBas_AOs) + double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs) + double precision,intent(in) :: X(nBas_AOs,nBas_MOs) + double precision,intent(in) :: ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs) + double precision,intent(inout):: ERI_MO(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs) + double precision,intent(in) :: dipole_int_AO(nBas_AOs,nBas_AOs,ncart) + double precision,intent(in) :: dipole_int_MO(nBas_MOs,nBas_MOs,ncart) ! Local variables logical :: dRPA = .false. integer :: nSCF - integer :: nBasSq + integer :: nBas_AOs_Sq integer :: ispin integer :: n_diis double precision :: ET @@ -113,7 +117,7 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Stuff - nBasSq = nBas*nBas + nBas_AOs_Sq = nBas_AOs*nBas_AOs ! TDA for T @@ -131,9 +135,29 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Memory allocation - allocate(Aph(nS,nS),Bph(nS,nS),eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), & - rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),err(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + allocate(Aph(nS,nS), Bph(nS,nS), Om(nS), XpY(nS,nS), XmY(nS,nS)) + + allocate(eGT(nBas_MOs)) + allocate(eOld(nBas_MOs)) + allocate(Z(nBas_MOs)) + + allocate(c(nBas_AOs,nBas_MOs)) + + allocate(cp(nBas_MOs,nBas_MOs)) + allocate(Fp(nBas_MOs,nBas_MOs)) + allocate(Sig(nBas_MOs,nBas_MOs)) + + allocate(P(nBas_AOs,nBas_AOs)) + allocate(F(nBas_AOs,nBas_AOs)) + allocate(J(nBas_AOs,nBas_AOs)) + allocate(K(nBas_AOs,nBas_AOs)) + allocate(Sigp(nBas_AOs,nBas_AOs)) + allocate(err(nBas_AOs,nBas_AOs)) + + allocate(err_diis(nBas_AOs_Sq,max_diis), F_diis(nBas_AOs_Sq,max_diis)) + + allocate(rhoL(nBas_MOs,nBas_MOs,nS), rhoR(nBas_MOs,nBas_MOs,nS)) + ! Initialization @@ -161,20 +185,20 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Buid Hartree matrix - call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J) + call Hartree_matrix_AO_basis(nBas_AOs,P,ERI_AO,J) ! Compute exchange part of the self-energy - call exchange_matrix_AO_basis(nBas,P,ERI_AO,K) + call exchange_matrix_AO_basis(nBas_AOs,P,ERI_AO,K) ! AO to MO transformation of two-electron integrals - call AOtoMO_ERI_RHF(nBas,nBas,c,ERI_AO,ERI_MO) + call AOtoMO_ERI_RHF(nBas_AOs, nBas_MOs, c, ERI_AO, ERI_MO) ! Compute linear response - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO,Aph) - if(.not.TDA_T) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) + call phLR_A(ispin,dRPA,nBas_MOs,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO,Aph) + if(.not.TDA_T) call phLR_B(ispin,dRPA,nBas_MOs,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY) @@ -182,17 +206,17 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Compute correlation part of the self-energy - call GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,XmY,rhoL,rhoR) + call GTeh_excitation_density(nBas_MOs,nC,nO,nR,nS,ERI_MO,XpY,XmY,rhoL,rhoR) - if(regularize) call GTeh_regularization(nBas,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR) + if(regularize) call GTeh_regularization(nBas_MOs,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR) - call GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR,EcGM,Sig,Z) + call GTeh_self_energy(eta,nBas_MOs,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR,EcGM,Sig,Z) ! Make correlation self-energy Hermitian and transform it back to AO basis Sig = 0.5d0*(Sig + transpose(Sig)) - call MOtoAO(nBas,nBas,S,c,Sig,Sigp) + call MOtoAO(nBas_AOs, nBas_MOs, S, c, Sig, Sigp) ! Solve the quasi-particle equation @@ -207,7 +231,7 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,F) + call DIIS_extrapolation(rcond,nBas_AOs_Sq,nBas_AOs_Sq,n_diis,err_diis,F_diis,err,F) end if @@ -215,9 +239,8 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d Fp = matmul(transpose(X),matmul(F,X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,eGT) + call diagonalize_matrix(nBas_MOs, cp, eGT) c = matmul(X,cp) - Sigp = matmul(transpose(c),matmul(Sigp,c)) ! Compute new density matrix in the AO basis @@ -234,19 +257,19 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Kinetic energy - ET = trace_matrix(nBas,matmul(P,T)) + ET = trace_matrix(nBas_AOs,matmul(P,T)) ! Potential energy - EV = trace_matrix(nBas,matmul(P,V)) + EV = trace_matrix(nBas_AOs,matmul(P,V)) ! Hartree energy - EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) + EJ = 0.5d0*trace_matrix(nBas_AOs,matmul(P,J)) ! Exchange energy - Ex = 0.25d0*trace_matrix(nBas,matmul(P,K)) + Ex = 0.25d0*trace_matrix(nBas_AOs,matmul(P,K)) ! Total energy @@ -254,8 +277,9 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Print results - call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,Sigp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) + call dipole_moment(nBas_AOs,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) + call print_qsRGTeh(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c, Sig, & + Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGT, dipole) end do !------------------------------------------------------------------------ @@ -272,19 +296,21 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) + deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, Om, XpY, XmY, rhoL, rhoR, err, err_diis, F_diis) + stop end if ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,Om,XpY,XmY,rhoL,rhoR,err,err_diis,F_diis) + deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, Om, XpY, XmY, rhoL, rhoR, err, err_diis, F_diis) ! Perform BSE calculation ! if(BSE) then -! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, & +! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas_AOs,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, & ! eGT,eGT,EcBSE) ! if(exchange_kernel) then @@ -319,7 +345,7 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! end if -! call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) +! call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas_AOs,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) ! write(*,*) ! write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GT/qsRGTpp.f90 b/src/GT/qsRGTpp.f90 index 9f4e78b..1f8e17e 100644 --- a/src/GT/qsRGTpp.f90 +++ b/src/GT/qsRGTpp.f90 @@ -1,6 +1,9 @@ -subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA, & - dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & - S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + +! --- + +subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, TDA_T, TDA, & + dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, & + nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! Perform a quasiparticle self-consistent GT calculation @@ -31,25 +34,26 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nC,nO,nV,nR,nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: PHF(nBas,nBas) - double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: T(nBas,nBas) - double precision,intent(in) :: V(nBas,nBas) - double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) - double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) - double precision,intent(in) :: dipole_int_MO(nBas,nBas,ncart) + double precision,intent(in) :: eHF(nBas_MOs) + double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs) + double precision,intent(in) :: PHF(nBas_AOs,nBas_AOs) + double precision,intent(in) :: S(nBas_AOs,nBas_AOs) + double precision,intent(in) :: T(nBas_AOs,nBas_AOs) + double precision,intent(in) :: V(nBas_AOs,nBas_AOs) + double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs) + double precision,intent(in) :: X(nBas_AOs,nBas_MOs) + double precision,intent(in) :: ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs) + double precision,intent(inout):: ERI_MO(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs) + double precision,intent(in) :: dipole_int_AO(nBas_AOs,nBas_AOs,ncart) + double precision,intent(in) :: dipole_int_MO(nBas_MOs,nBas_MOs,ncart) ! Local variables integer :: nSCF - integer :: nBasSq + integer :: nBas_AOs_Sq integer :: ispin integer :: iblock integer :: n_diis @@ -119,7 +123,7 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Stuff - nBasSq = nBas*nBas + nBas_AOs_Sq = nBas_AOs*nBas_AOs ! TDA for T @@ -137,16 +141,30 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Memory allocation - allocate(eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Z(nBas), & - error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + allocate(eGT(nBas_MOs)) + allocate(eOld(nBas_MOs)) + allocate(Z(nBas_MOs)) - allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & - Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & - rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & - Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & - Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & - rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt)) + allocate(c(nBas_AOs,nBas_MOs)) + + allocate(Fp(nBas_MOs,nBas_MOs)) + allocate(cp(nBas_MOs,nBas_MOs)) + allocate(Sig(nBas_MOs,nBas_MOs)) + + allocate(P(nBas_AOs,nBas_AOs)) + allocate(F(nBas_AOs,nBas_AOs)) + allocate(J(nBas_AOs,nBas_AOs)) + allocate(K(nBas_AOs,nBas_AOs)) + allocate(error(nBas_AOs,nBas_AOs)) + allocate(Sigp(nBas_AOs,nBas_AOs)) + + allocate(error_diis(nBas_AOs_Sq,max_diis)) + allocate(F_diis(nBas_AOs_Sq,max_diis)) + + allocate(Om1s(nVVs), X1s(nVVs,nVVs), Y1s(nOOs,nVVs), rho1s(nBas_MOs,nBas_MOs,nVVs)) + allocate(Om2s(nOOs), X2s(nVVs,nOOs), Y2s(nOOs,nOOs), rho2s(nBas_MOs,nBas_MOs,nOOs)) + allocate(Om1t(nVVt), X1t(nVVt,nVVt), Y1t(nOOt,nVVt), rho1t(nBas_MOs,nBas_MOs,nVVt)) + allocate(Om2t(nOOt), X2t(nVVt,nOOt), Y2t(nOOt,nOOt), rho2t(nBas_MOs,nBas_MOs,nOOt)) ! Initialization @@ -174,15 +192,15 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Buid Hartree matrix - call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J) + call Hartree_matrix_AO_basis(nBas_AOs,P,ERI_AO,J) ! Compute exchange part of the self-energy - call exchange_matrix_AO_basis(nBas,P,ERI_AO,K) + call exchange_matrix_AO_basis(nBas_AOs,P,ERI_AO,K) ! AO to MO transformation of two-electron integrals - call AOtoMO_ERI_RHF(nBas,nBas,c,ERI_AO,ERI_MO) + call AOtoMO_ERI_RHF(nBas_AOs, nBas_MOs, c, ERI_AO, ERI_MO) ! Compute linear response @@ -191,9 +209,9 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) + call ppLR_C(iblock,nBas_MOs,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nBas_MOs,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) + if(.not.TDA_T) call ppLR_B(iblock,nBas_MOs,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) @@ -204,9 +222,9 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) + call ppLR_C(iblock,nBas_MOs,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nBas_MOs,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) + if(.not.TDA_T) call ppLR_B(iblock,nBas_MOs,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) @@ -217,25 +235,25 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Compute correlation part of the self-energy - iblock = 3 - call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s) + iblock = 3 + call GTpp_excitation_density(iblock,nBas_MOs,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s) - iblock = 4 - call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t) + iblock = 4 + call GTpp_excitation_density(iblock,nBas_MOs,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t) if(regularize) then - call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s) - call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t) + call GTpp_regularization(eta,nBas_MOs,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s) + call GTpp_regularization(eta,nBas_MOs,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t) end if - call GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, & + call GTpp_self_energy(eta,nBas_MOs,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, & Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z) ! Make correlation self-energy Hermitian and transform it back to AO basis Sig = 0.5d0*(Sig + transpose(Sig)) - call MOtoAO(nBas,nBas,S,c,Sig,Sigp) + call MOtoAO(nBas_AOs, nBas_MOs, S, c, Sig, Sigp) ! Solve the quasi-particle equation @@ -249,22 +267,21 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d n_diis = min(n_diis+1,max_diis) if(abs(rcond) > 1d-7) then - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) + call DIIS_extrapolation(rcond,nBas_AOs_Sq,nBas_AOs_Sq,n_diis,error_diis,F_diis,error,F) else n_diis = 0 end if ! Diagonalize Hamiltonian in AO basis - Fp = matmul(transpose(X),matmul(F,X)) + Fp = matmul(transpose(X), matmul(F, X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,eGT) - c = matmul(X,cp) - Sigp = matmul(transpose(c),matmul(Sigp,c)) + call diagonalize_matrix(nBas_MOs, cp, eGT) + c = matmul(X, cp) ! Compute new density matrix in the AO basis - P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO))) ! Save quasiparticles energy for next cycle @@ -277,19 +294,19 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Kinetic energy - ET = trace_matrix(nBas,matmul(P,T)) + ET = trace_matrix(nBas_AOs,matmul(P,T)) ! Potential energy - EV = trace_matrix(nBas,matmul(P,V)) + EV = trace_matrix(nBas_AOs,matmul(P,V)) ! Hartree energy - EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) + EJ = 0.5d0*trace_matrix(nBas_AOs,matmul(P,J)) ! Exchange energy - Ex = 0.25d0*trace_matrix(nBas,matmul(P,K)) + Ex = 0.25d0*trace_matrix(nBas_AOs,matmul(P,K)) ! Total energy @@ -297,8 +314,10 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Print results - call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,Sigp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) + call dipole_moment(nBas_AOs,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) + call print_qsRGTpp(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, & + eGT, c, Sig, Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, & + EqsGT, dipole) end do !------------------------------------------------------------------------ @@ -315,19 +334,24 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) + deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, error, error_diis, F_diis) + deallocate(Om1s, X1s, Y1s, rho1s) + deallocate(Om2s, X2s, Y2s, rho2s) + deallocate(Om1t, X1t, Y1t, rho1t) + deallocate(Om2t, X2t, Y2t, rho2t) stop end if ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,error,error_diis,F_diis) + deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, error, error_diis, F_diis) ! Perform BSE calculation if(dophBSE) then - call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas_MOs,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & ERI_MO,dipole_int_MO,eGT,eGT,EcBSE) @@ -363,7 +387,7 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d end if - call GTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & + call GTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nBas_MOs,nC,nO,nV,nR,nS, & nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & Om2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcBSE) @@ -391,4 +415,9 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d end if + deallocate(Om1s, X1s, Y1s, rho1s) + deallocate(Om2s, X2s, Y2s, rho2s) + deallocate(Om1t, X1t, Y1t, rho1t) + deallocate(Om2t, X2t, Y2t, rho2t) + end subroutine diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 5796562..db16eab 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -330,11 +330,11 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doGT) then call wall_time(start_GT) - ! TODO - call RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & - maxSCF_GT,thresh_GT,max_diis_GT,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE, & - TDA_T,TDA,dBSE,dTDA,singlet,triplet,lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc, & - nBas_AOs,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + call RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevGTeh, doqsGTeh, & + maxSCF_GT, thresh_GT, max_diis_GT, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, doppBSE, & + TDA_T, TDA, dBSE, dTDA, singlet, triplet, lin_GT, eta_GT, reg_GT, nNuc, ZNuc, rNuc, ENuc, & + nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, & + dipole_int_MO, PHF, cHF, eHF) call wall_time(end_GT) t_GT = end_GT - start_GT