LCPQ
/
quack
mirror of https://github.com/pfloos/quack
4
1
Fork 0
Code Releases Activity

merge with T2

This commit is contained in:
Clotilde Marut 2022-02-07 10:45:18 +01:00
parent 5b436558ed 21c1ceb975
commit 82fd7c5fb4
49 changed files with 1249 additions and 437 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/dft
View File

@ -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)
1 0.0 0.0
0.0 0.0 0.0
# Ncentered ?
F
# Parameters for CC weight-dependent exchange functional

8
input/methods
View File

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

12
input/options
View File

@ -1,5 +1,5 @@
# HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess stability
1024 0.00001 T 2 1 1 F F
# HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess level_shift stability
256 0.0000001 T 5 2 1 T 1.0 F
# MP:
# CC: maxSCF thresh DIIS n_diis
@ -7,14 +7,14 @@
# spin: TDA singlet triplet spin_conserved spin_flip
F T T T T
# GF: maxSCF thresh DIIS n_diis lin eta renorm reg
256 0.00001 T 5 T 0.0 3 F
256 0.00001 T 5 T 0.0 3 F
# GW: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 reg
256 0.00001 T 5 T 0.0 F F F F F F
256 0.00001 T 5 T 0.0 F F F F F F
# GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg
256 0.00001 T 5 T 0.0 F F
# ACFDT: AC Kx XBS
F F F
F T T
# BSE: BSE dBSE dTDA evDyn
T T T F
T F T F
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

4
mol/h2.xyz
View File

@ -1,4 +1,4 @@
2
H 0.0 0.0 0.0
H 0.0 0.0 0.7
H 0. 0. 0.
H 0. 0. 0.741

2
src/CI/CISD.f90
View File

@ -86,7 +86,7 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI
write(*,*) 'nH = ',nH
write(*,*)
maxH = min(nH,41)
maxH = min(nH,51)
! Memory allocation

39
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)
! 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,17 @@ 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))
! 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 +74,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 +95,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,.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(.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 +114,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

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

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, &
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

8
src/GT/Bethe_Salpeter_Tmatrix.f90
View File

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

79
src/GT/Bethe_Salpeter_Tmatrix_so.f90 Normal file
View File

@ -0,0 +1,79 @@
subroutine Bethe_Salpeter_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,Omega1,X1,Y1,Omega2,X2,Y2,rho1,rho2, &
ERI,eT,eGT,EcBSE)
! Compute the Bethe-Salpeter excitation energies with the T-matrix kernel
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: eT(nBas)
double precision,intent(in) :: eGT(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega1(nVV)
double precision,intent(in) :: X1(nVV,nVV)
double precision,intent(in) :: Y1(nOO,nVV)
double precision,intent(in) :: Omega2(nOO)
double precision,intent(in) :: X2(nVV,nOO)
double precision,intent(in) :: Y2(nOO,nOO)
double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: rho2(nBas,nBas,nOO)
! Local variables
integer :: ispin
double precision :: EcRPA
double precision,allocatable :: TA(:,:),TB(:,:)
double precision,allocatable :: OmBSE(:)
double precision,allocatable :: XpY_BSE(:,:)
double precision,allocatable :: XmY_BSE(:,:)
! Output variables
double precision,intent(out) :: EcBSE
! Memory allocation
allocate(TA(nS,nS),TB(nS,nS),OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS))
!------------------!
! Compute T-matrix !
!------------------!
ispin = 4
call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eT,ERI, &
Omega1,X1,Y1,Omega2,X2,Y2,EcRPA)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,Omega1,rho1,Omega2,rho2,TA)
call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,Omega1,rho1,Omega2,rho2,TB)
!------------------!
! Singlet manifold !
!------------------!
ispin = 3
EcBSE = 0d0
! Compute BSE singlet excitation energies
call linear_response_BSE(ispin,.false.,.false.,.true.,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)
end subroutine Bethe_Salpeter_Tmatrix_so

8
src/GT/G0T0.f90
View File

@ -238,11 +238,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Free memory
deallocate(Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s, &
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,rho1t,rho2t)
! Compute the BSE correlation energy via the adiabatic connection
if(doACFDT) then
@ -260,7 +255,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
end if
call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
ERI_MO,eHF,eG0T0,EcAC)
nOOs,nVVs,nOOt,nVVt,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t, &
Omega2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eHF,eG0T0,EcAC)
if(exchange_kernel) then

3
src/GT/evGT.f90
View File

@ -265,7 +265,8 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
end if
call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
ERI_MO,eGT,eGT,EcAC)
nOOs,nVVs,nOOt,nVVt,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t, &
Omega2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcAC)
if(exchange_kernel) then

80
src/GT/excitation_density_Tmatrix_so.f90 Normal file
View File

@ -0,0 +1,80 @@
subroutine excitation_density_Tmatrix_so(nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
! Compute excitation densities for T-matrix self-energy
implicit none
! Input variables
integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: X1(nVV,nVV)
double precision,intent(in) :: Y1(nOO,nVV)
double precision,intent(in) :: X2(nVV,nOO)
double precision,intent(in) :: Y2(nOO,nOO)
! Local variables
integer :: j,k,l
integer :: b,c,d
integer :: p,q
integer :: ab,cd,ij,kl
double precision,external :: Kronecker_delta
! Output variables
double precision,intent(out) :: rho1(nBas,nBas,nVV)
double precision,intent(out) :: rho2(nBas,nBAs,nOO)
! Initialization
rho1(:,:,:) = 0d0
rho2(:,:,:) = 0d0
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do ab=1,nVV
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab)
end do
end do
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab)
end do
end do
end do
do ij=1,nOO
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij)
end do
end do
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij)
end do
end do
end do
end do
end do
end subroutine excitation_density_Tmatrix_so

3
src/GT/qsGT.f90
View File

@ -389,7 +389,8 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T
end if
call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
ERI_MO,eGT,eGT,EcAC)
nOOs,nVVs,nOOt,nVVt,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t, &
Omega2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

63
src/GT/renormalization_factor_Tmatrix_so.f90 Normal file
View File

