LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-07 06:33:55 +01:00
Code Releases Activity

clean linear response and static BSE2 kernel

Browse Source
This commit is contained in:
Pierre-Francois Loos 2022-01-31 13:19:42 +01:00
parent 23e520cea9
commit f7718c8ab1
25 changed files with 557 additions and 139 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
input/dft
View File

@ -19,11 +19,11 @@
# Number of states in ensemble (nEns) # Number of states in ensemble (nEns)
2 2
# occupation numbers # 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 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 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 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 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
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 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
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 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) # Ensemble weights: wEns(1),...,wEns(nEns-1)
0.5 0.0 0.0 0.1 0.0 0.0
# Ncentered ? # Ncentered ?
F F
# Parameters for CC weight-dependent exchange functional # Parameters for CC weight-dependent exchange functional

6
input/methods
View File

@ -1,5 +1,5 @@
# RHF UHF KS MOM # RHF UHF KS MOM
F T F F T F F F
# MP2* MP3 MP2-F12 # MP2* MP3 MP2-F12
F F F F F F
# CCD pCCD DCD CCSD CCSD(T) # CCD pCCD DCD CCSD CCSD(T)
@ -9,13 +9,13 @@
# CIS* CIS(D) CID CISD FCI # CIS* CIS(D) CID CISD FCI
F F F F F F F F F F
# RPA* RPAx* crRPA ppRPA # RPA* RPAx* crRPA ppRPA
F F F T F F F F
# G0F2* evGF2* qsGF2* G0F3 evGF3 # G0F2* evGF2* qsGF2* G0F3 evGF3
F F F F F F F F F F
# G0W0* evGW* qsGW* ufG0W0 ufGW # G0W0* evGW* qsGW* ufG0W0 ufGW
F F F F F F F F F F
# G0T0 evGT qsGT # G0T0 evGT qsGT
F F F T F F
# MCMP2 # MCMP2
F F
# * unrestricted version available # * unrestricted version available

2
input/options
View File

@ -15,6 +15,6 @@
# ACFDT: AC Kx XBS # ACFDT: AC Kx XBS
F F F F F F
# BSE: BSE dBSE dTDA evDyn # BSE: BSE dBSE dTDA evDyn
T F T F T T T F
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T 1000000 100000 10 0.3 10000 1234 T

42
src/GF/BSE2.f90
View File

@ -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) 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 implicit none
include 'parameters.h' 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 :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:) double precision,allocatable :: XmY(:,:,:)
double precision :: rho double precision :: rho
double precision,allocatable :: A_sta(:,:,:)
double precision,allocatable :: B_sta(:,:,:)
! Output variables ! Output variables
@ -40,7 +42,8 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,
! Memory allocation ! 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 ! Singlet manifold
@ -51,13 +54,20 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,
ispin = 1 ispin = 1
EcBSE(ispin) = 0d0 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 ! Compute BSE2 excitation energies
call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI, & call linear_response_BSE(ispin,.false.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,-A_sta(:,:,ispin),-B_sta(:,:,ispin), &
OmBSE(:,ispin),rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin)) call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin))
call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, & call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
! Compute dynamic correction for BSE via perturbation theory ! 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 if(evDyn) then
call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & 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 else
call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & 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 end if
@ -88,13 +98,17 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,
ispin = 2 ispin = 2
EcBSE(ispin) = 0d0 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 ! Compute BSE2 excitation energies
call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF(:),ERI(:,:,:,:), & call linear_response_BSE(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,-A_sta(:,:,ispin),-B_sta(:,:,ispin), &
OmBSE(:,ispin),rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin)) call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin))
call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
! Compute dynamic correction for BSE via perturbation theory ! 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 if(evDyn) then
call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & 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 else
call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & 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 end if

157
src/GF/BSE2_A_matrix_static.f90 Normal file
View File

@ -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

157
src/GF/BSE2_B_matrix_static.f90 Normal file
View File

@ -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

14
src/GF/BSE2_dynamic_perturbation.f90
View File

@ -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 ! 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) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGF(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) :: OmBSE(nS)
double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XpY(nS,nS)
double precision,intent(in) :: XmY(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 if(dTDA) then
ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) 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 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)) & ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) &
+ dot_product(Y,matmul(ZAm_dyn,Y)) + dot_product(Y,matmul(ZAm_dyn,Y))
OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) & OmDyn(ia) = dot_product(X,matmul(Ap_dyn - A_sta,X)) &
- dot_product(Y,matmul(Am_dyn,Y)) & - dot_product(Y,matmul(Am_dyn - A_sta,Y)) &
+ dot_product(X,matmul(B_dyn,Y)) & + dot_product(X,matmul(B_dyn - B_sta,Y)) &
- dot_product(Y,matmul(B_dyn,X)) - dot_product(Y,matmul(B_dyn - B_sta,X))
end if end if

