From 2a6f83dbf902c47904fc05fd3d7ebd4a7d5e223b Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 22 Sep 2020 22:13:51 +0200 Subject: [PATCH] UBSE --- input/options | 4 +- src/QuAcK/G0W0.f90 | 24 +++-- src/QuAcK/QuAcK.f90 | 22 ++--- src/QuAcK/UG0W0.f90 | 59 ++++++------ src/QuAcK/print_UG0W0.f90 | 57 ++++++++---- src/QuAcK/print_excitation.f90 | 10 +- src/QuAcK/renormalization_factor.f90 | 74 +++++++-------- src/QuAcK/unrestricted_excitation_density.f90 | 16 ++-- src/QuAcK/unrestricted_linear_response.f90 | 18 ++-- .../unrestricted_linear_response_A_matrix.f90 | 26 ++++-- .../unrestricted_linear_response_B_matrix.f90 | 26 ++++-- .../unrestricted_renormalization_factor.f90 | 93 +++++++++++++++++++ 12 files changed, 287 insertions(+), 142 deletions(-) create mode 100644 src/QuAcK/unrestricted_renormalization_factor.f90 diff --git a/input/options b/input/options index 52b183f..a5abf91 100644 --- a/input/options +++ b/input/options @@ -9,10 +9,10 @@ # 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 T 0.0 F F F F F + 256 0.00001 T 5 T 0.001 F F F F F # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn - F T F F + T F F F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index 1b32919..b6ffcbc 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -71,23 +71,31 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, ! SOSEX correction - if(SOSEX) write(*,*) 'SOSEX correction activated!' - write(*,*) + if(SOSEX) then + write(*,*) 'SOSEX correction activated!' + write(*,*) + end if ! COHSEX approximation - if(COHSEX) write(*,*) 'COHSEX approximation activated!' - write(*,*) + if(COHSEX) then + write(*,*) 'COHSEX approximation activated!' + write(*,*) + end if ! TDA for W - if(TDA_W) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' - write(*,*) + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) + end if ! TDA - if(TDA) write(*,*) 'Tamm-Dancoff approximation activated!' - write(*,*) + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if ! Spin manifold diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index c9cd58e..0bb56a2 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -57,9 +57,9 @@ program QuAcK double precision,allocatable :: ERI_MO(:,:,:,:) integer :: bra integer :: ket - double precision,allocatable :: ERI_MO_aa(:,:,:,:) - double precision,allocatable :: ERI_MO_ab(:,:,:,:) - double precision,allocatable :: ERI_MO_bb(:,:,:,:) + double precision,allocatable :: ERI_MO_aaaa(:,:,:,:) + double precision,allocatable :: ERI_MO_aabb(:,:,:,:) + double precision,allocatable :: ERI_MO_bbbb(:,:,:,:) double precision,allocatable :: ERI_ERF_AO(:,:,:,:) double precision,allocatable :: ERI_ERF_MO(:,:,:,:) double precision,allocatable :: F12(:,:,:,:),Yuk(:,:,:,:),FC(:,:,:,:,:,:) @@ -330,25 +330,25 @@ program QuAcK ! Memory allocation - allocate(ERI_MO_aa(nBas,nBas,nBas,nBas),ERI_MO_ab(nBas,nBas,nBas,nBas),ERI_MO_bb(nBas,nBas,nBas,nBas)) + allocate(ERI_MO_aaaa(nBas,nBas,nBas,nBas),ERI_MO_aabb(nBas,nBas,nBas,nBas),ERI_MO_bbbb(nBas,nBas,nBas,nBas)) ! 4-index transform for (aa|aa) block bra = 1 ket = 1 - call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_aa) + call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_aaaa) - ! 4-index transform for (bb|bb) block + ! 4-index transform for (aa|bb) block bra = 1 ket = 2 - call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_ab) + call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_aabb) - ! 4-index transform for (aa|bb) block + ! 4-index transform for (bb|bb) block bra = 2 ket = 2 - call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_bb) + call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_bbbb) else @@ -382,7 +382,7 @@ program QuAcK if(unrestricted) then - call UMP2(nBas,nC,nO,nV,nR,ERI_MO_aa,ERI_MO_ab,ERI_MO_bb,ENuc,EUHF,eHF,EcMP2) + call UMP2(nBas,nC,nO,nV,nR,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ENuc,EUHF,eHF,EcMP2) else @@ -749,7 +749,7 @@ program QuAcK call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, & singlet_manifold,triplet_manifold,linGW,eta_GW,nBas,nC,nO,nV,nR,nS, & - ENuc,EUHF,Hc,ERI_MO_aa,ERI_MO_ab,ERI_MO_bb,PHF,cHF,eHF,eG0W0) + ENuc,EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,PHF,cHF,eHF,eG0W0) else call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & diff --git a/src/QuAcK/UG0W0.f90 b/src/QuAcK/UG0W0.f90 index cdce17f..0e0b38e 100644 --- a/src/QuAcK/UG0W0.f90 +++ b/src/QuAcK/UG0W0.f90 @@ -1,6 +1,6 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, & - singlet_manifold,triplet_manifold,linearize,eta,nBas,nC,nO,nV,nR,nS, & - ENuc,EUHF,Hc,ERI_aa,ERI_ab,ERI_bb,PHF,cHF,eHF,eGW) + spin_conserved,spin_flip,linearize,eta,nBas,nC,nO,nV,nR,nS, & + ENuc,EUHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,PHF,cHF,eHF,eGW) ! Perform unrestricted G0W0 calculation @@ -20,8 +20,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev logical,intent(in) :: dBSE logical,intent(in) :: dTDA logical,intent(in) :: evDyn - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip logical,intent(in) :: linearize double precision,intent(in) :: eta @@ -37,9 +37,9 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev double precision,intent(in) :: cHF(nBas,nBas,nspin) double precision,intent(in) :: PHF(nBas,nBas,nspin) double precision,intent(in) :: Hc(nBas,nBas,nspin) - double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas) + 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) ! Local variables @@ -82,18 +82,24 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! COHSEX approximation - if(COHSEX) write(*,*) 'COHSEX approximation activated!' - write(*,*) + if(COHSEX) then + write(*,*) 'COHSEX approximation activated!' + write(*,*) + end if ! TDA for W - if(TDA_W) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' - write(*,*) + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) + end if ! TDA - if(TDA) write(*,*) 'Tamm-Dancoff approximation activated!' - write(*,*) + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if ! Memory allocation @@ -113,15 +119,15 @@ 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,nSa,nSb,nSt,1d0, & - eHF,ERI_aa,ERI_ab,ERI_bb,rho(:,:,:,ispin),EcRPA,Omega,XpY,XmY) + eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,rho,EcRPA,Omega,XpY,XmY) - if(print_W) call print_excitation('RPA@UHF',3,nSt,Omega) + if(print_W) call print_excitation('RPA@UHF',5,nSt,Omega) !----------------------! ! Excitation densities ! !----------------------! - call unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ab,ERI_bb,XpY,rho) + call unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho) !---------------------! ! Compute self-energy ! @@ -133,14 +139,12 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! Compute renormalization factor ! !--------------------------------! -! call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & -! Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:)) + call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,eHF,Omega,rho,Z) !-----------------------------------! ! Solve the quasi-particle equation ! !-----------------------------------! - Z(:,:) = 1d0 eGWlin(:,:) = eHF(:,:) + Z(:,:)*SigC(:,:) if(linearize) then @@ -163,14 +167,12 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! Dump results - do ispin=1,nspin - call print_G0W0(nBas,nO(ispin),eHF(:,ispin),ENuc,EUHF,SigC(:,ispin),Z(:,ispin),eGW(:,ispin),EcRPA) - end do + call print_UG0W0(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA) ! Compute the RPA correlation energy -! call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & -! rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho,EcRPA,Omega,XpY,XmY) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -181,10 +183,11 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! Perform BSE calculation -! if(BSE) then + if(BSE) then -! call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & -! nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) + call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, & + nBas,nC,nO,nV,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb, & + eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) ! if(exchange_kernel) then ! @@ -239,6 +242,6 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! end if -! end if + end if end subroutine UG0W0 diff --git a/src/QuAcK/print_UG0W0.f90 b/src/QuAcK/print_UG0W0.f90 index b73e886..3cd21cd 100644 --- a/src/QuAcK/print_UG0W0.f90 +++ b/src/QuAcK/print_UG0W0.f90 @@ -1,44 +1,61 @@ -subroutine print_UG0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA) +subroutine print_UG0W0(nBas,nO,e,ENuc,EHF,SigC,Z,eGW,EcRPA) ! Print one-electron energies and other stuff for G0W0 implicit none include 'parameters.h' - integer,intent(in) :: nBas,nO + integer,intent(in) :: nBas + integer,intent(in) :: nO(nspin) double precision,intent(in) :: ENuc double precision,intent(in) :: EHF double precision,intent(in) :: EcRPA - double precision,intent(in) :: e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas) + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: SigC(nBas,nspin) + double precision,intent(in) :: Z(nBas,nspin) + double precision,intent(in) :: eGW(nBas,nspin) - integer :: x,HOMO,LUMO + integer :: p + double precision :: HOMO + double precision :: LUMO double precision :: Gap ! HOMO and LUMO - HOMO = nO - LUMO = HOMO + 1 - Gap = eGW(LUMO)-eGW(HOMO) + HOMO = max(eGW(nO(1),1),eGW(nO(2),2)) + LUMO = min(eGW(nO(1)+1,1),eGW(nO(2)+1,2)) + Gap = LUMO - HOMO ! Dump results - write(*,*)'-------------------------------------------------------------------------------' - write(*,*)' One-shot G0W0 calculation' - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & - '|','#','|','e_HF (eV)','|','Sigma_c (eV)','|','Z','|','e_QP (eV)','|' - write(*,*)'-------------------------------------------------------------------------------' + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,*)' Unrestricted one-shot G0W0 calculation (eV)' + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(A1,A3,A1,2A15,A1,2A15,A1,2A15,A1,2A15,A1)') & + '|','#','|','e_HF up','e_HF dw','|','Sig_c up','Sig_c dw','|', & + 'Z up','Z dw','|','e_QP up','e_QP dw','|' + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' - do x=1,nBas - write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & - '|',x,'|',e(x)*HaToeV,'|',SigmaC(x)*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|' + do p=1,nBas + write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & + '|',p,'|',e(p,1)*HaToeV,e(p,2)*HaToeV,'|',SigC(p,1)*HaToeV,SigC(p,2)*HaToeV,'|', & + Z(p,1),Z(p,2),'|',eGW(p,1)*HaToeV,eGW(p,2)*HaToeV,'|' enddo - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A30,F15.6)') 'G0W0 HOMO energy (eV):',eGW(HOMO)*HaToeV - write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',eGW(LUMO)*HaToeV + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'G0W0 HOMO energy (eV):',HOMO*HaToeV + write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',LUMO*HaToeV write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV - write(*,*)'-------------------------------------------------------------------------------' + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'RPA@HF total energy =',ENuc + EHF + EcRPA + write(*,'(2X,A30,F15.6)') 'RPA@HF correlation energy =',EcRPA + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' write(*,*) end subroutine print_UG0W0 diff --git a/src/QuAcK/print_excitation.f90 b/src/QuAcK/print_excitation.f90 index ea8fa8b..e5ce17b 100644 --- a/src/QuAcK/print_excitation.f90 +++ b/src/QuAcK/print_excitation.f90 @@ -13,18 +13,20 @@ subroutine print_excitation(method,ispin,nS,Omega) ! Local variables - character*7 :: spin_manifold + character*14 :: spin_manifold integer,parameter :: maxS = 32 integer :: ia if(ispin == 1) spin_manifold = 'singlet' if(ispin == 2) spin_manifold = 'triplet' - if(ispin == 3) spin_manifold = 'alp-bet' - if(ispin == 4) spin_manifold = 'alp-alp' + if(ispin == 3) spin_manifold = 'alpha-beta' + if(ispin == 4) spin_manifold = 'alpha-alpha' + if(ispin == 5) spin_manifold = 'spin-conserved' + if(ispin == 6) spin_manifold = 'spin-flip' write(*,*) write(*,*)'-------------------------------------------------------------' - write(*,'(1X,A1,1X,A14,A14,A7,A9,A15)')'|',method,' calculation: ',spin_manifold,' manifold',' |' + write(*,'(1X,A1,1X,A14,A14,A14,A9,A13)')'|',method,' calculation: ',spin_manifold,' manifold',' |' write(*,*)'-------------------------------------------------------------' write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') & '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' diff --git a/src/QuAcK/renormalization_factor.f90 b/src/QuAcK/renormalization_factor.f90 index 070c1e5..15e3bf2 100644 --- a/src/QuAcK/renormalization_factor.f90 +++ b/src/QuAcK/renormalization_factor.f90 @@ -35,42 +35,10 @@ subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r Z(:) = 1d0 return - - end if - -! Occupied part of the correlation self-energy - - do x=nC+1,nBas-nR - do i=nC+1,nO - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - eps = e(x) - e(i) + Omega(jb) - Z(x) = Z(x) - 2d0*rho(x,i,jb)**2*(eps/(eps**2 + eta**2))**2 - end do - end do - end do - end do - -! Virtual part of the correlation self-energy - - do x=nC+1,nBas-nR - do a=nO+1,nBas-nR - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - eps = e(x) - e(a) - Omega(jb) - Z(x) = Z(x) - 2d0*rho(x,a,jb)**2*(eps/(eps**2 + eta**2))**2 - end do - end do - end do - end do - - ! SOSEX correction - - if(SOSEX) then + +! SOSEX correction + + elseif(SOSEX) then ! Occupied part of the correlation self-energy @@ -102,7 +70,39 @@ subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r end do end do - endif + else + + ! Occupied part of the correlation self-energy + + do x=nC+1,nBas-nR + do i=nC+1,nO + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + eps = e(x) - e(i) + Omega(jb) + Z(x) = Z(x) - 2d0*rho(x,i,jb)**2*(eps/(eps**2 + eta**2))**2 + end do + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do x=nC+1,nBas-nR + do a=nO+1,nBas-nR + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + eps = e(x) - e(a) - Omega(jb) + Z(x) = Z(x) - 2d0*rho(x,a,jb)**2*(eps/(eps**2 + eta**2))**2 + end do + end do + end do + end do + + end if ! Compute renormalization factor from derivative of SigC diff --git a/src/QuAcK/unrestricted_excitation_density.f90 b/src/QuAcK/unrestricted_excitation_density.f90 index 5ea602b..5822d65 100644 --- a/src/QuAcK/unrestricted_excitation_density.f90 +++ b/src/QuAcK/unrestricted_excitation_density.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ab,ERI_bb,XpY,rho) +subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho) ! Compute excitation densities for unrestricted reference @@ -14,9 +14,9 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ integer,intent(in) :: nSa integer,intent(in) :: nSb integer,intent(in) :: nSt - double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas) + 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) :: XpY(nSt,nSt) ! Local variables @@ -45,7 +45,7 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aa(p,j,q,b)*XpY(ia,jb) + rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aaaa(p,j,q,b)*XpY(ia,jb) enddo enddo @@ -58,7 +58,7 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_ab(p,j,q,b)*XpY(ia,jb) + rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aabb(p,j,q,b)*XpY(ia,jb) enddo enddo @@ -81,7 +81,7 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_ab(j,p,b,q)*XpY(ia,jb) + rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_aabb(j,p,b,q)*XpY(ia,jb) enddo enddo @@ -94,7 +94,7 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_bb(p,j,q,b)*XpY(ia,jb) + rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_bbbb(p,j,q,b)*XpY(ia,jb) enddo enddo diff --git a/src/QuAcK/unrestricted_linear_response.f90 b/src/QuAcK/unrestricted_linear_response.f90 index 0a547da..b8d8d52 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_aa,ERI_ab,ERI_bb,rho,EcRPA,Omega,XpY,XmY) + e,ERI_aaaa,ERI_aabb,ERI_bbbb,rho,EcRPA,Omega,XpY,XmY) ! Compute linear response for unrestricted formalism @@ -24,9 +24,9 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, 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_aa(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas) + 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) ! Local variables @@ -53,18 +53,20 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, ! Build A and B matrices - call unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,e,ERI_aa,ERI_ab,ERI_bb,A) + call unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,A) -! if(BSE) call Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,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,Omega,rho,A) ! Tamm-Dancoff approximation B = 0d0 if(.not. TDA) then - call unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aa,ERI_ab,ERI_bb,B) + call unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,B) -! if(BSE) call Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,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,Omega,rho,B) end if diff --git a/src/QuAcK/unrestricted_linear_response_A_matrix.f90 b/src/QuAcK/unrestricted_linear_response_A_matrix.f90 index 60d2d87..75c63b0 100644 --- a/src/QuAcK/unrestricted_linear_response_A_matrix.f90 +++ b/src/QuAcK/unrestricted_linear_response_A_matrix.f90 @@ -1,5 +1,5 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & - e,ERI_aa,ERI_ab,ERI_bb,A_lr) + e,ERI_aaaa,ERI_aabb,ERI_bbbb,A_lr) ! Compute linear response @@ -20,9 +20,9 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa integer,intent(in) :: nSt double precision,intent(in) :: lambda double precision,intent(in) :: e(nBas,nspin) - double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas) + 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) ! Local variables @@ -58,7 +58,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa jb = jb + 1 A_lr(ia,jb) = (e(a,1) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & - + lambda*ERI_aa(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_aa(i,b,j,a) + + lambda*ERI_aaaa(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_aaaa(i,b,j,a) end do end do @@ -76,7 +76,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - A_lr(ia,nSa+jb) = lambda*ERI_ab(i,b,a,j) + A_lr(ia,nSa+jb) = lambda*ERI_aabb(i,b,a,j) end do end do @@ -94,7 +94,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - A_lr(nSa+ia,jb) = lambda*ERI_ab(b,i,j,a) + A_lr(nSa+ia,jb) = lambda*ERI_aabb(b,i,j,a) end do end do @@ -113,7 +113,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,2) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & - + lambda*ERI_bb(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_bb(i,b,j,a) + + lambda*ERI_bbbb(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_bbbb(i,b,j,a) end do end do @@ -122,5 +122,15 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa end if +!----------------------------------------------- +! Build A matrix for spin-flip transitions +!----------------------------------------------- + + if(ispin == 2) then + + print*,'spin-flip transition NYI' + + end if + end subroutine unrestricted_linear_response_A_matrix diff --git a/src/QuAcK/unrestricted_linear_response_B_matrix.f90 b/src/QuAcK/unrestricted_linear_response_B_matrix.f90 index 22213df..400a498 100644 --- a/src/QuAcK/unrestricted_linear_response_B_matrix.f90 +++ b/src/QuAcK/unrestricted_linear_response_B_matrix.f90 @@ -1,5 +1,5 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & - ERI_aa,ERI_ab,ERI_bb,B_lr) + ERI_aaaa,ERI_aabb,ERI_bbbb,B_lr) ! Compute linear response @@ -19,9 +19,9 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa integer,intent(in) :: nSb integer,intent(in) :: nSt double precision,intent(in) :: lambda - double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas) + 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) ! Local variables @@ -56,7 +56,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - B_lr(ia,jb) = lambda*ERI_aa(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_aa(i,j,b,a) + B_lr(ia,jb) = lambda*ERI_aaaa(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_aaaa(i,j,b,a) end do end do @@ -74,7 +74,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - B_lr(ia,nSa+jb) = lambda*ERI_ab(i,j,a,b) + B_lr(ia,nSa+jb) = lambda*ERI_aabb(i,j,a,b) end do end do @@ -92,7 +92,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - B_lr(nSa+ia,jb) = lambda*ERI_ab(j,i,b,a) + B_lr(nSa+ia,jb) = lambda*ERI_aabb(j,i,b,a) end do end do @@ -110,7 +110,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - B_lr(nSa+ia,nSa+jb) = lambda*ERI_bb(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_bb(i,j,b,a) + B_lr(nSa+ia,nSa+jb) = lambda*ERI_bbbb(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_bbbb(i,j,b,a) end do end do @@ -119,5 +119,15 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa end if +!----------------------------------------------- +! Build A matrix for spin-flip transitions +!----------------------------------------------- + + if(ispin == 2) then + + print*,'spin-flip transition NYI' + + end if + end subroutine unrestricted_linear_response_B_matrix diff --git a/src/QuAcK/unrestricted_renormalization_factor.f90 b/src/QuAcK/unrestricted_renormalization_factor.f90 new file mode 100644 index 0000000..403060d --- /dev/null +++ b/src/QuAcK/unrestricted_renormalization_factor.f90 @@ -0,0 +1,93 @@ +subroutine unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSa,nSb,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) :: nSa + integer,intent(in) :: nSb + 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,j,a,b,p,q,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