@ -0,0 +1,63 @@
subroutine renormalization_factor_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,Z)
! Compute renormalization factor of the T-matrix self-energy
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: Omega2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO)
! Local variables
integer :: i,j,k,l,a,b,c,d,p,cd,kl
double precision :: eps
! Output variables
double precision,intent(out) :: Z(nBas)
! Initialize
Z(:) = 0d0
!----------------------------------------------
! T-matrix renormalization factor in the spinorbital basis
!----------------------------------------------
! Occupied part of the T-matrix self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do cd=1,nVV
eps = e(p) + e(i) - Omega1(cd)
Z(p) = Z(p) + (rho1(p,i,cd)/eps)**2
enddo
enddo
enddo
! Virtual part of the T-matrix self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do kl=1,nOO
eps = e(p) + e(a) - Omega2(kl)
Z(p) = Z(p) + (rho2(p,a,kl)/eps)**2
enddo
enddo
enddo
! Compute renormalization factor from derivative of SigT
Z(:) = 1d0/(1d0 + Z(:))
end subroutine renormalization_factor_Tmatrix_so

63
src/GT/self_energy_Tmatrix_diag_so.f90 Normal file
View File

@ -0,0 +1,63 @@
subroutine self_energy_Tmatrix_diag_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,SigT)
! Compute diagonal of the correlation part of the T-matrix self-energy
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: Omega2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO)
! Local variables
integer :: i,j,k,l,a,b,c,d,p,cd,kl
double precision :: eps
! Output variables
double precision,intent(out) :: SigT(nBas)
! Initialize
SigT(:) = 0d0
!----------------------------------------------
! T-matrix self-energy in the spinorbital basis
!----------------------------------------------
! Occupied part of the T-matrix self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do cd=1,nVV
eps = e(p) + e(i) - Omega1(cd)
SigT(p) = SigT(p) + rho1(p,i,cd)**2/eps
enddo
enddo
enddo
! Virtual part of the T-matrix self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do kl=1,nOO
eps = e(p) + e(a) - Omega2(kl)
SigT(p) = SigT(p) + rho2(p,a,kl)**2/eps
enddo
enddo
enddo
end subroutine self_energy_Tmatrix_diag_so

124
src/GT/soG0T0.f90 Normal file
View File

@ -0,0 +1,124 @@
subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
! Perform G0W0 calculation with a T-matrix self-energy (G0T0) in the spinorbital basis
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: ispin
integer :: nOO
integer :: nVV
double precision :: EcRPA
double precision :: EcGM
double precision :: EcBSE
integer :: nBas2,nC2,nO2,nV2,nR2,nS2
double precision,allocatable :: Omega1(:)
double precision,allocatable :: X1(:,:)
double precision,allocatable :: Y1(:,:)
double precision,allocatable :: rho1(:,:,:)
double precision,allocatable :: Omega2(:)
double precision,allocatable :: X2(:,:)
double precision,allocatable :: Y2(:,:)
double precision,allocatable :: rho2(:,:,:)
double precision,allocatable :: SigT(:)
double precision,allocatable :: Z(:)
double precision,allocatable :: eG0T0(:)
double precision,allocatable :: seHF(:)
double precision,allocatable :: sERI(:,:,:,:)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| One-shot soG0T0 calculation |'
write(*,*)'************************************************'
write(*,*)
! Define occupied and virtual spaces
nBas2 = 2*nBas
nO2 = 2*nO
nV2 = 2*nV
nC2 = 2*nC
nR2 = 2*nR
nS2 = nO2*nV2
! Spatial to spin orbitals
allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI)
! Dimensions of the rr-RPA linear reponse matrices
nOO = nO2*(nO2 - 1)/2
nVV = nV2*(nV2 - 1)/2
! Memory allocation
allocate(Omega1(nVV),X1(nVV,nVV),Y1(nOO,nVV), &
Omega2(nOO),X2(nVV,nOO),Y2(nOO,nOO), &
rho1(nBas2,nBas2,nVV),rho2(nBas2,nBas2,nOO), &
eG0T0(nBas2),SigT(nBas2),Z(nBas2))
!----------------------------------------------
! Spinorbital basis
!----------------------------------------------
ispin = 4
! Compute linear response
call linear_response_pp(ispin,.false.,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,1d0,seHF,sERI, &
Omega1,X1,Y1,Omega2,X2,Y2,EcRPA)
call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1)
call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2)
! Compute excitation densities for the T-matrix
call excitation_density_Tmatrix_so(nBas2,nC2,nO2,nV2,nR2,nOO,nVV,sERI,X1,Y1,rho1,X2,Y2,rho2)
!----------------------------------------------
! Compute T-matrix version of the self-energy
!----------------------------------------------
call self_energy_Tmatrix_diag_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF, &
Omega1,rho1,Omega2,rho2,SigT)
! Compute renormalization factor for T-matrix self-energy
call renormalization_factor_Tmatrix_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF, &
Omega1,rho1,Omega2,rho2,Z)
!----------------------------------------------
! Solve the quasi-particle equation
!----------------------------------------------
eG0T0(:) = seHF(:) + Z(:)*SigT(:)
!----------------------------------------------
! Dump results
!----------------------------------------------
call print_G0T0(nBas2,nO2,seHF,ENuc,ERHF,SigT,Z,eG0T0,EcGM,EcRPA)
call Bethe_Salpeter_Tmatrix_so(eta,nBas2,nC2,nO2,nV2,nR2,nS2,nOO,nVV,Omega1,X1,Y1,Omega2,X2,Y2,rho1,rho2, &
sERI,seHF,eG0T0,EcBSE)
end subroutine soG0T0

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 :: 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))

8
src/GW/Bethe_Salpeter_AB_matrix_dynamic.f90
View File

@ -90,16 +90,16 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
eps_Am = - OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Am/(eps_Am**2 + eta**2)
eps_Bp = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b))
eps_Bp = - OmRPA(kc) - (eGW(a) - eGW(b))
chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2)
eps_Bp = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i))
eps_Bp = - OmRPA(kc) - (eGW(j) - eGW(i))
chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2)
eps_Bm = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(b))
eps_Bm = - OmRPA(kc) - (eGW(a) - eGW(b))
chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2)
eps_Bm = - OmBSE - OmRPA(kc) - (eGW(j) - eGW(i))
eps_Bm = - OmRPA(kc) - (eGW(j) - eGW(i))
chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2)
enddo