14
src/GF/BSE2_dynamic_perturbation_iterative.f90
View File

@ -1,5 +1,5 @@
subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int, & 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 ! 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) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGF(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) :: OmBSE(nS)
double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XpY(nS,nS)
double precision,intent(in) :: XmY(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 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)) ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X))
else 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)) & ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) &
+ dot_product(Y,matmul(ZAm_dyn,Y)) + dot_product(Y,matmul(ZAm_dyn,Y))
OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) & OmDyn(ia) = dot_product(X,matmul(Ap_dyn - A_sta,X)) &
- dot_product(Y,matmul(Am_dyn,Y)) & - dot_product(Y,matmul(Am_dyn - A_sta,Y)) &
+ dot_product(X,matmul(B_dyn,Y)) & + dot_product(X,matmul(B_dyn - B_sta,Y)) &
- dot_product(Y,matmul(B_dyn,X)) - dot_product(Y,matmul(B_dyn - B_sta,X))
end if end if

4
src/GT/Bethe_Salpeter_Tmatrix.f90
View File

@ -124,7 +124,7 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
! Compute BSE singlet excitation energies ! 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, & 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)) EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin))
@ -163,7 +163,7 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
! Compute BSE triplet excitation energies ! 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, & 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)) EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin))
call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, &

2
src/GT/Bethe_Salpeter_Tmatrix_so.f90
View File

@ -71,7 +71,7 @@ subroutine Bethe_Salpeter_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,Omega1,X1,Y
! Compute BSE singlet excitation energies ! Compute BSE singlet excitation energies
call linear_response_Tmatrix(ispin,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TA,TB, & 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) EcBSE,OmBSE,XpY_BSE,XmY_BSE)
call print_excitation('BSE@GT ',ispin,nS,OmBSE) call print_excitation('BSE@GT ',ispin,nS,OmBSE)

22
src/GW/Bethe_Salpeter.f90
View File

@ -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 :: XpY_BSE(:,:,:)
double precision,allocatable :: XmY_BSE(:,:,:) double precision,allocatable :: XmY_BSE(:,:,:)
double precision,allocatable :: WA_sta(:,:)
double precision,allocatable :: WB_sta(:,:)
! Output variables ! Output variables
double precision,intent(out) :: EcBSE(nspin) double precision,intent(out) :: EcBSE(nspin)
! Memory allocation ! Memory allocation
allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) 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)) 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 ! 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 isp_W = 1
EcRPA = 0d0 EcRPA = 0d0
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, & call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_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 ! Singlet manifold
!------------------- !-------------------
@ -73,8 +79,8 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,
! Compute BSE excitation energies ! Compute BSE excitation energies
call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, & call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,WA_sta,WB_sta, &
rho_RPA,EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin)) call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin))
call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, & call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, &
OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) 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 ! Compute BSE excitation energies
call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, & call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,WA_sta,WB_sta, &
rho_RPA,EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin)) call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin))
call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, &
OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))

8
src/GW/G0W0.f90
View File

@ -109,8 +109,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
! Compute screening ! ! Compute screening !
!-------------------! !-------------------!
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0, & call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0, &
eHF,ERI_MO,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) eHF,ERI_MO,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA) 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 ! 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, & call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI_MO, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) EcRPA,OmRPA,XpY_RPA,XmY_RPA)
!--------------! !--------------!
! Dump results ! ! Dump results !

8
src/GW/G0W0_SOSEX.f90
View File

@ -93,8 +93,8 @@ subroutine G0W0_SOSEX(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDy
!-------------------! !-------------------!
do ispin=1,nspin do ispin=1,nspin
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI_MO, & call linear_response(ispin,.true.,TDA_W,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)) EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin))
if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA) if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA)
end do 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 ! Compute the RPA correlation energy
do ispin=1,nspin 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), & call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eSOSEX,ERI_MO, &
rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin))
end do end do
!--------------! !--------------!

4
src/GW/evGW.f90
View File

@ -151,8 +151,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
if(.not. GW0 .or. nSCF == 0) then 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, & call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) EcRPA,OmRPA,XpY_RPA,XmY_RPA)
end if end if

