LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-07-26 12:47:35 +02:00
Code Releases Activity

fix conflicts

Browse Source
This commit is contained in:
Clotilde Marut 2022-01-06 11:20:49 +01:00
commit 3c39e0aa3a
117 changed files with 2514 additions and 769 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
include/parameters.h
View File

@ -20,4 +20,3 @@
double precision,parameter :: CxLDA = - (3d0/4d0)*(3d0/pi)**(1d0/3d0)
double precision,parameter :: CxLSDA = - (3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)

2
input/methods
View File

@ -9,7 +9,7 @@
# CIS* CIS(D) CID CISD FCI
F F F F F
# RPA* RPAx* crRPA ppRPA
F F F F
F F F F
# G0F2* evGF2* qsGF2* G0F3 evGF3
F F F F F
# G0W0* evGW* qsGW* ufG0W0 ufGW

16
input/options
View File

@ -1,18 +1,20 @@
# HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess stability
1024 0.00001 T 5 1 1 T F
1024 0.00001 T 5 1 1 T F
# MP:
# CC: maxSCF thresh DIIS n_diis
64 0.00001 T 5
# spin: TDA singlet triplet spin_conserved spin_flip
F T T T T
# GF: maxSCF thresh DIIS n_diis lin eta renorm
256 0.00001 T 5 T 0.00367493 3
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0
256 0.00001 T 5 T 0.0 F F F F F
F T F T T
# GF: maxSCF thresh DIIS n_diis lin eta renorm reg
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
# 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
# BSE: BSE dBSE dTDA evDyn
F F F F
T T T F
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

2
mol/h2.xyz
View File

@ -1,4 +1,4 @@
2
H 0. 0. 0.
H 0. 0. 0.5
H 0. 0. 2.000000

2
src/CC/CCSD_Ec_nc.f90
View File

@ -1,6 +1,6 @@
subroutine CCSD_Ec_nc(nO,nV,t1,t2,Fov,OOVV,EcCCSD)
! Compute the CCSD correlatio energy in non-conanical form
! Compute the CCSD correlatio energy in non-canonical form
implicit none

218
src/CI/CID.f90 Normal file
View File