8
src/GW/G0W0.f90
View File

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

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

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

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

3
src/GW/self_energy_exchange_diag.f90
View File

@ -14,7 +14,8 @@ subroutine self_energy_exchange_diag(nBas,c,P,ERI,SigX)
! Local variables
integer :: q,mu,nu
integer :: mu,nu
integer :: q
double precision,allocatable :: Fx(:,:)
! Output variables

55
src/HF/RHF.f90
View File

@ -1,4 +1,5 @@
subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,ERHF,e,c,P,Vx)
subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,ERHF,e,c,P,Vx)
! Perform restricted Hartree-Fock calculation
@ -7,8 +8,11 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
! Input variables
integer,intent(in) :: maxSCF,max_diis,guess_type
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
double precision,intent(in) :: thresh
double precision,intent(in) :: level_shift
integer,intent(in) :: nBas
integer,intent(in) :: nO
@ -46,7 +50,6 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
double precision,allocatable :: K(:,:)
double precision,allocatable :: cp(:,:)
double precision,allocatable :: Fp(:,:)
double precision,allocatable :: ON(:)
! Output variables
@ -71,30 +74,22 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
! Memory allocation
allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas), &
cp(nBas,nBas),Fp(nBas,nBas),ON(nBas), &
allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas),cp(nBas,nBas),Fp(nBas,nBas), &
error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Guess coefficients and eigenvalues
! Guess coefficients and density matrix
call mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
! ON(:) = 0d0
! do i=1,nO
! ON(i) = 1d0
! ON(i) = dble(2*i-1)
! end do
! call density_matrix(nBas,ON,c,P)
call mo_guess(nBas,guess_type,S,Hc,X,c)
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! Initialization
n_diis = 0
F_diis(:,:) = 0d0
error_diis(:,:) = 0d0
Conv = 1d0
nSCF = 0
rcond = 0d0
Conv = 1d0
n_diis = 0
nSCF = 0
rcond = 0d0
!------------------------------------------------------------------------
! Main SCF loop
@ -123,17 +118,25 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
! Check convergence
error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
Conv = maxval(abs(error))
Conv = maxval(abs(error))
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
if(abs(rcond) > 1d-7) then
if(abs(rcond) > 1d-15) then
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
else
n_diis = 0
end if
! Level-shifting
if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nO,S,c,F)
! Diagonalize Fock matrix
Fp = matmul(transpose(X),matmul(F,X))
@ -144,14 +147,12 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
! Density matrix
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! call density_matrix(nBas,ON,c,P)
! Compute HF energy
ERHF = trace_matrix(nBas,matmul(P,Hc)) &
+ 0.5d0*trace_matrix(nBas,matmul(P,J)) &
+ 0.25d0*trace_matrix(nBas,matmul(P,K))
ERHF = trace_matrix(nBas,matmul(P,Hc)) &
+ 0.5d0*trace_matrix(nBas,matmul(P,J)) &
+ 0.25d0*trace_matrix(nBas,matmul(P,K))
! Compute HOMO-LUMO gap
@ -205,6 +206,6 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
! Compute Vx for post-HF calculations
call mo_fock_exchange_potential(nBas,c,K,Vx)
call mo_fock_exchange_potential(nBas,c,P,ERI,Vx)
end subroutine RHF

43
src/HF/UHF.f90
View File

@ -1,4 +1,5 @@
subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,e,c,P,Vx)
subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,e,c,P,Vx)
! Perform unrestricted Hartree-Fock calculation
@ -11,6 +12,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
logical,intent(in) :: mix
double precision,intent(in) :: level_shift
double precision,intent(in) :: thresh
integer,intent(in) :: nBas
@ -79,22 +81,12 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO
K(nBas,nBas,nspin),err(nBas,nBas,nspin),cp(nBas,nBas,nspin), &
err_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin))
! Guess coefficients and eigenvalues
! Guess coefficients and demsity matrices
if(guess_type == 1) then
F(:,:,:) = 0d0
do ispin=1,nspin
F(:,:,ispin) = Hc(:,:)
end do
else if(guess_type == 2) then
do ispin=1,nspin
call random_number(F(:,:,ispin))
end do
end if
do ispin=1,nspin
call mo_guess(nBas,guess_type,S,Hc,X,c(:,:,ispin))
P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
end do
! Initialization
@ -179,13 +171,28 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
if(minval(rcond(:)) > 1d-7) then
if(minval(rcond(:)) > 1d-15) then
do ispin=1,nspin
if(nO(ispin) > 1) call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,1:n_diis,ispin), &
F_diis(:,1:n_diis,ispin),err(:,:,ispin),F(:,:,ispin))
end do
else
n_diis = 0
end if
! Level-shifting
if(level_shift > 0d0 .and. Conv > thresh) then
do ispin=1,nspin
call level_shifting(level_shift,nBas,nO(ispin),S,c(:,:,ispin),F(:,:,ispin))
end do
end if
!------------------------------------------------------------------------
@ -253,7 +260,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO
! Compute Vx for post-HF calculations
do ispin=1,nspin
call mo_fock_exchange_potential(nBas,c(:,:,ispin),K(:,:,ispin),Vx(:,ispin))
call mo_fock_exchange_potential(nBas,c(:,:,ispin),P(:,:,ispin),ERI,Vx(:,ispin))
end do
end subroutine UHF

33
src/HF/core_guess.f90 Normal file
View File

@ -0,0 +1,33 @@
subroutine core_guess(nBas,Hc,X,c)
! Core guess of the molecular orbitals for HF calculation
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
! Local variables
double precision,allocatable :: cp(:,:)
double precision,allocatable :: e(:)
! Output variables
double precision,intent(out) :: c(nBas,nBas)
! Memory allocation
allocate(cp(nBas,nBas),e(nBas))
! Core guess
cp(:,:) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:)))
call diagonalize_matrix(nBas,cp,e)
c(:,:) = matmul(X(:,:),cp(:,:))
end subroutine

