From 08f3567d205d0935955e00bb937e26755e8d3682 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 23 Sep 2020 22:48:54 +0200 Subject: [PATCH] sf-BSE --- input/options | 4 +- src/QuAcK/UG0W0.f90 | 8 +- src/QuAcK/UHF.f90 | 3 +- src/QuAcK/URPAx.f90 | 8 +- src/QuAcK/UdRPA.f90 | 8 +- src/QuAcK/unrestricted_Bethe_Salpeter.f90 | 82 ++++---- .../unrestricted_Bethe_Salpeter_A_matrix.f90 | 180 +++++++++-------- .../unrestricted_Bethe_Salpeter_B_matrix.f90 | 185 ++++++++++-------- src/QuAcK/unrestricted_linear_response.f90 | 17 +- .../unrestricted_linear_response_A_matrix.f90 | 4 +- .../unrestricted_linear_response_B_matrix.f90 | 14 +- 11 files changed, 275 insertions(+), 238 deletions(-) diff --git a/input/options b/input/options index d55ed53..337b5e1 100644 --- a/input/options +++ b/input/options @@ -5,11 +5,11 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: singlet triplet spin_conserved spin_flip TDA - T T T F F + T T T T F # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.0 3 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 - 256 0.00001 T 5 F 0.0 F F F F F + 256 0.00001 T 5 T 0.0 F F F F F # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn diff --git a/src/QuAcK/UG0W0.f90 b/src/QuAcK/UG0W0.f90 index ac808cd..10b2b58 100644 --- a/src/QuAcK/UG0W0.f90 +++ b/src/QuAcK/UG0W0.f90 @@ -112,8 +112,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ispin = 1 - call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, & - eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) + call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & + eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) if(print_W) call print_excitation('RPA@UHF',5,nS_sc,Omega_sc) @@ -165,8 +165,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! Compute the RPA correlation energy - call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) + call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/QuAcK/UHF.f90 b/src/QuAcK/UHF.f90 index 28fc3c8..8559d03 100644 --- a/src/QuAcK/UHF.f90 +++ b/src/QuAcK/UHF.f90 @@ -232,6 +232,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUH ! Compute final UHF energy - call print_UHF(nBas,nO(:),e(:,:),c(:,:,:),ENuc,ET(:),EV(:),EJ(:),Ex(:),EUHF) + call matout(nBas,2,e) + call print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) end subroutine UHF diff --git a/src/QuAcK/URPAx.f90 b/src/QuAcK/URPAx.f90 index 4ae48cc..149c62d 100644 --- a/src/QuAcK/URPAx.f90 +++ b/src/QuAcK/URPAx.f90 @@ -73,8 +73,8 @@ subroutine URPAx(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) - call unrestricted_linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sc,EcRPAx(ispin),Omega_sc,XpY_sc,XmY_sc) + call unrestricted_linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPAx(ispin),Omega_sc,XpY_sc,XmY_sc) call print_excitation('URPAx ',5,nS_sc,Omega_sc) ! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) @@ -95,8 +95,8 @@ subroutine URPAx(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) - call unrestricted_linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,1d0,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sf,EcRPAx(ispin),Omega_sf,XpY_sf,XmY_sf) + call unrestricted_linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,1d0,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sf,rho_sf,EcRPAx(ispin),Omega_sf,XpY_sf,XmY_sf) call print_excitation('URPAx ',6,nS_sf,Omega_sf) ! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) diff --git a/src/QuAcK/UdRPA.f90 b/src/QuAcK/UdRPA.f90 index da3d019..9fff2be 100644 --- a/src/QuAcK/UdRPA.f90 +++ b/src/QuAcK/UdRPA.f90 @@ -73,8 +73,8 @@ subroutine UdRPA(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) - call unrestricted_linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) + call unrestricted_linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) call print_excitation('URPA ',5,nS_sc,Omega_sc) ! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) @@ -95,8 +95,8 @@ subroutine UdRPA(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) - call unrestricted_linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,1d0,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf) + call unrestricted_linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,1d0,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sf,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf) call print_excitation('URPA ',6,nS_sf,Omega_sf) ! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter.f90 index 42e7d6d..992ee9b 100644 --- a/src/QuAcK/unrestricted_Bethe_Salpeter.f90 +++ b/src/QuAcK/unrestricted_Bethe_Salpeter.f90 @@ -36,8 +36,8 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, integer :: ispin integer :: isp_W + integer :: nS_aa,nS_bb,nS_sc - integer :: nS_ab,nS_ba,nS_sf double precision,allocatable :: OmRPA_sc(:) double precision,allocatable :: XpY_RPA_sc(:,:) double precision,allocatable :: XmY_RPA_sc(:,:) @@ -46,6 +46,11 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, double precision,allocatable :: XpY_BSE_sc(:,:) double precision,allocatable :: XmY_BSE_sc(:,:) + integer :: nS_ab,nS_ba,nS_sf + double precision,allocatable :: OmBSE_sf(:) + double precision,allocatable :: XpY_BSE_sf(:,:) + double precision,allocatable :: XmY_BSE_sf(:,:) + ! Output variables double precision,intent(out) :: EcRPA(nspin) @@ -55,38 +60,41 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, ! Spin-conserved excitations ! !----------------------------! - if(spin_conserved) then + isp_W = 1 + + ! Memory allocation + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + allocate(OmRPA_sc(nS_sc),XpY_RPA_sc(nS_sc,nS_sc),XmY_RPA_sc(nS_sc,nS_sc),rho_RPA_sc(nBas,nBas,nS_sc,nspin)) + + ! Compute spin-conserved RPA screening + + call unrestricted_linear_response(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,ERI_abab,OmRPA_sc,rho_RPA_sc,EcRPA(isp_W), & + OmRPA_sc,XpY_RPA_sc,XmY_RPA_sc) + +! call print_excitation('RPA@UG0W0',5,nS_sc,OmRPA_sc) + + call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb, & + XpY_RPA_sc,rho_RPA_sc) + + if(spin_conserved) then ispin = 1 - isp_W = 1 + EcBSE(ispin) = 0d0 - ! Memory allocation - - nS_aa = nS(1) - nS_bb = nS(2) - nS_sc = nS_aa + nS_bb - - allocate(OmRPA_sc(nS_sc),XpY_RPA_sc(nS_sc,nS_sc),XmY_RPA_sc(nS_sc,nS_sc),rho_RPA_sc(nBas,nBas,nS_sc,nspin)) allocate(OmBSE_sc(nS_sc),XpY_BSE_sc(nS_sc,nS_sc),XmY_BSE_sc(nS_sc,nS_sc)) - ! Compute spin-conserved RPA screening - - call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, & - eW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_RPA_sc,EcRPA(ispin), & - OmRPA_sc,XpY_RPA_sc,XmY_RPA_sc) - -! call print_excitation('RPA@UG0W0',5,nS_sc,OmRPA_sc) - - call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb, & - XpY_RPA_sc,rho_RPA_sc) - ! Compute spin-conserved BSE excitation energies OmBSE_sc(:) = OmRPA_sc(:) - call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_RPA_sc,EcBSE(ispin), & + call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA_sc,rho_RPA_sc,EcBSE(ispin), & OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) call print_excitation('BSE@UG0W0',5,nS_sc,OmBSE_sc) @@ -117,25 +125,25 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, ! Spin-flip excitations ! !-----------------------! -!if(spin_flip) then + if(spin_flip) then -! ispin = 2 -! isp_W = 1 -! EcBSE(ispin) = 0d0 + ispin = 2 -! ! Compute (singlet) RPA screening + EcBSE(ispin) = 0d0 -! call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & -! rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) -! call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA(:,:,ispin),rho_RPA(:,:,:,ispin)) + ! Memory allocation -! ! Compute BSE excitation energies + nS_ab = (nO(1) - nC(1))*(nV(2) - nR(2)) + nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1)) + nS_sf = nS_ab + nS_ba + + allocate(OmBSE_sf(nS_sf),XpY_BSE_sf(nS_sf,nS_sf),XmY_BSE_sf(nS_sf,nS_sf)) -! OmBSE(:,ispin) = OmRPA(:,ispin) + call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,1d0, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA_sc,rho_RPA_sc,EcBSE(ispin), & + OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) -! call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & -! rho_RPA(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) -! call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin)) + call print_excitation('BSE@UG0W0',6,nS_sf,OmBSE_sf) !------------------------------------------------- ! Compute the dynamical screening at the BSE level @@ -157,6 +165,6 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, ! end if -! end if + end if end subroutine unrestricted_Bethe_Salpeter diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 index 175704a..55d9a25 100644 --- a/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 +++ b/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & +subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,A_lr) ! Compute the extra term for Bethe-Salpeter equation for linear response in the unrestricted formalism @@ -8,6 +8,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt ! Input variables + integer,intent(in) :: ispin integer,intent(in) :: nBas integer,intent(in) :: nC(nspin) integer,intent(in) :: nO(nspin) @@ -16,14 +17,15 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt integer,intent(in) :: nSa integer,intent(in) :: nSb integer,intent(in) :: nSt + integer,intent(in) :: nSsc double precision,intent(in) :: eta double precision,intent(in) :: lambda double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega(nSt) - double precision,intent(in) :: rho(nBas,nBas,nSt,nspin) + double precision,intent(in) :: Omega(nSsc) + double precision,intent(in) :: rho(nBas,nBas,nSsc,nspin) ! Local variables @@ -35,108 +37,116 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt double precision,intent(out) :: A_lr(nSt,nSt) -!--------------------------------! -! Build part A of the BSE matrix ! -!--------------------------------! +!--------------------------------------------------! +! Build BSE matrix for spin-conserving transitions ! +!--------------------------------------------------! - ! aaaa block + if(ispin == 1) then - ia = 0 - do i=nC(1)+1,nO(1) - do a=nO(1)+1,nBas-nR(1) - ia = ia + 1 - jb = 0 - do j=nC(1)+1,nO(1) - do b=nO(1)+1,nBas-nR(1) - jb = jb + 1 + ! aaaa block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + chi = 0d0 + do kc=1,nSsc + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps + enddo + + A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aaaa(i,b,j,a) + 4d0*lambda*chi - chi = 0d0 - do kc=1,nSt - eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps & - + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps enddo - - A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aaaa(i,b,j,a) + 2d0*lambda*chi - enddo enddo enddo - enddo - ! aabb block - - ia = 0 - do i=nC(1)+1,nO(1) - do a=nO(1)+1,nBas-nR(1) - ia = ia + 1 - jb = 0 - do j=nC(2)+1,nO(2) - do b=nO(2)+1,nBas-nR(2) - jb = jb + 1 + ! bbbb block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + chi = 0d0 + do kc=1,nSsc + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps + enddo + + A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,b,j,a) + 4d0*lambda*chi - chi = 0d0 - do kc=1,nSt - eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps & - + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps enddo - - A_lr(ia,nSa+jb) = A_lr(ia,nSa+jb) - lambda*ERI_aabb(i,b,j,a) + 2d0*lambda*chi - enddo enddo enddo - enddo - ! bbaa block + end if - ia = 0 - do i=nC(2)+1,nO(2) - do a=nO(2)+1,nBas-nR(2) - ia = ia + 1 - jb = 0 - do j=nC(1)+1,nO(1) - do b=nO(1)+1,nBas-nR(1) - jb = jb + 1 - - chi = 0d0 - do kc=1,nSt - eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps & - + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps - enddo +!--------------------------------------------! +! Build BSE matrix for spin-flip transitions ! +!--------------------------------------------! - A_lr(nSa+ia,jb) = A_lr(nSa+ia,jb) - lambda*ERI_aabb(b,i,a,j) + 2d0*lambda*chi + if(ispin == 2) then - enddo - enddo - enddo - enddo + ! abab block - ! bbbb block + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 - ia = 0 - do i=nC(2)+1,nO(2) - do a=nO(2)+1,nBas-nR(2) - ia = ia + 1 - jb = 0 - do j=nC(2)+1,nO(2) - do b=nO(2)+1,nBas-nR(2) - jb = jb + 1 - - chi = 0d0 - do kc=1,nSt - eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps & - + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps - enddo + chi = 0d0 + do kc=1,nSsc + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,j,kc,1)*rho(a,b,kc,2)*Omega(kc)/eps + enddo - A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,b,j,a) + 2d0*lambda*chi + A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_abab(i,b,j,a) + 4d0*lambda*chi - enddo - enddo - enddo - enddo + end do + end do + end do + end do + + ! baba block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + chi = 0d0 + do kc=1,nSsc + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,j,kc,2)*rho(a,b,kc,1)*Omega(kc)/eps + enddo + + A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_abab(b,i,a,j) + 4d0*lambda*chi + + end do + end do + end do + end do + + end if end subroutine unrestricted_Bethe_Salpeter_A_matrix diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 index 20dd14e..e954814 100644 --- a/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 +++ b/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 @@ -1,6 +1,5 @@ -subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & - ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, & - Omega,rho,B_lr) +subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,B_lr) ! Compute the extra term for Bethe-Salpeter equation for linear response @@ -9,6 +8,7 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt ! Input variables + integer,intent(in) :: ispin integer,intent(in) :: nBas integer,intent(in) :: nC(nspin) integer,intent(in) :: nO(nspin) @@ -17,14 +17,15 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt integer,intent(in) :: nSa integer,intent(in) :: nSb integer,intent(in) :: nSt + integer,intent(in) :: nSsc double precision,intent(in) :: eta double precision,intent(in) :: lambda double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega(nSt) - double precision,intent(in) :: rho(nBas,nBas,nSt,nspin) + double precision,intent(in) :: Omega(nSsc) + double precision,intent(in) :: rho(nBas,nBas,nSsc,nspin) ! Local variables @@ -36,104 +37,118 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt double precision,intent(out) :: B_lr(nSt,nSt) - ! alpha-alpha block +!--------------------------------------------------! +! Build BSE matrix for spin-conserving transitions ! +!--------------------------------------------------! - ia = 0 - do i=nC(1)+1,nO(1) - do a=nO(1)+1,nBas-nR(1) - ia = ia + 1 - jb = 0 - do j=nC(1)+1,nO(1) - do b=nO(1)+1,nBas-nR(1) - jb = jb + 1 + if(ispin == 1) then + + ! aaaa block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + chi = 0d0 + do kc=1,nSsc + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps + enddo + + B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI_aaaa(i,j,b,a) + 4d0*lambda*chi - chi = 0d0 - do kc=1,nSt - eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps & - + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps enddo - - B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI_aaaa(i,j,b,a) + 2d0*lambda*chi - enddo enddo enddo - enddo - - ! alpha-beta block - - ia = 0 - do i=nC(1)+1,nO(1) - do a=nO(1)+1,nBas-nR(1) - ia = ia + 1 - jb = 0 - do j=nC(2)+1,nO(2) - do b=nO(2)+1,nBas-nR(2) - jb = jb + 1 - chi = 0d0 - do kc=1,nSt - eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps & - + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps + + ! bbbb block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + chi = 0d0 + do kc=1,nSsc + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps + enddo + + B_lr(nSa+ia,nSa+jb) = B_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,j,b,a) + 4d0*lambda*chi + enddo - - B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_aabb(i,j,b,a) + 2d0*lambda*chi - enddo enddo enddo - enddo - ! beta-alpha block + end if - ia = 0 - do i=nC(2)+1,nO(2) - do a=nO(2)+1,nBas-nR(2) - ia = ia + 1 - jb = 0 - do j=nC(1)+1,nO(1) - do b=nO(1)+1,nBas-nR(1) - jb = jb + 1 - - chi = 0d0 - do kc=1,nSt - eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps & - + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps - enddo - B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_aabb(j,i,a,b) + 2d0*lambda*chi +!--------------------------------------------! +! Build BSE matrix for spin-flip transitions ! +!--------------------------------------------! - enddo - enddo - enddo - enddo + if(ispin == 2) then - ! beta-beta block + ! abba block - ia = 0 - do i=nC(2)+1,nO(2) - do a=nO(2)+1,nBas-nR(2) - ia = ia + 1 - jb = 0 - do j=nC(2)+1,nO(2) - do b=nO(2)+1,nBas-nR(2) - jb = jb + 1 - - chi = 0d0 - do kc=1,nSt - eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps & - + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps - enddo + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 - B_lr(nSa+ia,nSa+jb) = B_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,j,b,a) + 2d0*lambda*chi + chi = 0d0 + do kc=1,nSsc + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,b,kc,1)*rho(a,j,kc,2)*Omega(kc)/eps + enddo - enddo - enddo - enddo - enddo + B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_abab(i,a,b,j) + 4d0*lambda*chi + + end do + end do + end do + end do + + ! baab block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + chi = 0d0 + do kc=1,nSsc + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,b,kc,2)*rho(a,j,kc,1)*Omega(kc)/eps + enddo + + B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_abab(b,j,i,a) + 4d0*lambda*chi + + end do + end do + end do + end do + + end if end subroutine unrestricted_Bethe_Salpeter_B_matrix diff --git a/src/QuAcK/unrestricted_linear_response.f90 b/src/QuAcK/unrestricted_linear_response.f90 index 29a6dfe..b61cb38 100644 --- a/src/QuAcK/unrestricted_linear_response.f90 +++ b/src/QuAcK/unrestricted_linear_response.f90 @@ -1,5 +1,5 @@ -subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & - e,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho,EcRPA,Omega,XpY,XmY) +subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, & + e,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_RPA,rho_RPA,EcRPA,Omega,XpY,XmY) ! Compute linear response for unrestricted formalism @@ -21,13 +21,16 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, integer,intent(in) :: nSa integer,intent(in) :: nSb integer,intent(in) :: nSt + integer,intent(in) :: nSsc double precision,intent(in) :: lambda double precision,intent(in) :: e(nBas,nspin) - double precision,intent(in) :: rho(nBas,nBas,nSt,nspin) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas) + + double precision,intent(in) :: Omega_RPA(nSsc) + double precision,intent(in) :: rho_RPA(nBas,nBas,nSsc,nspin) ! Local variables @@ -58,8 +61,8 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A) if(BSE) & - call unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & - ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,A) + call unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_RPA,rho_RPA,A) ! Tamm-Dancoff approximation @@ -70,8 +73,8 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,B) if(BSE) & - call unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & - ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,B) + call unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_RPA,rho_RPA,B) end if diff --git a/src/QuAcK/unrestricted_linear_response_A_matrix.f90 b/src/QuAcK/unrestricted_linear_response_A_matrix.f90 index ded4c78..6c282e5 100644 --- a/src/QuAcK/unrestricted_linear_response_A_matrix.f90 +++ b/src/QuAcK/unrestricted_linear_response_A_matrix.f90 @@ -143,7 +143,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa jb = jb + 1 A_lr(ia,jb) = (e(a,2) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & - + lambda*ERI_abab(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_abab(i,b,j,a) + - (1d0 - delta_dRPA)*lambda*ERI_abab(i,b,j,a) end do end do @@ -162,7 +162,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa jb = jb + 1 A_lr(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & - + lambda*ERI_abab(b,i,j,a) - (1d0 - delta_dRPA)*lambda*ERI_abab(b,i,a,j) + - (1d0 - delta_dRPA)*lambda*ERI_abab(b,i,a,j) end do end do diff --git a/src/QuAcK/unrestricted_linear_response_B_matrix.f90 b/src/QuAcK/unrestricted_linear_response_B_matrix.f90 index 5ac615d..5cf001c 100644 --- a/src/QuAcK/unrestricted_linear_response_B_matrix.f90 +++ b/src/QuAcK/unrestricted_linear_response_B_matrix.f90 @@ -128,7 +128,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa B_lr(:,:) = 0d0 - ! abab block + ! abba block ia = 0 do i=nC(1)+1,nO(1) @@ -136,28 +136,28 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa ia = ia + 1 jb = 0 do j=nC(2)+1,nO(2) - do b=nO(2)+1,nBas-nR(2) + do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - B_lr(ia,jb) = lambda*ERI_abab(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_abab(i,j,b,a) + B_lr(ia,nSa+jb) = - (1d0 - delta_dRPA)*lambda*ERI_abab(i,a,b,j) end do end do end do end do - ! bbbb block + ! baab block ia = 0 do i=nC(2)+1,nO(2) do a=nO(1)+1,nBas-nR(1) ia = ia + 1 jb = 0 - do j=nC(2)+1,nO(2) - do b=nO(1)+1,nBas-nR(1) + do j=nC(1)+1,nO(1) + do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - B_lr(nSa+ia,nSa+jb) = lambda*ERI_abab(j,i,b,a) - (1d0 - delta_dRPA)*lambda*ERI_abab(j,i,a,b) + B_lr(nSa+ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI_abab(b,j,i,a) end do end do