diff --git a/src/GF/RGF.f90 b/src/GF/RGF.f90 index 3493ff4..ea52e82 100644 --- a/src/GF/RGF.f90 +++ b/src/GF/RGF.f90 @@ -1,7 +1,8 @@ -subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,maxSCF, & +subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2, & + renorm,maxSCF, & thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize, & eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, & - S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + S,X,T,V,Hc,ERI_AO,ERI_MO,CAP,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) ! Green's function module @@ -13,6 +14,7 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max logical,intent(in) :: dotest logical,intent(in) :: doG0F2 + logical,intent(in) :: docG0F2 logical,intent(in) :: doevGF2 logical,intent(in) :: doqsGF2 logical,intent(in) :: doufG0F02 @@ -52,6 +54,7 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max double precision,intent(in) :: cHF(nBas,nOrb) double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: CAP(nBas,nBas) double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) @@ -168,4 +171,21 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max end if +!------------------------------------------------------------------------ +! Compute complex G0F2 electronic binding energies +!------------------------------------------------------------------------ + + if(doG0F2) then + + call wall_time(start_GF) + call cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet, & + linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS, & + ENuc,ERHF,ERI_MO,CAP,dipole_int_MO,eHF) + call wall_time(end_GF) + + t_GF = end_GF - start_GF + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds' + write(*,*) + + end if end subroutine diff --git a/src/GW/cRG0W0.f90 b/src/GW/cRG0W0.f90 index b01ed30..74f84c7 100644 --- a/src/GW/cRG0W0.f90 +++ b/src/GW/cRG0W0.f90 @@ -135,10 +135,7 @@ subroutine cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TD ! Linearized or graphical solution? Re_eGWlin(:) = eHF(:) + Re_Z(:)*Re_SigC(:) - Im_Z(:)*Im_SigC(:) - Im_eGWlin(:) = Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:) - do p = 1, nOrb - Im_eGWlin(p) = Im_eGWlin(p) + CAP(p,p) - end do + Im_eGWlin(:) = e_cap(:) + Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:) if(linearize) then diff --git a/src/GW/cRGW_self_energy_diag.f90 b/src/GW/cRGW_self_energy_diag.f90 index eb0717c..53c8f34 100644 --- a/src/GW/cRGW_self_energy_diag.f90 +++ b/src/GW/cRGW_self_energy_diag.f90 @@ -25,8 +25,8 @@ subroutine cRGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Re_S integer :: i,a,p,m double precision :: num,eps double precision :: eta_tilde - double precision :: Re_DS(nBas) - double precision :: Im_DS(nBas) + double precision,allocatable :: Re_DS(:) + double precision,allocatable :: Im_DS(:) ! Output variables @@ -37,6 +37,7 @@ subroutine cRGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Re_S double precision,intent(out) :: EcGM ! Initialize + allocate(Re_DS(nBas),Im_DS(nBas)) Re_Sig(:) = 0d0 Im_Sig(:) = 0d0 Re_DS(:) = 0d0 @@ -97,5 +98,5 @@ subroutine cRGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Re_S ! Compute renormalization factor from derivative Re_Z(:) = (1d0-Re_DS(:))/((1d0 - Re_DS(:))**2 + Im_DS(:)**2) Im_Z(:) = Im_DS(:)/((1d0 - Re_DS(:))**2 + Im_DS(:)**2) - + deallocate(Re_DS,Im_DS) end subroutine diff --git a/src/HF/cRHF.f90 b/src/HF/cRHF.f90 index 04f84e8..7775731 100644 --- a/src/HF/cRHF.f90 +++ b/src/HF/cRHF.f90 @@ -146,7 +146,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,ENuc, & ! Total energy - ERHF = ET + EV+ EW + EJ + EK + ERHF = ET + EV + EW + EJ + EK ! DIIS extrapolation ! @@ -154,7 +154,6 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,ENuc, & n_diis = min(n_diis+1,max_diis) call complex_DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,F) - if (rcond<1.0d-10) write(*,*) "!!! DIIS system ill conditioned: rcond = ", rcond," !!!" end if ! Level shift if(level_shift > 0d0 .and. Conv > thresh) call complex_level_shifting(level_shift,nBas,nBas,nO,S,c,F) diff --git a/src/HF/print_cRHF.f90 b/src/HF/print_cRHF.f90 index 17bfab7..2295a7f 100644 --- a/src/HF/print_cRHF.f90 +++ b/src/HF/print_cRHF.f90 @@ -46,7 +46,7 @@ subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EW,EJ, EK, ERHF) write(*,'(A69)') '---------------------------------------------------------' write(*,'(A33)') ' Summary ' write(*,'(A69)') '---------------------------------------------------------' - write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' One-electron energy = ',real(ET + EV),'+',aimag(ET+EV),'i',' au' + write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' One-electron energy = ',real(ET + EV +EW),'+',aimag(ET+EV+EW),'i',' au' write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' Kinetic energy = ',real(ET),'+',aimag(ET),'i',' au' write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' Potential energy = ',real(EV),'+',aimag(EV),'i',' au' write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' CAP energy = ',real(EW),'+',aimag(EW),'i',' au' @@ -81,5 +81,10 @@ subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EW,EJ, EK, ERHF) write(*,'(A37)') '-------------------------------' call complex_vecout(nOrb, eHF) write(*,*) + write(*,'(A37)') '-------------------------------' + write(*,'(A37)') ' cRHF orbital energies (eV) ' + write(*,'(A37)') '-------------------------------' + call complex_vecout(nOrb, eHF*HaToeV) + write(*,*) end subroutine diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index fb27d49..f11986b 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -13,7 +13,7 @@ program QuAcK logical :: dophRPA,dophRPAx,docrRPA,doppRPA logical :: doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3 logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW - logical :: docG0W0 + logical :: docG0W0,docG0F2 logical :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh integer :: nNuc @@ -127,7 +127,7 @@ program QuAcK doG0W0,doevGW,doqsGW,doufG0W0,doufGW, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp, & doG0T0eh,doevGTeh,doqsGTeh, & - docG0W0, & + docG0W0,docG0F2, & doRtest,doUtest,doGtest) !--------------------------! @@ -258,7 +258,7 @@ program QuAcK dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & - docG0W0, & + docG0W0,docG0F2, & nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & S,T,V,Hc,CAP,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, & guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, & diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index ec63f8c..adad853 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -3,7 +3,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF, dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & - docG0W0, & + docG0W0,docG0F2, & nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & S,T,V,Hc,CAP_AO,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, & guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,singlet,triplet,TDA, & @@ -34,7 +34,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF, logical,intent(in) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW logical,intent(in) :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp logical,intent(in) :: doG0T0eh,doevGTeh,doqsGTeh - logical,intent(in) :: docG0W0 + logical,intent(in) :: docG0W0,docG0F2 integer,intent(in) :: nNuc,nBas,nOrb integer,intent(in) :: nC @@ -207,7 +207,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF, ! 4-index transform call AOtoMO_ERI_RHF(nBas,nOrb,cHF,ERI_AO,ERI_MO) - if (docG0W0) call AOtoMO(nBas,nOrb,cHF,CAP_AO,CAP_MO) ! Add the different cases for cGW methods when implemented + if (docG0W0 .or. docG0F2) call AOtoMO(nBas,nOrb,cHF,CAP_AO,CAP_MO) ! Add the different cases for cGW methods when implemented call wall_time(end_AOtoMO) t_AOtoMO = end_AOtoMO - start_AOtoMO @@ -327,15 +327,15 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF, ! Green's function module ! !-------------------------! - doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3 +doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3 .or. docG0F2 if(doGF) then call wall_time(start_GF) - call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm_GF,maxSCF_GF, & + call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2,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,nOrb,nC,nO,nV,nR,nS,ERHF, & - S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + S,X,T,V,Hc,ERI_AO,ERI_MO,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) call wall_time(end_GF) t_GF = end_GF - start_GF diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index 1e20afc..e3eef02 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -1,5 +1,5 @@ subroutine read_methods(working_dir, & - doRHF,doUHF,doGHF,doROHF,doHFB,docRHF, & + doRHF,doUHF,doGHF,doROHF,doHFB,docRHF, & doMP2,doMP3, & doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & do_drCCD,do_rCCD,do_crCCD,do_lCCD, & @@ -10,7 +10,7 @@ subroutine read_methods(working_dir, & doG0W0,doevGW,doqsGW,doufG0W0,doufGW, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp, & doG0T0eh,doevGTeh,doqsGTeh, & - docG0W0, & + docG0W0,docG0F2, & doRtest,doUtest,doGtest) ! Read desired methods @@ -33,7 +33,7 @@ subroutine read_methods(working_dir, & logical,intent(out) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW logical,intent(out) :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp logical,intent(out) :: doG0T0eh,doevGTeh,doqsGTeh - logical,intent(out) :: docG0W0 + logical,intent(out) :: docG0W0,docG0F2 logical,intent(out) :: doRtest,doUtest,doGtest @@ -205,12 +205,15 @@ subroutine read_methods(working_dir, & if(ans2 == 'T') doevGTeh = .true. if(ans3 == 'T') doqsGTeh = .true. - ! Read CAP-GW methods + ! Read Complex methods docG0W0 = .false. + doG0F2 = .false. + read(1,*) - read(1,*) ans1 + read(1,*) ans1,ans2 if(ans1 == 'T') docG0W0 = .true. + if(ans2 == 'T') docG0F2 = .true. ! Read test