diff --git a/src/GF/RG0F2.f90 b/src/GF/RG0F2.f90 index 7cb3c40..53aeb1e 100644 --- a/src/GF/RG0F2.f90 +++ b/src/GF/RG0F2.f90 @@ -51,7 +51,7 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize, ! Memory allocation - allocate(SigC(nBas),Z(nBas),eGFlin(nBas),eGF(nBas)) + allocate(SigC(nBas), Z(nBas), eGFlin(nBas), eGF(nBas)) ! Frequency-dependent second-order contribution @@ -133,4 +133,6 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize, end if + deallocate(SigC, Z, eGFlin, eGF) + end subroutine diff --git a/src/GF/RGF.f90 b/src/GF/RGF.f90 index 2fe50f2..72911e4 100644 --- a/src/GF/RGF.f90 +++ b/src/GF/RGF.f90 @@ -1,7 +1,10 @@ -subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,maxSCF,thresh,max_diis, & - dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,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) + +! --- + +subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, renorm, maxSCF, & + thresh, max_diis, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, linearize, & + eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, EHF, & + S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, epsHF) ! Green's function module @@ -39,7 +42,7 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max 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 @@ -47,18 +50,18 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max integer,intent(in) :: nS double precision,intent(in) :: EHF - double precision,intent(in) :: epsHF(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(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) :: epsHF(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 @@ -71,12 +74,13 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max if(doG0F2) then call wall_time(start_GF) - call RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & - nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) + call RG0F2(dotest, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, & + linearize, eta, regularize, nBas_MOs, nC, nO, nV, nR, nS, & + ENuc, EHF, ERI_MO, dipole_int_MO, epsHF) 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' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds' write(*,*) end if @@ -89,12 +93,12 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max call wall_time(start_GF) call evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, & - singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF, & - ERI,dipole_int,epsHF) + singlet,triplet,linearize,eta,regularize,nBas_MOs,nC,nO,nV,nR,nS,ENuc,EHF, & + ERI_MO,dipole_int_MO,epsHF) 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' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds' write(*,*) end if @@ -106,12 +110,14 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max if(doqsGF2) then call wall_time(start_GF) - call qsRGF2(dotest,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) + call qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & + dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, & + rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, EHF, S, & + X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, epsHF) 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' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGF2 = ',t_GF,' seconds' write(*,*) end if @@ -123,11 +129,11 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max if(doufG0F02) then call wall_time(start_GF) - call ufRG0F02(dotest,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,epsHF) + call ufRG0F02(dotest, nBas_MOs, nC, nO, nV, nR, nS, ENuc, EHF, ERI_MO, epsHF) call wall_time(end_GF) t_GF = end_GF - start_GF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0F02 = ',t_GF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ufG0F02 = ',t_GF,' seconds' write(*,*) end if @@ -139,11 +145,11 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max if(doG0F3) then call wall_time(start_GF) - call RG0F3(dotest,renorm,nBas,nC,nO,nV,nR,ERI,epsHF) + call RG0F3(dotest, renorm, nBas_MOs, nC, nO, nV, nR, ERI_MO, epsHF) 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' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF3 = ',t_GF,' seconds' write(*,*) end if @@ -155,11 +161,11 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max if(doevGF3) then call wall_time(start_GF) - call evRGF3(dotest,maxSCF,thresh,max_diis,renorm,nBas,nC,nO,nV,nR,ERI,epsHF) + call evRGF3(dotest, maxSCF, thresh, max_diis, renorm, nBas_MOs, nC, nO, nV, nR, ERI_MO, epsHF) 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' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF3 = ',t_GF,' seconds' write(*,*) end if diff --git a/src/GF/evRGF2.f90 b/src/GF/evRGF2.f90 index 4dc2610..c3c4637 100644 --- a/src/GF/evRGF2.f90 +++ b/src/GF/evRGF2.f90 @@ -62,7 +62,7 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si ! Memory allocation - allocate(SigC(nBas),Z(nBas),eGF(nBas),eOld(nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis)) + allocate(SigC(nBas), Z(nBas), eGF(nBas), eOld(nBas), error_diis(nBas,max_diis), e_diis(nBas,max_diis)) ! Initialization @@ -189,4 +189,6 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si end if + deallocate(SigC, Z, eGF, eOld, error_diis, e_diis) + end subroutine diff --git a/src/GF/print_qsRGF2.f90 b/src/GF/print_qsRGF2.f90 index 845c5a0..42132c0 100644 --- a/src/GF/print_qsRGF2.f90 +++ b/src/GF/print_qsRGF2.f90 @@ -1,4 +1,8 @@ -subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ,Ex,Ec,EqsGF,dipole) + +! --- + +subroutine print_qsRGF2(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGF, c, & + SigC, Z, ENuc, ET, EV, EJ, Ex, Ec, EqsGF, dipole) ! Print one-electron energies and other stuff for qsGF2 @@ -7,17 +11,17 @@ subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,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 double precision,intent(in) :: Conv double precision,intent(in) :: thresh - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: eGF(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) :: eGF(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) :: ET double precision,intent(in) :: EV double precision,intent(in) :: EJ @@ -53,7 +57,7 @@ subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ, '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' write(*,*)'-------------------------------------------------------------------------------' - do q=1,nBas + do q = 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)') & '|',q,'|',eHF(q)*HaToeV,'|',SigC(q,q)*HaToeV,'|',Z(q),'|',eGF(q)*HaToeV,'|' end do @@ -102,12 +106,12 @@ subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ, write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGF2 MO coefficients' write(*,'(A50)') '---------------------------------------' - call matout(nBas,nBas,c) + call matout(nBas_AOs, nBas_MOs, c) write(*,*) write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGF2 MO energies' write(*,'(A50)') '---------------------------------------' - call matout(nBas,1,eGF) + call matout(nBas_MOs, 1, eGF) write(*,*) end if diff --git a/src/GF/qsRGF2.f90 b/src/GF/qsRGF2.f90 index 9a2d71b..4287676 100644 --- a/src/GF/qsRGF2.f90 +++ b/src/GF/qsRGF2.f90 @@ -1,6 +1,10 @@ -subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,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 qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, 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 GF2 calculation @@ -29,25 +33,25 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si 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,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 :: n_diis double precision :: EqsGF2 @@ -94,7 +98,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si ! Stuff - nBasSq = nBas*nBas + nBas_AOs_Sq = nBas_AOs*nBas_AOs ! TDA @@ -105,9 +109,27 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si ! Memory allocation - allocate(eGF(nBas),eOld(nbas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),Z(nBas), & - error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + allocate(eGF(nBas_MOs)) + allocate(eOld(nBas_MOs)) + + allocate(c(nBas_AOs,nBas_MOs)) + + allocate(cp(nBas_MOs,nBas_MOs)) + allocate(Fp(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(Z(nBas_MOs)) + allocate(SigC(nBas_MOs,nBas_MOs)) + + allocate(SigCp(nBas_AOs,nBas_AOs)) + + allocate(error_diis(nBas_AOs_Sq,max_diis)) + allocate(F_diis(nBas_AOs_Sq,max_diis)) ! Initialization @@ -117,7 +139,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si Conv = 1d0 P(:,:) = PHF(:,:) eOld(:) = eHF(:) - eGF(:) = eHF(:) + eGF(:) = eHF(:) c(:,:) = cHF(:,:) F_diis(:,:) = 0d0 error_diis(:,:) = 0d0 @@ -135,25 +157,25 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si ! 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 self-energy and renormalization factor if(regularize) then - call GF2_reg_self_energy(eta,nBas,nC,nO,nV,nR,eGF,ERI_MO,SigC,Z) + call GF2_reg_self_energy(eta, nBas_MOs, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z) else - call GF2_self_energy(eta,nBas,nC,nO,nV,nR,eGF,ERI_MO,SigC,Z) + call GF2_self_energy(eta, nBas_MOs, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z) end if @@ -161,7 +183,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si SigC = 0.5d0*(SigC + transpose(SigC)) - call MOtoAO(nBas,nBas,S,c,SigC,SigCp) + call MOtoAO(nBas_AOs, nBas_MOs, S, c, SigC, SigCp) ! Solve the quasi-particle equation @@ -169,28 +191,27 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si ! Compute commutator and convergence criteria - error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) + error = matmul(F, matmul(P, S)) - matmul(matmul(S, P), F) ! DIIS extrapolation - n_diis = min(n_diis+1,max_diis) + 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,eGF) - c = matmul(X,cp) - SigCp = matmul(transpose(c),matmul(SigCp,c)) + call diagonalize_matrix(nBas_MOs, cp, eGF) + 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 @@ -203,23 +224,23 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si ! 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)) ! Correlation energy - call RMP2(.false.,regularize,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EqsGF2,eGF,Ec) + call RMP2(.false., regularize, nBas_MOs, nC, nO, nV, nR, ERI_MO, ENuc, EqsGF2, eGF, Ec) ! Total energy @@ -230,8 +251,9 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si ! Print results !------------------------------------------------------------------------ - call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigCp,Z,ENuc,ET,EV,EJ,Ex,Ec,EqsGF2,dipole) + call dipole_moment(nBas_AOs, P, nNuc, ZNuc, rNuc, dipole_int_AO, dipole) + call print_qsRGF2(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGF, & + c, SigC, Z, ENuc, ET, EV, EJ, Ex, Ec, EqsGF2, dipole) end do !------------------------------------------------------------------------ @@ -248,19 +270,21 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) + deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, error, error_diis, F_diis) stop end if ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,error,error_diis,F_diis) + deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, error, error_diis, F_diis) ! Perform BSE calculation if(dophBSE) then - call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGF,EcBSE) + call GF2_phBSE2(TDA, dBSE, dTDA, singlet, triplet, eta, nBas_MOs, nC, nO, & + nV, nR, nS, ERI_MO, dipole_int_MO, eGF, EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -278,7 +302,8 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si if(doppBSE) then - call GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI_MO,dipole_int_MO,eGF,EcBSE) + call GF2_ppBSE2(TDA, dBSE, dTDA, singlet, triplet, eta, nBas_MOs, & + nC, nO, nV, nR, ERI_MO, dipole_int_MO, eGF, EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index c9c3888..63a081c 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -269,9 +269,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doRPA) then call wall_time(start_RPA) - ! TODO - call RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, & - nBas_AOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF,cHF,S) + call RRPA(dotest, dophRPA, dophRPAx, docrRPA, doppRPA, TDA, doACFDT, exchange_kernel, singlet, triplet, & + nBas_MOs, nC, nO, nV, nR, nS, ENuc, ERHF, ERI_MO, dipole_int_MO, eHF) call wall_time(end_RPA) t_RPA = end_RPA - start_RPA @@ -289,11 +288,10 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doGF) then call wall_time(start_GF) - ! TODO - call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, & - dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF,eta_GF,reg_GF, & - 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 RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, renorm_GF, maxSCF_GF, & + thresh_GF, max_diis_GF, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, lin_GF, & + eta_GF, reg_GF, 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_GF) t_GF = end_GF - start_GF @@ -314,7 +312,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d ! TODO call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW, & doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, & - lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas_AOs,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, & + lin_GW,eta_GW,reg_GW,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 wall_time(end_GW) @@ -334,9 +332,9 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d 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, & + 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 wall_time(end_GT) diff --git a/src/RPA/RRPA.f90 b/src/RPA/RRPA.f90 index 93bf457..ab38932 100644 --- a/src/RPA/RRPA.f90 +++ b/src/RPA/RRPA.f90 @@ -1,5 +1,5 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF,cHF,S) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Random-phase approximation module @@ -29,8 +29,6 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) @@ -49,7 +47,7 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker 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' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RPA = ',t_RPA,' seconds' write(*,*) end if @@ -65,7 +63,7 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker 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' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RPAx = ',t_RPA,' seconds' write(*,*) end if @@ -81,7 +79,7 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker 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' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for pp-RPA = ',t_RPA,' seconds' write(*,*) end if @@ -97,7 +95,7 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker 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' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for pp-RPA = ',t_RPA,' seconds' write(*,*) end if