37
src/HF/huckel_guess.f90
View File

@ -1,47 +1,39 @@
subroutine huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
subroutine huckel_guess(nBas,S,Hc,X,c)
! Hickel guess of the molecular orbitals for HF calculation
! Hickel guess
implicit none
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nO
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(inout):: J(nBas,nBas)
double precision,intent(inout):: K(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(inout):: cp(nBas,nBas)
double precision,intent(inout):: F(nBas,nBas)
double precision,intent(inout):: Fp(nBas,nBas)
double precision,intent(inout):: e(nBas)
double precision,intent(inout):: P(nBas,nBas)
! Local variables
integer :: mu,nu
double precision :: a
double precision,allocatable :: F(:,:)
! Output variables
double precision,intent(out) :: c(nBas,nBas)
! Memory allocation
allocate(F(nBas,nBas))
! Extended Huckel parameter
a = 1.75d0
Fp(:,:) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:)))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c(:,:) = matmul(X(:,:),cp(:,:))
call Coulomb_matrix_AO_basis(nBas,P,ERI,J)
call exchange_matrix_AO_basis(nBas,P,ERI,K)
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
! GWH approximation
do mu=1,nBas
F(mu,mu) = Hc(mu,mu)
do nu=mu+1,nBas
F(mu,nu) = 0.5d0*a*S(mu,nu)*(Hc(mu,mu) + Hc(nu,nu))
@ -50,9 +42,6 @@ subroutine huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
enddo
enddo
Fp(:,:) = matmul(transpose(X(:,:)),matmul(F(:,:),X(:,:)))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c(:,:) = matmul(X(:,:),cp(:,:))
call core_guess(nBas,F,X,c)
end subroutine

17
src/HF/mo_fock_exchange_potential.f90
View File

@ -1,4 +1,4 @@
subroutine mo_fock_exchange_potential(nBas,c,Fx,Vx)
subroutine mo_fock_exchange_potential(nBas,c,P,ERI,Vx)
! Compute the exchange potential in the MO basis
@ -9,12 +9,14 @@ subroutine mo_fock_exchange_potential(nBas,c,Fx,Vx)
integer,intent(in) :: nBas
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: Fx(nBas,nBas)
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: mu,nu
integer :: p
integer :: q
double precision,allocatable :: Fx(:,:)
! Output variables
@ -22,13 +24,18 @@ subroutine mo_fock_exchange_potential(nBas,c,Fx,Vx)
! Compute Vx
allocate(Fx(nBas,nBas))
call exchange_matrix_AO_basis(nBas,P,ERI,Fx)
Vx(:) = 0d0
do p=1,nBas
do q=1,nBas
do mu=1,nBas
do nu=1,nBas
Vx(p) = Vx(p) + c(mu,p)*Fx(mu,nu)*c(nu,p)
Vx(q) = Vx(q) + c(mu,q)*Fx(mu,nu)*c(nu,q)
end do
end do
end do
deallocate(Fx)
end subroutine mo_fock_exchange_potential

20
src/HF/mo_guess.f90
View File

@ -1,4 +1,4 @@
subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
subroutine mo_guess(nBas,guess_type,S,Hc,X,c)
! Guess of the molecular orbitals for HF calculation
@ -7,19 +7,10 @@ subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nO
integer,intent(in) :: guess_type
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(inout):: J(nBas,nBas)
double precision,intent(inout):: K(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(inout):: cp(nBas,nBas)
double precision,intent(inout):: F(nBas,nBas)
double precision,intent(inout):: Fp(nBas,nBas)
double precision,intent(inout):: e(nBas)
double precision,intent(inout):: P(nBas,nBas)
! Local variables
@ -31,14 +22,11 @@ subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
if(guess_type == 1) then
Fp = matmul(transpose(X),matmul(Hc,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c = matmul(X,cp)
call core_guess(nBas,Hc,X,c)
elseif(guess_type == 2) then
call huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
call huckel_guess(nBas,S,Hc,X,c)
elseif(guess_type == 3) then
@ -51,6 +39,4 @@ subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
endif
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
end subroutine

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

42
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
@ -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
integer,intent(in) :: ispin
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) :: 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)
@ -42,12 +50,7 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda
call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A)
! print*,'A'
! call matout(nS,nS,A)
! print*,'TA'
! call matout(nS,nS,A_BSE)
A(:,:) = A(:,:) - A_BSE(:,:)
if(BSE) A(:,:) = A(:,:) - A_BSE(:,:)
! Tamm-Dancoff approximation
@ -63,12 +66,7 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda
call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B)
! print*,'B'
! call matout(nS,nS,B)
! 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
@ -82,10 +80,6 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda
if(minval(Omega) < 0d0) &
call print_warning('You may have instabilities in linear response: A-B is not positive definite!!')
! do ia=1,nS
! if(Omega(ia) < 0d0) Omega(ia) = 0d0
! end do
call ADAt(nS,AmB,1d0*sqrt(Omega),AmBSq)
call ADAt(nS,AmB,1d0/sqrt(Omega),AmBIv)
@ -95,10 +89,6 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda
if(minval(Omega) < 0d0) &
call print_warning('You may have instabilities in linear response: negative excitations!!')
! do ia=1,nS
! if(Omega(ia) < 0d0) Omega(ia) = 0d0
! end do
Omega = sqrt(Omega)
@ -112,6 +102,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

20
src/LR/unrestricted_linear_response_C_pp.f90
View File