@ -0,0 +1,218 @@
subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIin,Fin,E0)
! Perform configuration interaction with doubles
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
integer,intent(in) :: nBasin
integer,intent(in) :: nCin
integer,intent(in) :: nOin
integer,intent(in) :: nVin
integer,intent(in) :: nRin
double precision,intent(in) :: Fin(nBasin,nBasin)
double precision,intent(in) :: ERIin(nBasin,nBasin,nBasin,nBasin)
double precision,intent(in) :: E0
! Local variables
integer :: nBas
integer :: nC
integer :: nO
integer :: nV
integer :: nR
double precision,allocatable :: F(:,:)
double precision,allocatable :: sERI(:,:,:,:)
double precision,allocatable :: ERI(:,:,:,:)
logical :: dump_trans = .false.
integer :: i,j,k,l
integer :: a,b,c,d
integer :: ia,kc,iajb,kcld
integer :: ishift,jshift
integer :: ispin
integer :: nD
integer :: nH
integer :: maxH
double precision,external :: Kronecker_delta
double precision,allocatable :: H(:,:)
double precision,allocatable :: ECID(:)
double precision :: tmp
! Hello world
write(*,*)
write(*,*)'******************************************************'
write(*,*)'| Configuration Interaction with Singles and Doubles |'
write(*,*)'******************************************************'
write(*,*)
! Spatial to spin orbitals
nBas = 2*nBasin
nC = 2*nCin
nO = 2*nOin
nV = 2*nVin
nR = 2*nRin
allocate(F(nBas,nBas),sERI(nBas,nBas,nBas,nBas))
call spatial_to_spin_fock(nBasin,Fin,nBas,F)
call spatial_to_spin_ERI(nBasin,ERIin,nBas,sERI)
! Antysymmetrize ERIs
allocate(ERI(nBas,nBas,nBas,nBas))
call antisymmetrize_ERI(2,nBas,sERI,ERI)
deallocate(sERI)
! Compute CID matrix
nD = (nO - nC)*(nO - nC - 1)/2*(nV - nR)*(nV - nR - 1)/2
nH = 1 + nD
write(*,*) 'nD = ',nD
write(*,*) 'nH = ',nH
write(*,*)
maxH = min(nH,21)
! Memory allocation
allocate(H(nH,nH),ECID(nH))
! 00 block
ishift = 0
jshift = 0
H(ishift+1,jshift+1) = E0
print*,'00 block done...'
! 0D blocks
ishift = 0
jshift = 1
iajb = 0
do i=nC+1,nO
do a=1,nV-nR
do j=i+1,nO
do b=a+1,nV-nR
iajb = iajb + 1
tmp = ERI(i,j,nO+a,nO+b)
H(ishift+1,jshift+iajb) = tmp
H(jshift+iajb,ishift+1) = tmp
end do
end do
end do
end do
print*,'0D blocks done...'
! DD block
ishift = 1
jshift = 1
iajb = 0
do i=nC+1,nO
do a=1,nV-nR
do j=i+1,nO
do b=a+1,nV-nR
iajb = iajb + 1
kcld = 0
do k=nC+1,nO
do c=1,nV-nR
do l=k+1,nO
do d=c+1,nV-nR
kcld = kcld + 1
tmp = &
E0*Kronecker_delta(i,k)*Kronecker_delta(j,l)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
+ F(l,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,k) &
- F(l,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,k) &
- F(k,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,l) &
+ F(k,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,l) &
- F(l,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,k) &
+ F(l,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,k) &
+ F(k,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,l) &
- F(k,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,l) &
+ F(nO+a,nO+d)*Kronecker_delta(b,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
- F(nO+a,nO+c)*Kronecker_delta(b,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
- F(nO+a,nO+d)*Kronecker_delta(b,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
+ F(nO+a,nO+c)*Kronecker_delta(b,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
- F(nO+b,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+ F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+ F(nO+b,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
- F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
- ERI(k,l,i,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c) &
+ ERI(k,l,i,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
+ ERI(nO+a,l,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,k) &
- ERI(nO+a,l,nO+c,j)*Kronecker_delta(b,d)*Kronecker_delta(i,k) &
- ERI(nO+a,k,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,l) &
+ ERI(nO+a,k,nO+c,j)*Kronecker_delta(b,d)*Kronecker_delta(i,l) &
- ERI(nO+a,l,nO+d,i)*Kronecker_delta(b,c)*Kronecker_delta(j,k) &
+ ERI(nO+a,l,nO+c,i)*Kronecker_delta(b,d)*Kronecker_delta(j,k) &
+ ERI(nO+a,k,nO+d,i)*Kronecker_delta(b,c)*Kronecker_delta(j,l) &
- ERI(nO+a,k,nO+c,i)*Kronecker_delta(b,d)*Kronecker_delta(j,l) &
- ERI(nO+b,l,nO+d,j)*Kronecker_delta(a,c)*Kronecker_delta(i,k) &
+ ERI(nO+b,l,nO+c,j)*Kronecker_delta(a,d)*Kronecker_delta(i,k) &
+ ERI(nO+b,k,nO+d,j)*Kronecker_delta(a,c)*Kronecker_delta(i,l) &
- ERI(nO+b,k,nO+c,j)*Kronecker_delta(a,d)*Kronecker_delta(i,l) &
+ ERI(nO+b,l,nO+d,i)*Kronecker_delta(a,c)*Kronecker_delta(j,k) &
- ERI(nO+b,l,nO+c,i)*Kronecker_delta(a,d)*Kronecker_delta(j,k) &
- ERI(nO+b,k,nO+d,i)*Kronecker_delta(a,c)*Kronecker_delta(j,l) &
+ ERI(nO+b,k,nO+c,i)*Kronecker_delta(a,d)*Kronecker_delta(j,l) &
- ERI(nO+a,nO+b,nO+c,nO+d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+ ERI(nO+a,nO+b,nO+c,nO+d)*Kronecker_delta(i,k)*Kronecker_delta(j,l)
H(ishift+iajb,jshift+kcld) = tmp
end do
end do
end do
end do
end do
end do
end do
end do
print*,'DD block done...'
write(*,*)
write(*,*) 'Diagonalizing CID matrix...'
write(*,*)
call diagonalize_matrix(nH,H,ECID)
print*,'CID energies (au)'
call matout(maxH,1,ECID)
write(*,*)
print*,'CID excitation energies (eV)'
call matout(maxH-1,1,(ECID(2:maxH)-ECID(1))*HaToeV)
write(*,*)
if(dump_trans) then
print*,'Singlet CID transition vectors'
call matout(nH,nH,H)
write(*,*)
endif
end subroutine CID

342
src/CI/CISD.f90
View File

@ -1,4 +1,4 @@
subroutine CISD(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI,eHF)
subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIin,Fin,E0)
! Perform configuration interaction with singles and doubles
@ -9,23 +9,42 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI,eHF)
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
integer,intent(in) :: nBasin
integer,intent(in) :: nCin
integer,intent(in) :: nOin
integer,intent(in) :: nVin
integer,intent(in) :: nRin
double precision,intent(in) :: Fin(nBasin,nBasin)
double precision,intent(in) :: ERIin(nBasin,nBasin,nBasin,nBasin)
double precision,intent(in) :: E0
! Local variables
integer :: nBas
integer :: nC
integer :: nO
integer :: nV
integer :: nR
double precision,allocatable :: F(:,:)
double precision,allocatable :: sERI(:,:,:,:)
double precision,allocatable :: ERI(:,:,:,:)
logical :: dump_trans = .false.
integer :: i,j,k,l
integer :: a,b,c,d
integer :: ia,jb,iajb,kcld
integer :: ia,kc,iajb,kcld
integer :: ishift,jshift
integer :: ispin
integer :: nS
integer :: nD
integer :: nSD
integer :: nH
integer :: maxH
double precision,external :: Kronecker_delta
double precision,allocatable :: H(:,:),Omega(:)
double precision,allocatable :: H(:,:)
double precision,allocatable :: ECISD(:)
double precision :: tmp
! Hello world
@ -35,170 +54,257 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI,eHF)
write(*,*)'******************************************************'
write(*,*)
! Compute CIS matrix
! Spatial to spin orbitals
if(singlet_manifold) then
nBas = 2*nBasin
nC = 2*nCin
nO = 2*nOin
nV = 2*nVin
nR = 2*nRin
ispin = 1
allocate(F(nBas,nBas),sERI(nBas,nBas,nBas,nBas))
! Dimensions
call spatial_to_spin_fock(nBasin,Fin,nBas,F)
call spatial_to_spin_ERI(nBasin,ERIin,nBas,sERI)
nS = (nO - nC)*(nV - nR)
nD = (nO - nC)*(nO - nC + 1)/2*(nV - nR)*(nV - nR + 1)/2
nSD = 1 + nS + nD
! Antysymmetrize ERIs
print*,'nS = ',nS
print*,'nD = ',nD
print*,'nSD = ',nSD
allocate(ERI(nBas,nBas,nBas,nBas))
! Memory allocation
call antisymmetrize_ERI(2,nBas,sERI,ERI)
allocate(H(nSD,nSD),Omega(nSD))
! 0D block
deallocate(sERI)
ishift = 0
jshift = 1 + nS
! Compute CISD matrix
iajb = 0
nS = (nO - nC)*(nV - nR)
nD = (nO - nC)*(nO - nC - 1)/2*(nV - nR)*(nV - nR - 1)/2
nH = 1 + nS + nD
do i=nC+1,nO
do a=1,nV-nR
do j=i,nO
do b=a,nV-nR
write(*,*) 'nS = ',nS
write(*,*) 'nD = ',nD
write(*,*) 'nH = ',nH
write(*,*)
iajb = iajb + 1
H(ishift+1,jshift+iajb) = ERI(i,j,nO+a,nO+b)
maxH = min(nH,21)
! Memory allocation
allocate(H(nH,nH),ECISD(nH))
! 00 block
ishift = 0
jshift = 0
H(ishift+1,jshift+1) = E0
print*,'00 block done...'
! 0S blocks
ishift = 0
jshift = 1
ia = 0
do i=nC+1,nO
do a=1,nV-nR
ia = ia + 1
tmp = F(i,nO+a)
H(ishift+1,jshift+ia) = tmp
H(jshift+ia,ishift+1) = tmp
end do
end do
print*,'0S blocks done...'
! 0D blocks
ishift = 0
jshift = 1 + nS
iajb = 0
do i=nC+1,nO
do a=1,nV-nR
do j=i+1,nO
do b=a+1,nV-nR
iajb = iajb + 1
tmp = ERI(i,j,nO+a,nO+b)
H(ishift+1,jshift+iajb) = tmp
H(jshift+iajb,ishift+1) = tmp
end do
end do
end do
end do
! SS block
end do
print*,'0D blocks done...'
ishift = 1
jshift = 1
! SS block
ia = 0
jb = 0
ishift = 1
jshift = 1
do i=nC+1,nO
do a=1,nV-nR
ia = 0
do i=nC+1,nO
do a=1,nV-nR
ia = ia + 1
ia = ia + 1
kc = 0
do k=nC+1,nO
do c=1,nV-nR
do j=nC+1,nO
do b=1,nV-nR
kc = kc + 1
tmp = E0*Kronecker_delta(i,k)*Kronecker_delta(a,c) &
- F(i,k)*Kronecker_delta(a,c) &
+ F(nO+a,nO+c)*Kronecker_delta(i,k) &
- ERI(nO+a,k,nO+c,i)
jb = jb + 1
H(ishift+ia,jshift+kc) = tmp
H(ishift+ia,jshift+jb) &
= Kronecker_delta(i,j)*Kronecker_delta(a,b)*(eHF(nO+a) - eHF(i)) &
+ ERI(nO+a,j,i,nO+b) - ERI(nO+a,j,nO+b,i)
end do
end do
end do
end do
! SD block
end do
ishift = 1
jshift = 1 + nS
print*,'SS block done...'
ia = 0
kcld = 0
! SD blocks
do i=nC+1,nO
do a=1,nV-nR
ishift = 1
jshift = 1 + nS
ia = ia + 1
ia = 0
do i=nC+1,nO
do a=1,nV-nR
do k=nC+1,nO
do c=1,nV-nR
do l=k,nO
do d=c,nV-nR
ia = ia + 1
kcld = 0
kcld = kcld + 1
do k=nC+1,nO
do c=1,nV-nR
do l=k+1,nO
do d=c+1,nV-nR
H(ishift+ia,jshift+kcld) &
= Kronecker_delta(i,k)*(ERI(nO+a,l,nO+c,nO+d) - ERI(nO+a,l,nO+d,nO+c)) &
- Kronecker_delta(i,l)*(ERI(nO+a,k,nO+c,nO+d) - ERI(nO+a,k,nO+d,nO+c)) &
- Kronecker_delta(a,c)*(ERI(k,l,i,nO+d) - ERI(k,l,nO+d,i)) &
+ Kronecker_delta(a,d)*(ERI(k,l,i,nO+c) - ERI(k,l,nO+c,i))
kcld = kcld + 1
tmp = - F(l,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,k) &
+ F(l,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,k) &
- F(k,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,l) &
+ F(k,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,l) &
- ERI(k,l,nO+d,i)*Kronecker_delta(a,c) &
+ ERI(k,l,nO+c,i)*Kronecker_delta(a,d) &
- ERI(nO+a,l,nO+c,nO+d)*Kronecker_delta(i,k) &
+ ERI(nO+a,k,nO+c,nO+d)*Kronecker_delta(i,l)
H(ishift+ia,jshift+kcld) = tmp
H(jshift+kcld,ishift+ia) = tmp
end do
end do
end do
end do
end do
end do
end do
! DD block
print*,'SD blocks done...'
ishift = 1 + nS
jshift = 1 + nS
! DD block
iajb = 0
kcld = 0
ishift = 1 + nS
jshift = 1 + nS
do i=nC+1,nO
do a=1,nV-nR
do j=i,nO
do b=a,nV-nR
iajb = 0
do i=nC+1,nO
do a=1,nV-nR
do j=i+1,nO
do b=a+1,nV-nR
iajb = iajb + 1
iajb = iajb + 1
do k=nC+1,nO
do c=1,nV-nR
do l=k,nO
do d=c,nV-nR
kcld = kcld + 1
! H(ishift+iajb,jshift+kcld) &
! = Kronecker_delta(i,k)*(ERI(a,l,c,d) - ERI(a,l,d,c)) &
! - Kronecker_delta(i,l)*(ERI(a,k,c,d) - ERI(a,k,d,c)) &
! - Kronecker_delta(a,c)*(ERI(k,l,i,d) - ERI(k,l,d,i)) &
! + Kronecker_delta(a,d)*(ERI(k,l,i,c) - ERI(k,l,c,i))
end do
kcld = 0
do k=nC+1,nO
do c=1,nV-nR
do l=k+1,nO
do d=c+1,nV-nR
kcld = kcld + 1
tmp = &
E0*Kronecker_delta(i,k)*Kronecker_delta(j,l)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
+ F(l,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,k) &
- F(l,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,k) &
- F(k,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,l) &
+ F(k,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,l) &
- F(l,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,k) &
+ F(l,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,k) &
+ F(k,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,l) &
- F(k,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,l) &
+ F(nO+a,nO+d)*Kronecker_delta(b,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
- F(nO+a,nO+c)*Kronecker_delta(b,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
- F(nO+a,nO+d)*Kronecker_delta(b,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
+ F(nO+a,nO+c)*Kronecker_delta(b,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
- F(nO+b,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+ F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+ F(nO+b,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
- F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
- ERI(k,l,i,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c) &
+ ERI(k,l,i,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
+ ERI(nO+a,l,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,k) &
- ERI(nO+a,l,nO+c,j)*Kronecker_delta(b,d)*Kronecker_delta(i,k) &
- ERI(nO+a,k,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,l) &
+ ERI(nO+a,k,nO+c,j)*Kronecker_delta(b,d)*Kronecker_delta(i,l) &
- ERI(nO+a,l,nO+d,i)*Kronecker_delta(b,c)*Kronecker_delta(j,k) &
+ ERI(nO+a,l,nO+c,i)*Kronecker_delta(b,d)*Kronecker_delta(j,k) &
+ ERI(nO+a,k,nO+d,i)*Kronecker_delta(b,c)*Kronecker_delta(j,l) &
- ERI(nO+a,k,nO+c,i)*Kronecker_delta(b,d)*Kronecker_delta(j,l) &
- ERI(nO+b,l,nO+d,j)*Kronecker_delta(a,c)*Kronecker_delta(i,k) &
+ ERI(nO+b,l,nO+c,j)*Kronecker_delta(a,d)*Kronecker_delta(i,k) &
+ ERI(nO+b,k,nO+d,j)*Kronecker_delta(a,c)*Kronecker_delta(i,l) &
- ERI(nO+b,k,nO+c,j)*Kronecker_delta(a,d)*Kronecker_delta(i,l) &
+ ERI(nO+b,l,nO+d,i)*Kronecker_delta(a,c)*Kronecker_delta(j,k) &
- ERI(nO+b,l,nO+c,i)*Kronecker_delta(a,d)*Kronecker_delta(j,k) &
- ERI(nO+b,k,nO+d,i)*Kronecker_delta(a,c)*Kronecker_delta(j,l) &
+ ERI(nO+b,k,nO+c,i)*Kronecker_delta(a,d)*Kronecker_delta(j,l) &
- ERI(nO+a,nO+b,nO+c,nO+d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+ ERI(nO+a,nO+b,nO+c,nO+d)*Kronecker_delta(i,k)*Kronecker_delta(j,l)
H(ishift+iajb,jshift+kcld) = tmp
end do
end do
end do
end do
end do
end do
end do
end do
call diagonalize_matrix(nSD,H,Omega)
call print_excitation('CISD ',ispin,nS,Omega)
if(dump_trans) then
print*,'Singlet CISD transition vectors'
call matout(nSD,nSD,H)
write(*,*)
endif
print*,'DD block done...'
write(*,*)
write(*,*) 'Diagonalizing CISD matrix...'
write(*,*)
call diagonalize_matrix(nH,H,ECISD)
print*,'CISD energies (au)'
call matout(maxH,1,ECISD)
write(*,*)
print*,'CISD excitation energies (eV)'
call matout(maxH-1,1,(ECISD(2:maxH)-ECISD(1))*HaToeV)
write(*,*)
if(dump_trans) then
print*,'Singlet CISD transition vectors'
call matout(nH,nH,H)
write(*,*)
endif
! if(triplet_manifold) then
! ispin = 2
!
! call diagonalize_matrix(nSD,H,Omega)
! call print_excitation('CISD ',ispin,nSD,Omega)
! if(dump_trans) then
! print*,'Triplet CIS transition vectors'
! call matout(nSD,nSD,H)
! write(*,*)
! endif
! endif
end subroutine CISD

14
src/GF/G0F2.f90
View File

@ -1,4 +1,5 @@
subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,regularize, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform a one-shot second-order Green function calculation
@ -16,6 +17,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO
logical,intent(in) :: triplet
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nO
integer,intent(in) :: nC
@ -57,7 +59,15 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO
! Frequency-dependent second-order contribution
call self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,ERI,SigC,Z)
if(regularize) then
call regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,ERI,SigC,Z)
else
call self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,ERI,SigC,Z)
end if
if(linearize) then

13
src/GF/UG0F2.f90
View File

@ -1,4 +1,4 @@
subroutine UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
subroutine UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
S,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,eHF)
! Perform unrestricted G0W0 calculation
@ -18,6 +18,7 @@ subroutine UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linearize,eta,
logical,intent(in) :: spin_flip
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
@ -79,7 +80,15 @@ subroutine UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linearize,eta,
! Compute self-energy !
!---------------------!
call unrestricted_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eHF,SigC,Z)
if(regularize) then
call unrestricted_regularized_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eHF,SigC,Z)
else
call unrestricted_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eHF,SigC,Z)
end if
!-----------------------------------!
! Solve the quasi-particle equation !

14
src/GF/evGF2.f90
View File

@ -1,5 +1,5 @@
subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet, &
linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform eigenvalue self-consistent second-order Green function calculation
@ -20,6 +20,8 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet,
logical,intent(in) :: triplet
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nO
integer,intent(in) :: nC
@ -77,7 +79,15 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet,
! Frequency-dependent second-order contribution
call self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
if(regularize) then
call regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
else
call self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
end if
if(linearize) then

13
src/GF/evUGF2.f90
View File

@ -1,5 +1,5 @@
subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, &
eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, &
eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_aa,dipole_int_bb,cHF,eHF)
! Perform self-consistent eigenvalue-only GW calculation
@ -22,6 +22,7 @@ subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,
logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
@ -112,7 +113,15 @@ subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,
! Compute self-energy and renormalization factor !
!------------------------------------------------!
call unrestricted_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGF2,SigC,Z)
if(regularize) then
call unrestricted_regularized_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGF2,SigC,Z)
else
call unrestricted_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGF2,SigC,Z)
end if
!-----------------------------------!
! Solve the quasi-particle equation !

13
src/GF/qsGF2.f90
View File

@ -1,5 +1,5 @@
subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, &
eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! Perform a quasiparticle self-consistent GF2 calculation
@ -20,6 +20,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
@ -144,7 +145,15 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,
! Compute self-energy and renormalization factor
call self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI_MO,SigC,Z)
if(regularize) then
call regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI_MO,SigC,Z)
else
call self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI_MO,SigC,Z)
end if
! Make correlation self-energy Hermitian and transform it back to AO basis

