From f7718c8ab1fbd9fdeda2187c4fc85ff2fcf17d58 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 31 Jan 2022 13:19:42 +0100 Subject: [PATCH] clean linear response and static BSE2 kernel --- input/dft | 10 +- input/methods | 6 +- input/options | 2 +- src/GF/BSE2.f90 | 42 +++-- src/GF/BSE2_A_matrix_static.f90 | 157 ++++++++++++++++++ src/GF/BSE2_B_matrix_static.f90 | 157 ++++++++++++++++++ src/GF/BSE2_dynamic_perturbation.f90 | 14 +- .../BSE2_dynamic_perturbation_iterative.f90 | 14 +- src/GT/Bethe_Salpeter_Tmatrix.f90 | 8 +- src/GT/Bethe_Salpeter_Tmatrix_so.f90 | 4 +- src/GW/Bethe_Salpeter.f90 | 22 ++- src/GW/G0W0.f90 | 8 +- src/GW/G0W0_SOSEX.f90 | 8 +- src/GW/evGW.f90 | 4 +- src/GW/qsGW.f90 | 4 +- src/LR/linear_response.f90 | 14 +- ...se_Tmatrix.f90 => linear_response_BSE.f90} | 24 ++- src/RPA/ACFDT.f90 | 33 ++-- src/RPA/ACFDT_Tmatrix.f90 | 8 +- src/RPA/ACFDT_cr.f90 | 33 ++-- src/RPA/RPA.f90 | 5 +- src/RPA/RPAx.f90 | 5 +- src/RPA/crRPA.f90 | 5 +- src/eDFT/eDFT_UKS.f90 | 72 +++++--- src/utils/level_shifting.f90 | 37 +++++ 25 files changed, 557 insertions(+), 139 deletions(-) create mode 100644 src/GF/BSE2_A_matrix_static.f90 create mode 100644 src/GF/BSE2_B_matrix_static.f90 rename src/LR/{linear_response_Tmatrix.f90 => linear_response_BSE.f90} (78%) create mode 100644 src/utils/level_shifting.f90 diff --git a/input/dft b/input/dft index fba2c7e..24eda5d 100644 --- a/input/dft +++ b/input/dft @@ -19,11 +19,11 @@ # Number of states in ensemble (nEns) 2 # occupation numbers -1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 @@ -31,7 +31,7 @@ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 0.5 0.0 0.0 + 0.1 0.0 0.0 # Ncentered ? F # Parameters for CC weight-dependent exchange functional diff --git a/input/methods b/input/methods index df59e1c..07cc57d 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF KS MOM - F T F F + T F F F # MP2* MP3 MP2-F12 F F F # CCD pCCD DCD CCSD CCSD(T) @@ -9,13 +9,13 @@ # CIS* CIS(D) CID CISD FCI F F F F F # RPA* RPAx* crRPA ppRPA - F F F T + F F F F # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW F F F F F # G0T0 evGT qsGT - F F F + T F F # MCMP2 F # * unrestricted version available diff --git a/input/options b/input/options index 66bba38..a48a761 100644 --- a/input/options +++ b/input/options @@ -15,6 +15,6 @@ # ACFDT: AC Kx XBS F F F # BSE: BSE dBSE dTDA evDyn - T F T F + T T T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/GF/BSE2.f90 b/src/GF/BSE2.f90 index 25bc56c..2d365ed 100644 --- a/src/GF/BSE2.f90 +++ b/src/GF/BSE2.f90 @@ -1,6 +1,6 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,EcBSE) -! Compute the Bethe-Salpeter excitation energies +! Compute the second-order Bethe-Salpeter excitation energies implicit none include 'parameters.h' @@ -33,6 +33,8 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XmY(:,:,:) double precision :: rho + double precision,allocatable :: A_sta(:,:,:) + double precision,allocatable :: B_sta(:,:,:) ! Output variables @@ -40,7 +42,8 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, ! Memory allocation - allocate(OmBSE(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) + allocate(OmBSE(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin), & + A_sta(nS,nS,nspin),B_sta(nS,nS,nspin)) !------------------- ! Singlet manifold @@ -51,13 +54,20 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, ispin = 1 EcBSE(ispin) = 0d0 + ! Compute static kernel + + call BSE2_A_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,A_sta(:,:,ispin)) + if(.not.TDA) call BSE2_B_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_sta(:,:,ispin)) + + call matout(nS,nS,A_sta(:,:,ispin)) + call matout(nS,nS,B_sta(:,:,ispin)) + ! Compute BSE2 excitation energies - call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI, & - OmBSE(:,ispin),rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response_BSE(ispin,.false.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,-A_sta(:,:,ispin),-B_sta(:,:,ispin), & + EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin)) - call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, & - OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) ! Compute dynamic correction for BSE via perturbation theory @@ -67,11 +77,11 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, if(evDyn) then call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & - OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + A_sta(:,:,ispin),B_sta(:,:,ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) else call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & - OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + A_sta(:,:,ispin),B_sta(:,:,ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) end if @@ -88,13 +98,17 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, ispin = 2 EcBSE(ispin) = 0d0 + ! Compute static kernel + + call BSE2_A_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,A_sta(:,:,ispin)) + if(.not.TDA) call BSE2_B_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_sta(:,:,ispin)) + ! Compute BSE2 excitation energies - call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF(:),ERI(:,:,:,:), & - OmBSE(:,ispin),rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response_BSE(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,-A_sta(:,:,ispin),-B_sta(:,:,ispin), & + EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin)) - call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & - OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) ! Compute dynamic correction for BSE via perturbation theory @@ -103,11 +117,11 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, if(evDyn) then call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & - OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + A_sta(:,:,ispin),B_sta(:,:,ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) else call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & - OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + A_sta(:,:,ispin),B_sta(:,:,ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) end if diff --git a/src/GF/BSE2_A_matrix_static.f90 b/src/GF/BSE2_A_matrix_static.f90 new file mode 100644 index 0000000..dee7202 --- /dev/null +++ b/src/GF/BSE2_A_matrix_static.f90 @@ -0,0 +1,157 @@ +subroutine BSE2_A_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,A_sta) + +! Compute the resonant part of the static BSE2 matrix + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eGF(nBas) + +! Local variables + + double precision :: dem,num + integer :: i,j,k,l + integer :: a,b,c,d + integer :: ia,jb + +! Output variables + + double precision,intent(out) :: A_sta(nS,nS) + +! Initialization + + A_sta(:,:) = 0d0 + +! Second-order correlation kernel for the block A of the singlet manifold + + if(ispin == 1) then + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + do k=nC+1,nO + do c=nO+1,nBas-nR + + dem = - (eGF(c) - eGF(k)) + num = 2d0*ERI(j,k,i,c)*ERI(a,c,b,k) - ERI(j,k,i,c)*ERI(a,c,k,b) & + - ERI(j,k,c,i)*ERI(a,c,b,k) + 2d0*ERI(j,k,c,i)*ERI(a,c,k,b) + + A_sta(ia,jb) = A_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + dem = + (eGF(c) - eGF(k)) + num = 2d0*ERI(j,c,i,k)*ERI(a,k,b,c) - ERI(j,c,i,k)*ERI(a,k,c,b) & + - ERI(j,c,k,i)*ERI(a,k,b,c) + 2d0*ERI(j,c,k,i)*ERI(a,k,c,b) + + A_sta(ia,jb) = A_sta(ia,jb) + num*dem/(dem**2 + eta**2) + + end do + end do + + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + + dem = - (eGF(c) + eGF(d)) + num = 2d0*ERI(a,j,c,d)*ERI(c,d,i,b) - ERI(a,j,c,d)*ERI(c,d,b,i) & + - ERI(a,j,d,c)*ERI(c,d,i,b) + 2d0*ERI(a,j,d,c)*ERI(c,d,b,i) + + A_sta(ia,jb) = A_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + do k=nC+1,nO + do l=nC+1,nO + + dem = - (eGF(k) + eGF(l)) + num = 2d0*ERI(a,j,k,l)*ERI(k,l,i,b) - ERI(a,j,k,l)*ERI(k,l,b,i) & + - ERI(a,j,l,k)*ERI(k,l,i,b) + 2d0*ERI(a,j,l,k)*ERI(k,l,b,i) + + A_sta(ia,jb) = A_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + + end if + +! Second-order correlation kernel for the block A of the triplet manifold + + if(ispin == 2) then + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + do k=nC+1,nO + do c=nO+1,nBas-nR + + dem = - (eGF(c) - eGF(k)) + num = 2d0*ERI(j,k,i,c)*ERI(a,c,b,k) - ERI(j,k,i,c)*ERI(a,c,k,b) - ERI(j,k,c,i)*ERI(a,c,b,k) + + A_sta(ia,jb) = A_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + dem = + (eGF(c) - eGF(k)) + num = 2d0*ERI(j,c,i,k)*ERI(a,k,b,c) - ERI(j,c,i,k)*ERI(a,k,c,b) - ERI(j,c,k,i)*ERI(a,k,b,c) + + A_sta(ia,jb) = A_sta(ia,jb) + num*dem/(dem**2 + eta**2) + + end do + end do + + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + + dem = - (eGF(c) + eGF(d)) + num = ERI(a,j,c,d)*ERI(c,d,b,i) + ERI(a,j,d,c)*ERI(c,d,i,b) + + A_sta(ia,jb) = A_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + do k=nC+1,nO + do l=nC+1,nO + + dem = - (eGF(k) + eGF(l)) + num = ERI(a,j,k,l)*ERI(k,l,b,i) + ERI(a,j,l,k)*ERI(k,l,i,b) + + A_sta(ia,jb) = A_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + + end if + + +end subroutine BSE2_A_matrix_static diff --git a/src/GF/BSE2_B_matrix_static.f90 b/src/GF/BSE2_B_matrix_static.f90 new file mode 100644 index 0000000..e93abba --- /dev/null +++ b/src/GF/BSE2_B_matrix_static.f90 @@ -0,0 +1,157 @@ +subroutine BSE2_B_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_sta) + +! Compute the anti-resonant part of the static BSE2 matrix + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eGF(nBas) + +! Local variables + + double precision :: dem,num + integer :: i,j,k,l + integer :: a,b,c,d + integer :: ia,jb + +! Output variables + + double precision,intent(out) :: B_sta(nS,nS) + +! Initialization + + B_sta(:,:) = 0d0 + +! Second-order correlation kernel for the block A of the singlet manifold + + if(ispin == 1) then + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + do k=nC+1,nO + do c=nO+1,nBas-nR + + dem = + eGF(k) - eGF(c) + num = 2d0*ERI(b,k,i,c)*ERI(a,c,j,k) - ERI(b,k,i,c)*ERI(a,c,k,j) & + - ERI(b,k,c,i)*ERI(a,c,j,k) + 2d0*ERI(b,k,c,i)*ERI(a,c,k,j) + + B_sta(ia,jb) = B_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + dem = - eGF(c) + eGF(k) + num = 2d0*ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) & + - ERI(b,c,k,i)*ERI(a,k,j,c) + 2d0*ERI(b,c,k,i)*ERI(a,k,c,j) + + B_sta(ia,jb) = B_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + end do + end do + + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + + dem = - eGF(c) - eGF(d) + num = 2d0*ERI(a,b,c,d)*ERI(c,d,i,j) - ERI(a,b,c,d)*ERI(c,d,j,i) & + - ERI(a,b,d,c)*ERI(c,d,i,j) + 2d0*ERI(a,b,d,c)*ERI(c,d,j,i) + + B_sta(ia,jb) = B_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + do k=nC+1,nO + do l=nC+1,nO + + dem = + eGF(k) + eGF(l) + num = 2d0*ERI(a,b,k,l)*ERI(k,l,i,j) - ERI(a,b,k,l)*ERI(k,l,j,i) & + - ERI(a,b,l,k)*ERI(k,l,i,j) + 2d0*ERI(a,b,l,k)*ERI(k,l,j,i) + + B_sta(ia,jb) = B_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + + end if + +! Second-order correlation kernel for the block A of the triplet manifold + + if(ispin == 2) then + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + do k=nC+1,nO + do c=nO+1,nBas-nR + + dem = + eGF(k) - eGF(c) + num = 2d0*ERI(b,k,i,c)*ERI(a,c,j,k) - ERI(b,k,i,c)*ERI(a,c,k,j) - ERI(b,k,c,i)*ERI(a,c,j,k) + + B_sta(ia,jb) = B_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + dem = - eGF(c) + eGF(k) + num = 2d0*ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) - ERI(b,c,k,i)*ERI(a,k,j,c) + + B_sta(ia,jb) = B_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + end do + end do + + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + + dem = - eGF(c) - eGF(d) + num = ERI(a,b,c,d)*ERI(c,d,j,i) + ERI(a,b,d,c)*ERI(c,d,i,j) + + B_sta(ia,jb) = B_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + do k=nC+1,nO + do l=nC+1,nO + + dem = + eGF(k) + eGF(l) + num = ERI(a,b,k,l)*ERI(k,l,j,i) + ERI(a,b,l,k)*ERI(k,l,i,j) + + B_sta(ia,jb) = B_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + + end if + + +end subroutine BSE2_B_matrix_static diff --git a/src/GF/BSE2_dynamic_perturbation.f90 b/src/GF/BSE2_dynamic_perturbation.f90 index 94c7c59..50077f1 100644 --- a/src/GF/BSE2_dynamic_perturbation.f90 +++ b/src/GF/BSE2_dynamic_perturbation.f90 @@ -1,4 +1,4 @@ -subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,OmBSE,XpY,XmY) +subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,A_sta,B_sta,OmBSE,XpY,XmY) ! Compute dynamical effects via perturbation theory for BSE @@ -21,6 +21,8 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipo double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eGF(nBas) + double precision,intent(in) :: A_sta(nS,nS) + double precision,intent(in) :: B_sta(nS,nS) double precision,intent(in) :: OmBSE(nS) double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XmY(nS,nS) @@ -81,7 +83,7 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipo if(dTDA) then ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) - OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) + OmDyn(ia) = dot_product(X,matmul(Ap_dyn - A_sta,X)) else @@ -94,10 +96,10 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipo ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & + dot_product(Y,matmul(ZAm_dyn,Y)) - OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) & - - dot_product(Y,matmul(Am_dyn,Y)) & - + dot_product(X,matmul(B_dyn,Y)) & - - dot_product(Y,matmul(B_dyn,X)) + OmDyn(ia) = dot_product(X,matmul(Ap_dyn - A_sta,X)) & + - dot_product(Y,matmul(Am_dyn - A_sta,Y)) & + + dot_product(X,matmul(B_dyn - B_sta,Y)) & + - dot_product(Y,matmul(B_dyn - B_sta,X)) end if diff --git a/src/GF/BSE2_dynamic_perturbation_iterative.f90 b/src/GF/BSE2_dynamic_perturbation_iterative.f90 index eedfad1..9147bcd 100644 --- a/src/GF/BSE2_dynamic_perturbation_iterative.f90 +++ b/src/GF/BSE2_dynamic_perturbation_iterative.f90 @@ -1,5 +1,5 @@ subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int, & - eHF,eGF,OmBSE,XpY,XmY) + eHF,eGF,A_sta,B_sta,OmBSE,XpY,XmY) ! Compute self-consistently the dynamical effects via perturbation theory for BSE2 @@ -22,6 +22,8 @@ subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,n double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eGF(nBas) + double precision,intent(in) :: A_sta(nS,nS) + double precision,intent(in) :: B_sta(nS,nS) double precision,intent(in) :: OmBSE(nS) double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XmY(nS,nS) @@ -102,7 +104,7 @@ subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,n if(dTDA) then - OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) + OmDyn(ia) = dot_product(X,matmul(Ap_dyn - A_sta,X)) ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) else @@ -116,10 +118,10 @@ subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,n ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & + dot_product(Y,matmul(ZAm_dyn,Y)) - OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) & - - dot_product(Y,matmul(Am_dyn,Y)) & - + dot_product(X,matmul(B_dyn,Y)) & - - dot_product(Y,matmul(B_dyn,X)) + OmDyn(ia) = dot_product(X,matmul(Ap_dyn - A_sta,X)) & + - dot_product(Y,matmul(Am_dyn - A_sta,Y)) & + + dot_product(X,matmul(B_dyn - B_sta,Y)) & + - dot_product(Y,matmul(B_dyn - B_sta,X)) end if diff --git a/src/GT/Bethe_Salpeter_Tmatrix.f90 b/src/GT/Bethe_Salpeter_Tmatrix.f90 index b038740..ae8bd95 100644 --- a/src/GT/Bethe_Salpeter_Tmatrix.f90 +++ b/src/GT/Bethe_Salpeter_Tmatrix.f90 @@ -124,8 +124,8 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, ! Compute BSE singlet excitation energies - call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAt+TAs,TBt+TBs, & - EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call linear_response_BSE(ispin,.false.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAt+TAs,TBt+TBs, & + EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, & @@ -163,8 +163,8 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, ! Compute BSE triplet excitation energies - call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAt-TAs,TBt-TBs, & - EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call linear_response_BSE(ispin,.false.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAt-TAs,TBt-TBs, & + EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) diff --git a/src/GT/Bethe_Salpeter_Tmatrix_so.f90 b/src/GT/Bethe_Salpeter_Tmatrix_so.f90 index ded548e..8dffe26 100644 --- a/src/GT/Bethe_Salpeter_Tmatrix_so.f90 +++ b/src/GT/Bethe_Salpeter_Tmatrix_so.f90 @@ -71,8 +71,8 @@ subroutine Bethe_Salpeter_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,Omega1,X1,Y ! Compute BSE singlet excitation energies - call linear_response_Tmatrix(ispin,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TA,TB, & - EcBSE,OmBSE,XpY_BSE,XmY_BSE) + call linear_response_BSE(ispin,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TA,TB, & + EcBSE,OmBSE,XpY_BSE,XmY_BSE) call print_excitation('BSE@GT ',ispin,nS,OmBSE) diff --git a/src/GW/Bethe_Salpeter.f90 b/src/GW/Bethe_Salpeter.f90 index a71949a..2440585 100644 --- a/src/GW/Bethe_Salpeter.f90 +++ b/src/GW/Bethe_Salpeter.f90 @@ -42,14 +42,17 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, double precision,allocatable :: XpY_BSE(:,:,:) double precision,allocatable :: XmY_BSE(:,:,:) + double precision,allocatable :: WA_sta(:,:) + double precision,allocatable :: WB_sta(:,:) + ! Output variables double precision,intent(out) :: EcBSE(nspin) ! Memory allocation - allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) - allocate(OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin)) + allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & + WA_sta(nS,nS),WB_sta(nS,nS),OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin)) !--------------------------------- ! Compute (singlet) RPA screening @@ -58,10 +61,13 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, isp_W = 1 EcRPA = 0d0 - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA_sta) + call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB_sta) + !------------------- ! Singlet manifold !------------------- @@ -73,8 +79,8 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, ! Compute BSE excitation energies - call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, & - rho_RPA,EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,WA_sta,WB_sta, & + EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin)) call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, & OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) @@ -112,8 +118,8 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, ! Compute BSE excitation energies - call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, & - rho_RPA,EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,WA_sta,WB_sta, & + EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index aa35dd4..e6f49bc 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -109,8 +109,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & ! Compute screening ! !-------------------! - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0, & - eHF,ERI_MO,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0, & + eHF,ERI_MO,EcRPA,OmRPA,XpY_RPA,XmY_RPA) if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA) @@ -168,8 +168,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & ! Compute the RPA correlation energy - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI_MO,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI_MO, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) !--------------! ! Dump results ! diff --git a/src/GW/G0W0_SOSEX.f90 b/src/GW/G0W0_SOSEX.f90 index fa1eeb5..696fa3b 100644 --- a/src/GW/G0W0_SOSEX.f90 +++ b/src/GW/G0W0_SOSEX.f90 @@ -93,8 +93,8 @@ subroutine G0W0_SOSEX(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDy !-------------------! do ispin=1,nspin - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI_MO, & - OmRPA(:,ispin),rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) + call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI_MO, & + EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA) end do @@ -126,8 +126,8 @@ subroutine G0W0_SOSEX(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDy ! Compute the RPA correlation energy do ispin=1,nspin - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eSOSEX,ERI_MO,OmRPA(:,ispin), & - rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) + call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eSOSEX,ERI_MO, & + EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) end do !--------------! diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90 index f510975..c691616 100644 --- a/src/GW/evGW.f90 +++ b/src/GW/evGW.f90 @@ -151,8 +151,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, if(.not. GW0 .or. nSCF == 0) then - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) end if diff --git a/src/GW/qsGW.f90 b/src/GW/qsGW.f90 index 806ae6d..9f1d33b 100644 --- a/src/GW/qsGW.f90 +++ b/src/GW/qsGW.f90 @@ -181,8 +181,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, if(.not. GW0 .or. nSCF == 0) then - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO, & - OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) if(print_W) call print_excitation('RPA@qsGW ',ispin,nS,OmRPA) endif diff --git a/src/LR/linear_response.f90 b/src/LR/linear_response.f90 index ee80c72..e239f54 100644 --- a/src/LR/linear_response.f90 +++ b/src/LR/linear_response.f90 @@ -1,4 +1,4 @@ -subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Omega_RPA,rho_RPA,EcRPA,Omega,XpY,XmY) +subroutine linear_response(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Ec,Omega,XpY,XmY) ! Compute linear response @@ -9,7 +9,6 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E logical,intent(in) :: dRPA logical,intent(in) :: TDA - logical,intent(in) :: BSE double precision,intent(in) :: eta integer,intent(in) :: ispin integer,intent(in) :: nBas @@ -22,9 +21,6 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E double precision,intent(in) :: e(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega_RPA(nS) - double precision,intent(in) :: rho_RPA(nBas,nBas,nS) - ! Local variables double precision :: trace_matrix @@ -38,7 +34,7 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E ! Output variables - double precision,intent(out) :: EcRPA + double precision,intent(out) :: Ec double precision,intent(out) :: Omega(nS) double precision,intent(out) :: XpY(nS,nS) double precision,intent(out) :: XmY(nS,nS) @@ -51,8 +47,6 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A) - if(BSE) call Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega_RPA,rho_RPA,A) - ! Tamm-Dancoff approximation if(TDA) then @@ -67,8 +61,6 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B) - if(BSE) call Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega_RPA,rho_RPA,B) - ! Build A + B and A - B matrices ApB = A + B @@ -111,6 +103,6 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E ! Compute the RPA correlation energy - EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nS,A)) + Ec = 0.5d0*(sum(Omega) - trace_matrix(nS,A)) end subroutine linear_response diff --git a/src/LR/linear_response_Tmatrix.f90 b/src/LR/linear_response_BSE.f90 similarity index 78% rename from src/LR/linear_response_Tmatrix.f90 rename to src/LR/linear_response_BSE.f90 index 89b9207..b4078a0 100644 --- a/src/LR/linear_response_Tmatrix.f90 +++ b/src/LR/linear_response_BSE.f90 @@ -1,4 +1,4 @@ -subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A_BSE,B_BSE,EcRPA,Omega,XpY,XmY) +subroutine linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A_BSE,B_BSE,Ec,Omega,XpY,XmY) ! Compute linear response @@ -7,9 +7,17 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda ! Input variables - logical,intent(in) :: dRPA,TDA + logical,intent(in) :: dRPA + logical,intent(in) :: TDA + logical,intent(in) :: BSE double precision,intent(in) :: eta - integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS double precision,intent(in) :: lambda double precision,intent(in) :: e(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) @@ -29,7 +37,7 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda ! Output variables - double precision,intent(out) :: EcRPA + double precision,intent(out) :: Ec double precision,intent(out) :: Omega(nS) double precision,intent(out) :: XpY(nS,nS) double precision,intent(out) :: XmY(nS,nS) @@ -47,7 +55,7 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda ! print*,'TA' ! call matout(nS,nS,A_BSE) - A(:,:) = A(:,:) - A_BSE(:,:) + if(BSE) A(:,:) = A(:,:) - A_BSE(:,:) ! Tamm-Dancoff approximation @@ -68,7 +76,7 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda ! print*,'TB' ! call matout(nS,nS,B_BSE) - B(:,:) = B(:,:) - B_BSE(:,:) + if(BSE) B(:,:) = B(:,:) - B_BSE(:,:) ! Build A + B and A - B matrices @@ -112,6 +120,6 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda ! Compute the RPA correlation energy - EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nS,A)) + Ec = 0.5d0*(sum(Omega) - trace_matrix(nS,A)) -end subroutine linear_response_Tmatrix +end subroutine linear_response_BSE diff --git a/src/RPA/ACFDT.f90 b/src/RPA/ACFDT.f90 index 24d9de0..4a7f33a 100644 --- a/src/RPA/ACFDT.f90 +++ b/src/RPA/ACFDT.f90 @@ -32,6 +32,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB double precision,allocatable :: Ec(:,:) double precision :: EcRPA + double precision,allocatable :: WA(:,:) + double precision,allocatable :: WB(:,:) double precision,allocatable :: OmRPA(:) double precision,allocatable :: XpY_RPA(:,:) double precision,allocatable :: XmY_RPA(:,:) @@ -48,7 +50,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB ! Memory allocation allocate(Ec(nAC,nspin)) - allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) + allocate(WA(nS,nS),WB(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) ! Antisymmetrized kernel version @@ -69,10 +71,13 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB isp_W = 1 EcRPA = 0d0 - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA) + call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB) + ! Singlet manifold if(singlet) then @@ -94,15 +99,18 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB if(doXBS) then - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) ! call print_excitation('W^lambda: ',isp_W,nS,OmRPA) + call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA) + call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) + end if - call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & - rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,WA,WB, & + EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) @@ -143,14 +151,17 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB if(doXBS) then - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA) + call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) + end if - call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & - rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,WA,WB, & + EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) diff --git a/src/RPA/ACFDT_Tmatrix.f90 b/src/RPA/ACFDT_Tmatrix.f90 index 9b25e29..2009352 100644 --- a/src/RPA/ACFDT_Tmatrix.f90 +++ b/src/RPA/ACFDT_Tmatrix.f90 @@ -144,8 +144,8 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple end if - call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, & - EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response_BSE(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, & + EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) @@ -214,8 +214,8 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple end if - call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, & - EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response_BSE(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, & + EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) diff --git a/src/RPA/ACFDT_cr.f90 b/src/RPA/ACFDT_cr.f90 index d82c77e..7d128a6 100644 --- a/src/RPA/ACFDT_cr.f90 +++ b/src/RPA/ACFDT_cr.f90 @@ -33,6 +33,8 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta double precision,allocatable :: Ec(:,:) double precision :: EcRPA + double precision,allocatable :: WA(:,:) + double precision,allocatable :: WB(:,:) double precision,allocatable :: OmRPA(:) double precision,allocatable :: XpY_RPA(:,:) double precision,allocatable :: XmY_RPA(:,:) @@ -49,7 +51,7 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta ! Memory allocation allocate(Ec(nAC,nspin)) - allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) + allocate(WA(nS,nS),WB(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) ! Antisymmetrized kernel version @@ -70,10 +72,13 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta isp_W = 1 EcRPA = 0d0 - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA) + call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB) + ! Singlet manifold if(singlet) then @@ -95,15 +100,18 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta if(doXBS) then - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) ! call print_excitation('W^lambda: ',isp_W,nS,OmRPA) + call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA) + call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) + end if - call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & - rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,WA,WB, & + EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) @@ -144,14 +152,17 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta if(doXBS) then - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA) + call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) + end if - call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & - rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,WA,WB, & + EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) diff --git a/src/RPA/RPA.f90 b/src/RPA/RPA.f90 index 0dce1a7..e02423a 100644 --- a/src/RPA/RPA.f90 +++ b/src/RPA/RPA.f90 @@ -33,7 +33,6 @@ subroutine RPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR, double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XmY(:,:,:) - double precision :: rho double precision :: EcRPA(nspin) double precision :: EcAC(nspin) @@ -67,7 +66,7 @@ subroutine RPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR, ispin = 1 - call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), & + call linear_response(ispin,.true.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) @@ -80,7 +79,7 @@ subroutine RPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR, ispin = 2 - call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), & + call linear_response(ispin,.true.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) diff --git a/src/RPA/RPAx.f90 b/src/RPA/RPAx.f90 index f47351d..422d6ec 100644 --- a/src/RPA/RPAx.f90 +++ b/src/RPA/RPAx.f90 @@ -34,7 +34,6 @@ subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XmY(:,:,:) - double precision :: rho double precision :: EcRPAx(nspin) double precision :: EcAC(nspin) @@ -69,7 +68,7 @@ subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR ispin = 1 - call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Omega(:,ispin),rho, & + call linear_response(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('RPAx@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) @@ -82,7 +81,7 @@ subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR ispin = 2 - call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Omega(:,ispin),rho, & + call linear_response(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('RPAx@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) diff --git a/src/RPA/crRPA.f90 b/src/RPA/crRPA.f90 index 7671595..3312092 100644 --- a/src/RPA/crRPA.f90 +++ b/src/RPA/crRPA.f90 @@ -34,7 +34,6 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XmY(:,:,:) - double precision :: rho double precision :: EcRPAx(nspin) double precision :: EcAC(nspin) @@ -68,7 +67,7 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n ispin = 1 - call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Omega(:,ispin),rho, & + call linear_response(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI, & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('crRPA@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) @@ -81,7 +80,7 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n ispin = 2 - call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Omega(:,ispin),rho, & + call linear_response(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI, & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('crRPA@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90 index a9b527b..8888243 100644 --- a/src/eDFT/eDFT_UKS.f90 +++ b/src/eDFT/eDFT_UKS.f90 @@ -45,7 +45,9 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max integer :: xc_rung logical :: LDA_centered = .false. - integer :: nSCF,nBasSq + logical :: do_level_shift = .false. + integer :: nSCF + integer :: nBasSq integer :: n_diis integer :: nO(nspin) double precision :: conv @@ -121,31 +123,38 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max do ispin=1,nspin nO(ispin) = int(sum(occnum(:,ispin,1))) end do - - if(guess_type == 1) then + + do ispin=1,nspin + call mo_guess(nBas,nO(ispin),guess_type,S,Hc,ERI,J(:,:,ispin),Fx(:,:,ispin),X,cp(:,:,ispin),F(:,:,ispin), & + Fp(:,:,ispin),eKS(:,ispin),c(:,:,ispin),Pw(:,:,ispin)) + end do - do ispin=1,nspin - cp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:))) - call diagonalize_matrix(nBas,cp(:,:,ispin),eKS(:,ispin)) - c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin)) - end do + if(mix) call mix_guess(nBas,nO,c) - ! Mix guess to enforce symmetry breaking +! if(guess_type == 1) then - if(mix) call mix_guess(nBas,nO,c) +! do ispin=1,nspin +! cp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:))) +! call diagonalize_matrix(nBas,cp(:,:,ispin),eKS(:,ispin)) +! c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin)) +! end do - else if(guess_type == 2) then +! ! Mix guess to enforce symmetry breaking - do ispin=1,nspin - call random_number(F(:,:,ispin)) - end do +! if(mix) call mix_guess(nBas,nO,c) - else +! else if(guess_type == 2) then - print*,'Wrong guess option' - stop +! do ispin=1,nspin +! call random_number(F(:,:,ispin)) +! end do - end if +! else + +! print*,'Wrong guess option' +! stop + +! end if ! Initialization @@ -266,17 +275,32 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max if(nSCF > 1) conv = maxval(abs(err(:,:,:))) +! Level-shifting + + if(do_level_shift) then + + do ispin=1,nspin + call level_shifting(nBas,nO(ispin),S,c,F(:,:,ispin)) + end do + + end if + ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - if(minval(rcond(:)) > 1d-15) then - do ispin=1,nspin + do ispin=1,nspin + + if(rcond(ispin) > 1d-15) then + call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis, & err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin)) - end do - else - n_diis = 0 - end if + else + + n_diis = 0 + + end if + + end do ! Transform Fock matrix in orthogonal basis diff --git a/src/utils/level_shifting.f90 b/src/utils/level_shifting.f90 new file mode 100644 index 0000000..022b524 --- /dev/null +++ b/src/utils/level_shifting.f90 @@ -0,0 +1,37 @@ +subroutine level_shifting(nBas,nO,S,c,F) + +! Perform level-shifting on the Fock matrix + + implicit none + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nO + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: c(nBas,nBas) + +! Local variables + + double precision,allocatable :: F_MO(:,:) + double precision,allocatable :: Sc(:,:) + double precision,parameter :: LS = 0.1d0 + + integer :: a + +! Output variables + + double precision,intent(inout):: F(nBas,nBas) + + allocate(F_MO(nBas,nBas),Sc(nBas,nBas)) + + F_MO(:,:) = matmul(transpose(c),matmul(F,c)) + + do a=nO+1,nBas + F_MO(a,a) = F_MO(a,a) + LS + end do + + Sc(:,:) = matmul(S,c) + F(:,:) = matmul(Sc,matmul(F_MO,transpose(Sc))) + +end subroutine level_shifting