@ -1,5 +1,5 @@
subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,lambda,&
e,ERI_aaaa,ERI_aabb,ERI_bbbb,C_pp)
subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,&
lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,C_pp)
! Compute linear response
@ -54,7 +54,7 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP
do d=nO(2)+1,nBas-nR(2)
cd = cd + 1
C_pp(ab,cd) = (e(a,1) + e(b,2))*Kronecker_delta(a,c) &
*Kronecker_delta(b,d) + lambda*ERI_aabb(a,b,c,d)
*Kronecker_delta(b,d) + lambda*ERI_aabb(a,b,c,d)
end do
end do
end do
@ -79,9 +79,10 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP
do d=c+1,nBas-nR(1)
cd = cd + 1
C_pp(ab,cd) = (e(a,1) + e(b,1) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
+ lambda*(ERI_aaaa(a,b,c,d) - ERI_aaaa(a,b,d,c))
!write(*,*) C_pp(ab,cd)
C_pp(ab,cd) = (e(a,1) + e(b,1) - eF)*Kronecker_delta(a,c)&
*Kronecker_delta(b,d) + lambda*(ERI_aaaa(a,b,c,d) &
- ERI_aaaa(a,b,d,c))
end do
end do
end do
@ -103,12 +104,13 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP
cd = cd + 1
C_pp(ab,cd) = (e(a,2) + e(b,2) - eF)*Kronecker_delta(a,c) &
*Kronecker_delta(b,d) + lambda*(ERI_bbbb(a,b,c,d) - ERI_bbbb(a,b,d,c))
*Kronecker_delta(b,d) + lambda*(ERI_bbbb(a,b,c,d) &
- ERI_bbbb(a,b,d,c))
end do
end do
end do
end do
end do
end do
end if

16
src/LR/unrestricted_linear_response_D_pp.f90
View File

@ -1,5 +1,5 @@
subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,lambda,&
e,ERI_aaaa,ERI_aabb,ERI_bbbb,D_pp)
subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,&
lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,D_pp)
! Compute linear response
@ -55,7 +55,7 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH
do l=nC(2)+1,nO(2)
kl = kl + 1
D_pp(ij,kl) = -(e(i,1) + e(j,2))*Kronecker_delta(i,k)&
*Kronecker_delta(j,l) +lambda*ERI_aabb(i,j,k,l)
*Kronecker_delta(j,l) +lambda*ERI_aabb(i,j,k,l)
end do
end do
end do
@ -81,14 +81,15 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH
do l=k+1,nO(1)
kl = kl + 1
D_pp(ij,kl) = -(e(i,1) + e(j,1) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
+ lambda*(ERI_aaaa(i,j,k,l) - ERI_aaaa(i,j,l,k))
D_pp(ij,kl) = -(e(i,1) + e(j,1) - eF)*Kronecker_delta(i,k)&
*Kronecker_delta(j,l) + lambda*(ERI_aaaa(i,j,k,l) &
- ERI_aaaa(i,j,l,k))
end do
end do
end do
end do
end if
end if
if (ispin == 3) then
@ -104,7 +105,8 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH
kl = kl + 1
D_pp(ij,kl) = -(e(i,2) + e(j,2) - eF)*Kronecker_delta(i,k) &
*Kronecker_delta(j,l) + lambda*(ERI_bbbb(i,j,k,l) - ERI_bbbb(i,j,l,k))
*Kronecker_delta(j,l) + lambda*(ERI_bbbb(i,j,k,l) &
- ERI_bbbb(i,j,l,k))
end do
end do

53
src/LR/unrestricted_linear_response_pp.f90
View File

@ -1,6 +1,7 @@
subroutine unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt, &
nHaa,nHab,nHbb,nHt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega1,X1,Y1,Omega2,X2,Y2,&
EcRPA)
subroutine unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,&
nPt,nHaa,nHab,nHbb,nHt,lambda,e,ERI_aaaa,&
ERI_aabb,ERI_bbbb,Omega1,X1,Y1,Omega2,X2,Y2,&
EcRPA)
! Compute linear response for unrestricted formalism
@ -56,53 +57,47 @@ EcRPA)
! Memory allocation
allocate(C(nPt,nPt),B(nPt,nHt),D(nHt,nHt),M(nPt+nHt,nPt+nHt),Z(nPt+nHt,nPt+nHt)&
,Omega(nPt+nHt))
!write(*,*) 'Hello'
allocate(C(nPt,nPt),B(nPt,nHt),D(nHt,nHt),M(nPt+nHt,nPt+nHt),Z(nPt+nHt,nPt+nHt),&
Omega(nPt+nHt))
! Build C, B and D matrices for the pp channel
call unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,lambda,&
e,ERI_aaaa,ERI_aabb,ERI_bbbb,C)
call unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,&
lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,C)
call unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,nHaa,&
nHab,nHbb,nHt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,B)
call unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,nHaa,&
nHab,nHbb,nHt,lambda,ERI_aaaa,ERI_aabb,&
ERI_bbbb,B)
!call matout(nPt,nHt,B)
!write(*,*) 'Hello'
call unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,lambda,&
e,ERI_aaaa,ERI_aabb,ERI_bbbb,D)
call unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,&
lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,D)
!call matout(nHt,nHt,D)
!write(*,*) 'Hello'
! Diagonal blocks
M( 1:nPt , 1:nPt) = + C(1:nPt,1:nPt)
M(nPt+1:nPt+nHt,nPt+1:nPt+nHt) = - D(1:nHt,1:nHt)
M( 1:nPt , 1:nPt) = + C(1:nPt,1:nPt)
M(nPt+1:nPt+nHt,nPt+1:nPt+nHt) = - D(1:nHt,1:nHt)
! Off-diagonal blocks
! Off-diagonal blocks
M( 1:nPt ,nPt+1:nHt+nPt) = - B(1:nPt,1:nHt)
M(nPt+1:nHt+nPt, 1:nPt) = + transpose(B(1:nPt,1:nHt))
!call matout(nPt+nHt,nPt+nHt,M)
M( 1:nPt ,nPt+1:nHt+nPt) = - B(1:nPt,1:nHt)
M(nPt+1:nHt+nPt, 1:nPt) = + transpose(B(1:nPt,1:nHt))
! Diagonalize the p-h matrix
if(nHt+nPt > 0) call diagonalize_general_matrix(nHt+nPt,M,Omega,Z)
! Split the various quantities in p-p and h-h parts
! Split the various quantities in p-p and h-h parts
call sort_ppRPA(nHt,nPt,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),&
Y2(:,:))
call sort_ppRPA(nHt,nPt,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),&
Y2(:,:))
! Compute the RPA correlation energy
EcRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nPt,C(:,:)) - trace_matrix(nHt,D(:,:)) )
EcRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nPt,C(:,:)) &
- trace_matrix(nHt,D(:,:)) )
EcRPA1 = +sum(Omega1(:)) - trace_matrix(nPt,C(:,:))
EcRPA2 = -sum(Omega2(:)) - trace_matrix(nHt,D(:,:))
if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) &
print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!'
end subroutine unrestricted_linear_response_pp

