From 65ba14706c1a62e08d76dd09ff74d5b68051ad12 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 10 Sep 2024 09:34:46 +0200 Subject: [PATCH 1/2] 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(*,*)'-------------------------------------------------------------------------------' From 826d80d5941b0b24b1776685d58a317a0c05a4e5 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 10 Sep 2024 09:40:22 +0200 Subject: [PATCH 2/2] more cleanup --- src/GW/RG0W0.f90 | 41 +++++++++--------- src/GW/RGW.f90 | 8 ++-- src/GW/evRGW.f90 | 39 ++++++++--------- src/GW/ufBSE.f90 | 27 ++++++------ src/GW/ufG0W0.f90 | 45 ++++++++++---------- src/GW/ufGW.f90 | 105 +++++++++++++++++++++++----------------------- 6 files changed, 135 insertions(+), 130 deletions(-) diff --git a/src/GW/RG0W0.f90 b/src/GW/RG0W0.f90 index d54815f..34533b2 100644 --- a/src/GW/RG0W0.f90 +++ b/src/GW/RG0W0.f90 @@ -1,5 +1,5 @@ subroutine 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,dipole_int,eHF) + linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform G0W0 calculation @@ -28,6 +28,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA logical,intent(in) :: regularize integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -35,9 +36,9 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA integer,intent(in) :: nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) - double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) + double precision,intent(in) :: eHF(nOrb) ! Local variables @@ -88,15 +89,15 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! 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), & - eGW(nBas),eGWlin(nBas)) + allocate(Aph(nS,nS),Bph(nS,nS),SigC(nOrb),Z(nOrb),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS), & + eGW(nOrb),eGWlin(nOrb)) !-------------------! ! Compute screening ! !-------------------! - 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_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph) + if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,lambda,ERI,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) @@ -106,15 +107,15 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Compute spectral weights ! !--------------------------! - call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho) + call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho) !------------------------! ! Compute GW self-energy ! !------------------------! - if(regularize) call GW_regularization(nBas,nC,nO,nV,nR,nS,eHF,Om,rho) + if(regularize) call GW_regularization(nOrb,nC,nO,nV,nR,nS,eHF,Om,rho) - call RGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z) + call RGW_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z) !-----------------------------------! ! Solve the quasi-particle equation ! @@ -136,24 +137,24 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA write(*,*) ' *** Quasiparticle energies obtained by root search *** ' write(*,*) - call RGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eHF,eGW,Z) + call RGW_QP_graph(eta,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eHF,eGW,Z) end if ! Plot self-energy, renormalization factor, and spectral function - if(plot_self) call RGW_plot_self_energy(nBas,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho) + if(plot_self) call RGW_plot_self_energy(nOrb,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho) !--------------------! ! Cumulant expansion ! !--------------------! -! call RGWC(dotest,eta,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eHF,eGW,Z) +! call RGWC(dotest,eta,nOrb,nC,nO,nV,nR,nS,Om,rho,eHF,eHF,eGW,Z) ! Compute the RPA correlation energy - 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_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,lambda,eGW,ERI,Aph) + if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,lambda,ERI,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) @@ -161,14 +162,14 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Dump results ! !--------------! - call print_RG0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) + call print_RG0W0(nOrb,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) ! Perform BSE calculation if(dophBSE) then 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) + nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -188,7 +189,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA write(*,*) '-------------------------------------------------------------' write(*,*) - call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcBSE) + call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -205,7 +206,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA if(doppBSE) then - call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) + call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GW/RGW.f90 b/src/GW/RGW.f90 index 6970d9c..897969b 100644 --- a/src/GW/RGW.f90 +++ b/src/GW/RGW.f90 @@ -77,7 +77,7 @@ 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,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) + linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -94,7 +94,7 @@ 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_MO,eHF) + singlet,triplet,linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -150,7 +150,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre call wall_time(start_GW) ! TODO - call ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + call ufG0W0(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -167,7 +167,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre call wall_time(start_GW) ! TODO - call ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + call ufRGW(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW diff --git a/src/GW/evRGW.f90 b/src/GW/evRGW.f90 index 654f17f..f70d896 100644 --- a/src/GW/evRGW.f90 +++ b/src/GW/evRGW.f90 @@ -1,5 +1,5 @@ subroutine 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,dipole_int,eHF) + singlet,triplet,linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform self-consistent eigenvalue-only GW calculation @@ -32,14 +32,15 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop logical,intent(in) :: regularize integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) ! Local variables @@ -89,8 +90,8 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Memory allocation - allocate(Aph(nS,nS),Bph(nS,nS),eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas), & - Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis)) + allocate(Aph(nS,nS),Bph(nS,nS),eGW(nOrb),eOld(nOrb),Z(nOrb),SigC(nOrb), & + Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS),error_diis(nOrb,max_diis),e_diis(nOrb,max_diis)) ! Initialization @@ -113,20 +114,20 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Compute screening - 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(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph) + if(.not.TDA_W) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) ! Compute spectral weights - call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho) + call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho) ! Compute correlation part of the self-energy - if(regularize) call GW_regularization(nBas,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_diag(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + call RGW_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) ! Solve the quasi-particle equation @@ -142,7 +143,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop write(*,*) ' *** Quasiparticle energies obtained by root search *** ' write(*,*) - call RGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eOld,eOld,eGW,Z) + call RGW_QP_graph(eta,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,eOld,eOld,eGW,Z) end if @@ -152,7 +153,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Print results - call print_evRGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) + call print_evRGW(nOrb,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) ! Linear mixing or DIIS extrapolation @@ -164,7 +165,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop n_diis = min(n_diis+1,max_diis) if(abs(rcond) > 1d-7) then - call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGW-eOld,eGW) + call DIIS_extrapolation(rcond,nOrb,nOrb,n_diis,error_diis,e_diis,eGW-eOld,eGW) else n_diis = 0 end if @@ -203,8 +204,8 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop !--------------------! ! TODO - !call RGWC(dotest, eta, nBas, nC, nO, nV, nR, nS, Om, rho, eHF, eGW, eGW, Z) - call RGWC(dotest, eta, nBas, nC, nO, nV, nR, nS, Om, rho, eHF, eHF, eGW, Z) + !call RGWC(dotest, eta, nOrb, nC, nO, nV, nR, nS, Om, rho, eHF, eGW, eGW, Z) + call RGWC(dotest, eta, nOrb, nC, nO, nV, nR, nS, Om, rho, eHF, eHF, eGW, Z) ! Deallocate memory @@ -215,7 +216,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop if(dophBSE) then 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) + nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eGW,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -235,7 +236,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop write(*,*) '-----------------------------------------------------------' write(*,*) - call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcBSE) + call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -252,7 +253,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop if(doppBSE) then - call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) + call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GW/ufBSE.f90 b/src/GW/ufBSE.f90 index 24568f3..592a089 100644 --- a/src/GW/ufBSE.f90 +++ b/src/GW/ufBSE.f90 @@ -1,4 +1,4 @@ -subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) +subroutine ufBSE(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) ! Unfold BSE@GW equations @@ -8,6 +8,7 @@ subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) ! Input variables integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -15,9 +16,9 @@ subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) integer,intent(in) :: nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: eGW(nOrb) ! Local variables @@ -84,12 +85,12 @@ subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) ia = 0 do i=nC+1,nO - do a=nO+1,nBas-nR + do a=nO+1,nOrb-nR ia = ia + 1 jb = 0 do j=nC+1,nO - do b=nO+1,nBas-nR + do b=nO+1,nOrb-nR jb = jb + 1 H(ia,jb) = (eGW(a) - eGW(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & @@ -107,14 +108,14 @@ subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) iajb=0 do i=nC+1,nO - do a=nO+1,nBas-nR + do a=nO+1,nOrb-nR do j=nC+1,nO - do b=nO+1,nBas-nR + do b=nO+1,nOrb-nR iajb = iajb + 1 kc = 0 do k=nC+1,nO - do c=nO+1,nBas-nR + do c=nO+1,nOrb-nR kc = kc + 1 tmp = sqrt(2d0)*Kronecker_delta(k,j)*ERI(b,a,c,i) @@ -139,16 +140,16 @@ subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) iajb = 0 do i=nC+1,nO - do a=nO+1,nBas-nR + do a=nO+1,nOrb-nR do j=nC+1,nO - do b=nO+1,nBas-nR + do b=nO+1,nOrb-nR iajb = iajb + 1 kcld = 0 do k=nC+1,nO - do c=nO+1,nBas-nR + do c=nO+1,nOrb-nR do l=nC+1,nO - do d=nO+1,nBas-nR + do d=nO+1,nOrb-nR kcld = kcld + 1 tmp = ((eHF(a) + eGW(b) - eHF(i) - eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(a,c) & diff --git a/src/GW/ufG0W0.f90 b/src/GW/ufG0W0.f90 index d9b8063..6e9b2d0 100644 --- a/src/GW/ufG0W0.f90 +++ b/src/GW/ufG0W0.f90 @@ -1,4 +1,4 @@ -subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) +subroutine ufG0W0(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Unfold G0W0 equations @@ -11,6 +11,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) logical,intent(in) :: TDA_W integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -18,8 +19,8 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) integer,intent(in) :: nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: eHF(nOrb) ! Local variables @@ -93,10 +94,10 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Memory allocation - allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS)) + allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS)) - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) - call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) + call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) @@ -106,7 +107,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Compute spectral weights ! !--------------------------! - call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho) + call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho) deallocate(Aph,Bph,XpY,XmY) @@ -154,7 +155,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ija = 0 do i=nC+1,nO do j=nC+1,nO - do a=nO+1,nBas-nR + do a=nO+1,nOrb-nR ija = ija + 1 H(1 ,1+ija) = sqrt(2d0)*ERI(p,a,i,j) @@ -170,8 +171,8 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) iab = 0 do i=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR + do a=nO+1,nOrb-nR + do b=nO+1,nOrb-nR iab = iab + 1 H(1 ,1+n2h1p+iab) = sqrt(2d0)*ERI(p,i,b,a) @@ -188,13 +189,13 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ija = 0 do i=nC+1,nO do j=nC+1,nO - do a=nO+1,nBas-nR + do a=nO+1,nOrb-nR ija = ija + 1 klc = 0 do k=nC+1,nO do l=nC+1,nO - do c=nO+1,nBas-nR + do c=nO+1,nOrb-nR klc = klc + 1 H(1+ija,1+klc) & @@ -215,14 +216,14 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) iab = 0 do i=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR + do a=nO+1,nOrb-nR + do b=nO+1,nOrb-nR iab = iab + 1 kcd = 0 do k=nC+1,nO - do c=nO+1,nBas-nR - do d=nO+1,nBas-nR + do c=nO+1,nOrb-nR + do d=nO+1,nOrb-nR kcd = kcd + 1 H(1+n2h1p+iab,1+n2h1p+kcd) & @@ -304,7 +305,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) iab = 0 do ia=1,nS - do b=nO+1,nBas-nR + do b=nO+1,nOrb-nR iab = iab + 1 H(1+n2h1p+iab,1+n2h1p+iab) = eHF(b) + Om(ia) @@ -318,7 +319,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) iab = 0 do ia=1,nS - do b=nO+1,nBas-nR + do b=nO+1,nOrb-nR iab = iab + 1 H(1 ,1+n2h1p+iab) = sqrt(2d0)*rho(p,b,ia) @@ -409,7 +410,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ija = 0 do i=nC+1,nO do j=nC+1,nO - do a=nO+1,nBas-nR + do a=nO+1,nOrb-nR ija = ija + 1 if(abs(H(1+ija,s)) > cutoff2) & @@ -422,8 +423,8 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) iab = 0 do i=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR + do a=nO+1,nOrb-nR + do b=nO+1,nOrb-nR iab = iab + 1 if(abs(H(1+n2h1p+iab,s)) > cutoff2) & @@ -478,7 +479,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) iab = 0 do ia=1,nS - do b=nO+1,nBas-nR + do b=nO+1,nOrb-nR iab = iab + 1 if(abs(H(1+n2h1p+iab,s)) > cutoff2) & diff --git a/src/GW/ufGW.f90 b/src/GW/ufGW.f90 index 6e9d441..d2e8686 100644 --- a/src/GW/ufGW.f90 +++ b/src/GW/ufGW.f90 @@ -1,4 +1,4 @@ -subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) +subroutine ufRGW(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Unfold GW equations @@ -11,6 +11,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) logical,intent(in) :: TDA_W integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -18,8 +19,8 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) integer,intent(in) :: nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: eHF(nOrb) ! Local variables @@ -68,7 +69,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) n2h1p = nO*nO*nV n2p1h = nV*nV*nO - nH = nBas + n2h1p + n2p1h + nH = nOrb + n2h1p + n2p1h ! Memory allocation @@ -89,14 +90,14 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Memory allocation - allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS)) + allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS)) ! Spin manifold ispin = 1 - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) - call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) + call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) @@ -106,7 +107,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Compute spectral weights ! !--------------------------! - call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho) + call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho) deallocate(Aph,Bph,XpY,XmY) @@ -141,7 +142,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Block F ! !---------! - do p=nC+1,nBas-nR + do p=nC+1,nOrb-nR H(p,p) = eHF(p) end do @@ -149,16 +150,16 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Block V2h1p ! !-------------! - do p=nC+1,nBas-nR + do p=nC+1,nOrb-nR ija = 0 do i=nC+1,nO do j=nC+1,nO - do a=nO+1,nBas-nR + do a=nO+1,nOrb-nR ija = ija + 1 - H(p ,nBas+ija) = sqrt(2d0)*ERI(p,a,i,j) - H(nBas+ija,p ) = sqrt(2d0)*ERI(p,a,i,j) + H(p ,nOrb+ija) = sqrt(2d0)*ERI(p,a,i,j) + H(nOrb+ija,p ) = sqrt(2d0)*ERI(p,a,i,j) end do end do @@ -170,16 +171,16 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Block V2p1h ! !-------------! - do p=nC+1,nBas-nR + do p=nC+1,nOrb-nR iab = 0 do i=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR + do a=nO+1,nOrb-nR + do b=nO+1,nOrb-nR iab = iab + 1 - H(p ,nBas+n2h1p+iab) = sqrt(2d0)*ERI(p,i,b,a) - H(nBas+n2h1p+iab,p ) = sqrt(2d0)*ERI(p,i,b,a) + H(p ,nOrb+n2h1p+iab) = sqrt(2d0)*ERI(p,i,b,a) + H(nOrb+n2h1p+iab,p ) = sqrt(2d0)*ERI(p,i,b,a) end do end do @@ -194,16 +195,16 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ija = 0 do i=nC+1,nO do j=nC+1,nO - do a=nO+1,nBas-nR + do a=nO+1,nOrb-nR ija = ija + 1 klc = 0 do k=nC+1,nO do l=nC+1,nO - do c=nO+1,nBas-nR + do c=nO+1,nOrb-nR klc = klc + 1 - H(nBas+ija,nBas+klc) & + H(nOrb+ija,nOrb+klc) & = ((eHF(i) + eHF(j) - eHF(a))*Kronecker_delta(j,l)*Kronecker_delta(a,c) & - 2d0*ERI(j,c,a,l))*Kronecker_delta(i,k) @@ -221,17 +222,17 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) iab = 0 do i=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR + do a=nO+1,nOrb-nR + do b=nO+1,nOrb-nR iab = iab + 1 kcd = 0 do k=nC+1,nO - do c=nO+1,nBas-nR - do d=nO+1,nBas-nR + do c=nO+1,nOrb-nR + do d=nO+1,nOrb-nR kcd = kcd + 1 - H(nBas+n2h1p+iab,nBas+n2h1p+kcd) & + H(nOrb+n2h1p+iab,nOrb+n2h1p+kcd) & = ((eHF(a) + eHF(b) - eHF(i))*Kronecker_delta(i,k)*Kronecker_delta(a,c) & + 2d0*ERI(a,k,i,c))*Kronecker_delta(b,d) @@ -275,7 +276,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Block F ! !---------! - do p=nC+1,nBas-nR + do p=nC+1,nOrb-nR H(p,p) = eHF(p) end do @@ -283,15 +284,15 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Block W2h1p ! !-------------! - do p=nC+1,nBas-nR + do p=nC+1,nOrb-nR ija = 0 do i=nC+1,nO do ja=1,nS ija = ija + 1 - H(p ,nBas+ija) = sqrt(2d0)*rho(p,i,ja) - H(nBas+ija,p ) = sqrt(2d0)*rho(p,i,ja) + H(p ,nOrb+ija) = sqrt(2d0)*rho(p,i,ja) + H(nOrb+ija,p ) = sqrt(2d0)*rho(p,i,ja) end do end do @@ -302,15 +303,15 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Block W2p1h ! !-------------! - do p=nC+1,nBas-nR + do p=nC+1,nOrb-nR iab = 0 do ia=1,nS - do b=nO+1,nBas-nR + do b=nO+1,nOrb-nR iab = iab + 1 - H(p ,nBas+n2h1p+iab) = sqrt(2d0)*rho(p,b,ia) - H(nBas+n2h1p+iab,p ) = sqrt(2d0)*rho(p,b,ia) + H(p ,nOrb+n2h1p+iab) = sqrt(2d0)*rho(p,b,ia) + H(nOrb+n2h1p+iab,p ) = sqrt(2d0)*rho(p,b,ia) end do end do @@ -326,7 +327,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) do ja=1,nS ija = ija + 1 - H(nBas+ija,nBas+ija) = eHF(i) - Om(ja) + H(nOrb+ija,nOrb+ija) = eHF(i) - Om(ja) end do end do @@ -337,10 +338,10 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) iab = 0 do ia=1,nS - do b=nO+1,nBas-nR + do b=nO+1,nOrb-nR iab = iab + 1 - H(nBas+n2h1p+iab,nBas+n2h1p+iab) = eHF(b) + Om(ia) + H(nOrb+n2h1p+iab,nOrb+n2h1p+iab) = eHF(b) + Om(ia) end do end do @@ -375,7 +376,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) Z(:) = 0d0 do s=1,nH - do p=nC+1,nBas-nR + do p=nC+1,nOrb-nR Z(s) = Z(s) + H(p,s)**2 end do end do @@ -425,7 +426,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') & ' (',p,') ',H(p,s),H(p,s)**2 end do - do p=nO+1,nBas-nR + do p=nO+1,nOrb-nR if(abs(H(p,s)) > cutoff2) & write(*,'(1X,A16,I3,A7,1X,F15.6,1X,F15.6)') & ' (',p,') ',H(p,s),H(p,s)**2 @@ -434,12 +435,12 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ija = 0 do i=nC+1,nO do j=nC+1,nO - do a=nO+1,nBas-nR + do a=nO+1,nOrb-nR ija = ija + 1 - if(abs(H(nBas+ija,s)) > cutoff2) & + if(abs(H(nOrb+ija,s)) > cutoff2) & write(*,'(1X,A3,I3,A1,I3,A6,I3,A7,1X,F15.6,1X,F15.6)') & - ' (',i,',',j,') -> (',a,') ',H(nBas+ija,s),H(nBas+ija,s)**2 + ' (',i,',',j,') -> (',a,') ',H(nOrb+ija,s),H(nOrb+ija,s)**2 end do end do @@ -447,13 +448,13 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) iab = 0 do i=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR + do a=nO+1,nOrb-nR + do b=nO+1,nOrb-nR iab = iab + 1 - if(abs(H(nBas+n2h1p+iab,s)) > cutoff2) & + if(abs(H(nOrb+n2h1p+iab,s)) > cutoff2) & write(*,'(1X,A7,I3,A6,I3,A1,I3,A3,1X,F15.6,1X,F15.6)') & - ' (',i,') -> (',a,',',b,') ',H(nBas+n2h1p+iab,s),H(nBas+n2h1p+iab,s)**2 + ' (',i,') -> (',a,',',b,') ',H(nOrb+n2h1p+iab,s),H(nOrb+n2h1p+iab,s)**2 end do end do @@ -487,7 +488,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') & ' (',p,' ) ',H(p,s),H(p,s)**2 end do - do p=nO+1,nBas-nR + do p=nO+1,nOrb-nR if(abs(H(p,s)) > cutoff2) & write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') & ' (',p,' ) ',H(p,s),H(p,s)**2 @@ -498,21 +499,21 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) do ja=1,nS ija = ija + 1 - if(abs(H(nBas+ija,s)) > cutoff2) & + if(abs(H(nOrb+ija,s)) > cutoff2) & write(*,'(1X,A7,I3,A1,I3,A12,1X,F15.6,1X,F15.6)') & - ' (',i,',',ja,') ',H(nBas+ija,s),H(nBas+ija,s)**2 + ' (',i,',',ja,') ',H(nOrb+ija,s),H(nOrb+ija,s)**2 end do end do iab = 0 do ia=1,nS - do b=nO+1,nBas-nR + do b=nO+1,nOrb-nR iab = iab + 1 - if(abs(H(nBas+n2h1p+iab,s)) > cutoff2) & + if(abs(H(nOrb+n2h1p+iab,s)) > cutoff2) & write(*,'(1X,A7,I3,A1,I3,A12,1X,F15.6,1X,F15.6)') & - ' (',ia,',',b,') ',H(nBas+n2h1p+iab,s),H(nBas+n2h1p+iab,s)**2 + ' (',ia,',',b,') ',H(nOrb+n2h1p+iab,s),H(nOrb+n2h1p+iab,s)**2 end do end do