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