27
src/QuAcK/QuAcK.f90
View File

@ -111,7 +111,7 @@ program QuAcK
double precision :: start_Bas ,end_Bas ,t_Bas
integer :: maxSCF_HF,n_diis_HF
double precision :: thresh_HF
double precision :: thresh_HF,level_shift
logical :: DIIS_HF,guess_type,ortho_type,mix
integer :: maxSCF_CC,n_diis_CC
@ -180,15 +180,15 @@ program QuAcK
! Read options for methods
call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,dostab, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
TDA,singlet,triplet,spin_conserved,spin_flip, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW, &
COHSEX,SOSEX,TDA_W,G0W,GW0, &
maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
BSE,dBSE,dTDA,evDyn, &
call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,level_shift,dostab, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
TDA,singlet,triplet,spin_conserved,spin_flip, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW, &
COHSEX,SOSEX,TDA_W,G0W,GW0, &
maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
BSE,dBSE,dTDA,evDyn, &
nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift)
! Weird stuff
@ -292,7 +292,7 @@ program QuAcK
end if
call cpu_time(start_HF)
call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc, &
call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,F_AO,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc)
call cpu_time(end_HF)
@ -312,7 +312,7 @@ program QuAcK
unrestricted = .true.
call cpu_time(start_HF)
call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,nNuc,ZNuc,rNuc,ENuc, &
call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,eHF,cHF,PHF,Vxc)
call cpu_time(end_HF)
@ -332,7 +332,7 @@ program QuAcK
unrestricted = .true.
call cpu_time(start_KS)
call eDFT(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC, &
call eDFT(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC, &
nO,nV,nR,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, &
max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int_AO,EUHF,eHF,cHF,PHF,Vxc)
@ -1159,6 +1159,7 @@ program QuAcK
else
! call soG0T0(eta_GT,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF)
call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet, &
linGT,eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO, &
PHF,cHF,eHF,Vxc,eG0T0)

28
src/QuAcK/read_options.f90
View File

@ -1,12 +1,12 @@
subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,dostab, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
TDA,singlet,triplet,spin_conserved,spin_flip, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW, &
COHSEX,SOSEX,TDA_W,G0W,GW0, &
maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
BSE,dBSE,dTDA,evDyn, &
subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,level_shift,dostab, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
TDA,singlet,triplet,spin_conserved,spin_flip, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW, &
COHSEX,SOSEX,TDA_W,G0W,GW0, &
maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
BSE,dBSE,dTDA,evDyn, &
nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift)
! Read desired methods
@ -22,6 +22,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
integer,intent(out) :: guess_type
integer,intent(out) :: ortho_type
logical,intent(out) :: mix
double precision,intent(out) :: level_shift
logical,intent(out) :: dostab
integer,intent(out) :: maxSCF_CC
@ -100,14 +101,15 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
guess_type = 1
ortho_type = 1
mix = .false.
level_shift = 0d0
dostab = .false.
read(1,*)
read(1,*) maxSCF_HF,thresh_HF,answer1,n_diis_HF,guess_type,ortho_type,answer2,answer3
read(1,*) maxSCF_HF,thresh_HF,answer1,n_diis_HF,guess_type,ortho_type,answer2,level_shift,answer3
if(answer1 == 'T') DIIS_HF = .true.
if(answer2 == 'T') mix = .true.
if(answer3 == 'T') dostab = .true.
if(answer1 == 'T') DIIS_HF = .true.
if(answer2 == 'T') mix = .true.
if(answer3 == 'T') dostab = .true.
if(.not.DIIS_HF) n_diis_HF = 1

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 :: 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))

94
src/RPA/ACFDT_Tmatrix.f90
View File