13
src/GF/qsUGF2.f90
View File

@ -1,4 +1,4 @@
subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, &
subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,regularize, &
nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO, &
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
@ -20,6 +20,7 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,
logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
@ -178,7 +179,15 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,
! Compute self-energy and renormalization factor !
!------------------------------------------------!
call unrestricted_self_energy_GF2(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGF2,SigC,Z)
if(regularize) then
call unrestricted_regularized_self_energy_GF2(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGF2,SigC,Z)
else
call unrestricted_self_energy_GF2(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGF2,SigC,Z)
end if
! Make correlation self-energy Hermitian and transform it back to AO basis

92
src/GF/regularized_self_energy_GF2.f90 Normal file
View File

@ -0,0 +1,92 @@
subroutine regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
! Compute GF2 self-energy and its renormalization factor
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGF2(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: i,j,a,b
integer :: p,q
double precision :: eps
double precision :: num
double precision :: kappa
double precision :: fk,dfk
! Output variables
double precision,intent(out) :: SigC(nBas,nBas)
double precision,intent(out) :: Z(nBas)
! Initialize
SigC(:,:) = 0d0
Z(:) = 0d0
!-----------------------------------------!
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
!----------------------------------------------------!
! Compute GF2 self-energy and renormalization factor !
!----------------------------------------------------!
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(q,a,i,j)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,q) = SigC(p,q) + num*fk
if(p == q) Z(p) = Z(p) - num*dfk
end do
end do
end do
end do
end do
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(q,i,a,b)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,q) = SigC(p,q) + num*fk
if(p == q) Z(p) = Z(p) - num*dfk
end do
end do
end do
end do
end do
Z(:) = 1d0/(1d0 - Z(:))
end subroutine regularized_self_energy_GF2

88
src/GF/regularized_self_energy_GF2_diag.f90 Normal file
View File

@ -0,0 +1,88 @@
subroutine regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z)
! Compute diagonal part of the GF2 self-energy and its renormalization factor
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGF2(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: i,j,a,b
integer :: p
double precision :: eps
double precision :: num
double precision :: kappa
double precision :: fk,dfk
! Output variables
double precision,intent(out) :: SigC(nBas)
double precision,intent(out) :: Z(nBas)
! Initialize
SigC(:) = 0d0
Z(:) = 0d0
!-----------------------------------------!
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
!----------------------------------------------------!
! Compute GF2 self-energy and renormalization factor !
!----------------------------------------------------!
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p) = SigC(p) + num*fk
Z(p) = Z(p) - num*dfk
end do
end do
end do
end do
do p=nC+1,nBas-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p) = SigC(p) + num*fk
Z(p) = Z(p) - num*dfk
end do
end do
end do
end do
Z(:) = 1d0/(1d0 - Z(:))
end subroutine regularized_self_energy_GF2_diag

236
src/GF/unrestricted_regularized_self_energy_GF2.f90 Normal file
View File

