From 65ba14706c1a62e08d76dd09ff74d5b68051ad12 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 10 Sep 2024 09:34:46 +0200 Subject: [PATCH] cleanup in RG0W0 --- src/GW/RG0W0.f90 | 60 +++++++-------------- src/GW/RGW_phACFDT.f90 | 9 ++++ src/GW/RGW_phBSE.f90 | 29 +++++++++-- src/GW/RGW_ppBSE.f90 | 2 + src/GW/SRG_qsRGW.f90 | 29 +++-------- src/GW/evRGW.f90 | 32 ++---------- src/GW/qsRGW.f90 | 116 ++++++++++++++++------------------------- 7 files changed, 113 insertions(+), 164 deletions(-) diff --git a/src/GW/RG0W0.f90 b/src/GW/RG0W0.f90 index 37162b3..d54815f 100644 --- a/src/GW/RG0W0.f90 +++ b/src/GW/RG0W0.f90 @@ -41,9 +41,11 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Local variables - logical :: print_W = .true. - logical :: dRPA - integer :: ispin + logical :: print_W = .true. + logical :: plot_self = .false. + logical :: dRPA_W + integer :: isp_W + double precision :: lambda double precision :: EcRPA double precision :: EcBSE(nspin) double precision :: EcGM @@ -72,27 +74,18 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Initialization - dRPA = .true. - EcRPA = 0d0 + lambda = 1d0 -! TDA for W +! Spin manifold and TDA for dynamical screening + + isp_W = 1 + dRPA_W = .true. if(TDA_W) then - write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) 'Tamm-Dancoff approximation for dynamical screening!' write(*,*) end if -! TDA - - if(TDA) then - write(*,*) 'Tamm-Dancoff approximation activated!' - write(*,*) - end if - -! Spin manifold - - ispin = 1 - ! Memory allocation allocate(Aph(nS,nS),Bph(nS,nS),SigC(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), & @@ -102,8 +95,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Compute screening ! !-------------------! - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) - if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph) + if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) @@ -149,7 +142,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Plot self-energy, renormalization factor, and spectral function -! call RGW_plot_self_energy(nBas,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho) + if(plot_self) call RGW_plot_self_energy(nBas,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho) !--------------------! ! Cumulant expansion ! @@ -159,8 +152,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Compute the RPA correlation energy - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph) - if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,eGW,ERI,Aph) + if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) @@ -174,14 +167,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA if(dophBSE) then - call RGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) - - if(exchange_kernel) then - - EcBSE(1) = 0.5d0*EcBSE(1) - EcBSE(2) = 1.5d0*EcBSE(2) - - end if + call RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, & + nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -201,13 +188,6 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA write(*,*) '-------------------------------------------------------------' write(*,*) - if(doXBS) then - - write(*,*) '*** scaled screening version (XBS) ***' - write(*,*) - - end if - call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcBSE) write(*,*) @@ -227,8 +207,6 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) - EcBSE(2) = 3d0*EcBSE(2) - write(*,*) write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au' @@ -239,8 +217,6 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA write(*,*) end if - -! end if ! Testing zone diff --git a/src/GW/RGW_phACFDT.f90 b/src/GW/RGW_phACFDT.f90 index 5268c88..703b123 100644 --- a/src/GW/RGW_phACFDT.f90 +++ b/src/GW/RGW_phACFDT.f90 @@ -60,6 +60,15 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas, allocate(Aph(nS,nS),Bph(nS,nS),KA(nS,nS),KB(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & rho_RPA(nBas,nBas,nS),Om(nS),XpY(nS,nS),XmY(nS,nS)) +! eXtended BSE + + if(doXBS) then + + write(*,*) '*** scaled screening version (XBS) ***' + write(*,*) + + end if + ! Antisymmetrized kernel version if(exchange_kernel) then diff --git a/src/GW/RGW_phBSE.f90 b/src/GW/RGW_phBSE.f90 index 870c0d8..5953591 100644 --- a/src/GW/RGW_phBSE.f90 +++ b/src/GW/RGW_phBSE.f90 @@ -1,4 +1,5 @@ -subroutine RGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) +subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, & + nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -8,6 +9,7 @@ subroutine RGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO ! Input variables logical,intent(in) :: dophBSE2 + logical,intent(in) :: exchange_kernel logical,intent(in) :: TDA_W logical,intent(in) :: TDA logical,intent(in) :: dBSE @@ -63,6 +65,10 @@ subroutine RGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO Aph(nS,nS),Bph(nS,nS),KA_sta(nS,nS),KB_sta(nS,nS), & OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS)) +! Initialization + + EcBSE(:) = 0d0 + !--------------------------------- ! Compute (singlet) RPA screening !--------------------------------- @@ -79,6 +85,15 @@ subroutine RGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta) call RGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta) +!-----! +! TDA ! +!-----! + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated in phBSE!' + write(*,*) + end if + !------------------- ! Singlet manifold !------------------- @@ -86,7 +101,6 @@ subroutine RGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO if(singlet) then ispin = 1 - EcBSE(ispin) = 0d0 ! Compute BSE excitation energies @@ -143,7 +157,6 @@ subroutine RGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO if(triplet) then ispin = 2 - EcBSE(ispin) = 0d0 ! Compute BSE excitation energies @@ -168,4 +181,14 @@ subroutine RGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO end if +! Scale properly correlation energy if exchange is included in interaction kernel + + if(exchange_kernel) then + + EcBSE(1) = 0.5d0*EcBSE(1) + EcBSE(2) = 1.5d0*EcBSE(2) + + end if + + end subroutine diff --git a/src/GW/RGW_ppBSE.f90 b/src/GW/RGW_ppBSE.f90 index c4f6353..cfe44ed 100644 --- a/src/GW/RGW_ppBSE.f90 +++ b/src/GW/RGW_ppBSE.f90 @@ -180,6 +180,8 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) + EcBSE(ispin) = 3d0*EcBSE(ispin) + call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) !----------------------------------------------------! diff --git a/src/GW/SRG_qsRGW.f90 b/src/GW/SRG_qsRGW.f90 index 90cff7b..2552dbf 100644 --- a/src/GW/SRG_qsRGW.f90 +++ b/src/GW/SRG_qsRGW.f90 @@ -33,7 +33,8 @@ subroutine SRG_qsRGW(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, nOrb + integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -74,7 +75,7 @@ subroutine SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS double precision,external :: trace_matrix double precision :: dipole(ncart) - logical :: dRPA = .true. + logical :: dRPA_W = .true. logical :: print_W = .true. double precision,allocatable :: error_diis(:,:) double precision,allocatable :: F_diis(:,:) @@ -124,13 +125,6 @@ subroutine SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS write(*,*) end if -! TDA - - if(TDA) then - write(*,*) 'Tamm-Dancoff approximation activated!' - write(*,*) - end if - ! Memory allocation allocate(eGW(nOrb)) @@ -212,8 +206,8 @@ subroutine SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS call wall_time(tlr1) - call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph) - if(.not.TDA_W) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) + call phLR_A(ispin,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph) + if(.not.TDA_W) call phLR_B(ispin,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) @@ -344,21 +338,14 @@ subroutine SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS ! 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 RGW_phBSE(BSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, & - nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eGW, eGW, EcBSE) - - if(exchange_kernel) then - - EcBSE(1) = 0.5d0*EcBSE(1) - EcBSE(2) = 1.5d0*EcBSE(2) - - end if + call RGW_phBSE(BSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,& + nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GW/evRGW.f90 b/src/GW/evRGW.f90 index 6de0efc..654f17f 100644 --- a/src/GW/evRGW.f90 +++ b/src/GW/evRGW.f90 @@ -82,13 +82,6 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop write(*,*) end if -! TDA - - if(TDA) then - write(*,*) 'Tamm-Dancoff approximation activated!' - write(*,*) - end if - ! Linear mixing linear_mixing = .false. @@ -221,14 +214,8 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop if(dophBSE) then - call RGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGW,eGW,EcBSE) - - if(exchange_kernel) then - - EcBSE(1) = 0.5d0*EcBSE(1) - EcBSE(2) = 1.5d0*EcBSE(2) - - end if + call RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, & + nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGW,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -243,18 +230,11 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop if(doACFDT) then - write(*,*) '------------------------------------------------------' - write(*,*) 'Adiabatic connection version of BSE correlation energy' - write(*,*) '------------------------------------------------------' + write(*,*) '-----------------------------------------------------------' + write(*,*) 'Adiabatic connection version of BSE@evGW correlation energy' + write(*,*) '-----------------------------------------------------------' write(*,*) - if(doXBS) then - - write(*,*) '*** scaled screening version (XBS) ***' - write(*,*) - - end if - call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcBSE) write(*,*) @@ -274,8 +254,6 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) - EcBSE(2) = 3d0*EcBSE(2) - write(*,*) write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGW@RHF correlation energy (singlet) = ',EcBSE(1),' au' diff --git a/src/GW/qsRGW.f90 b/src/GW/qsRGW.f90 index 4990da1..e82caf8 100644 --- a/src/GW/qsRGW.f90 +++ b/src/GW/qsRGW.f90 @@ -1,10 +1,7 @@ - -! --- - -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, 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) +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,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) ! Perform a quasiparticle self-consistent GW calculation @@ -38,7 +35,8 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas, nOrb + integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -78,7 +76,7 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX double precision,external :: trace_matrix double precision :: dipole(ncart) - logical :: dRPA = .true. + logical :: dRPA_W = .true. logical :: print_W = .false. double precision,allocatable :: err_diis(:,:) double precision,allocatable :: F_diis(:,:) @@ -125,13 +123,6 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX write(*,*) end if -! TDA - - if(TDA) then - write(*,*) 'Tamm-Dancoff approximation activated!' - write(*,*) - end if - ! Memory allocation allocate(eGW(nOrb)) @@ -185,68 +176,68 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX ! 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,P,ERI_AO,J) + call exchange_matrix_AO_basis(nBas,P,ERI_AO,K) ! AO to MO transformation of two-electron integrals - do ixyz = 1, ncart - call AOtoMO(nBas, nOrb, c, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) + do ixyz=1,ncart + call AOtoMO(nBas,nOrb,c,dipole_int_AO(1,1,ixyz),dipole_int_MO(1,1,ixyz)) end do - call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) + call AOtoMO_ERI_RHF(nBas,nOrb,c,ERI_AO,ERI_MO) ! Compute linear response - call phLR_A(ispin, dRPA, nOrb, nC, nO, nV, nR, nS, 1d0, eGW, ERI_MO, Aph) - if(.not.TDA_W) call phLR_B(ispin, dRPA, nOrb, nC, nO, nV, nR, nS, 1d0, ERI_MO, Bph) + call phLR_A(ispin,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph) + if(.not.TDA_W) call phLR_B(ispin,dRPA_W,nOrb,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 RGW_excitation_density(nOrb, nC, nO, nR, nS, ERI_MO, XpY, rho) + call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,rho) - if(regularize) call GW_regularization(nOrb, nC, nO, nV, nR, nS, eGW, Om, rho) + if(regularize) call GW_regularization(nOrb,nC,nO,nV,nR,nS,eGW,Om,rho) - call RGW_self_energy(eta, nOrb, nC, nO, nV, nR, nS, eGW, Om, rho, EcGM, SigC, Z) + call RGW_self_energy(eta,nOrb,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, nOrb, S, c, SigC, SigCp) + call MOtoAO(nBas,nOrb,S,c,SigC,SigCp) ! Solve the quasi-particle equation F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:) if(nBas .ne. nOrb) then - call AOtoMO(nBas, nOrb, c(1,1), F(1,1), Fp(1,1)) - call MOtoAO(nBas, nOrb, S(1,1), c(1,1), Fp(1,1), F(1,1)) + call AOtoMO(nBas,nOrb,c(1,1),F(1,1),Fp(1,1)) + call MOtoAO(nBas,nOrb,S(1,1),c(1,1),Fp(1,1),F(1,1)) endif ! Compute commutator and convergence criteria - err = matmul(F, matmul(P, S)) - matmul(matmul(S, P), F) + err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) if(nSCF > 1) Conv = maxval(abs(err)) ! Kinetic energy - ET = trace_matrix(nBas, matmul(P, T)) + ET = trace_matrix(nBas,matmul(P,T)) ! Potential energy - EV = trace_matrix(nBas, matmul(P, V)) + EV = trace_matrix(nBas,matmul(P,V)) ! Hartree energy - EJ = 0.5d0*trace_matrix(nBas, matmul(P, J)) + EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) ! Exchange energy - EK = 0.25d0*trace_matrix(nBas, matmul(P, K)) + EK = 0.25d0*trace_matrix(nBas,matmul(P,K)) ! Total energy @@ -264,29 +255,29 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX ! Diagonalize Hamiltonian in AO basis if(nBas .eq. nOrb) then - Fp = matmul(transpose(X), matmul(F, X)) + Fp = matmul(transpose(X),matmul(F,X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nOrb, cp, eGW) - c = matmul(X, cp) + call diagonalize_matrix(nOrb,cp,eGW) + c = matmul(X,cp) else - Fp = matmul(transpose(c), matmul(F, c)) + Fp = matmul(transpose(c),matmul(F,c)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nOrb, cp, eGW) - c = matmul(c, cp) + call diagonalize_matrix(nOrb,cp,eGW) + c = matmul(c,cp) endif - call AOtoMO(nBas, nOrb, c, SigCp, SigC) + call AOtoMO(nBas,nOrb,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, nOrb, nO, nSCF, Conv, thresh, eHF, eGW, c, SigC, Z, & - ENuc, ET, EV, EJ, EK, EcGM, EcRPA, EqsGW, dipole) + call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) + call print_qsRGW(nBas,nOrb,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z, & + ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole) end do !------------------------------------------------------------------------ @@ -303,28 +294,21 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - 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) 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 RGW_phBSE(dophBSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, & - nOrb, nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eGW, eGW, EcBSE) - - if(exchange_kernel) then - - EcBSE(1) = 0.5d0*EcBSE(1) - EcBSE(2) = 1.5d0*EcBSE(2) - - end if + call RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, & + nOrb,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -339,18 +323,11 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX if(doACFDT) then - write(*,*) '------------------------------------------------------' - write(*,*) 'Adiabatic connection version of BSE correlation energy' - write(*,*) '------------------------------------------------------' + write(*,*) '-----------------------------------------------------------' + write(*,*) 'Adiabatic connection version of BSE@qsGW correlation energy' + write(*,*) '-----------------------------------------------------------' write(*,*) - if(doXBS) then - - write(*,*) '*** scaled screening version (XBS) ***' - write(*,*) - - end if - call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcBSE) write(*,*) @@ -368,10 +345,7 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX if(doppBSE) then - call RGW_ppBSE(TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, & - nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eHF, eGW, EcBSE) - - EcBSE(2) = 3d0*EcBSE(2) + call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------'