@ -1,5 +1,6 @@
subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
ERI,eT,eGT,EcAC)
nOOs,nVVs,nOOt,nVVt,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t, &
Omega2t,X2t,Y2t,rho1t,rho2t,ERI,eT,eGT,EcAC)
! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem for the T-matrix
@ -26,6 +27,28 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
integer,intent(in) :: nR
integer,intent(in) :: nS
integer,intent(in) :: nOOs
integer,intent(in) :: nOOt
integer,intent(in) :: nVVs
integer,intent(in) :: nVVt
double precision,intent(in) :: Omega1s(nVVs)
double precision,intent(in) :: X1s(nVVs,nVVs)
double precision,intent(in) :: Y1s(nOOs,nVVs)
double precision,intent(in) :: Omega2s(nOOs)
double precision,intent(in) :: X2s(nVVs,nOOs)
double precision,intent(in) :: Y2s(nOOs,nOOs)
double precision,intent(in) :: rho1s(nBas,nBas,nVVs)
double precision,intent(in) :: rho2s(nBas,nBas,nOOs)
double precision,intent(in) :: Omega1t(nVVt)
double precision,intent(in) :: X1t(nVVt,nVVt)
double precision,intent(in) :: Y1t(nOOt,nVVt)
double precision,intent(in) :: Omega2t(nOOt)
double precision,intent(in) :: X2t(nVVt,nOOt)
double precision,intent(in) :: Y2t(nOOt,nOOt)
double precision,intent(in) :: rho1t(nBas,nBas,nVVt)
double precision,intent(in) :: rho2t(nBas,nBas,nOOt)
double precision,intent(in) :: eT(nBas)
double precision,intent(in) :: eGT(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
@ -39,46 +62,23 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
double precision :: lambda
double precision,allocatable :: Ec(:,:)
integer :: nOOs,nOOt
integer :: nVVs,nVVt
double precision :: EcRPA(nspin)
double precision,allocatable :: TA(:,:)
double precision,allocatable :: TB(:,:)
double precision,allocatable :: TAs(:,:)
double precision,allocatable :: TBs(:,:)
double precision,allocatable :: TAt(:,:)
double precision,allocatable :: TBt(:,:)
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
double precision,allocatable :: Omega1s(:),Omega1t(:)
double precision,allocatable :: X1s(:,:),X1t(:,:)
double precision,allocatable :: Y1s(:,:),Y1t(:,:)
double precision,allocatable :: rho1s(:,:,:),rho1t(:,:,:)
double precision,allocatable :: Omega2s(:),Omega2t(:)
double precision,allocatable :: X2s(:,:),X2t(:,:)
double precision,allocatable :: Y2s(:,:),Y2t(:,:)
double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:)
! Output variables
double precision,intent(out) :: EcAC(nspin)
! Useful quantities
nOOs = nO*nO
nVVs = nV*nV
nOOt = nO*(nO-1)/2
nVVt = nV*(nV-1)/2
! Memory allocation
allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), &
Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), &
rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), &
Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), &
Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt))
allocate(TA(nS,nS),TB(nS,nS),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
allocate(TAs(nS,nS),TBs(nS,nS),TAt(nS,nS),TBt(nS,nS), &
Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
allocate(Ec(nAC,nspin))
! Antisymmetrized kernel version
@ -113,11 +113,6 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
lambda = rAC(iAC)
! Initialize T matrix
TA(:,:) = 0d0
TB(:,:) = 0d0
if(doXBS) then
isp_T = 1
@ -128,8 +123,8 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Omega1s,rho1s,Omega2s,rho2s,TAs)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Omega1s,rho1s,Omega2s,rho2s,TBs)
isp_T = 2
iblock = 4
@ -139,13 +134,13 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Omega1t,rho1t,Omega2t,rho2t,TAt)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Omega1t,rho1t,Omega2t,rho2t,TBt)
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,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TAt+TAs,TBt+TBs, &
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))
@ -183,11 +178,6 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
lambda = rAC(iAC)
! Initialize T matrix
TA(:,:) = 0d0
TB(:,:) = 0d0
if(doXBS) then
isp_T = 1
@ -198,8 +188,8 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Omega1s,rho1s,Omega2s,rho2s,TAs)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Omega1s,rho1s,Omega2s,rho2s,TBs)
isp_T = 2
iblock = 4
@ -209,13 +199,13 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Omega1t,rho1t,Omega2t,rho2t,TAt)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Omega1t,rho1t,Omega2t,rho2t,TBt)
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,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TAt-TAs,TBt-TBs, &
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))

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 :: 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))

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 :: 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))

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 :: 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))

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 :: 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))

78
src/RPA/ppURPA.f90
View File

@ -25,7 +25,7 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH
! Local variables
integer :: ispin
integer :: ispin,iblock
integer :: nPaa,nPbb,nPab,nP_sc,nP_sf
integer :: nHaa,nHbb,nHab,nH_sc,nH_sf
double precision,allocatable :: Omega1sc(:),Omega1sf(:)
@ -56,25 +56,28 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH
if(spin_conserved) then
ispin = 1
iblock = 1
!spin-conserved quantities
!Spin-conserved quantities
nPab = nV(1)*nV(2)
nHab = nO(1)*nO(2)
nPab = nV(1)*nV(2)
nHab = nO(1)*nO(2)
nP_sc = nPab
nH_sc = nHab
nP_sc = nPab
nH_sc = nHab
! Memory allocation
allocate(Omega1sc(nP_sc),X1sc(nP_sc,nP_sc),Y1sc(nH_sc,nP_sc), &
Omega2sc(nH_sc),X2sc(nP_sc,nH_sc),Y2sc(nH_sc,nH_sc))
allocate(Omega1sc(nP_sc),X1sc(nP_sc,nP_sc),Y1sc(nH_sc,nP_sc), &
Omega2sc(nH_sc),X2sc(nP_sc,nH_sc),Y2sc(nH_sc,nH_sc))
call unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nP_sc,nHaa,nHab,nHbb,nH_sc,1d0, &
e,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega1sc,X1sc,Y1sc,Omega2sc,X2sc,Y2sc,Ec_ppURPA(ispin))
call print_excitation('pp-RPA (N+2)',5,nP_sc,Omega1sc)
call print_excitation('pp-RPA (N-2)',5,nH_sc,Omega2sc)
call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, &
nP_sc,nHaa,nHab,nHbb,nH_sc,1d0,e,ERI_aaaa, &
ERI_aabb,ERI_bbbb,Omega1sc,X1sc,Y1sc, &
Omega2sc,X2sc,Y2sc,Ec_ppURPA(ispin))
call print_excitation('pp-RPA (N+2)',5,nP_sc,Omega1sc)
call print_excitation('pp-RPA (N-2)',5,nH_sc,Omega2sc)
endif
@ -83,40 +86,45 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH
if(spin_flip) then
ispin = 2
iblock = 2
!spin-flip quantities
!Spin-flip quantities
nPaa = nV(1)*(nV(1)-1)/2
nPbb = nV(2)*(nV(2)-1)/2
nPaa = nV(1)*(nV(1)-1)/2
nPbb = nV(2)*(nV(2)-1)/2
nP_sf = nPaa
nP_sf = nPaa
nHaa = nO(1)*(nO(1)-1)/2
nHbb = nO(2)*(nO(2)-1)/2
nHaa = nO(1)*(nO(1)-1)/2
nHbb = nO(2)*(nO(2)-1)/2
nH_sf = nHaa
nH_sf = nHaa
allocate(Omega1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), &
Omega2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf))
allocate(Omega1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), &
Omega2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf))
call unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nP_sf, &
nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega1sf,X1sf,Y1sf,Omega2sf,X2sf,Y2sf,&
Ec_ppURPA(ispin))
call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, &
nP_sf,nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa, &
ERI_aabb,ERI_bbbb,Omega1sf,X1sf,Y1sf, &
Omega2sf,X2sf,Y2sf,Ec_ppURPA(ispin))
ispin = 3
deallocate(Omega1sf,X1sf,Y1sf,Omega2sf,X2sf,Y2sf)
nP_sf = nPbb
nH_sf = nHbb
iblock = 3
!allocate(Omega1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), &
! Omega2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf))
nP_sf = nPbb
nH_sf = nHbb
call unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nP_sf, &
nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega1sf,X1sf,Y1sf,Omega2sf,X2sf,Y2sf,&
Ec_ppURPA(ispin))
allocate(Omega1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), &
Omega2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf))
call print_excitation('pp-RPA (N+2)',6,nP_sf,Omega1sf)
call print_excitation('pp-RPA (N-2)',6,nH_sf,Omega2sf)
call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,&
nP_sf,nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa,&
ERI_aabb,ERI_bbbb,Omega1sf,X1sf,Y1sf,&
Omega2sf,X2sf,Y2sf,Ec_ppURPA(ispin))
call print_excitation('pp-RPA (N+2)',6,nP_sf,Omega1sf)
call print_excitation('pp-RPA (N-2)',6,nH_sf,Omega2sf)
endif