@ -0,0 +1,236 @@
subroutine unrestricted_regularized_self_energy_GF2(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eGF2,SigC,Z)
! Perform unrestricted GF2 self-energy and its renormalization factor
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
double precision,intent(in) :: eta
double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: eGF2(nBas,nspin)
! Local variables
integer :: p,q
integer :: i,j,a,b
double precision :: eps,num
double precision :: kappa
double precision :: fk,dfk
! Output variables
double precision,intent(out) :: SigC(nBas,nBas,nspin)
double precision,intent(out) :: Z(nBas,nspin)
!---------------------!
! Compute self-energy |
!---------------------!
SigC(:,:,:) = 0d0
Z(:,:) = 0d0
!-----------------------------------------!
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
!----------------!
! Spin-up sector
!----------------!
do p=nC(1)+1,nBas-nR(1)
do q=nC(1)+1,nBas-nR(1)
! Addition part: aa
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
do b=nO(1)+1,nBas-nR(1)
eps = eGF2(p,1) + eHF(i,1) - eHF(a,1) - eHF(b,1)
num = ERI_aa(i,q,a,b)*ERI_aa(a,b,i,p) &
- ERI_aa(i,q,a,b)*ERI_aa(a,b,p,i)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,q,1) = SigC(p,q,1) + num*fk
if(p == q) Z(p,1) = Z(p,1) - num*dfk
enddo
enddo
enddo
! Addition part: ab
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
do b=nO(1)+1,nBas-nR(1)
eps = eGF2(p,1) + eHF(i,2) - eHF(a,2) - eHF(b,1)
num = ERI_ab(q,i,b,a)*ERI_ab(b,a,p,i)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,q,1) = SigC(p,q,1) + num*fk
if(p == q) Z(p,1) = Z(p,1) - num*dfk
enddo
enddo
enddo
! Removal part: aa
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
do j=nC(1)+1,nO(1)
eps = eGF2(p,1) + eHF(a,1) - eHF(i,1) - eHF(j,1)
num = ERI_aa(a,q,i,j)*ERI_aa(i,j,a,p) &
- ERI_aa(a,q,i,j)*ERI_aa(i,j,p,a)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,q,1) = SigC(p,q,1) + num*fk
if(p == q) Z(p,1) = Z(p,1) - num*dfk
enddo
enddo
enddo
! Removal part: ab
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
do j=nC(1)+1,nO(1)
eps = eGF2(p,1) + eHF(a,2) - eHF(i,2) - eHF(j,1)
num = ERI_ab(q,a,j,i)*ERI_ab(j,i,p,a)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,q,1) = SigC(p,q,1) + num*fk
if(p == q) Z(p,1) = Z(p,1) - num*dfk
enddo
enddo
enddo
enddo
enddo
!------------------!
! Spin-down sector !
!------------------!
do p=nC(2)+1,nBas-nR(2)
do q=nC(2)+1,nBas-nR(2)
! Addition part: bb
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
do b=nO(2)+1,nBas-nR(2)
eps = eGF2(p,2) + eHF(i,2) - eHF(a,2) - eHF(b,2)
num = ERI_bb(i,q,a,b)*ERI_bb(a,b,i,p) &
- ERI_bb(i,q,a,b)*ERI_bb(a,b,p,i)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,q,2) = SigC(p,q,2) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo
enddo
enddo
! Addition part: ab
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
do b=nO(2)+1,nBas-nR(2)
eps = eGF2(p,2) + eHF(i,1) - eHF(a,1) - eHF(b,2)
num = ERI_ab(i,q,a,b)*ERI_ab(a,b,i,p)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,q,2) = SigC(p,q,2) + num*fk
if(p == q) Z(p,2) = Z(p,2) - num*dfk
enddo
enddo
enddo
! Removal part: bb
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
do j=nC(2)+1,nO(2)
eps = eGF2(p,2) + eHF(a,2) - eHF(i,2) - eHF(j,2)
num = ERI_bb(a,q,i,j)*ERI_bb(i,j,a,p) &
- ERI_bb(a,q,i,j)*ERI_bb(i,j,p,a)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,q,2) = SigC(p,q,2) + num*fk
if(p == q) Z(p,2) = Z(p,2) - num*dfk
enddo
enddo
enddo
! Removal part: ab
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
do j=nC(2)+1,nO(2)
eps = eGF2(p,2) + eHF(a,1) - eHF(i,1) - eHF(j,2)
num = ERI_ab(a,q,i,j)*ERI_ab(i,j,a,p)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,q,2) = SigC(p,q,2) + num*fk
if(p == q) Z(p,2) = Z(p,2) - num*dfk
enddo
enddo
enddo
enddo
enddo
Z(:,:) = 1d0/(1d0 - Z(:,:))
end subroutine unrestricted_regularized_self_energy_GF2

231
src/GF/unrestricted_regularized_self_energy_GF2_diag.f90 Normal file
View File

@ -0,0 +1,231 @@
subroutine unrestricted_regularized_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eGF2,SigC,Z)
! Perform unrestricted GF2 self-energy and its renormalization factor
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
double precision,intent(in) :: eta
double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: eGF2(nBas,nspin)
! Local variables
integer :: p
integer :: i,j,a,b
double precision :: eps,num
double precision :: kappa
double precision :: fk,dfk
! Output variables
double precision,intent(out) :: SigC(nBas,nspin)
double precision,intent(out) :: Z(nBas,nspin)
!---------------------!
! Compute self-energy |
!---------------------!
SigC(:,:) = 0d0
!-----------------------------------------!
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
!----------------!
! Spin-up sector
!----------------!
do p=nC(1)+1,nBas-nR(1)
! Addition part: aa
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
do b=nO(1)+1,nBas-nR(1)
eps = eGF2(p,1) + eHF(i,1) - eHF(a,1) - eHF(b,1)
num = ERI_aa(i,p,a,b)*ERI_aa(a,b,i,p) &
- ERI_aa(i,p,a,b)*ERI_aa(a,b,p,i)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,1) = SigC(p,1) + num*fk
Z(p,1) = Z(p,1) - num*dfk
enddo
enddo
enddo
! Addition part: ab
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
do b=nO(1)+1,nBas-nR(1)
eps = eGF2(p,1) + eHF(i,2) - eHF(a,2) - eHF(b,1)
num = ERI_ab(p,i,b,a)*ERI_ab(b,a,p,i)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,1) = SigC(p,1) + num*fk
Z(p,1) = Z(p,1) - num*dfk
enddo
enddo
enddo
! Removal part: aa
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
do j=nC(1)+1,nO(1)
eps = eGF2(p,1) + eHF(a,1) - eHF(i,1) - eHF(j,1)
num = ERI_aa(a,p,i,j)*ERI_aa(i,j,a,p) &
- ERI_aa(a,p,i,j)*ERI_aa(i,j,p,a)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,1) = SigC(p,1) + num*fk
Z(p,1) = Z(p,1) - num*dfk
enddo
enddo
enddo
! Removal part: ab
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
do j=nC(1)+1,nO(1)
eps = eGF2(p,1) + eHF(a,2) - eHF(i,2) - eHF(j,1)
num = ERI_ab(p,a,j,i)*ERI_ab(j,i,p,a)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,1) = SigC(p,1) + num*fk
Z(p,1) = Z(p,1) - num*dfk
enddo
enddo
enddo
enddo
!------------------!
! Spin-down sector !
!------------------!
do p=nC(2)+1,nBas-nR(2)
! Addition part: bb
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
do b=nO(2)+1,nBas-nR(2)
eps = eGF2(p,2) + eHF(i,2) - eHF(a,2) - eHF(b,2)
num = ERI_bb(i,p,a,b)*ERI_bb(a,b,i,p) &
- ERI_bb(i,p,a,b)*ERI_bb(a,b,p,i)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,2) = SigC(p,2) + num*fk
Z(p,2) = Z(p,2) - num*dfk
enddo
enddo
enddo
! Addition part: ab
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
do b=nO(2)+1,nBas-nR(2)
eps = eGF2(p,2) + eHF(i,1) - eHF(a,1) - eHF(b,2)
num = ERI_ab(i,p,a,b)*ERI_ab(a,b,i,p)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,2) = SigC(p,2) + num*fk
Z(p,2) = Z(p,2) - num*dfk
enddo
enddo
enddo
! Removal part: bb
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
do j=nC(2)+1,nO(2)
eps = eGF2(p,2) + eHF(a,2) - eHF(i,2) - eHF(j,2)
num = ERI_bb(a,p,i,j)*ERI_bb(i,j,a,p) &
- ERI_bb(a,p,i,j)*ERI_bb(i,j,p,a)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,2) = SigC(p,2) + num*fk
Z(p,2) = Z(p,2) - num*dfk
enddo
enddo
enddo
! Removal part: ab
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
do j=nC(2)+1,nO(2)
eps = eGF2(p,2) + eHF(a,1) - eHF(i,1) - eHF(j,2)
num = ERI_ab(a,p,i,j)*ERI_ab(i,j,a,p)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
SigC(p,2) = SigC(p,2) + num*fk
Z(p,2) = Z(p,2) - num*dfk
enddo
enddo
enddo
enddo
Z(:,:) = 1d0/(1d0 - Z(:,:))
end subroutine unrestricted_regularized_self_energy_GF2_diag

