From eee5f6bac8dbdfbac539e5952f9624f40f1758ab Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 27 Jul 2023 21:59:18 +0200 Subject: [PATCH] rename UGW routines --- src/GW/UG0W0.f90 | 13 +- ...stricted_QP_graph.f90 => UGW_QP_graph.f90} | 4 +- ...density.f90 => UGW_excitation_density.f90} | 4 +- src/GW/UGW_phACFDT.f90 | 6 +- ...icted_Bethe_Salpeter.f90 => UGW_phBSE.f90} | 8 +- ...trix.f90 => UGW_phBSE_static_kernel_A.f90} | 6 +- ...trix.f90 => UGW_phBSE_static_kernel_B.f90} | 6 +- src/GW/UGW_self_energy.f90 | 134 ++++++++++++++++++ src/GW/UGW_self_energy_diag.f90 | 126 ++++++++++++++++ src/GW/evUGW.f90 | 9 +- src/GW/qsUGW.f90 | 9 +- .../unrestricted_renormalization_factor.f90 | 90 ------------ .../unrestricted_self_energy_correlation.f90 | 118 --------------- ...estricted_self_energy_correlation_diag.f90 | 111 --------------- src/LR/phULR.f90 | 8 +- src/RPA/phUACFDT.f90 | 6 +- 16 files changed, 298 insertions(+), 360 deletions(-) rename src/GW/{unrestricted_QP_graph.f90 => UGW_QP_graph.f90} (94%) rename src/GW/{unrestricted_excitation_density.f90 => UGW_excitation_density.f90} (93%) rename src/GW/{unrestricted_Bethe_Salpeter.f90 => UGW_phBSE.f90} (93%) rename src/GW/{unrestricted_Bethe_Salpeter_A_matrix.f90 => UGW_phBSE_static_kernel_A.f90} (93%) rename src/GW/{unrestricted_Bethe_Salpeter_B_matrix.f90 => UGW_phBSE_static_kernel_B.f90} (93%) create mode 100644 src/GW/UGW_self_energy.f90 create mode 100644 src/GW/UGW_self_energy_diag.f90 delete mode 100644 src/GW/unrestricted_renormalization_factor.f90 delete mode 100644 src/GW/unrestricted_self_energy_correlation.f90 delete mode 100644 src/GW/unrestricted_self_energy_correlation_diag.f90 diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index 9ae32e4..953f1bf 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -118,7 +118,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons ! Excitation densities ! !----------------------! - call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) !------------------------------------------------! ! Compute self-energy and renormalization factor ! @@ -131,8 +131,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons else - call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM) - call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z) + call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,Z,EcGM) end if @@ -154,8 +153,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons ! Find graphical solution of the QP equation do is=1,nspin - call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is), & - OmRPA,rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is)) + call UGW_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is), & + OmRPA,rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is)) end do end if @@ -177,8 +176,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons if(BSE) then - call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,S, & - ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eGW,EcBSE) + call UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,S, & + ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eGW,EcBSE) if(exchange_kernel) then diff --git a/src/GW/unrestricted_QP_graph.f90 b/src/GW/UGW_QP_graph.f90 similarity index 94% rename from src/GW/unrestricted_QP_graph.f90 rename to src/GW/UGW_QP_graph.f90 index 6162b7c..1a5899c 100644 --- a/src/GW/unrestricted_QP_graph.f90 +++ b/src/GW/UGW_QP_graph.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW) +subroutine UGW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW) ! Compute the graphical solution of the QP equation @@ -80,4 +80,4 @@ subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eG end do -end subroutine unrestricted_QP_graph +end subroutine diff --git a/src/GW/unrestricted_excitation_density.f90 b/src/GW/UGW_excitation_density.f90 similarity index 93% rename from src/GW/unrestricted_excitation_density.f90 rename to src/GW/UGW_excitation_density.f90 index 5822d65..571051b 100644 --- a/src/GW/unrestricted_excitation_density.f90 +++ b/src/GW/UGW_excitation_density.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho) +subroutine UGW_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho) ! Compute excitation densities for unrestricted reference @@ -103,4 +103,4 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ER enddo enddo -end subroutine unrestricted_excitation_density +end subroutine diff --git a/src/GW/UGW_phACFDT.f90 b/src/GW/UGW_phACFDT.f90 index 63fe4df..87e681d 100644 --- a/src/GW/UGW_phACFDT.f90 +++ b/src/GW/UGW_phACFDT.f90 @@ -95,7 +95,7 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,eW, & ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) ! Spin-conserved manifold @@ -122,7 +122,7 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) end if @@ -176,7 +176,7 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) end if diff --git a/src/GW/unrestricted_Bethe_Salpeter.f90 b/src/GW/UGW_phBSE.f90 similarity index 93% rename from src/GW/unrestricted_Bethe_Salpeter.f90 rename to src/GW/UGW_phBSE.f90 index 50b9801..eeba00d 100644 --- a/src/GW/unrestricted_Bethe_Salpeter.f90 +++ b/src/GW/UGW_phBSE.f90 @@ -1,6 +1,6 @@ -subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta, & - nBas,nC,nO,nV,nR,nS,S,ERI_aaaa,ERI_aabb,ERI_bbbb, & - dipole_int_aa,dipole_int_bb,cW,eW,eGW,EcBSE) +subroutine UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta, & + nBas,nC,nO,nV,nR,nS,S,ERI_aaaa,ERI_aabb,ERI_bbbb, & + dipole_int_aa,dipole_int_bb,cW,eW,eGW,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -82,7 +82,7 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_f call phULR(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & eW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) !----------------------------! diff --git a/src/GW/unrestricted_Bethe_Salpeter_A_matrix.f90 b/src/GW/UGW_phBSE_static_kernel_A.f90 similarity index 93% rename from src/GW/unrestricted_Bethe_Salpeter_A_matrix.f90 rename to src/GW/UGW_phBSE_static_kernel_A.f90 index 48752f8..b37bb4b 100644 --- a/src/GW/unrestricted_Bethe_Salpeter_A_matrix.f90 +++ b/src/GW/UGW_phBSE_static_kernel_A.f90 @@ -1,5 +1,5 @@ -subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,eGW, & - ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,A_lr) +subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,eGW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,A_lr) ! Compute the extra term for Bethe-Salpeter equation for linear response in the unrestricted formalism @@ -149,4 +149,4 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n end if -end subroutine unrestricted_Bethe_Salpeter_A_matrix +end subroutine diff --git a/src/GW/unrestricted_Bethe_Salpeter_B_matrix.f90 b/src/GW/UGW_phBSE_static_kernel_B.f90 similarity index 93% rename from src/GW/unrestricted_Bethe_Salpeter_B_matrix.f90 rename to src/GW/UGW_phBSE_static_kernel_B.f90 index c4a45f5..81cb378 100644 --- a/src/GW/unrestricted_Bethe_Salpeter_B_matrix.f90 +++ b/src/GW/UGW_phBSE_static_kernel_B.f90 @@ -1,5 +1,5 @@ -subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & - ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,B_lr) +subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & + ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,B_lr) ! Compute the extra term for Bethe-Salpeter equation for linear response @@ -150,4 +150,4 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n end if -end subroutine unrestricted_Bethe_Salpeter_B_matrix +end subroutine diff --git a/src/GW/UGW_self_energy.f90 b/src/GW/UGW_self_energy.f90 new file mode 100644 index 0000000..cc02b14 --- /dev/null +++ b/src/GW/UGW_self_energy.f90 @@ -0,0 +1,134 @@ +subroutine UGW_self_energy(eta,nBas,nC,nO,nV,nR,nSt,e,Om,rho,Sig,Z,EcGM) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nSt + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: Om(nSt) + double precision,intent(in) :: rho(nBas,nBas,nSt,nspin) + +! Local variables + + integer :: i,a,p,q,m + double precision :: num,eps + +! Output variables + + double precision,intent(out) :: Sig(nBas,nBas,nspin) + double precision,intent(out) :: Z(nBas,nspin) + double precision :: EcGM(nspin) + +! Initialize + + Sig(:,:,:) = 0d0 + Z(:,:) = 0d0 + EcGM(:) = 0d0 + +!--------------! +! Spin-up part ! +!--------------! + + ! Occupied part of the correlation self-energy + + do p=nC(1)+1,nBas-nR(1) + do q=nC(1)+1,nBas-nR(1) + do i=nC(1)+1,nO(1) + do m=1,nSt + eps = e(p,1) - e(i,1) + Om(m) + num = rho(p,i,m,1)*rho(q,i,m,1) + Sig(p,q,1) = Sig(p,q,1) + num*eps/(eps**2 + eta**2) + if(p == q) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC(1)+1,nBas-nR(1) + do q=nC(1)+1,nBas-nR(1) + do a=nO(1)+1,nBas-nR(1) + do m=1,nSt + eps = e(p,1) - e(a,1) - Om(m) + num = rho(p,a,m,1)*rho(q,a,m,1) + Sig(p,q,1) = Sig(p,q,1) + num*eps/(eps**2 + eta**2) + if(p == q) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + end do + end do + end do + + ! GM correlation energy + + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + do m=1,nSt + eps = e(a,1) - e(i,1) + Om(m) + num = rho(a,i,m,1)**2 + EcGM(1) = EcGM(1) - num*eps/(eps**2 + eta**2) + end do + end do + end do + +!----------------! +! Spin-down part ! +!----------------! + + ! Occupied part of the correlation self-energy + + do p=nC(2)+1,nBas-nR(2) + do q=nC(2)+1,nBas-nR(2) + do i=nC(2)+1,nO(2) + do m=1,nSt + eps = e(p,2) - e(i,2) + Om(m) + num = rho(p,i,m,2)*rho(q,i,m,2) + Sig(p,q,2) = Sig(p,q,2) + num*eps/(eps**2 + eta**2) + if(p == q) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC(2)+1,nBas-nR(2) + do q=nC(2)+1,nBas-nR(2) + do a=nO(2)+1,nBas-nR(2) + do m=1,nSt + eps = e(p,2) - e(a,2) - Om(m) + num = rho(p,a,m,2)*rho(q,a,m,2) + Sig(p,q,2) = Sig(p,q,2) + num*eps/(eps**2 + eta**2) + if(p == q) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + end do + end do + end do + + ! GM correlation energy + + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + do m=1,nSt + eps = e(a,2) - e(i,2) + Om(m) + num = rho(a,i,m,2)**2 + EcGM(2) = EcGM(2) - num*eps/(eps**2 + eta**2) + end do + end do + end do + +! Compute renormalization factor from derivative + + Z(:,:) = 1d0/(1d0 - Z(:,:)) + +end subroutine diff --git a/src/GW/UGW_self_energy_diag.f90 b/src/GW/UGW_self_energy_diag.f90 new file mode 100644 index 0000000..a895878 --- /dev/null +++ b/src/GW/UGW_self_energy_diag.f90 @@ -0,0 +1,126 @@ +subroutine UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nSt,e,Om,rho,Sig,Z,EcGM) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nSt + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: Om(nSt) + double precision,intent(in) :: rho(nBas,nBas,nSt,nspin) + +! Local variables + + integer :: i,a,p,m + double precision :: num,eps + +! Output variables + + double precision,intent(out) :: Sig(nBas,nspin) + double precision,intent(out) :: Z(nBas,nspin) + double precision :: EcGM(nspin) + +! Initialize + + Sig(:,:) = 0d0 + Z(:,:) = 0d0 + EcGM(:) = 0d0 + +!--------------! +! Spin-up part ! +!--------------! + + ! Occupied part of the correlation self-energy + + do p=nC(1)+1,nBas-nR(1) + do i=nC(1)+1,nO(1) + do m=1,nSt + eps = e(p,1) - e(i,1) + Om(m) + num = rho(p,i,m,1)**2 + Sig(p,1) = Sig(p,1) + num*eps/(eps**2 + eta**2) + Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC(1)+1,nBas-nR(1) + do a=nO(1)+1,nBas-nR(1) + do m=1,nSt + eps = e(p,1) - e(a,1) - Om(m) + num = rho(p,a,m,1)**2 + Sig(p,1) = Sig(p,1) + num*eps/(eps**2 + eta**2) + Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + end do + end do + + ! GM correlation energy + + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + do m=1,nSt + eps = e(a,1) - e(i,1) + Om(m) + num = rho(a,i,m,1)**2 + EcGM(1) = EcGM(1) - num*eps/(eps**2 + eta**2) + end do + end do + end do + +!----------------! +! Spin-down part ! +!----------------! + + ! Occupied part of the correlation self-energy + + do p=nC(2)+1,nBas-nR(2) + do i=nC(2)+1,nO(2) + do m=1,nSt + eps = e(p,2) - e(i,2) + Om(m) + num = rho(p,i,m,2)**2 + Sig(p,2) = Sig(p,2) + num*eps/(eps**2 + eta**2) + Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC(2)+1,nBas-nR(2) + do a=nO(2)+1,nBas-nR(2) + do m=1,nSt + eps = e(p,2) - e(a,2) - Om(m) + num = rho(p,a,m,2)**2 + Sig(p,2) = Sig(p,2) + num*eps/(eps**2 + eta**2) + Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + end do + end do + + ! GM correlation energy + + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + do m=1,nSt + eps = e(a,2) - e(i,2) + Om(m) + num = rho(a,i,m,2)**2 + EcGM(2) = EcGM(2) - num*eps/(eps**2 + eta**2) + end do + end do + end do + +! Compute renormalization factor from derivative + + Z(:,:) = 1d0/(1d0 - Z(:,:)) + +end subroutine diff --git a/src/GW/evUGW.f90 b/src/GW/evUGW.f90 index 0030bfc..e6a00f3 100644 --- a/src/GW/evUGW.f90 +++ b/src/GW/evUGW.f90 @@ -135,7 +135,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, ! Excitation densities ! !----------------------! - call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) !------------------------------------------------! ! Compute self-energy and renormalization factor ! @@ -148,8 +148,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, else - call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,EcGM) - call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z) + call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,Z,EcGM) end if @@ -222,8 +221,8 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, if(BSE) then - call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & - S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eGW,eGW,EcBSE) + call UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & + S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eGW,eGW,EcBSE) if(exchange_kernel) then diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90 index a835f52..dacdaaa 100644 --- a/src/GW/qsUGW.f90 +++ b/src/GW/qsUGW.f90 @@ -204,7 +204,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, ! Excitation densities ! !----------------------! - call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) !------------------------------------------------! ! Compute self-energy and renormalization factor ! @@ -217,8 +217,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, else - call unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,EcGM) - call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z) + call UGW_self_energy(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,Z,EcGM) end if @@ -362,8 +361,8 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, if(BSE) then - call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & - S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,c,eGW,eGW,EcBSE) + call UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & + S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,c,eGW,eGW,EcBSE) if(exchange_kernel) then diff --git a/src/GW/unrestricted_renormalization_factor.f90 b/src/GW/unrestricted_renormalization_factor.f90 deleted file mode 100644 index 147b483..0000000 --- a/src/GW/unrestricted_renormalization_factor.f90 +++ /dev/null @@ -1,90 +0,0 @@ -subroutine unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,Z) - -! Compute the renormalization factor in the unrestricted formalism - - implicit none - include 'parameters.h' - -! Input variables - - double precision,intent(in) :: eta - integer,intent(in) :: nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) - integer,intent(in) :: nSt - double precision,intent(in) :: e(nBas,nspin) - double precision,intent(in) :: Omega(nSt) - double precision,intent(in) :: rho(nBas,nBas,nSt,nspin) - -! Local variables - - integer :: i,a,p,jb - double precision :: eps - -! Output variables - - double precision,intent(out) :: Z(nBas,nspin) - -! Initialize - - Z(:,:) = 0d0 - -!--------------! -! Spin-up part ! -!--------------! - - ! Occupied part of the renormalization factor - - do p=nC(1)+1,nBas-nR(1) - do i=nC(1)+1,nO(1) - do jb=1,nSt - eps = e(p,1) - e(i,1) + Omega(jb) - Z(p,1) = Z(p,1) + rho(p,i,jb,1)**2*(eps/(eps**2 + eta**2))**2 - end do - end do - end do - - ! Virtual part of the correlation self-energy - - do p=nC(1)+1,nBas-nR(1) - do a=nO(1)+1,nBas-nR(1) - do jb=1,nSt - eps = e(p,1) - e(a,1) - Omega(jb) - Z(p,1) = Z(p,1) + rho(p,a,jb,1)**2*(eps/(eps**2 + eta**2))**2 - end do - end do - end do - -!----------------! -! Spin-down part ! -!----------------! - - ! Occupied part of the correlation self-energy - - do p=nC(2)+1,nBas-nR(2) - do i=nC(2)+1,nO(2) - do jb=1,nSt - eps = e(p,2) - e(i,2) + Omega(jb) - Z(p,2) = Z(p,2) + rho(p,i,jb,2)**2*(eps/(eps**2 + eta**2))**2 - end do - end do - end do - - ! Virtual part of the correlation self-energy - - do p=nC(2)+1,nBas-nR(2) - do a=nO(2)+1,nBas-nR(2) - do jb=1,nSt - eps = e(p,2) - e(a,2) - Omega(jb) - Z(p,2) = Z(p,2) + rho(p,a,jb,2)**2*(eps/(eps**2 + eta**2))**2 - end do - end do - end do - -! Final rescaling - - Z(:,:) = 1d0/(1d0 + Z(:,:)) - -end subroutine unrestricted_renormalization_factor diff --git a/src/GW/unrestricted_self_energy_correlation.f90 b/src/GW/unrestricted_self_energy_correlation.f90 deleted file mode 100644 index de1ae28..0000000 --- a/src/GW/unrestricted_self_energy_correlation.f90 +++ /dev/null @@ -1,118 +0,0 @@ -subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC,EcGM) - -! Compute diagonal of the correlation part of the self-energy - - implicit none - include 'parameters.h' - -! Input variables - - double precision,intent(in) :: eta - integer,intent(in) :: nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) - integer,intent(in) :: nSt - double precision,intent(in) :: e(nBas,nspin) - double precision,intent(in) :: Omega(nSt) - double precision,intent(in) :: rho(nBas,nBas,nSt,nspin) - -! Local variables - - integer :: i,a,p,q,jb - double precision :: eps - -! Output variables - - double precision,intent(out) :: SigC(nBas,nBas,nspin) - double precision :: EcGM(nspin) - -! Initialize - - SigC(:,:,:) = 0d0 - EcGM(:) = 0d0 - -!--------------! -! Spin-up part ! -!--------------! - - ! Occupied part of the correlation self-energy - - do p=nC(1)+1,nBas-nR(1) - do q=nC(1)+1,nBas-nR(1) - do i=nC(1)+1,nO(1) - do jb=1,nSt - eps = e(p,1) - e(i,1) + Omega(jb) - SigC(p,q,1) = SigC(p,q,1) + rho(p,i,jb,1)*rho(q,i,jb,1)*eps/(eps**2 + eta**2) - end do - end do - end do - end do - - ! Virtual part of the correlation self-energy - - do p=nC(1)+1,nBas-nR(1) - do q=nC(1)+1,nBas-nR(1) - do a=nO(1)+1,nBas-nR(1) - do jb=1,nSt - eps = e(p,1) - e(a,1) - Omega(jb) - SigC(p,q,1) = SigC(p,q,1) + rho(p,a,jb,1)*rho(q,a,jb,1)*eps/(eps**2 + eta**2) - end do - end do - end do - end do - - ! GM correlation energy - - do i=nC(1)+1,nO(1) - do a=nO(1)+1,nBas-nR(1) - do jb=1,nSt - eps = e(a,1) - e(i,1) + Omega(jb) - EcGM(1) = EcGM(1) - rho(a,i,jb,1)**2*eps/(eps**2 + eta**2) - end do - end do - end do - -!----------------! -! Spin-down part ! -!----------------! - - ! Occupied part of the correlation self-energy - - do p=nC(2)+1,nBas-nR(2) - do q=nC(2)+1,nBas-nR(2) - do i=nC(2)+1,nO(2) - do jb=1,nSt - eps = e(p,2) - e(i,2) + Omega(jb) - SigC(p,q,2) = SigC(p,q,2) + rho(p,i,jb,2)*rho(q,i,jb,2)*eps/(eps**2 + eta**2) - end do - end do - end do - end do - - ! Virtual part of the correlation self-energy - - do p=nC(2)+1,nBas-nR(2) - do q=nC(2)+1,nBas-nR(2) - do a=nO(2)+1,nBas-nR(2) - do jb=1,nSt - eps = e(p,2) - e(a,2) - Omega(jb) - SigC(p,q,2) = SigC(p,q,2) + rho(p,a,jb,2)*rho(q,a,jb,2)*eps/(eps**2 + eta**2) - end do - end do - end do - end do - - ! GM correlation energy - - do i=nC(2)+1,nO(2) - do a=nO(2)+1,nBas-nR(2) - do jb=1,nSt - eps = e(a,2) - e(i,2) + Omega(jb) - EcGM(2) = EcGM(2) - rho(a,i,jb,2)**2*eps/(eps**2 + eta**2) - end do - end do - end do - -end subroutine unrestricted_self_energy_correlation diff --git a/src/GW/unrestricted_self_energy_correlation_diag.f90 b/src/GW/unrestricted_self_energy_correlation_diag.f90 deleted file mode 100644 index 8290d9d..0000000 --- a/src/GW/unrestricted_self_energy_correlation_diag.f90 +++ /dev/null @@ -1,111 +0,0 @@ -subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC,EcGM) - -! Compute diagonal of the correlation part of the self-energy - - implicit none - include 'parameters.h' - -! Input variables - - double precision,intent(in) :: eta - integer,intent(in) :: nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) - integer,intent(in) :: nSt - double precision,intent(in) :: e(nBas,nspin) - double precision,intent(in) :: Omega(nSt) - double precision,intent(in) :: rho(nBas,nBas,nSt,nspin) - -! Local variables - - integer :: i,a,p,q,jb - double precision :: eps - -! Output variables - - double precision,intent(out) :: SigC(nBas,nspin) - double precision :: EcGM(nspin) - -! Initialize - - SigC(:,:) = 0d0 - EcGM(:) = 0d0 - -!--------------! -! Spin-up part ! -!--------------! - - ! Occupied part of the correlation self-energy - - do p=nC(1)+1,nBas-nR(1) - do i=nC(1)+1,nO(1) - do jb=1,nSt - eps = e(p,1) - e(i,1) + Omega(jb) - SigC(p,1) = SigC(p,1) + rho(p,i,jb,1)**2*eps/(eps**2 + eta**2) - end do - end do - end do - - ! Virtual part of the correlation self-energy - - do p=nC(1)+1,nBas-nR(1) - do a=nO(1)+1,nBas-nR(1) - do jb=1,nSt - eps = e(p,1) - e(a,1) - Omega(jb) - SigC(p,1) = SigC(p,1) + rho(p,a,jb,1)**2*eps/(eps**2 + eta**2) - end do - end do - end do - - ! GM correlation energy - - do i=nC(1)+1,nO(1) - do a=nO(1)+1,nBas-nR(1) - do jb=1,nSt - eps = e(a,1) - e(i,1) + Omega(jb) - EcGM(1) = EcGM(1) - rho(a,i,jb,1)**2*eps/(eps**2 + eta**2) - end do - end do - end do - -!----------------! -! Spin-down part ! -!----------------! - - ! Occupied part of the correlation self-energy - - do p=nC(2)+1,nBas-nR(2) - do i=nC(2)+1,nO(2) - do jb=1,nSt - eps = e(p,2) - e(i,2) + Omega(jb) - SigC(p,2) = SigC(p,2) + rho(p,i,jb,2)**2*eps/(eps**2 + eta**2) - end do - end do - end do - - ! Virtual part of the correlation self-energy - - do p=nC(2)+1,nBas-nR(2) - do a=nO(2)+1,nBas-nR(2) - do jb=1,nSt - eps = e(p,2) - e(a,2) - Omega(jb) - SigC(p,2) = SigC(p,2) + rho(p,a,jb,2)**2*eps/(eps**2 + eta**2) - end do - end do - end do - - ! GM correlation energy - - do i=nC(2)+1,nO(2) - do a=nO(2)+1,nBas-nR(2) - do jb=1,nSt - eps = e(a,2) - e(i,2) + Omega(jb) - EcGM(2) = EcGM(2) - rho(a,i,jb,2)**2*eps/(eps**2 + eta**2) - end do - end do - end do - - -end subroutine unrestricted_self_energy_correlation_diag diff --git a/src/LR/phULR.f90 b/src/LR/phULR.f90 index 53080b1..9a0018a 100644 --- a/src/LR/phULR.f90 +++ b/src/LR/phULR.f90 @@ -57,8 +57,8 @@ subroutine phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) if(BSE) & - call unrestricted_Bethe_Salpeter_A_matrix(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,Aph) + call UGW_phBSE_static_kernel_A(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,Aph) ! Tamm-Dancoff approximation @@ -75,8 +75,8 @@ subroutine phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) if(BSE) & - call unrestricted_Bethe_Salpeter_B_matrix(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,Bph) + call UGW_phBSE_static_kernel_B(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,Bph) ! Build A + B and A - B matrices diff --git a/src/RPA/phUACFDT.f90 b/src/RPA/phUACFDT.f90 index e913d55..2174029 100644 --- a/src/RPA/phUACFDT.f90 +++ b/src/RPA/phUACFDT.f90 @@ -94,7 +94,7 @@ subroutine phUACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,eW, & ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) ! Spin-conserved manifold @@ -121,7 +121,7 @@ subroutine phUACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) end if @@ -175,7 +175,7 @@ subroutine phUACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) end if