4
src/GW/qsGW.f90
View File

@ -181,8 +181,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
if(.not. GW0 .or. nSCF == 0) then 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, & call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO, &
OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) EcRPA,OmRPA,XpY_RPA,XmY_RPA)
if(print_W) call print_excitation('RPA@qsGW ',ispin,nS,OmRPA) if(print_W) call print_excitation('RPA@qsGW ',ispin,nS,OmRPA)
endif endif

14
src/LR/linear_response.f90
View File

@ -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 ! 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) :: dRPA
logical,intent(in) :: TDA logical,intent(in) :: TDA
logical,intent(in) :: BSE
double precision,intent(in) :: eta double precision,intent(in) :: eta
integer,intent(in) :: ispin integer,intent(in) :: ispin
integer,intent(in) :: nBas 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) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,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 ! Local variables
double precision :: trace_matrix 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 ! Output variables
double precision,intent(out) :: EcRPA double precision,intent(out) :: Ec
double precision,intent(out) :: Omega(nS) double precision,intent(out) :: Omega(nS)
double precision,intent(out) :: XpY(nS,nS) double precision,intent(out) :: XpY(nS,nS)
double precision,intent(out) :: XmY(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) 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 ! Tamm-Dancoff approximation
if(TDA) then 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) 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 ! Build A + B and A - B matrices
ApB = A + B 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 ! 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 end subroutine linear_response

24
src/LR/linear_response_Tmatrix.f90 → src/LR/linear_response_BSE.f90
View File

@ -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 ! Compute linear response
@ -7,9 +7,17 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda
! Input variables ! Input variables
logical,intent(in) :: dRPA,TDA logical,intent(in) :: dRPA
logical,intent(in) :: TDA
logical,intent(in) :: BSE
double precision,intent(in) :: eta 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) :: lambda
double precision,intent(in) :: e(nBas) double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,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 ! Output variables
double precision,intent(out) :: EcRPA double precision,intent(out) :: Ec
double precision,intent(out) :: Omega(nS) double precision,intent(out) :: Omega(nS)
double precision,intent(out) :: XpY(nS,nS) double precision,intent(out) :: XpY(nS,nS)
double precision,intent(out) :: XmY(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' ! print*,'TA'
! call matout(nS,nS,A_BSE) ! call matout(nS,nS,A_BSE)
A(:,:) = A(:,:) - A_BSE(:,:) if(BSE) A(:,:) = A(:,:) - A_BSE(:,:)
! Tamm-Dancoff approximation ! Tamm-Dancoff approximation
@ -68,7 +76,7 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda
! print*,'TB' ! print*,'TB'
! call matout(nS,nS,B_BSE) ! call matout(nS,nS,B_BSE)
B(:,:) = B(:,:) - B_BSE(:,:) if(BSE) B(:,:) = B(:,:) - B_BSE(:,:)
! Build A + B and A - B matrices ! 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 ! 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

33
src/RPA/ACFDT.f90
View File

@ -32,6 +32,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB
double precision,allocatable :: Ec(:,:) double precision,allocatable :: Ec(:,:)
double precision :: EcRPA double precision :: EcRPA
double precision,allocatable :: WA(:,:)
double precision,allocatable :: WB(:,:)
double precision,allocatable :: OmRPA(:) double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:) double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_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 ! Memory allocation
allocate(Ec(nAC,nspin)) 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)) allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
! Antisymmetrized kernel version ! Antisymmetrized kernel version
@ -69,10 +71,13 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB
isp_W = 1 isp_W = 1
EcRPA = 0d0 EcRPA = 0d0
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, & call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_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 ! Singlet manifold
if(singlet) then if(singlet) then
@ -94,15 +99,18 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB
if(doXBS) then if(doXBS) then
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
! call print_excitation('W^lambda: ',isp_W,nS,OmRPA) ! 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 end if
call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,WA,WB, &
rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, &
ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) 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 if(doXBS) then
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_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 end if
call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,WA,WB, &
rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) 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)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))

4
src/RPA/ACFDT_Tmatrix.f90
View File

@ -144,7 +144,7 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
end if end if
call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, & 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)) 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)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
@ -214,7 +214,7 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
end if end if
call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, & 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)) 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)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))

33
src/RPA/ACFDT_cr.f90
View File