19
src/MBPT/Bethe_Salpeter_Tmatrix.f90 → src/GT/Bethe_Salpeter_Tmatrix.f90
View File

@ -88,8 +88,13 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
call static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
if(.not.TDA) call static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
! print*,'aa block of TA'
! call matout(nS,nS,TA)
! print*,'aa block of TB'
! call matout(nS,nS,TB)
!----------------------------------------------
! Compute T-matrix for alpha-alpha block
@ -103,8 +108,13 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
call static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
if(.not.TDA) call static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
! print*,'aa+ab block of TA'
! call matout(nS,nS,TA)
! print*,'aa+ab block of TB'
! call matout(nS,nS,TB)
!-------------------
! Singlet manifold
@ -119,6 +129,7 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TA,TB, &
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, &
OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))

0
src/MBPT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90 → src/GT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90
View File

7
src/MBPT/G0T0.f90 → src/GT/G0T0.f90
View File

@ -1,5 +1,5 @@
subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet, &
linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0T0)
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0T0)
! Perform one-shot calculation with a T-matrix self-energy (G0T0)
@ -21,6 +21,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
logical,intent(in) :: triplet
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nC
@ -189,6 +190,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
iblock = 4
call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eG0T0,ERI_MO, &
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
EcRPA(1) = EcRPA(1) - EcRPA(2)
EcRPA(2) = 3d0*EcRPA(2)
@ -198,9 +200,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
if(BSE) then
! eG0T0(1) = -0.5507952119d0
! eG0T0(2) = +1.540259769d0
call Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,rho1t,rho2t, &
ERI_MO,dipole_int,eHF,eG0T0,EcBSE)

42
src/MBPT/dynamic_Tmatrix_A.f90 → src/GT/dynamic_Tmatrix_A.f90
View File

@ -58,11 +58,27 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
chi = 0d0
do cd=1,nVV
chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*Omega1(cd)/(Omega1(cd)**2 + eta**2)
eps = - Omega1(cd)
chi = chi + rho1(i,b,cd)*rho1(j,a,cd)*eps/(eps**2 + eta**2)
end do
do kl=1,nOO
chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2)
do kl=1,nOO
eps = + Omega2(kl)
chi = chi + rho2(i,b,kl)*rho2(j,a,kl)*eps/(eps**2 + eta**2)
end do
A_dyn(ia,jb) = A_dyn(ia,jb) - lambda*chi
chi = 0d0
do cd=1,nVV
eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j))
chi = chi + rho1(i,b,cd)*rho1(j,a,cd)*eps/(eps**2 + eta**2)
end do
do kl=1,nOO
eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b))
chi = chi + rho2(i,b,kl)*rho2(j,a,kl)*eps/(eps**2 + eta**2)
end do
A_dyn(ia,jb) = A_dyn(ia,jb) + 1d0*lambda*chi
@ -71,29 +87,15 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
do cd=1,nVV
eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j))
chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2)
end do
do kl=1,nOO
eps = + OmBSE - Omega2(kl) - (eGT(a) + eGT(b))
chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
end do
A_dyn(ia,jb) = A_dyn(ia,jb) - 1d0*lambda*chi
chi = 0d0
do cd=1,nVV
eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j))
chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
chi = chi + rho1(i,b,cd)*rho1(j,a,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
do kl=1,nOO
eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b))
chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
chi = chi + rho2(i,b,kl)*rho2(j,a,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + 1d0*lambda*chi
ZA_dyn(ia,jb) = ZA_dyn(ia,jb) - 1d0*lambda*chi
end do
end do

3
src/MBPT/evGT.f90 → src/GT/evGT.f90
View File

@ -1,5 +1,5 @@
subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas, &
BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nBas, &
nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0T0)
! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT)
@ -24,6 +24,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nC

0
src/MBPT/excitation_density_Tmatrix.f90 → src/GT/excitation_density_Tmatrix.f90
View File

0
src/MBPT/print_G0T0.f90 → src/GT/print_G0T0.f90
View File

0
src/MBPT/print_evGT.f90 → src/GT/print_evGT.f90
View File

0
src/MBPT/print_qsGT.f90 → src/GT/print_qsGT.f90
View File

3
src/MBPT/qsGT.f90 → src/GT/qsGT.f90
View File

@ -1,5 +1,5 @@
subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA, &
dBSE,dTDA,evDyn,singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! Perform a quasiparticle self-consistent GT calculation
@ -24,6 +24,7 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)

74
src/GT/regularized_renormalization_factor_Tmatrix.f90 Normal file
View File

@ -0,0 +1,74 @@
subroutine regularized_renormalization_factor_Tmatrix(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,a,p,cd,kl
double precision :: eps
double precision :: kappa
double precision :: fk,dfk
! Output variables
double precision,intent(out) :: Z(nBas)
!-----------------------------------------!
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
! 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)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
Z(p) = Z(p) - rho1(p,i,cd)**2*dfk
enddo
enddo
enddo
! Virtual part of the T-matrix self-energy
do p=nC+1,nBas-nR
do a=1,nV-nR
do kl=1,nOO
eps = e(p) + e(nO+a) - Omega2(kl)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
Z(p) = Z(p) - rho2(p,nO+a,kl)**2*dfk
enddo
enddo
enddo
end subroutine regularized_renormalization_factor_Tmatrix

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

@ -0,0 +1,80 @@
subroutine regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,SigT)
! Compute 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,a,p,q,cd,kl
double precision :: eps
double precision :: kappa
double precision :: fk
! Output variables
double precision,intent(inout) :: SigT(nBas,nBas)
!-----------------------------------------!
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
!----------------------------------------------
! Occupied part of the T-matrix self-energy
!----------------------------------------------
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
do cd=1,nVV
eps = e(p) + e(i) - Omega1(cd)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
SigT(p,q) = SigT(p,q) + rho1(p,i,cd)*rho1(q,i,cd)*eps/(eps**2 + eta**2)
enddo
enddo
enddo
enddo
!----------------------------------------------
! Virtual part of the T-matrix self-energy
!----------------------------------------------
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do a=nO+1,nBas-nR
do kl=1,nOO
eps = e(p) + e(a) - Omega2(kl)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
SigT(p,q) = SigT(p,q) + rho2(p,a,kl)*rho2(q,a,kl)*eps/(eps**2 + eta**2)
enddo
enddo
enddo
enddo
end subroutine regularized_self_energy_Tmatrix

76
src/GT/regularized_self_energy_Tmatrix_diag.f90 Normal file
View File

@ -0,0 +1,76 @@
subroutine regularized_self_energy_Tmatrix_diag(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,a,p,cd,kl
double precision :: eps
double precision :: kappa
double precision :: fk
! Output variables
double precision,intent(inout) :: SigT(nBas)
!-----------------------------------------!
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
!----------------------------------------------
! 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)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
SigT(p) = SigT(p) + rho1(p,i,cd)**2*fk
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)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
SigT(p) = SigT(p) + rho2(p,a,kl)**2*fk
enddo
enddo
enddo
end subroutine regularized_self_energy_Tmatrix_diag

2
src/MBPT/renormalization_factor_Tmatrix.f90 → src/GT/renormalization_factor_Tmatrix.f90
View File

@ -19,7 +19,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,
! Local variables
integer :: i,j,k,l,a,b,c,d,p,cd,kl
integer :: i,a,p,cd,kl
double precision :: eps
! Output variables

6
src/MBPT/self_energy_Tmatrix.f90 → src/GT/self_energy_Tmatrix.f90
View File

@ -23,7 +23,7 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2
! Local variables
integer :: i,j,k,l,a,b,c,d,p,q,cd,kl
integer :: i,a,p,q,cd,kl
double precision :: eps
! Output variables
@ -38,7 +38,7 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2
do q=nC+1,nBas-nR
do i=nC+1,nO
do cd=1,nVV
eps = e(p) + e(i) - Omega1(cd)
eps = e(p) + e(i) - Omega1(cd)
SigT(p,q) = SigT(p,q) + rho1(p,i,cd)*rho1(q,i,cd)*eps/(eps**2 + eta**2)
enddo
enddo
@ -53,7 +53,7 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2
do q=nC+1,nBas-nR
do a=nO+1,nBas-nR
do kl=1,nOO
eps = e(p) + e(a) - Omega2(kl)
eps = e(p) + e(a) - Omega2(kl)
SigT(p,q) = SigT(p,q) + rho2(p,a,kl)*rho2(q,a,kl)*eps/(eps**2 + eta**2)
enddo
enddo

6
src/MBPT/self_energy_Tmatrix_diag.f90 → src/GT/self_energy_Tmatrix_diag.f90
View File

