From c06871a0ffb646d5cf5b6aa994c38807db36fc45 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Thu, 29 Aug 2024 00:00:41 +0200 Subject: [PATCH] introduce nBas_MOs in RGW --- src/GW/RGW.f90 | 74 ++++++++++++----------- src/GW/SRG_qsGW.f90 | 124 ++++++++++++++++++++++++--------------- src/GW/print_qsRGW.f90 | 24 ++++---- src/GW/qsRGW.f90 | 130 +++++++++++++++++++++++++---------------- src/QuAcK/RQuAcK.f90 | 9 ++- 5 files changed, 217 insertions(+), 144 deletions(-) diff --git a/src/GW/RGW.f90 b/src/GW/RGW.f90 index 696039d..8498485 100644 --- a/src/GW/RGW.f90 +++ b/src/GW/RGW.f90 @@ -1,7 +1,10 @@ -subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thresh,max_diis,doACFDT, & - exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,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 RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxSCF, thresh, max_diis, doACFDT, & + exchange_kernel, doXBS, dophBSE, dophBSE2, doppBSE, TDA_W, 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) ! Restricted GW module @@ -43,7 +46,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre 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 @@ -51,18 +54,18 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre 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 @@ -76,11 +79,11 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre call wall_time(start_GW) call RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,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_GW) t_GW = end_GW - start_GW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_GW,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for G0W0 = ',t_GW,' seconds' write(*,*) end if @@ -93,11 +96,11 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre call wall_time(start_GW) call evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,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_AOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGW = ',t_GW,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for evGW = ',t_GW,' seconds' write(*,*) end if @@ -109,13 +112,14 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre if(doqsGW) then call wall_time(start_GW) - call qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & - 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 qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, & + TDA_W, TDA, dBSE, dTDA, doppBSE, 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_GW) t_GW = end_GW - start_GW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_GW,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GW,' seconds' write(*,*) end if @@ -127,13 +131,15 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre if(doSRGqsGW) then call wall_time(start_GW) - call SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA, & - singlet,triplet,eta,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 SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, & + dophBSE, dophBSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, & + 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_GW) t_GW = end_GW - start_GW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_GW,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GW,' seconds' write(*,*) end if @@ -145,11 +151,12 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre if(doufG0W0) then call wall_time(start_GW) - call ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + ! TODO + call ufG0W0(dotest,TDA_W,nBas_AOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0W0 = ',t_GW,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ufG0W0 = ',t_GW,' seconds' write(*,*) end if @@ -161,11 +168,12 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre if(doufGW) then call wall_time(start_GW) - call ufGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + ! TODO + call ufGW(dotest,TDA_W,nBas_AOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufGW = ',t_GW,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ufGW = ',t_GW,' seconds' write(*,*) end if diff --git a/src/GW/SRG_qsGW.f90 b/src/GW/SRG_qsGW.f90 index da87210..32c2a19 100644 --- a/src/GW/SRG_qsGW.f90 +++ b/src/GW/SRG_qsGW.f90 @@ -1,6 +1,10 @@ -subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_W,TDA, & - dBSE,dTDA,singlet,triplet,eta,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 SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, & + BSE, BSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, 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 GW calculation @@ -32,30 +36,30 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, 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(inout):: 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(inout):: dipole_int_MO(nBas_MOs,nBas_MOs,ncart) ! Local variables integer :: nSCF - integer :: nBasSq + integer :: nBas_AOs_Sq integer :: ispin integer :: ixyz integer :: n_diis @@ -114,7 +118,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! Stuff - nBasSq = nBas*nBas + nBas_AOs_Sq = nBas_AOs*nBas_AOs ! TDA for W @@ -132,9 +136,32 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! Memory allocation - allocate(eGW(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),Aph(nS,nS),Bph(nS,nS), & - Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS),error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + allocate(eGW(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(SigC(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(SigCp(nBas_AOs,nBas_AOs)) + + allocate(Aph(nS,nS)) + allocate(Bph(nS,nS)) + allocate(Om(nS)) + allocate(XpY(nS,nS)) + allocate(XmY(nS,nS)) + allocate(rho(nBas_MOs,nBas_MOs,nS)) + + allocate(error_diis(nBas_AOs_Sq,max_diis)) + allocate(F_diis(nBas_AOs_Sq,max_diis)) ! Initialization @@ -162,11 +189,11 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! Buid Hartree matrix call wall_time(t1) - 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) call wall_time(t2) tt=tt+t2-t1 @@ -174,11 +201,11 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, call wall_time(tao1) - do ixyz=1,ncart - call AOtoMO(nBas,nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) + do ixyz = 1, ncart + call AOtoMO(nBas_AOs, nBas_MOs, cHF, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) end do - call AOtoMO_ERI_RHF(nBas,nBas,c,ERI_AO,ERI_MO) + call AOtoMO_ERI_RHF(nBas_AOs, nBas_MOs, c, ERI_AO, ERI_MO) call wall_time(tao2) @@ -188,8 +215,8 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, call wall_time(tlr1) - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph) - if(.not.TDA_W) 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,eGW,ERI_MO,Aph) + if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas_MOs,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) @@ -203,13 +230,13 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, call wall_time(tex1) - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho) + call GW_excitation_density(nBas_MOs,nC,nO,nR,nS,ERI_MO,XpY,rho) call wall_time(tex2) tex=tex+tex2-tex1 call wall_time(tsrg1) - call SRG_self_energy(flow,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + call SRG_self_energy(flow,nBas_MOs,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) call wall_time(tsrg2) @@ -218,7 +245,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! Make correlation self-energy Hermitian and transform it back to AO basis call wall_time(tmo1) - call MOtoAO(nBas,nBas,S,c,SigC,SigCp) + call MOtoAO(nBas_AOs, nBas_MOs, S, c, SigC, SigCp) call wall_time(tmo2) tmo = tmo + tmo2 - tmo1 ! Solve the quasi-particle equation @@ -234,18 +261,18 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - 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) 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,eGW) - c = matmul(X,cp) + call diagonalize_matrix(nBas_MOs, cp, eGW) + c = matmul(X, cp) - call AOtoMO(nBas,nBas,c,SigCp,SigC) + call AOtoMO(nBas_AOs, nBas_MOs, c, SigCp, SigC) ! Compute new density matrix in the AO basis @@ -262,19 +289,19 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! 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 @@ -282,8 +309,9 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! Print results - call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole) + call dipole_moment(nBas_AOs,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) + call print_qsRGW(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGW, c, & + SigC, Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGW, dipole) end do !------------------------------------------------------------------------ @@ -300,6 +328,8 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) + deallocate(c, cp, P, F, Fp, J, K, SigC, Z, Om, XpY, XmY, rho, error, error_diis, F_diis) + stop end if @@ -313,17 +343,18 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! Cumulant expansion - call RGWC(dotest,eta,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,eGW,Z) + call RGWC(dotest,eta,nBas_MOs,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,eGW,Z) ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,SigC,Z,Om,XpY,XmY,rho,error,error_diis,F_diis) + deallocate(c, cp, P, F, Fp, J, K, SigC, Z, Om, XpY, XmY, rho, error, error_diis, F_diis) ! Perform BSE calculation if(BSE) then - call GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE) + call GW_phBSE(BSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nBas_MOs, & + nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eGW, eGW, EcBSE) if(exchange_kernel) then @@ -357,7 +388,8 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, end if - call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcBSE) + call GW_phACFDT(exchange_kernel, doXBS, .true., TDA_W, TDA, BSE, singlet, triplet, & + eta, nBas_MOs, nC, nO, nV, nR, nS, ERI_MO, eGW, eGW, EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GW/print_qsRGW.f90 b/src/GW/print_qsRGW.f90 index d09a670..7e90ce7 100644 --- a/src/GW/print_qsRGW.f90 +++ b/src/GW/print_qsRGW.f90 @@ -1,4 +1,8 @@ -subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole) + +! --- + +subroutine print_qsRGW(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGW, c, SigC, & + Z, ENuc, ET, EV, EJ, EK, EcGM, EcRPA, EqsGW, dipole) ! Print useful information about qsRGW calculation @@ -7,7 +11,7 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E ! 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_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E double precision,intent(in) :: EcRPA double precision,intent(in) :: Conv double precision,intent(in) :: thresh - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: eGW(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) :: eGW(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) :: EqsGW double precision,intent(in) :: dipole(ncart) @@ -59,7 +63,7 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E '|','#','|','e_HF (eV)','|','Sig_GW (eV)','|','Z','|','e_GW (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),'|',eGW(p)*HaToeV,'|' end do @@ -110,13 +114,13 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' Restricted qsGW orbital coefficients' write(*,'(A50)') '---------------------------------------' - call matout(nBas,nBas,c) + call matout(nBas_AOs, nBas_MOs, c) write(*,*) end if write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' Restricted qsGW orbital energies (au) ' write(*,'(A50)') '---------------------------------------' - call vecout(nBas,eGW) + call vecout(nBas_MOs, eGW) write(*,*) end if diff --git a/src/GW/qsRGW.f90 b/src/GW/qsRGW.f90 index 96af23b..395c30e 100644 --- a/src/GW/qsRGW.f90 +++ b/src/GW/qsRGW.f90 @@ -1,6 +1,10 @@ -subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & - 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 qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, & + TDA_W, TDA, dBSE, dTDA, doppBSE, 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 GW calculation @@ -34,30 +38,30 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop 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(inout):: 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_AOs) + 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(inout):: dipole_int_MO(nBas_MOs,nBas_MOs,ncart) ! Local variables integer :: nSCF - integer :: nBasSq + integer :: nBas_AOs_Sq integer :: ispin integer :: ixyz integer :: n_diis @@ -112,7 +116,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Stuff - nBasSq = nBas*nBas + nBas_AOs_Sq = nBas_AOs*nBas_AOs ! TDA for W @@ -130,10 +134,31 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Memory allocation - allocate(eGW(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), & - Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), & - err(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + allocate(eGW(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(SigC(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(err(nBas_AOs,nBas_AOs)) + allocate(SigCp(nBas_AOs,nBas_AOs)) + + allocate(Aph(nS,nS)) + allocate(Bph(nS,nS)) + allocate(Om(nS)) + allocate(XpY(nS,nS)) + allocate(XmY(nS,nS)) + allocate(rho(nBas_MOs,nBas_MOs,nS)) + + allocate(err_diis(nBas_AOs_Sq,max_diis)) + allocate(F_diis(nBas_AOs_Sq,max_diis)) ! Initialization @@ -160,38 +185,38 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Build Hartree-exchange matrix - call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J) - call exchange_matrix_AO_basis(nBas,P,ERI_AO,K) + call Hartree_matrix_AO_basis(nBas_AOs, P, ERI_AO, J) + call exchange_matrix_AO_basis(nBas_AOs, P, ERI_AO, K) ! AO to MO transformation of two-electron integrals - do ixyz=1,ncart - call AOtoMO(nBas,nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) + do ixyz = 1, ncart + call AOtoMO(nBas_AOs, nBas_MOs, c, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) end do - 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,eGW,ERI_MO,Aph) - if(.not.TDA_W) 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, eGW, ERI_MO, Aph) + if(.not.TDA_W) call phLR_B(ispin, dRPA, nBas_MOs, nC, nO, nV, nR, nS, 1d0, ERI_MO, Bph) - call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + call phLR(TDA_W, nS, Aph, Bph, EcRPA, Om, XpY, XmY) if(print_W) call print_excitation_energies('phRPA@GW@RHF','singlet',nS,Om) ! Compute correlation part of the self-energy - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho) + call GW_excitation_density(nBas_MOs, nC, nO, nR, nS, ERI_MO, XpY, rho) - if(regularize) call GW_regularization(nBas,nC,nO,nV,nR,nS,eGW,Om,rho) + if(regularize) call GW_regularization(nBas_MOs, nC, nO, nV, nR, nS, eGW, Om, rho) - call GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + call GW_self_energy(eta, nBas_MOs, nC, nO, nV, nR, nS, eGW, Om, rho, EcGM, SigC, Z) ! Make correlation self-energy Hermitian and transform it back to AO basis 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 @@ -205,19 +230,19 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! 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 - EK = 0.25d0*trace_matrix(nBas,matmul(P,K)) + EK = 0.25d0*trace_matrix(nBas_AOs, matmul(P, K)) ! Total energy @@ -228,26 +253,27 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop 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 ! 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,eGW) - c = matmul(X,cp) - call AOtoMO(nBas,nBas,c,SigCp,SigC) + call diagonalize_matrix(nBas_MOs, cp, eGW) + c = matmul(X, cp) + call AOtoMO(nBas_AOs, nBas_MOs, c, SigCp, SigC) ! Density matrix - P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO))) ! Print results - call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigCp,Z,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole) + call dipole_moment(nBas_AOs, P, nNuc, ZNuc, rNuc, dipole_int_AO, dipole) + call print_qsRGW(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGW, c, SigC, Z, & + ENuc, ET, EV, EJ, EK, EcGM, EcRPA, EqsGW, dipole) end do !------------------------------------------------------------------------ @@ -264,19 +290,21 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) + deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, Om, XpY, XmY, rho, err, err_diis, F_diis) stop end if ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,Om,XpY,XmY,rho,err,err_diis,F_diis) + deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, Om, XpY, XmY, rho, err, err_diis, F_diis) ! Perform BSE calculation if(dophBSE) then - call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE) + call GW_phBSE(dophBSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, & + nBas_MOs, nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eGW, eGW, EcBSE) if(exchange_kernel) then @@ -310,7 +338,8 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop end if - call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcBSE) + call GW_phACFDT(exchange_kernel, doXBS, .true., TDA_W, TDA, dophBSE, singlet, triplet, & + eta, nBas_MOs, nC, nO, nV, nR, nS, ERI_MO, eGW, eGW, EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -327,7 +356,8 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop if(doppBSE) then - call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGW,EcBSE) + call GW_ppBSE(TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nBas_MOs, & + nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eHF, eGW, EcBSE) EcBSE(2) = 3d0*EcBSE(2) diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 63a081c..5796562 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -309,11 +309,10 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doGW) then call wall_time(start_GW) - ! 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, & - ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + 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, 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_GW) t_GW = end_GW - start_GW