@ -33,6 +33,8 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta
double precision,allocatable :: Ec(:,:) double precision,allocatable :: Ec(:,:)
double precision :: EcRPA double precision :: EcRPA
double precision,allocatable :: WA(:,:)
double precision,allocatable :: WB(:,:)
double precision,allocatable :: OmRPA(:) double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:) double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_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 ! Memory allocation
allocate(Ec(nAC,nspin)) 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)) allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
! Antisymmetrized kernel version ! Antisymmetrized kernel version
@ -70,10 +72,13 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta
isp_W = 1 isp_W = 1
EcRPA = 0d0 EcRPA = 0d0
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, & call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_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 ! Singlet manifold
if(singlet) then if(singlet) then
@ -95,15 +100,18 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta
if(doXBS) then if(doXBS) then
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
! call print_excitation('W^lambda: ',isp_W,nS,OmRPA) ! 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 end if
call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,WA,WB, &
rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, &
ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) 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 if(doXBS) then
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_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 end if
call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,WA,WB, &
rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) 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)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))

5
src/RPA/RPA.f90
View File

@ -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 :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:) double precision,allocatable :: XmY(:,:,:)
double precision :: rho
double precision :: EcRPA(nspin) double precision :: EcRPA(nspin)
double precision :: EcAC(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 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)) EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPA@HF ',ispin,nS,Omega(:,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)) 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 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)) EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPA@HF ',ispin,nS,Omega(:,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)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))

5
src/RPA/RPAx.f90
View File

@ -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 :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:) double precision,allocatable :: XmY(:,:,:)
double precision :: rho
double precision :: EcRPAx(nspin) double precision :: EcRPAx(nspin)
double precision :: EcAC(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 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)) EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPAx@HF ',ispin,nS,Omega(:,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)) 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 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)) EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPAx@HF ',ispin,nS,Omega(:,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)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))

5
src/RPA/crRPA.f90
View File

@ -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 :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:) double precision,allocatable :: XmY(:,:,:)
double precision :: rho
double precision :: EcRPAx(nspin) double precision :: EcRPAx(nspin)
double precision :: EcAC(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 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)) EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('crRPA@HF ',ispin,nS,Omega(:,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)) 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 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)) EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('crRPA@HF ',ispin,nS,Omega(:,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)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))

60
src/eDFT/eDFT_UKS.f90
View File

@ -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 integer :: xc_rung
logical :: LDA_centered = .false. logical :: LDA_centered = .false.
integer :: nSCF,nBasSq logical :: do_level_shift = .false.
integer :: nSCF
integer :: nBasSq
integer :: n_diis integer :: n_diis
integer :: nO(nspin) integer :: nO(nspin)
double precision :: conv double precision :: conv
@ -122,30 +124,37 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
nO(ispin) = int(sum(occnum(:,ispin,1))) nO(ispin) = int(sum(occnum(:,ispin,1)))
end do end do
if(guess_type == 1) then
do ispin=1,nspin do ispin=1,nspin
cp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:))) call mo_guess(nBas,nO(ispin),guess_type,S,Hc,ERI,J(:,:,ispin),Fx(:,:,ispin),X,cp(:,:,ispin),F(:,:,ispin), &
call diagonalize_matrix(nBas,cp(:,:,ispin),eKS(:,ispin)) Fp(:,:,ispin),eKS(:,ispin),c(:,:,ispin),Pw(:,:,ispin))
c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
end do end do
! Mix guess to enforce symmetry breaking
if(mix) call mix_guess(nBas,nO,c) if(mix) call mix_guess(nBas,nO,c)
else if(guess_type == 2) then ! if(guess_type == 1) then
do ispin=1,nspin ! do ispin=1,nspin
call random_number(F(:,:,ispin)) ! cp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:)))
end do ! call diagonalize_matrix(nBas,cp(:,:,ispin),eKS(:,ispin))
! c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
! end do
else ! ! Mix guess to enforce symmetry breaking
print*,'Wrong guess option' ! if(mix) call mix_guess(nBas,nO,c)
stop
end if ! else if(guess_type == 2) then
! do ispin=1,nspin
! call random_number(F(:,:,ispin))
! end do
! else
! print*,'Wrong guess option'
! stop
! end if
! Initialization ! Initialization
@ -266,18 +275,33 @@ 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(:,:,:))) 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 ! DIIS extrapolation
n_diis = min(n_diis+1,max_diis) 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, & call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis, &
err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin)) err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin))
end do
else else
n_diis = 0 n_diis = 0
end if end if
end do
! Transform Fock matrix in orthogonal basis ! Transform Fock matrix in orthogonal basis
do ispin=1,nspin do ispin=1,nspin

37
src/utils/level_shifting.f90 Normal file
View File

@ -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