@ -23,7 +23,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O
! Local variables
integer :: i,j,k,l,a,b,c,d,p,cd,kl
integer :: i,a,p,cd,kl
double precision :: eps
! Output variables
@ -37,7 +37,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O
do p=nC+1,nBas-nR
do i=nC+1,nO
do cd=1,nVV
eps = e(p) + e(i) - Omega1(cd)
eps = e(p) + e(i) - Omega1(cd)
SigT(p) = SigT(p) + rho1(p,i,cd)**2*eps/(eps**2 + eta**2)
enddo
enddo
@ -50,7 +50,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do kl=1,nOO
eps = e(p) + e(a) - Omega2(kl)
eps = e(p) + e(a) - Omega2(kl)
SigT(p) = SigT(p) + rho2(p,a,kl)**2*eps/(eps**2 + eta**2)
enddo
enddo

20
src/MBPT/static_Tmatrix_TA.f90 → src/GT/static_Tmatrix_A.f90
View File

@ -1,4 +1,4 @@
subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rho1,Omega2,rho2,TA)
subroutine static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rho1,Omega2,rho2,TA)
! Compute the OOVV block of the static T-matrix for the resonant block
@ -7,6 +7,8 @@ subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
! Input variables
integer,intent(in) :: ispin
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -15,7 +17,6 @@ subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
integer,intent(in) :: nS
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: eta
double precision,intent(in) :: lambda
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega1(nVV)
@ -26,6 +27,7 @@ subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
! Local variables
double precision :: chi
double precision :: eps
integer :: i,j,a,b,ia,jb,kl,cd
! Output variables
@ -44,20 +46,22 @@ subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
chi = 0d0
do cd=1,nVV
! chi = chi + lambda*rho1(i,j,cd)*rho1(a,b,cd)*Omega1(cd)/(Omega1(cd)**2 + eta**2)
chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*Omega1(cd)/(Omega1(cd)**2 + eta**2)
eps = - Omega1(cd)
! chi = chi + lambda*rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2)
chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*eps/(eps**2 + eta**2)
enddo
do kl=1,nOO
! chi = chi + lambda*rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2)
chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2)
eps = + Omega2(kl)
! chi = chi - lambda*rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*eps/(eps**2 + eta**2)
enddo
TA(ia,jb) = TA(ia,jb) - 1d0*lambda*chi
TA(ia,jb) = TA(ia,jb) + lambda*chi
enddo
enddo
enddo
enddo
end subroutine static_Tmatrix_TA
end subroutine static_Tmatrix_A

16
src/MBPT/static_Tmatrix_TB.f90 → src/GT/static_Tmatrix_B.f90
View File

@ -1,4 +1,4 @@
subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rho1,Omega2,rho2,TB)
subroutine static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rho1,Omega2,rho2,TB)
! Compute the OVVO block of the static T-matrix for the coupling block
@ -7,6 +7,8 @@ subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
! Input variables
integer,intent(in) :: ispin
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -15,7 +17,6 @@ subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
integer,intent(in) :: nS
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: eta
double precision,intent(in) :: lambda
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega1(nVV)
@ -26,6 +27,7 @@ subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
! Local variables
double precision :: chi
double precision :: eps
integer :: i,j,a,b,ia,jb,kl,cd
! Output variables
@ -44,20 +46,22 @@ subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
chi = 0d0
do cd=1,nVV
eps = - Omega1(cd)
! chi = chi + lambda*rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/Omega1(cd)**2 + eta**2
chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/Omega1(cd)**2 + eta**2
chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2)
enddo
do kl=1,nOO
eps = + Omega2(kl)
! chi = chi + lambda*rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2
chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2
chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
enddo
TB(ia,jb) = TB(ia,jb) - 1d0*lambda*chi
TB(ia,jb) = TB(ia,jb) + lambda*chi
enddo
enddo
enddo
enddo
end subroutine static_Tmatrix_TB
end subroutine static_Tmatrix_B

0
src/MBPT/Bethe_Salpeter.f90 → src/GW/Bethe_Salpeter.f90
View File

0
src/MBPT/Bethe_Salpeter_AB_matrix_dynamic.f90 → src/GW/Bethe_Salpeter_AB_matrix_dynamic.f90
View File

0
src/MBPT/Bethe_Salpeter_A_matrix.f90 → src/GW/Bethe_Salpeter_A_matrix.f90
View File

7
src/MBPT/Bethe_Salpeter_A_matrix_dynamic.f90 → src/GW/Bethe_Salpeter_A_matrix_dynamic.f90
View File

@ -48,15 +48,16 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om
jb = jb + 1
chi = 0d0
do kc=1,maxS
chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2)
eps = OmRPA(kc)
chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps/(eps**2 + eta**2)
enddo
A_dyn(ia,jb) = A_dyn(ia,jb) - 4d0*lambda*chi
chi = 0d0
do kc=1,maxS
eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))

0
src/MBPT/Bethe_Salpeter_B_matrix.f90 → src/GW/Bethe_Salpeter_B_matrix.f90
View File

0
src/MBPT/Bethe_Salpeter_ZAB_matrix_dynamic.f90 → src/GW/Bethe_Salpeter_ZAB_matrix_dynamic.f90
View File

0
src/MBPT/Bethe_Salpeter_ZA_matrix_dynamic.f90 → src/GW/Bethe_Salpeter_ZA_matrix_dynamic.f90
View File

0
src/MBPT/Bethe_Salpeter_dynamic_perturbation.f90 → src/GW/Bethe_Salpeter_dynamic_perturbation.f90
View File

0
src/MBPT/Bethe_Salpeter_dynamic_perturbation_iterative.f90 → src/GW/Bethe_Salpeter_dynamic_perturbation_iterative.f90
View File

18
src/MBPT/G0W0.f90 → src/GW/G0W0.f90
View File

@ -1,5 +1,5 @@
subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet,triplet,linearize,eta, &
dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,regularize, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0W0)
! Perform G0W0 calculation
@ -24,6 +24,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
logical,intent(in) :: triplet
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nC
@ -124,13 +125,18 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
!------------------------!
call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX)
call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
!--------------------------------!
! Compute renormalization factor !
!--------------------------------!
if(regularize) then
call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
else
call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
end if
!-----------------------------------!
! Solve the quasi-particle equation !

0
src/MBPT/G0W0_SOSEX.f90 → src/GW/G0W0_SOSEX.f90
View File

0
src/MBPT/QP_graph.f90 → src/GW/QP_graph.f90
View File

0
src/MBPT/QP_roots.f90 → src/GW/QP_roots.f90
View File

0
src/MBPT/Sangalli_dynamic_perturbation.f90 → src/GW/Sangalli_dynamic_perturbation.f90
View File

0
src/MBPT/SigmaC.f90 → src/GW/SigmaC.f90
View File

3
src/MBPT/UG0W0.f90 → src/GW/UG0W0.f90
View File

@ -1,5 +1,5 @@
subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, &
linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eGW)
! Perform unrestricted G0W0 calculation
@ -24,6 +24,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
logical,intent(in) :: spin_flip
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)

0
src/MBPT/USigmaC.f90 → src/GW/USigmaC.f90
View File

0
src/MBPT/dSigmaC.f90 → src/GW/dSigmaC.f90
View File

0
src/MBPT/dUSigmaC.f90 → src/GW/dUSigmaC.f90
View File

46
src/MBPT/evGW.f90 → src/GW/evGW.f90
View File

@ -1,5 +1,5 @@
subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, &
G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, &
ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0W0)
! Perform self-consistent eigenvalue-only GW calculation
@ -29,6 +29,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -151,7 +153,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
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)
endif
end if
! Compute spectral weights
@ -161,17 +163,33 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
if(G0W) then
! call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
if(regularize) then
call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
else
call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
end if
else
! call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
if(regularize) then
endif
call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
else
call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
end if
end if
! Solve the quasi-particle equation
@ -198,9 +216,9 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGW-eOld,eGW)
else
n_diis = 0
endif
end if
endif
end if
! Save quasiparticles energy for next cycle
@ -210,7 +228,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
nSCF = nSCF + 1
enddo
end do
!------------------------------------------------------------------------
! End main loop
!------------------------------------------------------------------------
@ -231,7 +249,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
stop
endif
end if
! Deallocate memory
@ -288,6 +306,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
end if
endif
end if
end subroutine evGW

3
src/MBPT/evUGW.f90 → src/GW/evUGW.f90
View File

@ -1,5 +1,5 @@
subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc, &
G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc, &
EUHF,S,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0W0)
! Perform self-consistent eigenvalue-only GW calculation
@ -29,6 +29,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)

0
src/MBPT/exchange_matrix_MO_basis.f90 → src/GW/exchange_matrix_MO_basis.f90
View File

0
src/MBPT/excitation_density.f90 → src/GW/excitation_density.f90
View File

0
src/MBPT/excitation_density_SOSEX.f90 → src/GW/excitation_density_SOSEX.f90
View File