81
src/eDFT/eDFT_UKS.f90 → src/eDFT/UKS.f90
View File

@ -1,5 +1,5 @@
subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, &
nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,c,Pw,Vxc)
subroutine UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix,level_shift, &
nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,c,Pw,Vxc)
! Perform unrestricted Kohn-Sham calculation for ensembles
@ -20,6 +20,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
logical,intent(in) :: mix
double precision,intent(in) :: level_shift
double precision,intent(in) :: thresh
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
@ -45,10 +46,11 @@ 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
integer :: nSCF
integer :: nBasSq
integer :: n_diis
integer :: nO(nspin)
double precision :: conv
integer :: nO(nspin,nEns)
double precision :: Conv
double precision :: rcond(nspin)
double precision :: ET(nspin)
double precision :: EV(nspin)
@ -117,34 +119,22 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
! Guess coefficients and eigenvalues
nO(:) = 0
do ispin=1,nspin
nO(ispin) = int(sum(occnum(:,ispin,1)))
nO(:,:) = 0
do iEns=1,nEns
do ispin=1,nspin
nO(ispin,iEns) = int(sum(occnum(:,ispin,iEns)))
end do
end do
if(guess_type == 1) then
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
do ispin=1,nspin
call mo_guess(nBas,guess_type,S,Hc,X,c(:,:,ispin))
end do
! Mix guess to enforce symmetry breaking
if(mix) call mix_guess(nBas,nO,c)
else if(guess_type == 2) then
do ispin=1,nspin
call random_number(F(:,:,ispin))
end do
else
print*,'Wrong guess option'
stop
! Mix guess for UHF solution in singlet states
if(mix) then
write(*,*) '!! guess mixing disabled in UKS !!'
write(*,*)
end if
! Initialization
@ -175,7 +165,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
'|','#','|','E(KS)','|','Ex(KS)','|','Ec(KS)','|','Conv','|','nEl','|'
write(*,*)'------------------------------------------------------------------------------------------'
do while(conv > thresh .and. nSCF < maxSCF)
do while(Conv > thresh .and. nSCF < maxSCF)
! Increment
@ -264,19 +254,34 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
err(:,:,ispin) = matmul(F(:,:,ispin),matmul(Pw(:,:,ispin),S(:,:))) - matmul(matmul(S(:,:),Pw(:,:,ispin)),F(:,:,ispin))
end do
if(nSCF > 1) conv = maxval(abs(err(:,:,:)))
if(nSCF > 1) Conv = maxval(abs(err(:,:,:)))
! 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))
else
n_diis = 0
end if
end do
! Level-shifting
if(level_shift > 0d0 .and. Conv > thresh) then
do ispin=1,nspin
call level_shifting(level_shift,nBas,maxval(nO(ispin,:)),S,c,F(:,:,ispin))
end do
else
n_diis = 0
end if
end if
! Transform Fock matrix in orthogonal basis
@ -342,7 +347,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
! Dump results
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
'|',nSCF,'|',Ew + ENuc,'|',sum(Ex(:)),'|',sum(Ec(:)),'|',conv,'|',sum(nEl(:)),'|'
'|',nSCF,'|',Ew + ENuc,'|',sum(Ex(:)),'|',sum(Ec(:)),'|',Conv,'|',sum(nEl(:)),'|'
end do
write(*,*)'------------------------------------------------------------------------------------------'
@ -384,4 +389,4 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
call individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, &
AO,dAO,T,V,ERI,ENuc,eKS,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,occnum,Cx_choice,doNcentered,Ew)
end subroutine eDFT_UKS
end subroutine UKS

13
src/eDFT/eDFT.f90
View File

@ -1,4 +1,4 @@
subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC,nO,nV,nR, &
subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC,nO,nV,nR, &
nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, &
max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI,dipole_int,Ew,eKS,cKS,PKS,Vxc)
@ -15,6 +15,7 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
logical,intent(in) :: mix
logical,intent(in) :: level_shift
double precision,intent(in) :: thresh
integer,intent(in) :: nNuc
@ -165,8 +166,9 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n
end do
call cpu_time(start_KS)
call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC(1:nCC,1:nEns-1),nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, &
nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc)
call UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC(1:nCC,1:nEns-1),nGrid,weight, &
maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc)
call cpu_time(end_KS)
t_KS = end_KS - start_KS
@ -182,8 +184,9 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n
if(method == 'eDFT-UKS') then
call cpu_time(start_KS)
call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC(1:nCC,1:nEns-1),nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, &
nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc)
call UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC(1:nCC,1:nEns-1),nGrid,weight, &
maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc)
call cpu_time(end_KS)
t_KS = end_KS - start_KS

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

@ -0,0 +1,37 @@
subroutine level_shifting(level_shift,nBas,nO,S,c,F)
! Perform level-shifting on the Fock matrix
implicit none
! Input variables
double precision,intent(in) :: level_shift
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(:,:)
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) + level_shift
end do
Sc(:,:) = matmul(S,c)
F(:,:) = matmul(Sc,matmul(F_MO,transpose(Sc)))
end subroutine level_shifting