From 282cbcb517ddfe6fc04401dc407839a54c2aab0c Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 22 Sep 2020 15:32:26 +0200 Subject: [PATCH] UG0W0 OK --- input/options | 2 +- src/QuAcK/G0W0.f90 | 3 +- src/QuAcK/QuAcK.f90 | 2 +- src/QuAcK/UG0W0.f90 | 74 +++++----- src/QuAcK/print_G0W0.f90 | 10 +- src/QuAcK/self_energy_correlation_diag.f90 | 4 +- src/QuAcK/unrestricted_excitation_density.f90 | 57 ++++---- src/QuAcK/unrestricted_linear_response.f90 | 113 ++++++++++++++++ .../unrestricted_linear_response_A_matrix.f90 | 126 ++++++++++++++++++ .../unrestricted_linear_response_B_matrix.f90 | 123 +++++++++++++++++ ...estricted_self_energy_correlation_diag.f90 | 12 +- 11 files changed, 437 insertions(+), 89 deletions(-) create mode 100644 src/QuAcK/unrestricted_linear_response.f90 create mode 100644 src/QuAcK/unrestricted_linear_response_A_matrix.f90 create mode 100644 src/QuAcK/unrestricted_linear_response_B_matrix.f90 diff --git a/input/options b/input/options index daaa7f4..52b183f 100644 --- a/input/options +++ b/input/options @@ -13,6 +13,6 @@ # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn - T T F F + F T 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 bd76850..1b32919 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -1,6 +1,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linearize,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eGW) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI,PHF,cHF,eHF,eGW) ! Perform G0W0 calculation @@ -33,7 +33,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: H(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index d7f4fbc..c9cd58e 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -4,7 +4,7 @@ program QuAcK include 'parameters.h' logical :: doSph - logical :: unrestricted + logical :: unrestricted = .false. logical :: doRHF,doUHF,doMOM logical :: doMP2,doMP3,doMP2F12 logical :: doCCD,doCCSD,doCCSDT diff --git a/src/QuAcK/UG0W0.f90 b/src/QuAcK/UG0W0.f90 index b71122b..cdce17f 100644 --- a/src/QuAcK/UG0W0.f90 +++ b/src/QuAcK/UG0W0.f90 @@ -51,16 +51,14 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev integer :: nSa integer :: nSb integer :: nSt - double precision :: EcRPA(nspin) - double precision :: EcBSE(nspin) + double precision :: EcRPA + double precision :: EcBSE double precision :: EcAC(nspin) double precision,allocatable :: SigC(:,:) double precision,allocatable :: Z(:,:) double precision,allocatable :: Omega(:) - double precision,allocatable :: XpY_a(:,:) - double precision,allocatable :: XpY_b(:,:) - double precision,allocatable :: XmY_a(:,:) - double precision,allocatable :: XmY_b(:,:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) double precision,allocatable :: rho(:,:,:,:) double precision,allocatable :: eGWlin(:,:) @@ -80,7 +78,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! Initialization - EcRPA(:) = 0d0 + EcRPA = 0d0 ! COHSEX approximation @@ -103,53 +101,45 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev nSb = nS(2) nSt = nSa + nSb - allocate(SigC(nBas,nspin),Z(nBas,nspin),Omega(nSt),XpY_a(nSa,nSa),XpY_b(nSb,nSb),XmY_a(nSa,nSa),XmY_b(nSb,nSb), & + allocate(SigC(nBas,nspin),Z(nBas,nspin),Omega(nSt),XpY(nSt,nSt),XmY(nSt,nSt), & rho(nBas,nBas,nSt,nspin),eGWlin(nBas,nspin)) -! Compute linear response +!-------------------! +! Compute screening ! +!-------------------! -!---------------------------------------------- -! alpha-alpha block -!---------------------------------------------- +! Spin-conserving transition - ispin = 1 - iblock = 3 + ispin = 1 - call linear_response(iblock,.true.,TDA_W,.false.,eta,nBas,nC(ispin),nO(ispin),nV(ispin),nR(ispin),nSa,1d0, & - eHF(:,ispin),ERI_aa,rho(:,:,1:nSa,ispin),EcRPA(ispin),Omega(1:nSa),XpY_a,XmY_a) + 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) - if(print_W) call print_excitation('RPA@HF alpha',iblock,nSa,Omega(1:nSa)) + if(print_W) call print_excitation('RPA@UHF',3,nSt,Omega) -!---------------------------------------------- -! alpha-beta block -!---------------------------------------------- - - ispin = 2 - iblock = 3 - - call linear_response(iblock,.true.,TDA_W,.false.,eta,nBas,nC(ispin),nO(ispin),nV(ispin),nR(ispin),nSb,1d0, & - eHF(:,ispin),ERI_bb,rho(:,:,nSa+1:nSt,ispin),EcRPA(ispin),Omega(nSa+1:nSt),XpY_b,XmY_b) - - if(print_W) call print_excitation('RPA@HF beta ',iblock,nSb,Omega(nSa+1:nSt)) - -!---------------------------------------------- -! Excitation densities for alpha self-energy -!---------------------------------------------- +!----------------------! +! Excitation densities ! +!----------------------! - call unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ab,ERI_bb,XpY_a,XpY_b,rho) + call unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ab,ERI_bb,XpY,rho) -!---------------------- -! Compute self-energy -!---------------------- +!---------------------! +! Compute self-energy ! +!---------------------! call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,eHF,Omega,rho,SigC) -! Compute renormalization factor +!--------------------------------! +! Compute renormalization factor ! +!--------------------------------! ! call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & ! Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:)) -! Solve the quasi-particle equation +!-----------------------------------! +! Solve the quasi-particle equation ! +!-----------------------------------! + Z(:,:) = 1d0 eGWlin(:,:) = eHF(:,:) + Z(:,:)*SigC(:,:) @@ -174,7 +164,7 @@ 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(ispin)) + call print_G0W0(nBas,nO(ispin),eHF(:,ispin),ENuc,EUHF,SigC(:,ispin),Z(:,ispin),eGW(:,ispin),EcRPA) end do ! Compute the RPA correlation energy @@ -184,10 +174,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (singlet) =',EcRPA(1) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (triplet) =',EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy =',EcRPA(1) + EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy =',EcRPA + write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + EUHF + EcRPA write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/QuAcK/print_G0W0.f90 b/src/QuAcK/print_G0W0.f90 index 95d1051..d1ee3f4 100644 --- a/src/QuAcK/print_G0W0.f90 +++ b/src/QuAcK/print_G0W0.f90 @@ -40,11 +40,11 @@ subroutine print_G0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM) write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',eGW(LUMO)*HaToeV write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A30,F15.6)') 'RPA@G0W0 total energy =',ENuc + EHF + EcRPA -! write(*,'(2X,A30,F15.6)') 'RPA@G0W0 correlation energy =',EcRPA -! write(*,'(2X,A30,F15.6)') 'GM@G0W0 total energy =',ENuc + EHF + EcGM -! write(*,'(2X,A30,F15.6)') 'GM@G0W0 correlation energy =',EcGM -! write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'RPA@HF total energy =',ENuc + EHF + EcRPA + write(*,'(2X,A30,F15.6)') 'RPA@HF correlation energy =',EcRPA + write(*,'(2X,A30,F15.6)') 'GM@G0W0 total energy =',ENuc + EHF + EcGM + write(*,'(2X,A30,F15.6)') 'GM@G0W0 correlation energy =',EcGM + write(*,*)'-------------------------------------------------------------------------------' write(*,*) end subroutine print_G0W0 diff --git a/src/QuAcK/self_energy_correlation_diag.f90 b/src/QuAcK/self_energy_correlation_diag.f90 index 03afb7a..441806d 100644 --- a/src/QuAcK/self_energy_correlation_diag.f90 +++ b/src/QuAcK/self_energy_correlation_diag.f90 @@ -66,7 +66,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O EcGM = 0d0 do i=nC+1,nO - EcGM = EcGM + 0.5d0*SigC(i) + EcGM = EcGM - SigC(i) end do !----------------------------- @@ -143,7 +143,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O do a=nO+1,nBas-nR do jb=1,nS eps = e(a) - e(i) + Omega(jb) - EcGM = EcGM - 2d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2) + EcGM = EcGM - 4d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2) end do end do end do diff --git a/src/QuAcK/unrestricted_excitation_density.f90 b/src/QuAcK/unrestricted_excitation_density.f90 index a705278..5ea602b 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_a,XpY_b,rho) +subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ab,ERI_bb,XpY,rho) ! Compute excitation densities for unrestricted reference @@ -17,8 +17,7 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ 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) :: XpY_a(nSa,nSa) - double precision,intent(in) :: XpY_b(nSb,nSb) + double precision,intent(in) :: XpY(nSt,nSt) ! Local variables @@ -32,34 +31,34 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ rho(:,:,:,:) = 0d0 -!------------- -! alpha block -!------------- +!-------------! +! alpha block ! +!-------------! do p=nC(1)+1,nBas-nR(1) do q=nC(1)+1,nBas-nR(1) ! Same-spin contribution - do ia=1,nSa + do ia=1,nSt jb = 0 do j=nC(1)+1,nO(1) 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_a(ia,jb) + rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aa(p,j,q,b)*XpY(ia,jb) enddo enddo enddo ! Opposite-spin contribution - do ia=1,nSb - jb = 0 + do ia=1,nSt + jb = nSa do j=nC(2)+1,nO(2) do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - rho(p,q,nSa+ia,1) = rho(p,q,nSa+ia,1) + ERI_ab(p,j,q,b)*XpY_b(ia,jb) + rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_ab(p,j,q,b)*XpY(ia,jb) enddo enddo @@ -68,34 +67,34 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ enddo enddo -!------------ -! Beta block -!------------ +!------------! +! Beta block ! +!------------! do p=nC(2)+1,nBas-nR(2) do q=nC(2)+1,nBas-nR(2) - ! Same-spin contribution - do ia=1,nSb - jb = 0 - do j=nC(2)+1,nO(2) - 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_b(ia,jb) - - enddo - enddo - enddo - ! Opposite-spin contribution - do ia=1,nSa + do ia=1,nSt jb = 0 do j=nC(1)+1,nO(1) do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - rho(p,q,nSb+ia,2) = rho(p,q,nSb+ia,2) + ERI_ab(j,p,b,q)*XpY_a(ia,jb) + rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_ab(j,p,b,q)*XpY(ia,jb) + + enddo + enddo + enddo + + ! Same-spin contribution + do ia=1,nSt + jb = nSa + do j=nC(2)+1,nO(2) + 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) enddo enddo diff --git a/src/QuAcK/unrestricted_linear_response.f90 b/src/QuAcK/unrestricted_linear_response.f90 new file mode 100644 index 0000000..0a547da --- /dev/null +++ b/src/QuAcK/unrestricted_linear_response.f90 @@ -0,0 +1,113 @@ +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) + +! Compute linear response for unrestricted formalism + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + logical,intent(in) :: dRPA + logical,intent(in) :: TDA + logical,intent(in) :: BSE + 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) :: 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) + +! Local variables + + integer :: ia + double precision :: trace_matrix + double precision,allocatable :: A(:,:) + double precision,allocatable :: B(:,:) + double precision,allocatable :: ApB(:,:) + double precision,allocatable :: AmB(:,:) + double precision,allocatable :: AmBSq(:,:) + double precision,allocatable :: AmBIv(:,:) + double precision,allocatable :: Z(:,:) + +! Output variables + + double precision,intent(out) :: EcRPA + double precision,intent(out) :: Omega(nSt) + double precision,intent(out) :: XpY(nSt,nSt) + double precision,intent(out) :: XmY(nSt,nSt) + +! Memory allocation + + allocate(A(nSt,nSt),B(nSt,nSt),ApB(nSt,nSt),AmB(nSt,nSt),AmBSq(nSt,nSt),AmBIv(nSt,nSt),Z(nSt,nSt)) + +! 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) + +! if(BSE) call Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,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) + +! if(BSE) call Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B) + + end if + +! Build A + B and A - B matrices + + ApB = A + B + AmB = A - B + +! Diagonalize linear response matrix + + call diagonalize_matrix(nSt,AmB,Omega) + + if(minval(Omega) < 0d0) & + call print_warning('You may have instabilities in linear response: A-B is not positive definite!!') + + do ia=1,nSt + if(Omega(ia) < 0d0) Omega(ia) = 0d0 + end do + + call ADAt(nSt,AmB,1d0*sqrt(Omega),AmBSq) + call ADAt(nSt,AmB,1d0/sqrt(Omega),AmBIv) + + Z = matmul(AmBSq,matmul(ApB,AmBSq)) + + call diagonalize_matrix(nSt,Z,Omega) + + if(minval(Omega) < 0d0) & + call print_warning('You may have instabilities in linear response: negative excitations!!') + + do ia=1,nSt + if(Omega(ia) < 0d0) Omega(ia) = 0d0 + end do + + Omega = sqrt(Omega) + + XpY = matmul(transpose(Z),AmBSq) + call DA(nSt,1d0/sqrt(Omega),XpY) + + XmY = matmul(transpose(Z),AmBIv) + call DA(nSt,1d0*sqrt(Omega),XmY) + +! Compute the RPA correlation energy + + EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nSt,A)) + +end subroutine unrestricted_linear_response diff --git a/src/QuAcK/unrestricted_linear_response_A_matrix.f90 b/src/QuAcK/unrestricted_linear_response_A_matrix.f90 new file mode 100644 index 0000000..60d2d87 --- /dev/null +++ b/src/QuAcK/unrestricted_linear_response_A_matrix.f90 @@ -0,0 +1,126 @@ +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) + +! Compute linear response + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dRPA + integer,intent(in) :: ispin + 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) :: 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) + +! Local variables + + double precision :: delta_dRPA + double precision,external :: Kronecker_delta + + integer :: i,j,a,b,ia,jb + +! Output variables + + double precision,intent(out) :: A_lr(nSt,nSt) + +! Direct RPA + + delta_dRPA = 0d0 + if(dRPA) delta_dRPA = 1d0 + +!----------------------------------------------- +! Build A matrix for spin-conserving transitions +!----------------------------------------------- + + if(ispin == 1) then + + ! alpha-alpha 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 + + 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) + + end do + end do + end do + end do + + ! 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 + + A_lr(ia,nSa+jb) = lambda*ERI_ab(i,b,a,j) + + end do + end do + end do + end do + + ! beta-alpha 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(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + A_lr(nSa+ia,jb) = lambda*ERI_ab(b,i,j,a) + + end do + end do + end do + end do + + ! beta-beta 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 + + 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) + + end do + end do + end do + end do + + 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 new file mode 100644 index 0000000..22213df --- /dev/null +++ b/src/QuAcK/unrestricted_linear_response_B_matrix.f90 @@ -0,0 +1,123 @@ +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) + +! Compute linear response + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dRPA + integer,intent(in) :: ispin + 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) :: 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) + +! Local variables + + double precision :: delta_dRPA + double precision,external :: Kronecker_delta + + integer :: i,j,a,b,ia,jb + +! Output variables + + double precision,intent(out) :: B_lr(nSt,nSt) + +! Direct RPA + + delta_dRPA = 0d0 + if(dRPA) delta_dRPA = 1d0 + +!----------------------------------------------- +! Build A matrix for spin-conserving transitions +!----------------------------------------------- + + if(ispin == 1) then + + ! alpha-alpha 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 + + B_lr(ia,jb) = lambda*ERI_aa(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_aa(i,j,b,a) + + end do + end do + end do + end do + + ! 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 + + B_lr(ia,nSa+jb) = lambda*ERI_ab(i,j,a,b) + + end do + end do + end do + end do + + ! beta-alpha 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(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + B_lr(nSa+ia,jb) = lambda*ERI_ab(j,i,b,a) + + end do + end do + end do + end do + + ! beta-beta 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 + + B_lr(nSa+ia,nSa+jb) = lambda*ERI_bb(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_bb(i,j,b,a) + + end do + end do + end do + end do + + end if + + +end subroutine unrestricted_linear_response_B_matrix diff --git a/src/QuAcK/unrestricted_self_energy_correlation_diag.f90 b/src/QuAcK/unrestricted_self_energy_correlation_diag.f90 index 3c91767..123ec23 100644 --- a/src/QuAcK/unrestricted_self_energy_correlation_diag.f90 +++ b/src/QuAcK/unrestricted_self_energy_correlation_diag.f90 @@ -33,9 +33,9 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nS SigC(:,:) = 0d0 -!-------------- -! Spin-up part -!-------------- +!--------------! +! Spin-up part ! +!--------------! ! Occupied part of the correlation self-energy @@ -59,9 +59,9 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nS end do end do -!---------------- -! Spin-down part -!---------------- +!----------------! +! Spin-down part ! +!----------------! ! Occupied part of the correlation self-energy