0
src/MBPT/plot_GW.f90 → src/GW/plot_GW.f90
View File

0
src/MBPT/print_G0W0.f90 → src/GW/print_G0W0.f90
View File

0
src/MBPT/print_SOSEX.f90 → src/GW/print_SOSEX.f90
View File

0
src/MBPT/print_UG0W0.f90 → src/GW/print_UG0W0.f90
View File

0
src/MBPT/print_evGW.f90 → src/GW/print_evGW.f90
View File

0
src/MBPT/print_evUGW.f90 → src/GW/print_evUGW.f90
View File

0
src/MBPT/print_qsGW.f90 → src/GW/print_qsGW.f90
View File

0
src/MBPT/print_qsUGW.f90 → src/GW/print_qsUGW.f90
View File

31
src/MBPT/qsGW.f90 → src/GW/qsGW.f90
View File

@ -1,5 +1,5 @@
subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! Perform a quasiparticle self-consistent GW calculation
@ -27,6 +27,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
@ -192,15 +193,31 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
if(G0W) then
! call regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
if(regularize) then
call regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
else
call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
end if
else
! call regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
if(regularize) then
call regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
else
call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
endif
endif

0
src/MBPT/qsGW_PT.f90 → src/GW/qsGW_PT.f90
View File

3
src/MBPT/qsUGW.f90 → src/GW/qsUGW.f90
View File

@ -1,5 +1,5 @@
subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, &
G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, &
nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa, &
dipole_int_bb,PHF,cHF,eHF)
@ -28,6 +28,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)

87
src/GW/regularized_renormalization_factor.f90 Normal file
View File

@ -0,0 +1,87 @@
subroutine regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z)
! Compute the regularized version of the GW renormalization factor
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: COHSEX
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
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
integer :: i,a,p,jb
double precision :: eps
double precision :: kappa
double precision :: fk,dfk
! Output variables
double precision,intent(out) :: Z(nBas)
! Initialize
Z(:) = 0d0
!-----------------------------------------!
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
! static COHSEX approximation
if(COHSEX) then
Z(:) = 1d0
return
else
! Occupied part of the correlation self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
eps = e(p) - e(i) + Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
Z(p) = Z(p) - 2d0*rho(p,i,jb)**2*dfk
end do
end do
end do
! Virtual part of the correlation self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(p) - e(a) - Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
Z(p) = Z(p) - 2d0*rho(p,a,jb)**2*dfk
end do
end do
end do
end if
! Compute renormalization factor from derivative of SigC
Z(:) = 1d0/(1d0 - Z(:))
end subroutine regularized_renormalization_factor

7
src/MBPT/regularized_self_energy_correlation.f90 → src/GW/regularized_self_energy_correlation.f90
View File

@ -9,7 +9,12 @@ subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,
logical,intent(in) :: COHSEX
double precision,intent(in) :: eta
integer,intent(in) :: 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) :: e(nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)

7
src/MBPT/regularized_self_energy_correlation_diag.f90 → src/GW/regularized_self_energy_correlation_diag.f90
View File

@ -23,7 +23,6 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,
integer :: i,a,p,q,jb
double precision :: eps
double precision,external :: SigC_dcgw
double precision :: kappa
double precision :: fk
@ -37,9 +36,9 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,
SigC(:) = 0d0
!---------------------------------------------!
! Parameters for regularized MP2 calculations !
!---------------------------------------------!
!-----------------------------------------!
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0

33
src/MBPT/renormalization_factor.f90 → src/GW/renormalization_factor.f90
View File

@ -9,14 +9,19 @@ subroutine renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z)
logical,intent(in) :: COHSEX
double precision,intent(in) :: eta
integer,intent(in) :: 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) :: e(nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
integer :: i,j,a,b,x,jb
integer :: p,i,a,jb
double precision :: eps
! Output variables
@ -38,30 +43,22 @@ subroutine renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z)
! Occupied part of the correlation self-energy
do x=nC+1,nBas-nR
do p=nC+1,nBas-nR
do i=nC+1,nO
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
eps = e(x) - e(i) + Omega(jb)
Z(x) = Z(x) - 2d0*rho(x,i,jb)**2*(eps/(eps**2 + eta**2))**2
end do
do jb=1,nS
eps = e(p) - e(i) + Omega(jb)
Z(p) = Z(p) - 2d0*rho(p,i,jb)**2*(eps/(eps**2 + eta**2))**2
end do
end do
end do
! Virtual part of the correlation self-energy
do x=nC+1,nBas-nR
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
eps = e(x) - e(a) - Omega(jb)
Z(x) = Z(x) - 2d0*rho(x,a,jb)**2*(eps/(eps**2 + eta**2))**2
end do
do jb=1,nS
eps = e(p) - e(a) - Omega(jb)
Z(p) = Z(p) - 2d0*rho(p,a,jb)**2*(eps/(eps**2 + eta**2))**2
end do
end do
end do

0
src/MBPT/renormalization_factor_SOSEX.f90 → src/GW/renormalization_factor_SOSEX.f90
View File

7
src/MBPT/self_energy_correlation.f90 → src/GW/self_energy_correlation.f90
View File

@ -9,7 +9,12 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec
logical,intent(in) :: COHSEX
double precision,intent(in) :: eta
integer,intent(in) :: 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) :: e(nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)

0
src/MBPT/self_energy_correlation_SOSEX_diag.f90 → src/GW/self_energy_correlation_SOSEX_diag.f90
View File

1
src/MBPT/self_energy_correlation_diag.f90 → src/GW/self_energy_correlation_diag.f90
View File

@ -23,7 +23,6 @@ subroutine self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r
integer :: i,a,p,q,jb
double precision :: eps
double precision,external :: SigC_dcgw
! Output variables

0
src/MBPT/self_energy_exchange.f90 → src/GW/self_energy_exchange.f90
View File

0
src/MBPT/self_energy_exchange_diag.f90 → src/GW/self_energy_exchange_diag.f90
View File

0
src/MBPT/static_screening_WA.f90 → src/GW/static_screening_WA.f90
View File

0
src/MBPT/static_screening_WB.f90 → src/GW/static_screening_WB.f90
View File

3
src/MBPT/ufBSE.f90 → src/GW/ufBSE.f90
View File

@ -1,4 +1,4 @@
subroutine ufBSE(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
! Unfold BSE@GW equations
@ -7,7 +7,6 @@ subroutine ufBSE(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO

3
src/MBPT/ufG0W0.f90 → src/GW/ufG0W0.f90
View File

@ -1,4 +1,4 @@
subroutine ufG0W0(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
subroutine ufG0W0(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Unfold G0W0 equations
@ -7,7 +7,6 @@ subroutine ufG0W0(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO

3
src/MBPT/ufGW.f90 → src/GW/ufGW.f90
View File

@ -1,4 +1,4 @@
subroutine ufGW(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
subroutine ufGW(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Unfold GW equations
@ -7,7 +7,6 @@ subroutine ufGW(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO

0
src/MBPT/unrestricted_Bethe_Salpeter.f90 → src/GW/unrestricted_Bethe_Salpeter.f90
View File

0
src/MBPT/unrestricted_Bethe_Salpeter_A_matrix.f90 → src/GW/unrestricted_Bethe_Salpeter_A_matrix.f90
View File

0
src/MBPT/unrestricted_Bethe_Salpeter_A_matrix_dynamic.f90 → src/GW/unrestricted_Bethe_Salpeter_A_matrix_dynamic.f90
View File

0
src/MBPT/unrestricted_Bethe_Salpeter_B_matrix.f90 → src/GW/unrestricted_Bethe_Salpeter_B_matrix.f90
View File

0
src/MBPT/unrestricted_Bethe_Salpeter_ZA_matrix_dynamic.f90 → src/GW/unrestricted_Bethe_Salpeter_ZA_matrix_dynamic.f90
View File

0
src/MBPT/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90 → src/GW/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90
View File

0
src/MBPT/unrestricted_QP_graph.f90 → src/GW/unrestricted_QP_graph.f90
View File

0
src/MBPT/unrestricted_excitation_density.f90 → src/GW/unrestricted_excitation_density.f90
View File

111
src/GW/unrestricted_regularized_renormalization_factor.f90 Normal file
View File

@ -0,0 +1,111 @@
subroutine unrestricted_regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,Z)
! Compute the renormalization factor in the unrestricted formalism
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nSt
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: Omega(nSt)
double precision,intent(in) :: rho(nBas,nBas,nSt,nspin)
! Local variables
integer :: i,a,p,jb
double precision :: eps
double precision :: kappa
double precision :: fk,dfk
! Output variables
double precision,intent(out) :: Z(nBas,nspin)
! Initialize
Z(:,:) = 0d0
!-----------------------------------------!
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
!--------------!
! Spin-up part !
!--------------!
! Occupied part of the renormalization factor
do p=nC(1)+1,nBas-nR(1)
do i=nC(1)+1,nO(1)
do jb=1,nSt
eps = e(p,1) - e(i,1) + Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
Z(p,1) = Z(p,1) + rho(p,i,jb,1)**2*dfk
end do
end do
end do
! Virtual part of the correlation self-energy
do p=nC(1)+1,nBas-nR(1)
do a=nO(1)+1,nBas-nR(1)
do jb=1,nSt
eps = e(p,1) - e(a,1) - Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
Z(p,1) = Z(p,1) + rho(p,a,jb,1)**2*dfk
end do
end do
end do
!----------------!
! Spin-down part !
!----------------!
! Occupied part of the correlation self-energy
do p=nC(2)+1,nBas-nR(2)
do i=nC(2)+1,nO(2)
do jb=1,nSt
eps = e(p,2) - e(i,2) + Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
Z(p,2) = Z(p,2) + rho(p,i,jb,2)**2*dfk
end do
end do
end do
! Virtual part of the correlation self-energy
do p=nC(2)+1,nBas-nR(2)
do a=nO(2)+1,nBas-nR(2)
do jb=1,nSt
eps = e(p,2) - e(a,2) - Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
Z(p,2) = Z(p,2) + rho(p,a,jb,2)**2*dfk
end do
end do
end do
! Final rescaling
Z(:,:) = 1d0/(1d0 + Z(:,:))
end subroutine unrestricted_regularized_renormalization_factor

133
src/GW/unrestricted_regularized_self_energy_correlation.f90 Normal file
View File

@ -0,0 +1,133 @@
subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC,EcGM)
! Compute diagonal of the correlation part of the self-energy
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nSt
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: Omega(nSt)
double precision,intent(in) :: rho(nBas,nBas,nSt,nspin)
! Local variables
integer :: i,a,p,q,jb
double precision :: eps
double precision :: kappa
double precision :: fk
! Output variables
double precision,intent(out) :: SigC(nBas,nBas,nspin)
double precision :: EcGM(nspin)
! Initialize
SigC(:,:,:) = 0d0
EcGM(:) = 0d0
!-----------------------------------------!
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
!--------------!
! Spin-up part !
!--------------!
! Occupied part of the correlation self-energy
do p=nC(1)+1,nBas-nR(1)
do q=nC(1)+1,nBas-nR(1)
do i=nC(1)+1,nO(1)
do jb=1,nSt
eps = e(p,1) - e(i,1) + Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
SigC(p,q,1) = SigC(p,q,1) + rho(p,i,jb,1)*rho(q,i,jb,1)*eps/(eps**2 + eta**2)
end do
end do
end do
end do
! Virtual part of the correlation self-energy
do p=nC(1)+1,nBas-nR(1)
do q=nC(1)+1,nBas-nR(1)
do a=nO(1)+1,nBas-nR(1)
do jb=1,nSt
eps = e(p,1) - e(a,1) - Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
SigC(p,q,1) = SigC(p,q,1) + rho(p,a,jb,1)*rho(q,a,jb,1)*eps/(eps**2 + eta**2)
end do
end do
end do
end do
! GM correlation energy
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
do jb=1,nSt
eps = e(a,1) - e(i,1) + Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
EcGM(1) = EcGM(1) - rho(a,i,jb,1)**2*eps/(eps**2 + eta**2)
end do
end do
end do
!----------------!
! Spin-down part !
!----------------!
! Occupied part of the correlation self-energy
do p=nC(2)+1,nBas-nR(2)
do q=nC(2)+1,nBas-nR(2)
do i=nC(2)+1,nO(2)
do jb=1,nSt
eps = e(p,2) - e(i,2) + Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
SigC(p,q,2) = SigC(p,q,2) + rho(p,i,jb,2)*rho(q,i,jb,2)*eps/(eps**2 + eta**2)
end do
end do
end do
end do
! Virtual part of the correlation self-energy
do p=nC(2)+1,nBas-nR(2)
do q=nC(2)+1,nBas-nR(2)
do a=nO(2)+1,nBas-nR(2)
do jb=1,nSt
eps = e(p,2) - e(a,2) - Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
SigC(p,q,2) = SigC(p,q,2) + rho(p,a,jb,2)*rho(q,a,jb,2)*eps/(eps**2 + eta**2)
end do
end do
end do
end do
! GM correlation energy
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
do jb=1,nSt
eps = e(a,2) - e(i,2) + Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
EcGM(2) = EcGM(2) - rho(a,i,jb,2)**2*eps/(eps**2 + eta**2)
end do
end do
end do
end subroutine unrestricted_self_energy_correlation

126
src/GW/unrestricted_regularized_self_energy_correlation_diag.f90 Normal file
View File

@ -0,0 +1,126 @@
subroutine unrestricted_regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC,EcGM)
! Compute diagonal of the correlation part of the self-energy
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nSt
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: Omega(nSt)
double precision,intent(in) :: rho(nBas,nBas,nSt,nspin)
! Local variables
integer :: i,a,p,q,jb
double precision :: eps
double precision :: kappa
double precision :: fk
! Output variables
double precision,intent(out) :: SigC(nBas,nspin)
double precision :: EcGM(nspin)
! Initialize
SigC(:,:) = 0d0
EcGM(:) = 0d0
!-----------------------------------------!
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
!--------------!
! Spin-up part !
!--------------!
! Occupied part of the correlation self-energy
do p=nC(1)+1,nBas-nR(1)
do i=nC(1)+1,nO(1)
do jb=1,nSt
eps = e(p,1) - e(i,1) + Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
SigC(p,1) = SigC(p,1) + rho(p,i,jb,1)**2*fk
end do
end do
end do
! Virtual part of the correlation self-energy
do p=nC(1)+1,nBas-nR(1)
do a=nO(1)+1,nBas-nR(1)
do jb=1,nSt
eps = e(p,1) - e(a,1) - Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
SigC(p,1) = SigC(p,1) + rho(p,a,jb,1)**2*fk
end do
end do
end do
! GM correlation energy
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
do jb=1,nSt
eps = e(a,1) - e(i,1) + Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
EcGM(1) = EcGM(1) - rho(a,i,jb,1)**2*fk
end do
end do
end do
!----------------!
! Spin-down part !
!----------------!
! Occupied part of the correlation self-energy
do p=nC(2)+1,nBas-nR(2)
do i=nC(2)+1,nO(2)
do jb=1,nSt
eps = e(p,2) - e(i,2) + Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
SigC(p,2) = SigC(p,2) + rho(p,i,jb,2)**2*fk
end do
end do
end do
! Virtual part of the correlation self-energy
do p=nC(2)+1,nBas-nR(2)
do a=nO(2)+1,nBas-nR(2)
do jb=1,nSt
eps = e(p,2) - e(a,2) - Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
SigC(p,2) = SigC(p,2) + rho(p,a,jb,2)**2*fk
end do
end do
end do
! GM correlation energy
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
do jb=1,nSt
eps = e(a,2) - e(i,2) + Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
EcGM(2) = EcGM(2) - rho(a,i,jb,2)**2*fk
end do
end do
end do
end subroutine unrestricted_regularized_self_energy_correlation_diag

2
src/MBPT/unrestricted_renormalization_factor.f90 → src/GW/unrestricted_renormalization_factor.f90
View File

@ -20,7 +20,7 @@ subroutine unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,
! Local variables
integer :: i,j,a,b,p,q,jb
integer :: i,a,p,jb
double precision :: eps
! Output variables

2
src/MBPT/unrestricted_self_energy_correlation.f90 → src/GW/unrestricted_self_energy_correlation.f90
View File

@ -20,7 +20,7 @@ subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega
! Local variables
integer :: i,j,a,b,p,q,jb
integer :: i,a,p,q,jb
double precision :: eps
! Output variables

2
src/MBPT/unrestricted_self_energy_correlation_diag.f90 → src/GW/unrestricted_self_energy_correlation_diag.f90
View File

@ -20,7 +20,7 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e,
! Local variables
integer :: i,j,a,b,p,q,jb
integer :: i,a,p,q,jb
double precision :: eps
! Output variables

6
src/HF/RHF.f90
View File

@ -1,4 +1,4 @@
subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,e,c,P,Vx)
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)
! Perform restricted Hartree-Fock calculation
@ -45,7 +45,6 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
double precision,allocatable :: J(:,:)
double precision,allocatable :: K(:,:)
double precision,allocatable :: cp(:,:)
double precision,allocatable :: F(:,:)
double precision,allocatable :: Fp(:,:)
double precision,allocatable :: ON(:)
@ -56,6 +55,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
double precision,intent(out) :: c(nBas,nBas)
double precision,intent(out) :: P(nBas,nBas)
double precision,intent(out) :: Vx(nBas)
double precision,intent(out) :: F(nBas,nBas)
! Hello world
@ -72,7 +72,7 @@ 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),F(nBas,nBas),ON(nBas), &
cp(nBas,nBas),Fp(nBas,nBas),ON(nBas), &
error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Guess coefficients and eigenvalues

Some files were not shown because too many files have changed in this diff Show More