LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-05-07 07:44:43 +02:00
Code Releases Activity

fix conflict

Browse Source
This commit is contained in:
Clotilde Marut 2020-09-29 11:07:36 +02:00
commit 9074a23700
68 changed files with 3167 additions and 865 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
examples/basis.Ne.cc-pvdz
View File

@ -1,30 +1,29 @@
1 6
S 8
1 17880.0000000 0.0007380
2 2683.0000000 0.0056770
3 611.5000000 0.0288830
4 173.5000000 0.1085400
5 56.6400000 0.2909070
6 20.4200000 0.4483240
7 7.8100000 0.2580260
8 1.6530000 0.0150630
1 17880.0000000 0.0007380
2 2683.0000000 0.0056770
3 611.5000000 0.0288830
4 173.5000000 0.1085400
5 56.6400000 0.2909070
6 20.4200000 0.4483240
7 7.8100000 0.2580260
8 1.6530000 0.0150630
S 8
1 17880.0000000 -0.0001720
2 2683.0000000 -0.0013570
3 611.5000000 -0.0067370
4 173.5000000 -0.0276630
5 56.6400000 -0.0762080
6 20.4200000 -0.1752270
7 7.8100000 -0.1070380
8 1.6530000 0.5670500
1 17880.0000000 -0.0001720
2 2683.0000000 -0.0013570
3 611.5000000 -0.0067370
4 173.5000000 -0.0276630
5 56.6400000 -0.0762080
6 20.4200000 -0.1752270
7 7.8100000 -0.1070380
8 1.6530000 0.5670500
S 1
1 0.4869000 1.0000000
1 0.4869000 1.0000000
P 3
1 28.3900000 0.0460870
2 6.2700000 0.2401810
3 1.6950000 0.5087440
1 28.3900000 0.0460870
2 6.2700000 0.2401810
3 1.6950000 0.5087440
P 1
1 0.4317000 1.0000000
1 0.4317000 1.0000000
D 1
1 2.2020000 1.0000000

4
examples/molecule.H2
View File

@ -1,5 +1,5 @@
# nAt nEla nElb nCore nRyd
2 1 1 0 0
# Znuc x y z
H 0. 0. 0.
H 0. 0. 1.4
H 0. 0. -0.7
H 0. 0. 0.7

9
input/basis
View File

@ -7,12 +7,3 @@ S 1
1 0.1220000 1.0000000
P 1
1 0.7270000 1.0000000
2 3
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
S 1
1 0.1220000 1.0000000
P 1
1 0.7270000 1.0000000

9
input/methods
View File

@ -1,20 +1,21 @@
# RHF UHF MOM
F T F
# MP2 MP3 MP2-F12
# MP2* MP3 MP2-F12
F F F
# CCD CCSD CCSD(T)
F F F
# drCCD rCCD lCCD pCCD
F F F F
# CIS CID CISD
# CIS* CID CISD
F F F
# RPA RPAx ppRPA
# RPA* RPAx* ppRPA
F F F
# G0F2 evGF2 G0F3 evGF3
F F F F
# G0W0 evGW qsGW
# G0W0* evGW* qsGW
T F F
# G0T0 evGT qsGT
F F F
# MCMP2
F
# * unrestricted version available

6
input/molecule
View File

@ -1,5 +1,5 @@
# nAt nEla nElb nCore nRyd
2 1 1 0 0
2 4 3 0 0
# Znuc x y z
H 0. 0. 0.
H 0. 0. 1.4
C 0. 0. -0.16245872
H 0. 0. 1.93436816

4
input/molecule.xyz
View File

@ -1,4 +1,4 @@
2
H 0.0000000000 0.0000000000 0.0000000000
H 0.0000000000 0.0000000000 0.7408481486
C 0.0000000000 0.0000000000 -0.0859694585
H 0.0000000000 0.0000000000 1.0236236215

10
input/options
View File

@ -1,11 +1,11 @@
# RHF: maxSCF thresh DIIS n_diis guess_type ortho_type
64 0.0000001 T 5 1 1
# HF: maxSCF thresh DIIS n_diis guess_type ortho_type
64 0.00001 T 5 1 1
# MP:
# CC: maxSCF thresh DIIS n_diis
64 0.0000001 T 5
# spin: singlet triplet TDA
T T F
# spin: singlet triplet spin_conserved spin_flip TDA
T T T F F
# GF: maxSCF thresh DIIS n_diis lin eta renorm
256 0.00001 T 5 T 0.0 3
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0
@ -13,6 +13,6 @@
# ACFDT: AC Kx XBS
F F T
# BSE: BSE dBSE dTDA evDyn
T T F F
T F T F
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

67
src/QuAcK/ACFDT.f90
View File

@ -1,5 +1,4 @@
subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eW,e,Omega,XpY,XmY,rho,EcAC)
subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,e,EcAC)
! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem
@ -15,8 +14,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
logical,intent(in) :: TDA_W
logical,intent(in) :: TDA
logical,intent(in) :: BSE
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
@ -24,11 +23,6 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision :: Omega(nS,nspin)
double precision :: XpY(nS,nS,nspin)
double precision :: XmY(nS,nS,nspin)
double precision :: rho(nBas,nBas,nS,nspin)
! Local variables
integer :: ispin
@ -37,6 +31,16 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
double precision :: lambda
double precision,allocatable :: Ec(:,:)
double precision :: EcRPA
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
! Output variables
double precision,intent(out) :: EcAC(nspin)
@ -44,6 +48,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
! Memory allocation
allocate(Ec(nAC,nspin))
allocate(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
@ -58,12 +64,20 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
EcAC(:) = 0d0
Ec(:,:) = 0d0
! Compute (singlet) RPA screening
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 excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
! Singlet manifold
if(singlet_manifold) then
if(singlet) then
ispin = 1
isp_W = 1
write(*,*) '--------------'
write(*,*) 'Singlet states'
@ -80,17 +94,17 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
if(doXBS) then
call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, &
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
end if
call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
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 ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, &
ERI(:,:,:,:),XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
@ -98,6 +112,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin))
if(exchange_kernel) EcAC(ispin) = 0.5d0*EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
@ -107,7 +123,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
! Triplet manifold
if(triplet_manifold) then
if(triplet) then
ispin = 2
isp_W = 1
@ -127,17 +143,16 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
if(doXBS) then
call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, &
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
end if
call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
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 ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, &
ERI(:,:,:,:),XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
@ -145,6 +160,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin))
if(exchange_kernel) EcAC(ispin) = 1.5d0*EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'

16
src/QuAcK/AOtoMO_integral_transform.f90
View File

@ -1,15 +1,15 @@
subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis)
subroutine AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,c,ERI_AO_basis,ERI_MO_basis)
! AO to MO transformation of two-electron integrals via the semi-direct O(N^5) algorithm
! bra and ket are the spin of (bra|ket)
! bra and ket are the spin of (bra1 bra2|ket1 ket2)
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: bra
integer,intent(in) :: ket
integer,intent(in) :: bra1,bra2
integer,intent(in) :: ket1,ket2
integer,intent(in) :: nBas
double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas),c(nBas,nBas,nspin)
@ -35,7 +35,7 @@ subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis)
do la=1,nBas
do nu=1,nBas
do mu=1,nBas
scr(mu,nu,la,l) = scr(mu,nu,la,l) + ERI_AO_basis(mu,nu,la,si)*c(si,l,ket)
scr(mu,nu,la,l) = scr(mu,nu,la,l) + ERI_AO_basis(mu,nu,la,si)*c(si,l,ket2)
enddo
enddo
enddo
@ -49,7 +49,7 @@ subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis)
do nu=1,nBas
do i=1,nBas
do mu=1,nBas
ERI_MO_basis(i,nu,la,l) = ERI_MO_basis(i,nu,la,l) + c(mu,i,bra)*scr(mu,nu,la,l)
ERI_MO_basis(i,nu,la,l) = ERI_MO_basis(i,nu,la,l) + c(mu,i,bra1)*scr(mu,nu,la,l)
enddo
enddo
enddo
@ -63,7 +63,7 @@ subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis)
do la=1,nBas
do nu=1,nBas
do i=1,nBas
scr(i,nu,k,l) = scr(i,nu,k,l) + ERI_MO_basis(i,nu,la,l)*c(la,k,bra)
scr(i,nu,k,l) = scr(i,nu,k,l) + ERI_MO_basis(i,nu,la,l)*c(la,k,bra2)
enddo
enddo
enddo
@ -77,7 +77,7 @@ subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis)
do j=1,nBas
do i=1,nBas
do nu=1,nBas
ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j,ket)*scr(i,nu,k,l)
ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j,ket1)*scr(i,nu,k,l)
enddo
! print*,i,k,j,l,ERI_MO_basis(i,j,k,l)
enddo

34
src/QuAcK/BSE2.f90
View File

@ -1,5 +1,4 @@
subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF,EcBSE)
subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,EcBSE)
! Compute the Bethe-Salpeter excitation energies
@ -12,8 +11,8 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: evDyn
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
integer,intent(in) :: nBas
@ -25,6 +24,7 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
@ -46,15 +46,15 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
! Singlet manifold
!-------------------
if(singlet_manifold) then
if(singlet) then
ispin = 1
EcBSE(ispin) = 0d0
! Compute BSE2 excitation energies
call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF(:),ERI(:,:,:,:), &
rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
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 print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin))
! Compute dynamic correction for BSE via perturbation theory
@ -63,12 +63,12 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
if(evDyn) then
call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, &
ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
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))
else
call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, &
ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, &
OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
end if
@ -80,7 +80,7 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
! Triplet manifold
!-------------------
if(triplet_manifold) then
if(triplet) then
ispin = 2
EcBSE(ispin) = 0d0
@ -88,7 +88,7 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
! Compute BSE2 excitation energies
call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF(:),ERI(:,:,:,:), &
rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
OmBSE(:,ispin),rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin))
! Compute dynamic correction for BSE via perturbation theory
@ -97,12 +97,12 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
if(evDyn) then
call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, &
ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
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))
else
call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, &
ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, &
OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
end if

5
src/QuAcK/BSE2_dynamic_perturbation.f90
View File

@ -1,4 +1,4 @@
subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF,OmBSE,XpY,XmY)
subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,OmBSE,XpY,XmY)
! Compute dynamical effects via perturbation theory for BSE
@ -18,6 +18,7 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,
integer,intent(in) :: nS
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
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) :: OmBSE(nS)
@ -56,7 +57,7 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,
! Print main components of transition vectors
call print_transition_vectors(nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY)
call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)
gapGF = eGF(nO+1) - eGF(nO)

6
src/QuAcK/BSE2_dynamic_perturbation_iterative.f90
View File

@ -1,4 +1,5 @@
subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF,OmBSE,XpY,XmY)
subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int, &
eHF,eGF,OmBSE,XpY,XmY)
! Compute self-consistently the dynamical effects via perturbation theory for BSE2
@ -18,6 +19,7 @@ subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,n
integer,intent(in) :: nS
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
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) :: OmBSE(nS)
@ -58,7 +60,7 @@ subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,n
! Print main components of transition vectors
call print_transition_vectors(nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY)
call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY)
if(dTDA) then
write(*,*)

87
src/QuAcK/Bethe_Salpeter.f90
View File

@ -1,5 +1,4 @@
subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,EcRPA,EcBSE)
subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE)
! Compute the Bethe-Salpeter excitation energies
@ -13,60 +12,67 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: evDyn
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eW(nBas)
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision :: OmRPA(nS,nspin)
double precision :: XpY_RPA(nS,nS,nspin)
double precision :: XmY_RPA(nS,nS,nspin)
double precision :: rho_RPA(nBas,nBas,nS,nspin)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
integer :: ispin
integer :: isp_W
double precision :: EcRPA
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: OmBSE(:,:)
double precision,allocatable :: XpY_BSE(:,:,:)
double precision,allocatable :: XmY_BSE(:,:,:)
! Output variables
double precision,intent(out) :: EcRPA(nspin)
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))
!---------------------------------
! Compute (singlet) RPA screening
!---------------------------------
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 excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
!-------------------
! Singlet manifold
!-------------------
if(singlet_manifold) then
if(singlet) then
ispin = 1
isp_W = 1
EcBSE(ispin) = 0d0
! Compute (singlet) RPA screening
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA(:,:,ispin),rho_RPA(:,:,:,ispin))
! Compute BSE excitation energies
OmBSE(:,ispin) = OmRPA(:,ispin)
call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho_RPA(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin))
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 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))
!-------------------------------------------------
! Compute the dynamical screening at the BSE level
@ -78,12 +84,12 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man
if(evDyn) then
call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), &
XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin))
call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, &
OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
else
call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), &
XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin))
call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, &
OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
end if
end if
@ -94,25 +100,18 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man
! Triplet manifold
!-------------------
if(triplet_manifold) then
if(triplet) then
ispin = 2
isp_W = 1
EcBSE(ispin) = 0d0
! Compute (singlet) RPA screening
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA(:,:,ispin),rho_RPA(:,:,:,ispin))
! Compute BSE excitation energies
OmBSE(:,ispin) = OmRPA(:,ispin)
call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho_RPA(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin))
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 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))
!-------------------------------------------------
! Compute the dynamical screening at the BSE level
@ -124,12 +123,12 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man
if(evDyn) then
call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), &
XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin))
call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, &
OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
else
call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), &
XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin))
call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, &
OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
end if
end if

25
src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90
View File

@ -1,5 +1,4 @@
subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho, &
Ap,Am,Bp,Bm)
subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rhO_RPA,OmBSE,Ap,Am,Bp,Bm)
! Compute the dynamic part of the Bethe-Salpeter equation matrices
@ -13,8 +12,8 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
double precision,intent(in) :: lambda
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS)
double precision,intent(in) :: OmBSE
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
@ -60,8 +59,8 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
do kc=1,maxS
chi_A = chi_A + rho(i,j,kc)*rho(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2)
chi_B = chi_B + rho(i,b,kc)*rho(a,j,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2)
chi_A = chi_A + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2)
chi_B = chi_B + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2)
enddo
@ -80,28 +79,28 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
do kc=1,maxS
eps_Ap = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2)
chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2)
eps_Ap = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2)
chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2)
eps_Am = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*eps_Am/(eps_Am**2 + eta**2)
chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Am/(eps_Am**2 + eta**2)
eps_Am = - OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*eps_Am/(eps_Am**2 + eta**2)
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))
chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2)
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))
chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2)
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))
chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2)
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))
chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2)
chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2)
enddo

10
src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,A_dyn)
subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,A_dyn)
! Compute the dynamic part of the Bethe-Salpeter equation matrices
@ -12,8 +12,8 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om
double precision,intent(in) :: lambda
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS)
double precision,intent(in) :: OmBSE
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
@ -48,7 +48,7 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om
chi = 0d0
do kc=1,maxS
chi = chi + rho(i,j,kc)*rho(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2)
chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2)
enddo
@ -58,10 +58,10 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om
do kc=1,maxS
eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
chi = chi + rho(i,j,kc)*rho(a,b,kc)*eps/(eps**2 + eta**2)
chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps/(eps**2 + eta**2)
eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
chi = chi + rho(i,j,kc)*rho(a,b,kc)*eps/(eps**2 + eta**2)
chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps/(eps**2 + eta**2)
enddo

10
src/QuAcK/Bethe_Salpeter_B_matrix.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B_lr)
subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,B_lr)
! Compute the extra term for Bethe-Salpeter equation for linear response
@ -11,8 +11,8 @@ subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,
double precision,intent(in) :: eta
double precision,intent(in) :: lambda
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS)
! Local variables
@ -35,8 +35,8 @@ subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,
chi = 0d0
do kc=1,nS
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,b,kc)*rho(a,j,kc)*Omega(kc)/eps
eps = OmRPA(kc)**2 + eta**2
chi = chi + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*OmRPA(kc)/eps
enddo
B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI(i,j,b,a) + 4d0*lambda*chi

20
src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho, &
subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE, &
ZAp,ZAm,ZBp,ZBm)
! Compute the dynamic part of the renormalization for the Bethe-Salpeter equation matrices
@ -13,8 +13,8 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,
double precision,intent(in) :: lambda
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS)
double precision,intent(in) :: OmBSE
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
@ -63,28 +63,28 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,
do kc=1,maxS
eps_Ap = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2
chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2
eps_Ap = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2
chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2
eps_Am = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2
chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2
eps_Am = - OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2
chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2
eps_Bp = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b))
chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2
chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2
eps_Bp = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i))
chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2
chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2
eps_Bm = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(b))
chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2
chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2
eps_Bm = - OmBSE - OmRPA(kc) - (eGW(j) - eGW(i))
chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2
chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2
enddo

8
src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,ZA_dyn)
subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,ZA_dyn)
! Compute the dynamic part of the Bethe-Salpeter equation matrices
@ -12,8 +12,8 @@ subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
double precision,intent(in) :: lambda
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS)
double precision,intent(in) :: OmBSE
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
@ -49,10 +49,10 @@ subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
do kc=1,maxS
eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
chi = chi + rho(i,j,kc)*rho(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
chi = chi + rho(i,j,kc)*rho(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo

16
src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90
View File

@ -1,4 +1,5 @@
subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho)
subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int, &
OmRPA,rho_RPA,OmBSE,XpY,XmY)
! Compute dynamical effects via perturbation theory for BSE
@ -17,11 +18,12 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,
integer,intent(in) :: nS
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS)
double precision,intent(in) :: OmBSE(nS)
double precision,intent(in) :: XpY(nS,nS)
double precision,intent(in) :: XmY(nS,nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
@ -55,7 +57,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,
! Print main components of transition vectors
call print_transition_vectors(nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY)
call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)
if(dTDA) then
write(*,*)
@ -84,11 +86,11 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,
! Resonant part of the BSE correction for dynamical TDA
call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho,Ap_dyn)
call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),Ap_dyn)
! Renormalization factor of the resonant parts for dynamical TDA
call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho,ZAp_dyn)
call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),ZAp_dyn)
ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X))
OmDyn(ia) = dot_product(X,matmul( Ap_dyn,X))
@ -97,12 +99,12 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,
! Resonant and anti-resonant part of the BSE correction
call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho, &
call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia), &
Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn)
! Renormalization factor of the resonant and anti-resonant parts
call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho, &
call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia), &
ZAp_dyn,ZAm_dyn,ZBp_dyn,ZBm_dyn)
ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) &

15
src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90
View File

@ -1,4 +1,5 @@
subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho)
subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int, &
OmRPA,rho_RPA,OmBSE,XpY,XmY)
! Compute self-consistently the dynamical effects via perturbation theory for BSE
@ -17,11 +18,12 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
integer,intent(in) :: nS
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS)
double precision,intent(in) :: OmBSE(nS)
double precision,intent(in) :: XpY(nS,nS)
double precision,intent(in) :: XmY(nS,nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
@ -54,7 +56,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
! Print main components of transition vectors
call print_transition_vectors(nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY)
call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)
if(dTDA) then
write(*,*)
@ -97,8 +99,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
! Resonant part of the BSE correction
call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmOld(ia),rho(:,:,:), &
Ap_dyn(:,:))
call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia),Ap_dyn)
OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:)))
@ -106,8 +107,8 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
! Anti-resonant part of the BSE correction
call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmOld(ia),rho(:,:,:), &
Ap_dyn(:,:),Am_dyn(:,:),Bp_dyn(:,:),Bm_dyn(:,:))
call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia), &
Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn)
OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) &
- dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) &

3
src/QuAcK/CIS.f90
View File

@ -1,5 +1,5 @@
subroutine CIS(singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,eHF)
nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
! Perform configuration interaction single calculation`
@ -13,6 +13,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables

5
src/QuAcK/G0F2.f90
View File

@ -1,5 +1,5 @@
subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform a one-shot second-order Green function calculation
@ -27,6 +27,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
@ -115,7 +116,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
if(BSE) then
call BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF2,EcBSE)
call BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE)
end if

21
src/QuAcK/G0T0.f90
View File

@ -1,5 +1,5 @@
subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eG0T0)
subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, &
linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF,eG0T0)
! Perform one-shot calculation with a T-matrix self-energy (G0T0)
@ -17,8 +17,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,sing
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: evDyn
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: linearize
double precision,intent(in) :: eta
@ -32,6 +32,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,sing
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
@ -212,10 +213,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,sing
if(BSE) then
allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin))
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eHF,eG0T0,Omega,XpY,XmY,rho,EcRPA,EcBSE)
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eG0T0,EcBSE)
if(exchange_kernel) then
@ -249,8 +248,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,sing
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eHF,eG0T0,Omega,XpY,XmY,rho,EcAC)
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eHF,eG0T0,EcAC)
if(exchange_kernel) then
@ -270,6 +269,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,sing
end if
allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin))
end if
end subroutine G0T0

147
src/QuAcK/G0W0.f90
View File

@ -1,6 +1,6 @@
subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linearize,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eGW)
subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet,triplet,linearize,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI,dipole_int,PHF,cHF,eHF,eGW)
! Perform G0W0 calculation
@ -21,8 +21,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: evDyn
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: linearize
double precision,intent(in) :: eta
@ -33,24 +33,23 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: H(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
logical :: print_W = .true.
integer :: ispin
double precision :: EcRPA(nspin)
double precision :: EcRPA
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision :: EcGM
double precision,allocatable :: SigC(:)
double precision,allocatable :: Z(:)
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
double precision,allocatable :: rho(:,:,:,:)
double precision,allocatable :: rhox(:,:,:,:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: eGWlin(:)
@ -68,27 +67,35 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
! Initialization
EcRPA(:) = 0d0
EcRPA = 0d0
! SOSEX correction
if(SOSEX) write(*,*) 'SOSEX correction activated!'
write(*,*)
if(SOSEX) then
write(*,*) 'SOSEX correction activated but BUG!'
stop
end if
! COHSEX approximation
if(COHSEX) write(*,*) 'COHSEX approximation activated!'
write(*,*)
if(COHSEX) then
write(*,*) 'COHSEX approximation activated!'
write(*,*)
end if
! TDA for W
if(TDA_W) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
write(*,*)
if(TDA_W) then
write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
write(*,*)
end if
! TDA
if(TDA) write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
end if
! Spin manifold
@ -96,34 +103,43 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
! Memory allocation
allocate(SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin), &
rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin),eGWlin(nBas))
allocate(SigC(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS),eGWlin(nBas))
! Compute linear response
!-------------------!
! Compute screening !
!-------------------!
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0, &
eHF,ERI,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
if(print_W) call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin))
if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA)
! Compute correlation part of the self-energy
!--------------------------!
! Compute spectral weights !
!--------------------------!
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rhox(:,:,:,ispin))
!------------------------!
! Compute GW self-energy !
!------------------------!
call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC)
call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
! Compute renormalization factor
!--------------------------------!
! Compute renormalization factor !
!--------------------------------!
call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:))
call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
! Solve the quasi-particle equation
!-----------------------------------!
! Solve the quasi-particle equation !
!-----------------------------------!
eGWlin(:) = eHF(:) + Z(:)*SigC(:)
! Linearized or graphical solution?
if(linearize) then
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
@ -131,51 +147,48 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
eGW(:) = eGWlin(:)
! Find all the roots of the QP equation if necessary
! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin)
else
! Find graphical solution of the QP equation
call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,OmRPA,rho_RPA,eGWlin,eGW)
! Find all the roots of the QP equation if necessary
! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin)
end if
! Dump results
call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM)
! Compute the RPA correlation energy
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
!--------------!
! Dump results !
!--------------!
call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
! Deallocate memory
deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eGWlin)
! Plot stuff
! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin))
! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,OmRPA,rho_RPA)
! Perform BSE calculation
if(BSE) then
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
if(exchange_kernel) then
EcRPA(1) = 0.5d0*EcRPA(1)
EcRPA(2) = 1.5d0*EcRPA(1)
EcBSE(1) = 0.5d0*EcBSE(1)
EcBSE(2) = 1.5d0*EcBSE(2)
end if
@ -204,15 +217,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcAC)
if(exchange_kernel) then
EcAC(1) = 0.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(1)
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

219
src/QuAcK/QuAcK.f90
View File

@ -4,7 +4,7 @@ program QuAcK
include 'parameters.h'
logical :: doSph
logical :: unrestricted
logical :: unrestricted = .false.
logical :: doRHF,doUHF,doMOM
logical :: doMP2,doMP3,doMP2F12
logical :: doCCD,doCCSD,doCCSDT
@ -52,14 +52,25 @@ program QuAcK
integer :: TrialType
double precision,allocatable :: cTrial(:),gradient(:),hessian(:,:)
double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),H(:,:),X(:,:)
double precision,allocatable :: S(:,:)
double precision,allocatable :: T(:,:)
double precision,allocatable :: V(:,:)
double precision,allocatable :: Hc(:,:)
double precision,allocatable :: H(:,:)
double precision,allocatable :: X(:,:)
double precision,allocatable :: dipole_int(:,:,:)
double precision,allocatable :: dipole_int_aa(:,:,:)
double precision,allocatable :: dipole_int_bb(:,:,:)
double precision,allocatable :: ERI_AO(:,:,:,:)
double precision,allocatable :: ERI_MO(:,:,:,:)
integer :: bra
integer :: ket
double precision,allocatable :: ERI_MO_aa(:,:,:,:)
double precision,allocatable :: ERI_MO_ab(:,:,:,:)
double precision,allocatable :: ERI_MO_bb(:,:,:,:)
integer :: ixyz
integer :: ispin
integer :: bra1,bra2
integer :: ket1,ket2
double precision,allocatable :: ERI_MO_aaaa(:,:,:,:)
double precision,allocatable :: ERI_MO_aabb(:,:,:,:)
double precision,allocatable :: ERI_MO_bbbb(:,:,:,:)
double precision,allocatable :: ERI_MO_abab(:,:,:,:)
double precision,allocatable :: ERI_ERF_AO(:,:,:,:)
double precision,allocatable :: ERI_ERF_MO(:,:,:,:)
double precision,allocatable :: F12(:,:,:,:),Yuk(:,:,:,:),FC(:,:,:,:,:,:)
@ -101,8 +112,10 @@ program QuAcK
double precision :: thresh_CC
logical :: DIIS_CC
logical :: singlet_manifold
logical :: triplet_manifold
logical :: singlet
logical :: triplet
logical :: spin_conserved
logical :: spin_flip
logical :: TDA
integer :: maxSCF_GF,n_diis_GF,renormGF
@ -156,7 +169,7 @@ program QuAcK
call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
singlet_manifold,triplet_manifold,TDA, &
singlet,triplet,spin_conserved,spin_flip,TDA, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW, &
COHSEX,SOSEX,TDA_W,G0W,GW0, &
@ -316,7 +329,6 @@ program QuAcK
write(*,*) 'AO to MO transformation... Please be patient'
write(*,*)
if(doSph) then
allocate(ERI_MO(nBas,nBas,nBas,nBas))
@ -328,39 +340,81 @@ program QuAcK
if(unrestricted) then
! Read and transform dipole-related integrals
allocate(dipole_int_aa(nBas,nBas,ncart),dipole_int_bb(nBas,nBas,ncart))
call read_dipole_integrals(nBas,dipole_int_aa)
call read_dipole_integrals(nBas,dipole_int_bb)
do ixyz=1,ncart
call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_aa(:,:,ixyz))
call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_bb(:,:,ixyz))
end do
! Memory allocation
allocate(ERI_MO_aa(nBas,nBas,nBas,nBas),ERI_MO_ab(nBas,nBas,nBas,nBas),ERI_MO_bb(nBas,nBas,nBas,nBas))
allocate(ERI_MO_aaaa(nBas,nBas,nBas,nBas),ERI_MO_aabb(nBas,nBas,nBas,nBas),ERI_MO_bbbb(nBas,nBas,nBas,nBas))
! 4-index transform for (aa|aa) block
bra = 1
ket = 1
call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_aa)
bra1 = 1
bra2 = 1
ket1 = 1
ket2 = 1
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_aaaa)
! 4-index transform for (bb|bb) block
bra = 1
ket = 2
call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_ab)
! 4-index transform for (aa|bb) block
bra = 2
ket = 2
call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_bb)
bra1 = 1
bra2 = 1
ket1 = 2
ket2 = 2
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_aabb)
! 4-index transform for (bb|bb) block
bra1 = 2
bra2 = 2
ket1 = 2
ket2 = 2
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_bbbb)
if(spin_flip) then
allocate(ERI_MO_abab(nBas,nBas,nBas,nBas))
! 4-index transform for (ab|ab) block
bra1 = 1
bra2 = 2
ket1 = 1
ket2 = 2
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_abab)
end if
else
! Memory allocation
allocate(ERI_MO(nBas,nBas,nBas,nBas))
! Read and transform dipole-related integrals
allocate(dipole_int(nBas,nBas,ncart))
call read_dipole_integrals(nBas,dipole_int)
do ixyz=1,ncart
call AOtoMO_transform(nBas,cHF,dipole_int(:,:,ixyz))
end do
! 4-index transform
bra = 1
ket = 1
call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO)
bra1 = 1
bra2 = 1
ket1 = 1
ket2 = 1
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO)
end if
@ -382,7 +436,7 @@ program QuAcK
if(unrestricted) then
call UMP2(nBas,nC,nO,nV,nR,ERI_MO_aa,ERI_MO_ab,ERI_MO_bb,ENuc,EUHF,eHF,EcMP2)
call UMP2(nBas,nC,nO,nV,nR,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ENuc,EUHF,eHF,EcMP2)
else
@ -560,7 +614,16 @@ program QuAcK
if(doCIS) then
call cpu_time(start_CIS)
call CIS(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF)
if(unrestricted) then
call UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_MO_aaaa,ERI_MO_aabb, &
ERI_MO_bbbb,ERI_MO_abab,dipole_int,eHF)
else
call CIS(singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF)
end if
call cpu_time(end_CIS)
t_CIS = end_CIS - start_CIS
@ -576,7 +639,7 @@ program QuAcK
if(doCID) then
call cpu_time(start_CID)
! call CID(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI_MO,eHF)
! call CID(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,eHF)
call cpu_time(end_CID)
t_CID = end_CID - start_CID
@ -592,7 +655,7 @@ program QuAcK
if(doCISD) then
call cpu_time(start_CISD)
call CISD(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI_MO,eHF)
! call CISD(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,eHF)
call cpu_time(end_CISD)
t_CISD = end_CISD - start_CISD
@ -608,8 +671,16 @@ program QuAcK
if(doRPA) then
call cpu_time(start_RPA)
call RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,0d0, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
if(unrestricted) then
call UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int,eHF)
else
call dRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF)
end if
call cpu_time(end_RPA)
t_RPA = end_RPA - start_RPA
@ -625,8 +696,16 @@ program QuAcK
if(doRPAx) then
call cpu_time(start_RPAx)
call RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,0d0, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
if(unrestricted) then
call URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF)
else
call RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF)
end if
call cpu_time(end_RPAx)
t_RPAx = end_RPAx - start_RPAx
@ -642,8 +721,7 @@ program QuAcK
if(doppRPA) then
call cpu_time(start_ppRPA)
call ppRPA(singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF)
call ppRPA(singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_ppRPA)
t_ppRPA = end_ppRPA - start_ppRPA
@ -659,7 +737,7 @@ program QuAcK
! if(doADC) then
! call cpu_time(start_ADC)
! call ADC(singlet_manifold,triplet_manifold,maxSCF_GF,thresh_GF,n_diis_GF, &
! call ADC(singlet,triplet,maxSCF_GF,thresh_GF,n_diis_GF, &
! nBas,nC,nO,nV,nR,eHF,ERI_MO)
! call cpu_time(end_ADC)
@ -676,8 +754,8 @@ program QuAcK
if(doG0F2) then
call cpu_time(start_GF2)
call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGF, &
eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linGF, &
eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF)
call cpu_time(end_GF2)
t_GF2 = end_GF2 - start_GF2
@ -694,8 +772,8 @@ program QuAcK
call cpu_time(start_GF2)
call evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF_GF,thresh_GF,n_diis_GF, &
singlet_manifold,triplet_manifold,linGF, &
eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
singlet,triplet,linGF,eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, &
ERI_MO,dipole_int,eHF)
call cpu_time(end_GF2)
t_GF2 = end_GF2 - start_GF2
@ -747,14 +825,13 @@ program QuAcK
call cpu_time(start_G0W0)
if(unrestricted) then
call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, &
singlet_manifold,triplet_manifold,linGW,eta_GW,nBas,nC,nO,nV,nR,nS, &
ENuc,EUHF,Hc,ERI_MO_aa,ERI_MO_ab,ERI_MO_bb,PHF,cHF,eHF,eG0W0)
call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, &
linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab, &
dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eG0W0)
else
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,PHF,cHF,eHF,eG0W0)
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, &
linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,dipole_int,PHF,cHF,eHF,eG0W0)
end if
@ -773,9 +850,19 @@ program QuAcK
if(doevGW) then
call cpu_time(start_evGW)
call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, &
BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0)
if(unrestricted) then
call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc, &
EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb, &
PHF,cHF,eHF,eG0W0)
else
call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, &
BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,dipole_int,PHF,cHF,eHF,eG0W0)
end if
call cpu_time(end_evGW)
t_evGW = end_evGW - start_evGW
@ -792,8 +879,8 @@ program QuAcK
call cpu_time(start_qsGW)
call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, &
BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,PHF,cHF,eHF)
BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF)
call cpu_time(end_qsGW)
t_qsGW = end_qsGW - start_qsGW
@ -811,9 +898,9 @@ program QuAcK
if(doG0T0) then
call cpu_time(start_G0T0)
call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0)
call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF,eG0T0)
call cpu_time(end_G0T0)
t_G0T0 = end_G0T0 - start_G0T0
@ -829,9 +916,9 @@ program QuAcK
if(doevGT) then
call cpu_time(start_evGT)
call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, &
BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0)
call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, &
BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF,eG0T0)
call cpu_time(end_evGT)
t_evGT = end_evGT - start_evGT
@ -840,10 +927,6 @@ program QuAcK
end if
!------------------------------------------------------------------------
! Perform evGT calculatiom
!------------------------------------------------------------------------
!------------------------------------------------------------------------
! Information for Monte Carlo calculations
!------------------------------------------------------------------------
@ -947,8 +1030,8 @@ program QuAcK
call cpu_time(start_G0W0)
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_ERF_MO,PHF,cHF,eHF,eG0W0)
dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_ERF_MO,dipole_int,PHF,cHF,eHF,eG0W0)
call cpu_time(end_G0W0)
t_G0W0 = end_G0W0 - start_G0W0
@ -961,8 +1044,8 @@ program QuAcK
call cpu_time(start_G0T0)
call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, &
singlet_manifold,triplet_manifold,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,eHF,eG0T0)
singlet,triplet,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,dipole_int,eHF,eG0T0)
call cpu_time(end_G0T0)
t_G0T0 = end_G0T0 - start_G0T0

40
src/QuAcK/RPAx.f90
View File

@ -1,5 +1,5 @@
subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e)
subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, &
ERI,dipole_int,eHF)
! Perform random phase approximation calculation with exchange (aka TDHF)
@ -9,11 +9,12 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
! Input variables
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: singlet_manifold
logical,intent(in) :: singlet
double precision,intent(in) :: eta
logical,intent(in) :: triplet_manifold
logical,intent(in) :: triplet
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -22,8 +23,9 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
@ -55,28 +57,27 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
! Singlet manifold
if(singlet_manifold) then
if(singlet) then
ispin = 1
call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Omega(:,ispin),rho, &
EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPAx ',ispin,nS,Omega(:,ispin))
call print_transition_vectors(nBas,nC,nO,nV,nR,nS,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))
endif
! Triplet manifold
if(triplet_manifold) then
if(triplet) then
ispin = 2
call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), &
EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPAx ',ispin,nS,Omega(:,ispin))
call print_transition_vectors(nBas,nC,nO,nV,nR,nS,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))
endif
@ -105,15 +106,8 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
write(*,*) '-------------------------------------------------------'
write(*,*)
call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC)
if(exchange_kernel) then
EcAC(1) = 0.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(2)
end if
call ACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

26
src/QuAcK/SigmaC.f90
View File

@ -1,4 +1,4 @@
double precision function SigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
double precision function SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
! Compute diagonal of the correlation part of the self-energy
@ -7,7 +7,7 @@ double precision function SigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
! Input variables
integer,intent(in) :: x
integer,intent(in) :: p
double precision,intent(in) :: w
double precision,intent(in) :: eta
integer,intent(in) :: nBas
@ -22,7 +22,7 @@ double precision function SigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
! Local variables
integer :: i,j,a,b,p,jb
integer :: i,a,jb
double precision :: eps
! Initialize
@ -32,26 +32,18 @@ double precision function SigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
! Occupied part of the correlation self-energy
do i=nC+1,nO
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
eps = w - e(i) + Omega(jb)
SigmaC = SigmaC + 2d0*rho(x,i,jb)**2*eps/(eps**2 + eta**2)
enddo
do jb=1,nS
eps = w - e(i) + Omega(jb)
SigmaC = SigmaC + 2d0*rho(p,i,jb)**2*eps/(eps**2 + eta**2)
enddo
enddo
! Virtual part of the correlation self-energy
do a=nO+1,nBas-nR
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
eps = w - e(a) - Omega(jb)
SigmaC = SigmaC + 2d0*rho(x,a,jb)**2*eps/(eps**2 + eta**2)
enddo
do jb=1,nS
eps = w - e(a) - Omega(jb)
SigmaC = SigmaC + 2d0*rho(p,a,jb)**2*eps/(eps**2 + eta**2)
enddo
enddo

128
src/QuAcK/UCIS.f90 Normal file
View File

@ -0,0 +1,128 @@
subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int,eHF)
! Perform configuration interaction single calculation`
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip
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) :: nS(nspin)
double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart,nspin)
! Local variables
logical :: dump_matrix = .false.
logical :: dump_trans = .false.
integer :: ispin
double precision :: lambda
integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: A_sc(:,:)
double precision,allocatable :: Omega_sc(:)
integer :: nS_ab,nS_ba,nS_sf
double precision,allocatable :: A_sf(:,:)
double precision,allocatable :: Omega_sf(:)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Configuration Interaction Singles |'
write(*,*)'************************************************'
write(*,*)
! Adiabatic connection scaling
lambda = 1d0
!----------------------------!
! Spin-conserved transitions !
!----------------------------!
if(spin_conserved) then
ispin = 1
! Memory allocation
nS_aa = nS(1)
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
allocate(A_sc(nS_sc,nS_sc),Omega_sc(nS_sc))
call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eHF, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_sc)
if(dump_matrix) then
print*,'CIS matrix (spin-conserved transitions)'
call matout(nS_sc,nS_sc,A_sc)
write(*,*)
endif
call diagonalize_matrix(nS_sc,A_sc,Omega_sc)
call print_excitation('UCIS ',5,nS_sc,Omega_sc)
if(dump_trans) then
print*,'Spin-conserved CIS transition vectors'
call matout(nS_sc,nS_sc,A_sc)
write(*,*)
endif
deallocate(A_sc,Omega_sc)
endif
!-----------------------!
! Spin-flip transitions !
!-----------------------!
if(spin_flip) then
ispin = 2
! Memory allocation
nS_ab = (nO(1) - nC(1))*(nV(2) - nR(2))
nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1))
nS_sf = nS_ab + nS_ba
allocate(A_sf(nS_sf,nS_sf),Omega_sf(nS_sf))
call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,lambda,eHF, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_sf)
if(dump_matrix) then
print*,'CIS matrix (spin-conserved transitions)'
call matout(nS_sf,nS_sf,A_sf)
write(*,*)
endif
call diagonalize_matrix(nS_sf,A_sf,Omega_sf)
call print_excitation('UCIS ',6,nS_sf,Omega_sf)
if(dump_trans) then
print*,'Spin-flip CIS transition vectors'
call matout(nS_sf,nS_sf,A_sf)
write(*,*)
endif
deallocate(A_sf,Omega_sf)
endif
end subroutine UCIS

166
src/QuAcK/UG0W0.f90
View File

@ -1,6 +1,6 @@
subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, &
singlet_manifold,triplet_manifold,linearize,eta,nBas,nC,nO,nV,nR,nS, &
ENuc,EUHF,Hc,ERI_aa,ERI_ab,ERI_bb,PHF,cHF,eHF,eGW)
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,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, &
dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eGW)
! Perform unrestricted G0W0 calculation
@ -20,8 +20,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: evDyn
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip
logical,intent(in) :: linearize
double precision,intent(in) :: eta
@ -37,31 +37,28 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
double precision,intent(in) :: cHF(nBas,nBas,nspin)
double precision,intent(in) :: PHF(nBas,nBas,nspin)
double precision,intent(in) :: Hc(nBas,nBas,nspin)
double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart)
! Local variables
logical :: print_W = .true.
integer :: is
integer :: ispin
integer :: iblock
integer :: bra
integer :: ket
integer :: nSa
integer :: nSb
integer :: nSt
double precision :: EcRPA(nspin)
double precision :: EcRPA
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision,allocatable :: SigC(:,:)
double precision,allocatable :: Z(:,:)
double precision,allocatable :: Omega(:)
double precision,allocatable :: XpY_a(:,:)
double precision,allocatable :: XpY_b(:,:)
double precision,allocatable :: XmY_a(:,:)
double precision,allocatable :: XmY_b(:,:)
double precision,allocatable :: rho(:,:,:,:)
integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:,:)
double precision,allocatable :: eGWlin(:,:)
@ -80,77 +77,73 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
! Initialization
EcRPA(:) = 0d0
EcRPA = 0d0
! COHSEX approximation
if(COHSEX) write(*,*) 'COHSEX approximation activated!'
write(*,*)
if(COHSEX) then
write(*,*) 'COHSEX approximation activated!'
write(*,*)
end if
! TDA for W
if(TDA_W) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
write(*,*)
if(TDA_W) then
write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
write(*,*)
end if
! TDA
if(TDA) write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
end if
! Memory allocation
nSa = nS(1)
nSb = nS(2)
nSt = nSa + nSb
nS_aa = nS(1)
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
allocate(SigC(nBas,nspin),Z(nBas,nspin),Omega(nSt),XpY_a(nSa,nSa),XpY_b(nSb,nSb),XmY_a(nSa,nSa),XmY_b(nSb,nSb), &
rho(nBas,nBas,nSt,nspin),eGWlin(nBas,nspin))
allocate(SigC(nBas,nspin),Z(nBas,nspin),OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc), &
rho_RPA(nBas,nBas,nS_sc,nspin),eGWlin(nBas,nspin))
! Compute linear response
!-------------------!
! Compute screening !
!-------------------!
!----------------------------------------------
! alpha-alpha block
!----------------------------------------------
! Spin-conserving transition
ispin = 1
iblock = 3
ispin = 1
call linear_response(iblock,.true.,TDA_W,.false.,eta,nBas,nC(ispin),nO(ispin),nV(ispin),nR(ispin),nSa,1d0, &
eHF(:,ispin),ERI_aa,rho(:,:,1:nSa,ispin),EcRPA(ispin),Omega(1:nSa),XpY_a,XmY_a)
call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
if(print_W) call print_excitation('RPA@HF alpha',iblock,nSa,Omega(1:nSa))
if(print_W) call print_excitation('RPA@UHF ',5,nS_sc,OmRPA)
!----------------------------------------------
! alpha-beta block
!----------------------------------------------
ispin = 2
iblock = 3
call linear_response(iblock,.true.,TDA_W,.false.,eta,nBas,nC(ispin),nO(ispin),nV(ispin),nR(ispin),nSb,1d0, &
eHF(:,ispin),ERI_bb,rho(:,:,nSa+1:nSt,ispin),EcRPA(ispin),Omega(nSa+1:nSt),XpY_b,XmY_b)
if(print_W) call print_excitation('RPA@HF beta ',iblock,nSb,Omega(nSa+1:nSt))
!----------------------------------------------
! Excitation densities for alpha self-energy
!----------------------------------------------
!----------------------!
! Excitation densities !
!----------------------!
call unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ab,ERI_bb,XpY_a,XpY_b,rho)
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
!----------------------
! Compute self-energy
!----------------------
!---------------------!
! Compute self-energy !
!---------------------!
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,eHF,Omega,rho,SigC)
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC)
! Compute renormalization factor
!--------------------------------!
! Compute renormalization factor !
!--------------------------------!
! call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, &
! Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:))
call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z)
!-----------------------------------!
! Solve the quasi-particle equation !
!-----------------------------------!
! Solve the quasi-particle equation
Z(:,:) = 1d0
eGWlin(:,:) = eHF(:,:) + Z(:,:)*SigC(:,:)
if(linearize) then
@ -164,39 +157,32 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
! Find graphical solution of the QP equation
! do is=1,nspin
! call QP_graph(nBas,nC(:,is),nO(:,is),nV(:,is),nR(:,is),nS(:,is),eta,eHF(:,is),Omega(:,is), &
! rho(:,:,:,ispin),eGWlin(:,is),eGW(:,is))
! end do
do is=1,nspin
call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is),OmRPA, &
rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is))
end do
end if
! Compute RPA correlation energy
call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
! Dump results
do ispin=1,nspin
call print_G0W0(nBas,nO(ispin),eHF(:,ispin),ENuc,EUHF,SigC(:,ispin),Z(:,ispin),eGW(:,ispin),EcRPA(ispin))
end do
call print_UG0W0(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA)
! Compute the RPA correlation energy
! Free memory
! call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
! rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
deallocate(OmRPA,XpY_RPA,XmY_RPA,rho_RPA)
! Perform BSE calculation
! if(BSE) then
if(BSE) then
! call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, &
! nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,eHF,eGW,EcBSE)
! if(exchange_kernel) then
!
@ -251,6 +237,6 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
! end if
! end if
end if
end subroutine UG0W0

5
src/QuAcK/UHF.f90
View File

@ -166,7 +166,8 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUH
n_diis = min(n_diis+1,max_diis)
do ispin=1,nspin
call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin))
if(nO(ispin) > 1) call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,:,ispin),F_diis(:,:,ispin), &
err(:,:,ispin),F(:,:,ispin))
end do
! Reset DIIS if required
@ -232,6 +233,6 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUH
! Compute final UHF energy
call print_UHF(nBas,nO(:),e(:,:),c(:,:,:),ENuc,ET(:),EV(:),EJ(:),Ex(:),EUHF)
call print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF)
end subroutine UHF

158
src/QuAcK/URPAx.f90 Normal file
View File

@ -0,0 +1,158 @@
subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,e)
! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism
implicit none
include 'parameters.h'
include 'quadrature.h'
! Input variables
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip
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) :: nS(nspin)
double precision,intent(in) :: eta
double precision,intent(in) :: ENuc
double precision,intent(in) :: EUHF
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart)
! Local variables
integer :: ispin
integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: Omega_sc(:)
double precision,allocatable :: XpY_sc(:,:)
double precision,allocatable :: XmY_sc(:,:)
integer :: nS_ab,nS_ba,nS_sf
double precision,allocatable :: Omega_sf(:)
double precision,allocatable :: XpY_sf(:,:)
double precision,allocatable :: XmY_sf(:,:)
double precision :: rho_sc,rho_sf
double precision :: EcRPAx(nspin)
double precision :: EcAC(nspin)
! Hello world
write(*,*)
write(*,*)'*********************************************************************'
write(*,*)'| Unrestricted random phase approximation calculation with exchange |'
write(*,*)'*********************************************************************'
write(*,*)
! Initialization
EcRPAx(:) = 0d0
EcAC(:) = 0d0
! Spin-conserved transitions
if(spin_conserved) then
ispin = 1
! Memory allocation
nS_aa = nS(1)
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc))
call unrestricted_linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPAx(ispin),Omega_sc,XpY_sc,XmY_sc)
call print_excitation('URPAx ',5,nS_sc,Omega_sc)
call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
Omega_sc,XpY_sc,XmY_sc)
deallocate(Omega_sc,XpY_sc,XmY_sc)
endif
! Spin-flip transitions
if(spin_flip) then
ispin = 2
! Memory allocation
nS_ab = (nO(1) - nC(1))*(nV(2) - nR(2))
nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1))
nS_sf = nS_ab + nS_ba
allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf))
call unrestricted_linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,1d0,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sf,rho_sf,EcRPAx(ispin),Omega_sf,XpY_sf,XmY_sf)
call print_excitation('URPAx ',6,nS_sf,Omega_sf)
! call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
deallocate(Omega_sf,XpY_sf,XmY_sf)
endif
! if(exchange_kernel) then
! EcRPAx(1) = 0.5d0*EcRPAx(1)
! EcRPAx(2) = 1.5d0*EcRPAx(2)
! end if
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy (spin-conserved) =',EcRPAx(1)
write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy (spin-flip) =',EcRPAx(2)
write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy =',EcRPAx(1) + EcRPAx(2)
write(*,'(2X,A50,F20.10)') 'Tr@URPAx total energy =',ENuc + EUHF + EcRPAx(1) + EcRPAx(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Compute the correlation energy via the adiabatic connection
! if(doACFDT) then
! write(*,*) '-------------------------------------------------------'
! write(*,*) 'Adiabatic connection version of RPAx correlation energy'
! write(*,*) '-------------------------------------------------------'
! write(*,*)
! call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,singlet,triplet,eta, &
! nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC)
! if(exchange_kernel) then
! EcAC(1) = 0.5d0*EcAC(1)
! EcAC(2) = 1.5d0*EcAC(2)
! end if
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (singlet) =',EcAC(1)
! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (triplet) =',EcAC(2)
! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy =',EcAC(1) + EcAC(2)
! write(*,'(2X,A50,F20.10)') 'AC@RPAx total energy =',ENuc + EUHF + EcAC(1) + EcAC(2)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
! end if
end subroutine URPAx

48
src/QuAcK/USigmaC.f90 Normal file
View File

@ -0,0 +1,48 @@
double precision function USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
! Compute diagonal of the correlation part of the self-energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: p
double precision,intent(in) :: w
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,jb
double precision :: eps
! Initialize
USigmaC = 0d0
! Occupied part of the correlation self-energy
do i=nC+1,nO
do jb=1,nS
eps = w - e(i) + Omega(jb)
USigmaC = uSigmaC + rho(p,i,jb)**2*eps/(eps**2 + eta**2)
end do
end do
do a=nO+1,nBas-nR
do jb=1,nS
eps = w - e(a) - Omega(jb)
USigmaC = USigmaC + rho(p,a,jb)**2*eps/(eps**2 + eta**2)
end do
end do
end function USigmaC

156
src/QuAcK/UdRPA.f90 Normal file
View File

@ -0,0 +1,156 @@
subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,e)
! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism
implicit none
include 'parameters.h'
include 'quadrature.h'
! Input variables
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip
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) :: nS(nspin)
double precision,intent(in) :: eta
double precision,intent(in) :: ENuc
double precision,intent(in) :: EUHF
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart)
! Local variables
integer :: ispin
integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: Omega_sc(:)
double precision,allocatable :: XpY_sc(:,:)
double precision,allocatable :: XmY_sc(:,:)
integer :: nS_ab,nS_ba,nS_sf
double precision,allocatable :: Omega_sf(:)
double precision,allocatable :: XpY_sf(:,:)
double precision,allocatable :: XmY_sf(:,:)
double precision :: rho_sc,rho_sf
double precision :: EcRPA(nspin)
double precision :: EcAC(nspin)
! Hello world
write(*,*)
write(*,*)'**************************************************************'
write(*,*)'| Unrestricted direct random phase approximation calculation |'
write(*,*)'**************************************************************'
write(*,*)
! Initialization
EcRPA(:) = 0d0
EcAC(:) = 0d0
! Spin-conserved transitions
if(spin_conserved) then
ispin = 1
! Memory allocation
nS_aa = nS(1)
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc))
call unrestricted_linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc)
call print_excitation('URPA ',5,nS_sc,Omega_sc)
call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
Omega_sc,XpY_sc,XmY_sc)
endif
! Spin-flip transitions
if(spin_flip) then
ispin = 2
! Memory allocation
nS_ab = (nO(1) - nC(1))*(nV(2) - nR(2))
nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1))
nS_sf = nS_ab + nS_ba
allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf))
call unrestricted_linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,1d0,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sf,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf)
call print_excitation('URPA ',6,nS_sf,Omega_sf)
! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif
! if(exchange_kernel) then
! EcRPA(1) = 0.5d0*EcRPA(1)
! EcRPA(2) = 1.5d0*EcRPA(2)
! end if
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@URPA correlation energy (spin-conserved) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@URPA correlation energy (spin-flip) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@URPA correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@URPA total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Compute the correlation energy via the adiabatic connection
! if(doACFDT) then
! write(*,*) '-------------------------------------------------------'
! write(*,*) 'Adiabatic connection version of RPA correlation energy'
! write(*,*) '-------------------------------------------------------'
! write(*,*)
! call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,singlet,triplet,eta, &
! nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC)
! if(exchange_kernel) then
! EcAC(1) = 0.5d0*EcAC(1)
! EcAC(2) = 1.5d0*EcAC(2)
! end if
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy (singlet) =',EcAC(1)
! write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy (triplet) =',EcAC(2)
! write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy =',EcAC(1) + EcAC(2)
! write(*,'(2X,A50,F20.10)') 'AC@RPA total energy =',ENuc + EUHF + EcAC(1) + EcAC(2)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
! end if
end subroutine UdRPA

37
src/QuAcK/RPA.f90 → src/QuAcK/dRPA.f90
View File

@ -1,5 +1,4 @@
subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e)
subroutine dRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform a direct random phase approximation calculation
@ -9,10 +8,11 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
! Input variables
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
@ -22,8 +22,9 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
@ -40,7 +41,7 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
write(*,*)
write(*,*)'***********************************************'
write(*,*)'| random-phase approximation calculation |'
write(*,*)'| Random-phase approximation calculation |'
write(*,*)'***********************************************'
write(*,*)
@ -55,32 +56,34 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
! Singlet manifold
if(singlet_manifold) then
if(singlet) then
ispin = 1
call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin))
call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif
! Triplet manifold
if(triplet_manifold) then
if(triplet) then
ispin = 2
call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin))
call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif
if(exchange_kernel) then
EcRPA(1) = 0.5d0*EcRPA(1)
EcRPA(2) = 1.5d0*EcRPA(1)
EcRPA(2) = 1.5d0*EcRPA(2)
end if
@ -103,13 +106,13 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
write(*,*) '------------------------------------------------------'
write(*,*)
call ACFDT(exchange_kernel,.false.,.true.,.false.,.false.,.false.,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC)
call ACFDT(exchange_kernel,.false.,.true.,.false.,TDA,.false.,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC)
if(exchange_kernel) then
EcAC(1) = 0.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(2)
end if
@ -125,4 +128,4 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
end if
end subroutine RPA
end subroutine dRPA

50
src/QuAcK/dUSigmaC.f90 Normal file
View File

@ -0,0 +1,50 @@
double precision function dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
! Compute the derivative of the correlation part of the self-energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: p
double precision,intent(in) :: w
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,jb
double precision :: eps
! Initialize
dUSigmaC = 0d0
! Occupied part of the correlation self-energy
do i=nC+1,nO
do jb=1,nS
eps = w - e(i) + Omega(jb)
dUSigmaC = dUSigmaC + rho(p,i,jb)**2*(eps/(eps**2 + eta**2))**2
end do
end do
! Virtual part of the correlation self-energy
do a=nO+1,nBas-nR
do jb=1,nS
eps = w - e(a) - Omega(jb)
dUSigmaC = dUSigmaC + rho(p,a,jb)**2*(eps/(eps**2 + eta**2))**2
end do
end do
end function dUSigmaC

5
src/QuAcK/evGF2.f90
View File

@ -1,5 +1,5 @@
subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifold, &
linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform eigenvalue self-consistent second-order Green function calculation
@ -30,6 +30,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet_manifold
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
@ -168,7 +169,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet_manifold
if(BSE) then
call BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF2,EcBSE)
call BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE)
end if

19
src/QuAcK/evGT.f90
View File

@ -1,6 +1,6 @@
subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eG0T0)
BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, &
eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF,eG0T0)
! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT)
@ -21,8 +21,8 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: evDyn
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
integer,intent(in) :: nBas
@ -35,6 +35,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: eG0T0(nBas)
@ -260,10 +261,8 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
if(BSE) then
allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin))
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eGT,eGT,Omega,XpY,XmY,rho,EcRPA,EcBSE)
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGT,eGT,EcRPA,EcBSE)
if(exchange_kernel) then
@ -297,8 +296,8 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eGT,eGT,Omega,XpY,XmY,rho,EcAC)
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eGT,eGT,EcAC)
if(exchange_kernel) then

121
src/QuAcK/evGW.f90
View File

@ -1,5 +1,6 @@
subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,G0W,GW0, &
singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0W0)
subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI, &
dipole_int,PHF,cHF,eHF,eG0W0)
! Perform self-consistent eigenvalue-only GW calculation
@ -26,8 +27,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
logical,intent(in) :: evDyn
logical,intent(in) :: G0W
logical,intent(in) :: GW0
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eHF(nBas)
@ -35,8 +36,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: eG0W0(nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: H(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
@ -46,7 +47,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
integer :: n_diis
double precision :: rcond
double precision :: Conv
double precision :: EcRPA(nspin)
double precision :: EcRPA
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision :: EcGM
@ -57,11 +58,10 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
double precision,allocatable :: eOld(:)
double precision,allocatable :: Z(:)
double precision,allocatable :: SigC(:)
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
double precision,allocatable :: rho(:,:,:,:)
double precision,allocatable :: rhox(:,:,:,:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
! Hello world
@ -73,23 +73,45 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
! SOSEX correction
if(SOSEX) write(*,*) 'SOSEX correction activated!'
write(*,*)
if(SOSEX) then
write(*,*) 'SOSEX correction activated but BUG!'
stop
end if
! COHSEX approximation
if(COHSEX) write(*,*) 'COHSEX approximation activated!'
write(*,*)
if(COHSEX) then
write(*,*) 'COHSEX approximation activated!'
write(*,*)
end if
! TDA for W
if(TDA_W) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
write(*,*)
if(TDA_W) then
write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
write(*,*)
end if
! TDA
if(TDA) write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
end if
! GW0
if(GW0) then
write(*,*) 'GW0 scheme activated!'
write(*,*)
end if
! G0W
if(G0W) then
write(*,*) 'G0W scheme activated!'
write(*,*)
end if
! Linear mixing
@ -98,10 +120,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
! Memory allocation
allocate(eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas),Omega(nS,nspin), &
XpY(nS,nS,nspin),XmY(nS,nS,nspin), &
rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), &
error_diis(nBas,max_diis),e_diis(nBas,max_diis))
allocate(eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), &
rho_RPA(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
! Initialization
@ -121,36 +141,30 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
do while(Conv > thresh .and. nSCF <= maxSCF)
! Compute linear response
! Compute screening
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, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
endif
! Compute spectral weights
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
! Compute correlation part of the self-energy
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rhox(:,:,:,ispin))
! Correlation self-energy
if(G0W) then
call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC)
call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:))
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)
else
call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC)
call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:))
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)
endif
@ -164,8 +178,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
! Print results
! call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW)
call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
! Linear mixing or DIIS extrapolation
@ -215,23 +228,22 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
endif
! Dump the RPA correlation energy
! Deallocate memory
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
deallocate(eOld,Z,SigC,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error_diis,e_diis)
! Perform BSE calculation
if(BSE) then
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGW,eGW,EcBSE)
if(exchange_kernel) then
EcBSE(1) = 0.5d0*EcBSE(1)
EcBSE(2) = 1.5d0*EcBSE(2)
end if
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
@ -258,8 +270,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcAC)
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

308
src/QuAcK/evUGW.f90 Normal file
View File

@ -0,0 +1,308 @@
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, &
ERHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eG0W0)
! Perform self-consistent eigenvalue-only GW calculation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
double precision,intent(in) :: thresh
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: doXBS
logical,intent(in) :: COHSEX
logical,intent(in) :: BSE
logical,intent(in) :: TDA_W
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: evDyn
logical,intent(in) :: G0W
logical,intent(in) :: GW0
logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip
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) :: nS(nspin)
double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: cHF(nBas,nBas,nspin)
double precision,intent(in) :: PHF(nBas,nBas,nspin)
double precision,intent(in) :: eG0W0(nBas,nspin)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart)
! Local variables
logical :: linear_mixing
integer :: is
integer :: ispin
integer :: nSCF
integer :: n_diis
double precision :: rcond(nspin)
double precision :: Conv
double precision :: EcRPA
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision :: EcGM
double precision :: alpha
double precision,allocatable :: error_diis(:,:,:)
double precision,allocatable :: e_diis(:,:,:)
double precision,allocatable :: eGW(:,:)
double precision,allocatable :: eOld(:,:)
double precision,allocatable :: Z(:,:)
integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: SigC(:,:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:,:)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Self-consistent evGW calculation |'
write(*,*)'************************************************'
write(*,*)
! COHSEX approximation
if(COHSEX) then
write(*,*) 'COHSEX approximation activated!'
write(*,*)
end if
! TDA for W
if(TDA_W) then
write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
write(*,*)
end if
! TDA
if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
end if
! GW0
if(GW0) then
write(*,*) 'GW0 scheme activated!'
write(*,*)
end if
! G0W
if(G0W) then
write(*,*) 'G0W scheme activated!'
write(*,*)
end if
! Linear mixing
linear_mixing = .false.
alpha = 0.2d0
! Memory allocation
nS_aa = nS(1)
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
allocate(eGW(nBas,nspin),eOld(nBas,nspin),Z(nBas,nspin),SigC(nBas,nspin),OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc), &
XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin),error_diis(nBas,max_diis,nspin),e_diis(nBas,max_diis,nspin))
! Initialization
nSCF = 0
ispin = 1
n_diis = 0
Conv = 1d0
e_diis(:,:,:) = 0d0
error_diis(:,:,:) = 0d0
eGW(:,:) = eG0W0(:,:)
eOld(:,:) = eGW(:,:)
Z(:,:) = 1d0
!------------------------------------------------------------------------
! Main loop
!------------------------------------------------------------------------
do while(Conv > thresh .and. nSCF <= maxSCF)
! Compute screening
if(.not. GW0 .or. nSCF == 0) then
call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
endif
!----------------------!
! Excitation densities !
!----------------------!
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
!------------------------------------------------!
! Compute self-energy and renormalization factor !
!------------------------------------------------!
if(G0W) then
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC)
call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z)
else
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC)
call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z)
endif
!-----------------------------------!
! Solve the quasi-particle equation !
!-----------------------------------!
eGW(:,:) = eHF(:,:) + SigC(:,:)
! Convergence criteria
Conv = maxval(abs(eGW(:,:) - eOld(:,:)))
! Print results
call print_evUGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA)
! Linear mixing or DIIS extrapolation
if(linear_mixing) then
eGW(:,:) = alpha*eGW(:,:) + (1d0 - alpha)*eOld(:,:)
else
n_diis = min(n_diis+1,max_diis)
do is=1,nspin
call DIIS_extrapolation(rcond(ispin),nBas,nBas,n_diis,error_diis(:,:,is), &
e_diis(:,:,is),eGW(:,is)-eOld(:,is),eGW(:,is))
end do
! Reset DIIS if required
if(minval(rcond(:)) < 1d-15) n_diis = 0
endif
! Save quasiparticles energy for next cycle
eOld(:,:) = eGW(:,:)
! Increment
nSCF = nSCF + 1
enddo
!------------------------------------------------------------------------
! End main loop
!------------------------------------------------------------------------
! Plot stuff
! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin))
! Did it actually converge?
if(nSCF == maxSCF+1) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
stop
endif
! Deallocate memory
deallocate(eOld,Z,SigC,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error_diis,e_diis)
! Perform BSE calculation
if(BSE) then
call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,eGW,eGW,EcBSE)
! if(exchange_kernel) then
! EcBSE(1) = 0.5d0*EcBSE(1)
! EcBSE(2) = 1.5d0*EcBSE(2)
! end if
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy (singlet) =',EcBSE(1)
! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy (triplet) =',EcBSE(2)
! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2)
! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
! Compute the BSE correlation energy via the adiabatic connection
! if(doACFDT) then
! write(*,*) '------------------------------------------------------'
! write(*,*) 'Adiabatic connection version of BSE correlation energy'
! write(*,*) '------------------------------------------------------'
! write(*,*)
! if(doXBS) then
! write(*,*) '*** scaled screening version (XBS) ***'
! write(*,*)
! end if
! call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcAC)
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (singlet) =',EcAC(1)
! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (triplet) =',EcAC(2)
! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy =',EcAC(1) + EcAC(2)
! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
! end if
endif
end subroutine evUGW

10
src/QuAcK/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,rho,EcRPA,Omega,XpY,XmY)
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)
! Compute linear response
@ -12,8 +12,10 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: lambda
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: rho(nBas,nBas,nS)
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
@ -42,7 +44,7 @@ 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,rho,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
@ -51,7 +53,7 @@ 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,rho,B)
if(BSE) call Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega_RPA,rho_RPA,B)
end if

4
src/QuAcK/plot_GW.f90
View File

@ -1,4 +1,4 @@
subroutine plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega,rho,rhox)
subroutine plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega,rho)
! Dump several GW quantities for external plotting
@ -8,7 +8,7 @@ subroutine plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega,rho,rhox)
! Input variables
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eHF(nBas),eGW(nBas),Omega(nS),rho(nBas,nBas,nS),rhox(nBas,nBas,nS)
double precision,intent(in) :: eHF(nBas),eGW(nBas),Omega(nS),rho(nBas,nBas,nS)
! Local variables

10
src/QuAcK/print_G0W0.f90
View File

@ -40,11 +40,11 @@ subroutine print_G0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM)
write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',eGW(LUMO)*HaToeV
write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A30,F15.6)') 'RPA@G0W0 total energy =',ENuc + EHF + EcRPA
! write(*,'(2X,A30,F15.6)') 'RPA@G0W0 correlation energy =',EcRPA
! write(*,'(2X,A30,F15.6)') 'GM@G0W0 total energy =',ENuc + EHF + EcGM
! write(*,'(2X,A30,F15.6)') 'GM@G0W0 correlation energy =',EcGM
! write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A30,F15.6)') 'RPA@G0W0 total energy =',ENuc + EHF + EcRPA
write(*,'(2X,A30,F15.6)') 'RPA@G0W0 correlation energy =',EcRPA
write(*,'(2X,A30,F15.6)') 'GM@G0W0 total energy =',ENuc + EHF + EcGM
write(*,'(2X,A30,F15.6)') 'GM@G0W0 correlation energy =',EcGM
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end subroutine print_G0W0

71
src/QuAcK/print_UG0W0.f90
View File

@ -1,44 +1,71 @@
subroutine print_UG0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA)
subroutine print_UG0W0(nBas,nO,e,ENuc,EHF,SigC,Z,eGW,EcRPA)
! Print one-electron energies and other stuff for G0W0
implicit none
include 'parameters.h'
integer,intent(in) :: nBas,nO
integer,intent(in) :: nBas
integer,intent(in) :: nO(nspin)
double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF
double precision,intent(in) :: EcRPA
double precision,intent(in) :: e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas)
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: SigC(nBas,nspin)
double precision,intent(in) :: Z(nBas,nspin)
double precision,intent(in) :: eGW(nBas,nspin)
integer :: x,HOMO,LUMO
double precision :: Gap
integer :: p
integer :: ispin
double precision :: HOMO(nspin)
double precision :: LUMO(nspin)
double precision :: Gap(nspin)
! HOMO and LUMO
HOMO = nO
LUMO = HOMO + 1
Gap = eGW(LUMO)-eGW(HOMO)
do ispin=1,nspin
if(nO(ispin) > 0) then
HOMO(ispin) = eGW(nO(ispin),ispin)
LUMO(ispin) = eGW(nO(ispin)+1,ispin)
Gap(ispin) = LUMO(ispin) - HOMO(ispin)
else
HOMO(ispin) = 0d0
LUMO(ispin) = e(1,ispin)
Gap(ispin) = 0d0
end if
end do
! Dump results
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)' One-shot G0W0 calculation'
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
'|','#','|','e_HF (eV)','|','Sigma_c (eV)','|','Z','|','e_QP (eV)','|'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,*)' Unrestricted one-shot G0W0 calculation (eV)'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(A1,A3,A1,A30,A1,A30,A1,A30,A1,A30,A1)') &
'|',' ','|','e_HF ','|','Sig_c ','|','Z ','|','e_QP ','|'
write(*,'(A1,A3,A1,2A15,A1,2A15,A1,2A15,A1,2A15,A1)') &
'|','#','|','up ','dw ','|','up ','dw ','|','up ','dw ','|','up ','dw ','|'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
do x=1,nBas
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',x,'|',e(x)*HaToeV,'|',SigmaC(x)*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|'
do p=1,nBas
write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') &
'|',p,'|',e(p,1)*HaToeV,e(p,2)*HaToeV,'|',SigC(p,1)*HaToeV,SigC(p,2)*HaToeV,'|', &
Z(p,1),Z(p,2),'|',eGW(p,1)*HaToeV,eGW(p,2)*HaToeV,'|'
enddo
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A30,F15.6)') 'G0W0 HOMO energy (eV):',eGW(HOMO)*HaToeV
write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',eGW(LUMO)*HaToeV
write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(2X,A30,F15.6)') 'G0W0 HOMO energy (eV):',maxval(HOMO(:))*HaToeV
write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',minval(LUMO(:))*HaToeV
write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(2X,A30,F15.6)') 'RPA@G0W0 total energy =',ENuc + EHF + EcRPA
write(*,'(2X,A30,F15.6)') 'RPA@G0W0 correlation energy =',EcRPA
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,*)
end subroutine print_UG0W0

31
src/QuAcK/print_UHF.f90
View File

@ -16,18 +16,24 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF)
double precision,intent(in) :: Ex(nspin)
double precision,intent(in) :: EUHF
integer :: HOMO(nspin)
integer :: LUMO(nspin)
integer :: ispin
double precision :: HOMO(nspin)
double precision :: LUMO(nspin)
double precision :: Gap(nspin)
! HOMO and LUMO
HOMO(:) = nO(:)
LUMO(:) = HOMO(:) + 1
Gap(1) = e(LUMO(1),1) - e(HOMO(1),1)
Gap(2) = e(LUMO(2),2) - e(HOMO(2),2)
do ispin=1,nspin
if(nO(ispin) > 0) then
HOMO(ispin) = e(nO(ispin),ispin)
LUMO(ispin) = e(nO(ispin)+1,ispin)
Gap(ispin) = LUMO(ispin) - HOMO(ispin)
else
HOMO(ispin) = 0d0
LUMO(ispin) = e(1,ispin)
Gap(ispin) = 0d0
end if
end do
! Dump results
@ -62,12 +68,12 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF)
write(*,'(A40,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au'
write(*,'(A40,1X,F16.10,A3)') ' UHF energy: ',EUHF + ENuc,' au'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,F13.6,A3)') ' UHF HOMO a energy:',e(HOMO(1),1)*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' UHF LUMO a energy:',e(LUMO(1),1)*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' UHF HOMO a energy:',HOMO(1)*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' UHF LUMO a energy:',LUMO(1)*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' UHF HOMOa-LUMOa gap:',Gap(1)*HatoeV,' eV'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,F13.6,A3)') ' UHF HOMO b energy:',e(HOMO(2),2)*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' UHF LUMO b energy:',e(LUMO(2),2)*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' UHF HOMO b energy:',HOMO(2)*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' UHF LUMO b energy:',LUMO(2)*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' UHF HOMOb-LUMOb gap :',Gap(2)*HatoeV,' eV'
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
@ -78,6 +84,7 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF)
write(*,'(A50)') 'UHF spin-up orbital coefficients '
write(*,'(A50)') '-----------------------------------------'
call matout(nBas,nBas,c(:,:,1))
write(*,*)
write(*,'(A50)') '-----------------------------------------'
write(*,'(A50)') 'UHF spin-down orbital coefficients '
write(*,'(A50)') '-----------------------------------------'

26
src/QuAcK/print_evGW.f90
View File

@ -1,4 +1,4 @@
subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW)
subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigC,Z,eGW,EcRPA,EcGM)
! Print one-electron energies and other stuff for evGW
@ -8,9 +8,15 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW)
integer,intent(in) :: nBas,nO,nSCF
double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF
double precision,intent(in) :: Conv,e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas)
double precision,intent(in) :: Conv
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: SigC(nBas)
double precision,intent(in) :: Z(nBas)
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: EcRPA
double precision,intent(in) :: EcGM
integer :: x,HOMO,LUMO
integer :: p,HOMO,LUMO
double precision :: Gap
! HOMO and LUMO
@ -32,9 +38,9 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW)
'|','#','|','e_HF (eV)','|','Sigma_c (eV)','|','Z','|','e_QP (eV)','|'
write(*,*)'-------------------------------------------------------------------------------'
do x=1,nBas
do p=1,nBas
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',x,'|',e(x)*HaToeV,'|',SigmaC(x)*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|'
'|',p,'|',e(p)*HaToeV,'|',SigC(p)*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|'
enddo
write(*,*)'-------------------------------------------------------------------------------'
@ -45,11 +51,11 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW)
write(*,'(2X,A30,F15.6)') 'evGW LUMO energy (eV):',eGW(LUMO)*HaToeV
write(*,'(2X,A30,F15.6)') 'evGW HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A30,F15.6)') 'RPA@evGW total energy =',ENuc + EHF + EcRPA
! write(*,'(2X,A30,F15.6)') 'RPA@evGW correlation energy =',EcRPA
! write(*,'(2X,A30,F15.6)') 'GM@evGW total energy =',ENuc + EHF + EcGM
! write(*,'(2X,A30,F15.6)') 'GM@evGW correlation energy =',EcGM
! write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A30,F15.6)') 'RPA@evGW total energy =',ENuc + EHF + EcRPA
write(*,'(2X,A30,F15.6)') 'RPA@evGW correlation energy =',EcRPA
write(*,'(2X,A30,F15.6)') 'GM@evGW total energy =',ENuc + EHF + EcGM
write(*,'(2X,A30,F15.6)') 'GM@evGW correlation energy =',EcGM
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end subroutine print_evGW

81
src/QuAcK/print_evUGW.f90 Normal file
View File

@ -0,0 +1,81 @@
subroutine print_evUGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigC,Z,eGW,EcRPA)
! Print one-electron energies and other stuff for evGW
implicit none
include 'parameters.h'
integer,intent(in) :: nBas
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nSCF
double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF
double precision,intent(in) :: EcRPA
double precision,intent(in) :: Conv
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: SigC(nBas,nspin)
double precision,intent(in) :: Z(nBas,nspin)
double precision,intent(in) :: eGW(nBas,nspin)
integer :: p
integer :: ispin
double precision :: HOMO(nspin)
double precision :: LUMO(nspin)
double precision :: Gap(nspin)
! HOMO and LUMO
do ispin=1,nspin
if(nO(ispin) > 0) then
HOMO(ispin) = eGW(nO(ispin),ispin)
LUMO(ispin) = eGW(nO(ispin)+1,ispin)
Gap(ispin) = LUMO(ispin) - HOMO(ispin)
else
HOMO(ispin) = 0d0
LUMO(ispin) = e(1,ispin)
Gap(ispin) = 0d0
end if
end do
! Dump results
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
if(nSCF < 10) then
write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent evG',nSCF,'W',nSCF,' calculation'
else
write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent evG',nSCF,'W',nSCF,' calculation'
endif
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(A1,A3,A1,A30,A1,A30,A1,A30,A1,A30,A1)') &
'|',' ','|','e_HF ','|','Sig_c ','|','Z ','|','e_QP ','|'
write(*,'(A1,A3,A1,2A15,A1,2A15,A1,2A15,A1,2A15,A1)') &
'|','#','|','up ','dw ','|','up ','dw ','|','up ','dw ','|','up ','dw ','|'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
do p=1,nBas
write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') &
'|',p,'|',e(p,1)*HaToeV,e(p,2)*HaToeV,'|',SigC(p,1)*HaToeV,SigC(p,2)*HaToeV,'|', &
Z(p,1),Z(p,2),'|',eGW(p,1)*HaToeV,eGW(p,2)*HaToeV,'|'
enddo
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(2X,A10,I3)') 'Iteration ',nSCF
write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(2X,A30,F15.6)') 'evGW HOMO energy (eV):',maxval(HOMO(:))*HaToeV
write(*,'(2X,A30,F15.6)') 'evGW LUMO energy (eV):',minval(LUMO(:))*HaToeV
write(*,'(2X,A30,F15.6)') 'evGW HOMO-LUMO gap (eV):',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(2X,A30,F15.6)') 'RPA@evGW total energy =',ENuc + EHF + EcRPA
write(*,'(2X,A30,F15.6)') 'RPA@evGW correlation energy =',EcRPA
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,*)
end subroutine print_evUGW

12
src/QuAcK/print_excitation.f90
View File

@ -7,24 +7,26 @@ subroutine print_excitation(method,ispin,nS,Omega)
! Input variables
character*12,intent(in) :: method
character*12,intent(in) :: method
integer,intent(in) :: ispin,nS
double precision,intent(in) :: Omega(nS)
! Local variables
character*7 :: spin_manifold
character*14 :: spin_manifold
integer,parameter :: maxS = 32
integer :: ia
if(ispin == 1) spin_manifold = 'singlet'
if(ispin == 2) spin_manifold = 'triplet'
if(ispin == 3) spin_manifold = 'alp-bet'
if(ispin == 4) spin_manifold = 'alp-alp'
if(ispin == 3) spin_manifold = 'alpha-beta'
if(ispin == 4) spin_manifold = 'alpha-alpha'
if(ispin == 5) spin_manifold = 'spin-conserved'
if(ispin == 6) spin_manifold = 'spin-flip'
write(*,*)
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A1,1X,A14,A14,A7,A9,A15)')'|',method,' calculation: ',spin_manifold,' manifold',' |'
write(*,'(1X,A14,A14,A14,A9)') method,' calculation: ',spin_manifold,' manifold'
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
'|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'

5
src/QuAcK/print_qsGW.f90
View File

@ -70,8 +70,8 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig
write(*,'(2X,A30,F15.6)') 'qsGW GM total energy =',EqsGW + ENuc + EcGM
write(*,'(2X,A30,F15.6)') 'qsGW exchange energy =',Ex
write(*,'(2X,A30,F15.6)') 'qsGW correlation energy =',Ec
! write(*,'(2X,A30,F15.6)') 'RPA@qsGW correlation energy =',EcRPA
! write(*,'(2X,A30,F15.6)') 'GM@qsGW correlation energy =',EcGM
write(*,'(2X,A30,F15.6)') 'RPA@qsGW correlation energy =',EcRPA
write(*,'(2X,A30,F15.6)') 'GM@qsGW correlation energy =',EcGM
write(*,*)'-------------------------------------------'
write(*,*)
@ -95,7 +95,6 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig
write(*,'(A32,1X,F16.10)') ' Electronic energy ',EqsGW
write(*,'(A32,1X,F16.10)') ' Nuclear repulsion ',ENuc
write(*,'(A32,1X,F16.10)') ' qsGW energy ',ENuc + EqsGW
! write(*,'(A32,1X,F16.10)') ' RPA corr. energy ',EcRPA
write(*,'(A50)') '---------------------------------------'
write(*,*)

77
src/QuAcK/print_transition_vectors.f90
View File

@ -1,4 +1,4 @@
subroutine print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega,XpY,XmY)
subroutine print_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY)
! Print transition vectors for linear response calculation
@ -7,44 +7,86 @@ subroutine print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega,XpY,XmY)
! Input variables
logical,intent(in) :: spin_allowed
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 :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: XpY(nS,nS)
double precision,intent(in) :: XmY(nS,nS)
! Local variables
logical :: debug = .false.
integer :: ia,jb,i,j,a,b
integer :: ixyz
integer,parameter :: maxS = 10
double precision,parameter :: thres_vec = 0.1d0
double precision :: norm
double precision,parameter :: thres_vec = 0.1d0
double precision,allocatable :: X(:)
double precision,allocatable :: Y(:)
double precision,allocatable :: f(:,:)
double precision,allocatable :: os(:)
! Memory allocation
allocate(X(nS),Y(nS))
allocate(X(nS),Y(nS),f(nS,ncart),os(nS))
! Compute dipole moments and oscillator strengths
f(:,:) = 0d0
if(spin_allowed) then
do ia=1,nS
do ixyz=1,ncart
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
f(ia,ixyz) = f(ia,ixyz) + dipole_int(j,b,ixyz)*XpY(ia,jb)
end do
end do
end do
end do
f(:,:) = sqrt(2d0)*f(:,:)
do ia=1,nS
os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2)
end do
if(debug) then
write(*,*) '------------------------'
write(*,*) ' Dipole moments (X Y Z) '
write(*,*) '------------------------'
call matout(nS,ncart,f)
write(*,*)
write(*,*) '----------------------'
write(*,*) ' Oscillator strengths '
write(*,*) '----------------------'
call matout(nS,1,os)
write(*,*)
end if
end if
! Print details about excitations
write(*,*)
do ia=1,min(nS,maxS)
X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
norm = 0d0
do jb=1,nS
norm = norm + X(jb)*X(jb)
end do
norm = sqrt(norm)
print*,'--------------------------------'
write(*,'(A15,I3,A2,F10.6,A3)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV'
print*,'--------------------------------'
print*,'---------------------------------------------'
write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A1)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' (f = ',os(ia),')'
print*,'---------------------------------------------'
jb = 0
do j=nC+1,nO
@ -54,12 +96,6 @@ subroutine print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega,XpY,XmY)
end do
end do
norm = 0d0
do jb=1,nS
norm = norm + Y(jb)*Y(jb)
end do
norm = sqrt(norm)
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
@ -71,4 +107,7 @@ subroutine print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega,XpY,XmY)
end do
write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:))
write(*,*)
end subroutine print_transition_vectors

148
src/QuAcK/print_unrestricted_transition_vectors.f90 Normal file
View File

@ -0,0 +1,148 @@
subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb, &
Omega,XpY,XmY)
! Print transition vectors for linear response calculation
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: spin_allowed
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) :: nS(nspin)
integer,intent(in) :: nSa
integer,intent(in) :: nSb
integer,intent(in) :: nSt
double precision :: dipole_int_aa(nBas,nBas,ncart)
double precision :: dipole_int_bb(nBas,nBas,ncart)
double precision,intent(in) :: Omega(nSt)
double precision,intent(in) :: XpY(nSt,nSt)
double precision,intent(in) :: XmY(nSt,nSt)
! Local variables
logical :: debug = .false.
integer :: ia,jb,i,j,a,b
integer :: ixyz
integer :: ispin
integer,parameter :: maxS = 10
double precision :: norm
double precision,parameter :: thres_vec = 0.1d0
double precision,allocatable :: X(:)
double precision,allocatable :: Y(:)
double precision,allocatable :: f(:,:)
double precision,allocatable :: os(:)
! Memory allocation
allocate(X(nSt),Y(nSt),f(nSt,ncart),os(nSt))
! Compute dipole moments and oscillator strengths
f(:,:) = 0d0
if(spin_allowed) then
do ia=1,nSt
do ixyz=1,ncart
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
f(ia,ixyz) = f(ia,ixyz) + dipole_int_aa(j,b,ixyz)*XpY(ia,jb)
end do
end do
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
f(ia,ixyz) = f(ia,ixyz) + dipole_int_bb(j,b,ixyz)*XpY(ia,nSa+jb)
end do
end do
end do
end do
do ia=1,nSt
os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2)
end do
if(debug) then
write(*,*) '----------------'
write(*,*) ' Dipole moments '
write(*,*) '----------------'
call matout(nSt,ncart,f(:,:))
write(*,*)
write(*,*) '----------------------'
write(*,*) ' Oscillator strengths '
write(*,*) '----------------------'
call matout(nSt,1,os(:))
write(*,*)
end if
end if
! Print details about excitations
do ia=1,min(nSt,maxS)
X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
print*,'---------------------------------------------'
write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A1)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' (f = ',os(ia),')'
print*,'---------------------------------------------'
! Spin-up transitions
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'A = ',X(jb)/sqrt(2d0)
end do
end do
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'A = ',Y(jb)/sqrt(2d0)
end do
end do
! Spin-down transitions
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B -> ',b,'B = ',X(jb)/sqrt(2d0)
end do
end do
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B <- ',b,'B = ',Y(jb)/sqrt(2d0)
end do
end do
write(*,*)
end do
write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:))
write(*,*)
end subroutine print_unrestricted_transition_vectors

95
src/QuAcK/qsGW.f90
View File

@ -1,6 +1,6 @@
subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,G0W,GW0,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF)
subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V, &
Hc,ERI_AO_basis,ERI_MO_basis,dipole_int,PHF,cHF,eHF)
! Perform a quasiparticle self-consistent GW calculation
@ -25,8 +25,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
logical,intent(in) :: evDyn
logical,intent(in) :: G0W
logical,intent(in) :: GW0
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ENuc
@ -41,6 +41,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO_basis(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
@ -49,7 +50,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
integer :: ispin
integer :: n_diis
double precision :: EqsGW
double precision :: EcRPA(nspin)
double precision :: EcRPA
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision :: EcGM
@ -58,11 +59,10 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
double precision,external :: trace_matrix
double precision,allocatable :: error_diis(:,:)
double precision,allocatable :: F_diis(:,:)
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
double precision,allocatable :: rho(:,:,:,:)
double precision,allocatable :: rhox(:,:,:,:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: c(:,:)
double precision,allocatable :: cp(:,:)
double precision,allocatable :: eGW(:)
@ -96,29 +96,37 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! SOSEX correction
if(SOSEX) write(*,*) 'SOSEX correction activated!'
write(*,*)
if(SOSEX) then
write(*,*) 'SOSEX correction activated but BUG!'
stop
end if
! COHSEX approximation
if(COHSEX) write(*,*) 'COHSEX approximation activated!'
write(*,*)
if(COHSEX) then
write(*,*) 'COHSEX approximation activated!'
write(*,*)
end if
! TDA for W
if(TDA_W) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
write(*,*)
if(TDA_W) then
write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
write(*,*)
end if
! TDA
if(TDA) write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
end if
! Memory allocation
allocate(eGW(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), &
Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), &
OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), &
error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Initialization
@ -149,36 +157,30 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! AO to MO transformation of two-electron integrals
call AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis)
call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO_basis,ERI_MO_basis)
! Compute linear response
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_basis, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
endif
! Compute correlation part of the self-energy
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rhox(:,:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY_RPA,rho_RPA)
if(G0W) then
call self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC)
call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z)
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)
else
call self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC)
call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z)
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
@ -221,7 +223,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! Print results
! call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA(ispin),EcGM,EqsGW)
call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA,EcGM,EqsGW)
! Increment
@ -258,23 +260,23 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
endif
! Dump RPA correlation energy
! Deallocate memory
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW total energy =',ENuc + EqsGW + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error,error_diis,F_diis)
! Perform BSE calculation
if(BSE) then
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO_basis,dipole_int, &
eGW,eGW,EcBSE)
if(exchange_kernel) then
EcBSE(1) = 0.5d0*EcBSE(1)
EcBSE(2) = 1.5d0*EcBSE(2)
end if
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
@ -301,8 +303,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcAC)
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

26
src/QuAcK/read_options.f90
View File

@ -1,6 +1,6 @@
subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
singlet_manifold,triplet_manifold,TDA, &
singlet,triplet,spin_conserved,spin_flip,TDA, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW, &
COHSEX,SOSEX,TDA_W,G0W,GW0, &
@ -26,8 +26,10 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
logical,intent(out) :: DIIS_CC
integer,intent(out) :: n_diis_CC
logical,intent(out) :: singlet_manifold
logical,intent(out) :: triplet_manifold
logical,intent(out) :: singlet
logical,intent(out) :: triplet
logical,intent(out) :: spin_conserved
logical,intent(out) :: spin_flip
logical,intent(out) :: TDA
integer,intent(out) :: maxSCF_GF
@ -113,16 +115,20 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
! Read excited state options
singlet_manifold = .false.
triplet_manifold = .false.
TDA = .false.
singlet = .false.
triplet = .false.
spin_conserved = .false.
spin_flip = .false.
TDA = .false.
read(1,*)
read(1,*) answer1,answer2,answer3
read(1,*) answer1,answer2,answer3,answer4,answer5
if(answer1 == 'T') singlet_manifold = .true.
if(answer2 == 'T') triplet_manifold = .true.
if(answer3 == 'T') TDA = .true.
if(answer1 == 'T') singlet = .true.
if(answer2 == 'T') triplet = .true.
if(answer3 == 'T') spin_conserved = .true.
if(answer4 == 'T') spin_flip = .true.
if(answer5 == 'T') TDA = .true.
! Read Green function options

48
src/QuAcK/renormalization_factor.f90
View File

@ -1,4 +1,4 @@
subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,Z)
subroutine renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z)
! Compute renormalization factor for GW
@ -8,13 +8,11 @@ subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r
! Input variables
logical,intent(in) :: COHSEX
logical,intent(in) :: SOSEX
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
double precision,intent(in) :: rhox(nBas,nBas,nS)
! Local variables
@ -35,42 +33,8 @@ subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r
Z(:) = 1d0
return
end if
! Occupied part of the correlation self-energy
do x=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
end do
end do
end do
! Virtual part of the correlation self-energy
do x=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
end do
end do
end do
! SOSEX correction
if(SOSEX) then
else
! Occupied part of the correlation self-energy
@ -81,7 +45,7 @@ subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r
do b=nO+1,nBas-nR
jb = jb + 1
eps = e(x) - e(i) + Omega(jb)
Z(x) = Z(x) - (rho(x,i,jb)/eps)*(rhox(x,i,jb)/eps)
Z(x) = Z(x) - 2d0*rho(x,i,jb)**2*(eps/(eps**2 + eta**2))**2
end do
end do
end do
@ -96,13 +60,13 @@ subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r
do b=nO+1,nBas-nR
jb = jb + 1
eps = e(x) - e(a) - Omega(jb)
Z(x) = Z(x) - (rho(x,a,jb)/eps)*(rhox(x,a,jb)/eps)
Z(x) = Z(x) - 2d0*rho(x,a,jb)**2*(eps/(eps**2 + eta**2))**2
end do
end do
end do
end do
endif
end if
! Compute renormalization factor from derivative of SigC

42
src/QuAcK/self_energy_correlation.f90
View File

@ -1,4 +1,4 @@
subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,EcGM,SigC)
subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC)
! Compute correlation part of the self-energy
@ -8,13 +8,11 @@ subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,
! Input variables
logical,intent(in) :: COHSEX
logical,intent(in) :: SOSEX
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
double precision,intent(in) :: rhox(nBas,nBas,nS)
! Local variables
@ -30,7 +28,9 @@ subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,
SigC = 0d0
! COHSEX static approximation
!-----------------------------!
! COHSEX static approximation !
!-----------------------------!
if(COHSEX) then
@ -65,6 +65,10 @@ subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,
else
!----------------!
! GW self-energy !
!----------------!
! Occupied part of the correlation self-energy
do x=nC+1,nBas-nR
@ -91,36 +95,6 @@ subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,
enddo
enddo
if(SOSEX) then
! SOSEX: occupied part of the correlation self-energy
do x=nC+1,nBas-nR
do y=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
eps = e(x) - e(i) + Omega(jb)
SigC(x,y) = SigC(x,y) - rho(x,i,jb)*rhox(y,i,jb)*eps/(eps**2 + eta**2)
enddo
enddo
enddo
enddo
! SOSEX: virtual part of the correlation self-energy
do x=nC+1,nBas-nR
do y=nC+1,nBas-nR
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(x) - e(a) - Omega(jb)
SigC(x,y) = SigC(x,y) - rho(x,a,jb)*rhox(y,a,jb)*eps/(eps**2 + eta**2)
enddo
enddo
enddo
enddo
endif
endif
end subroutine self_energy_correlation

47
src/QuAcK/self_energy_correlation_diag.f90
View File

@ -1,4 +1,4 @@
subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,EcGM,SigC)
subroutine self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC)
! Compute diagonal of the correlation part of the self-energy
@ -8,7 +8,6 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
! Input variables
logical,intent(in) :: COHSEX
logical,intent(in) :: SOSEX
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
@ -19,7 +18,6 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
double precision,intent(in) :: rhox(nBas,nBas,nS)
! Local variables
@ -66,46 +64,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
EcGM = 0d0
do i=nC+1,nO
EcGM = EcGM + 0.5d0*SigC(i)
end do
!-----------------------------
! SOSEX self-energy *BUG*
!-----------------------------
elseif(SOSEX) then
! SOSEX: 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)
SigC(p) = SigC(p) - rho(p,i,jb)*rhox(p,i,jb)*eps/(eps**2 + eta**2)
end do
end do
end do
! SOSEX: 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)
SigC(p) = SigC(p) - rho(p,a,jb)*rhox(p,a,jb)*eps/(eps**2 + eta**2)
end do
end do
end do
! GM correlation energy
do i=nC+1,nO
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(a) - e(i) + Omega(jb)
EcGM = EcGM + rho(a,i,jb)*rhox(a,i,jb)*eps/(eps**2 + eta**2)
end do
end do
EcGM = EcGM - SigC(i)
end do
!-----------------------------
@ -143,7 +102,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(a) - e(i) + Omega(jb)
EcGM = EcGM - 2d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2)
EcGM = EcGM - 4d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2)
end do
end do
end do

2
src/QuAcK/spatial_to_spin_MO_energy.f90
View File

@ -1,6 +1,6 @@
subroutine spatial_to_spin_MO_energy(nBas,e,nBas2,se)
! Convert ERIs from spatial to spin orbitals
! Convert MO energies from spatial to spin orbitals
implicit none

171
src/QuAcK/unrestricted_Bethe_Salpeter.f90 Normal file
View File

@ -0,0 +1,171 @@
subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, &
nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, &
dipole_int_aa,dipole_int_bb,eW,eGW,EcBSE)
! Compute the Bethe-Salpeter excitation energies
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: TDA_W
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: evDyn
logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip
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) :: nS(nspin)
double precision,intent(in) :: eW(nBas,nspin)
double precision,intent(in) :: eGW(nBas,nspin)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart)
! Local variables
integer :: ispin
integer :: isp_W
integer :: nS_aa,nS_bb,nS_sc
double precision :: EcRPA
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:,:)
double precision,allocatable :: OmBSE_sc(:)
double precision,allocatable :: XpY_BSE_sc(:,:)
double precision,allocatable :: XmY_BSE_sc(:,:)
integer :: nS_ab,nS_ba,nS_sf
double precision,allocatable :: OmBSE_sf(:)
double precision,allocatable :: XpY_BSE_sf(:,:)
double precision,allocatable :: XmY_BSE_sf(:,:)
! Output variables
double precision,intent(out) :: EcBSE(nspin)
! Memory allocation
nS_aa = nS(1)
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin))
!--------------------------!
! Spin-conserved screening !
!--------------------------!
isp_W = 1
EcRPA = 0d0
! Compute spin-conserved RPA screening
call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
eW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
!----------------------------!
! Spin-conserved excitations !
!----------------------------!
if(spin_conserved) then
ispin = 1
EcBSE(ispin) = 0d0
allocate(OmBSE_sc(nS_sc),XpY_BSE_sc(nS_sc,nS_sc),XmY_BSE_sc(nS_sc,nS_sc))
! Compute spin-conserved BSE excitation energies
call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcBSE(ispin), &
OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
call print_excitation('BSE@UG0W0',5,nS_sc,OmBSE_sc)
call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
!-------------------------------------------------
! Compute the dynamical screening at the BSE level
!-------------------------------------------------
! if(dBSE) then
! ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized)
!
! if(evDyn) then
!
! call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), &
! XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin))
! else
!
! call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), &
! XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin))
! end if
! end if
end if
!-----------------------!
! Spin-flip excitations !
!-----------------------!
if(spin_flip) then
ispin = 2
EcBSE(ispin) = 0d0
! Memory allocation
nS_ab = (nO(1) - nC(1))*(nV(2) - nR(2))
nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1))
nS_sf = nS_ab + nS_ba
allocate(OmBSE_sf(nS_sf),XpY_BSE_sf(nS_sf,nS_sf),XmY_BSE_sf(nS_sf,nS_sf))
call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcBSE(ispin), &
OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf)
call print_excitation('BSE@UG0W0',6,nS_sf,OmBSE_sf)
!-------------------------------------------------
! Compute the dynamical screening at the BSE level
!-------------------------------------------------
! if(dBSE) then
! ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized)
! if(evDyn) then
!
! call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), &
! XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin))
! else
!
! call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), &
! XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin))
! end if
! end if
end if
end subroutine unrestricted_Bethe_Salpeter

152
src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 Normal file
View File

@ -0,0 +1,152 @@
subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,A_lr)
! Compute the extra term for Bethe-Salpeter equation for linear response in the unrestricted formalism
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: ispin
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nSa
integer,intent(in) :: nSb
integer,intent(in) :: nSt
integer,intent(in) :: nS_sc
double precision,intent(in) :: eta
double precision,intent(in) :: lambda
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega(nS_sc)
double precision,intent(in) :: rho(nBas,nBas,nS_sc,nspin)
! Local variables
double precision :: chi
double precision :: eps
integer :: i,j,a,b,ia,jb,kc
! Output variables
double precision,intent(out) :: A_lr(nSt,nSt)
!--------------------------------------------------!
! Build BSE matrix for spin-conserving transitions !
!--------------------------------------------------!
if(ispin == 1) then
! aaaa block
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
chi = 0d0
do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps
enddo
A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aaaa(i,b,j,a) + 2d0*lambda*chi
enddo
enddo
enddo
enddo
! bbbb block
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
chi = 0d0
do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps
enddo
A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,b,j,a) + 2d0*lambda*chi
enddo
enddo
enddo
enddo
end if
!--------------------------------------------!
! Build BSE matrix for spin-flip transitions !
!--------------------------------------------!
if(ispin == 2) then
! abab block
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
chi = 0d0
do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,j,kc,1)*rho(a,b,kc,2)*Omega(kc)/eps
enddo
A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_abab(i,b,j,a) + 2d0*lambda*chi
end do
end do
end do
end do
! baba block
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
chi = 0d0
do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,j,kc,2)*rho(a,b,kc,1)*Omega(kc)/eps
enddo
A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_abab(b,i,a,j) + 2d0*lambda*chi
end do
end do
end do
end do
end if
end subroutine unrestricted_Bethe_Salpeter_A_matrix

154
src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 Normal file
View File

@ -0,0 +1,154 @@
subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,B_lr)
! Compute the extra term for Bethe-Salpeter equation for linear response
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: ispin
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nSa
integer,intent(in) :: nSb
integer,intent(in) :: nSt
integer,intent(in) :: nS_sc
double precision,intent(in) :: eta
double precision,intent(in) :: lambda
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega(nS_sc)
double precision,intent(in) :: rho(nBas,nBas,nS_sc,nspin)
! Local variables
double precision :: chi
double precision :: eps
integer :: i,j,a,b,ia,jb,kc
! Output variables
double precision,intent(out) :: B_lr(nSt,nSt)
!--------------------------------------------------!
! Build BSE matrix for spin-conserving transitions !
!--------------------------------------------------!
if(ispin == 1) then
! aaaa block
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
chi = 0d0
do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps
enddo
B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI_aaaa(i,j,b,a) + 2d0*lambda*chi
enddo
enddo
enddo
enddo
! bbbb block
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
chi = 0d0
do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps
enddo
B_lr(nSa+ia,nSa+jb) = B_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,j,b,a) + 2d0*lambda*chi
enddo
enddo
enddo
enddo
end if
!--------------------------------------------!
! Build BSE matrix for spin-flip transitions !
!--------------------------------------------!
if(ispin == 2) then
! abba block
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
chi = 0d0
do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,b,kc,1)*rho(a,j,kc,2)*Omega(kc)/eps
enddo
B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_abab(i,a,b,j) + 2d0*lambda*chi
end do
end do
end do
end do
! baab block
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
chi = 0d0
do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,b,kc,2)*rho(a,j,kc,1)*Omega(kc)/eps
enddo
B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_abab(b,j,i,a) + 2d0*lambda*chi
end do
end do
end do
end do
end if
end subroutine unrestricted_Bethe_Salpeter_B_matrix

83
src/QuAcK/unrestricted_QP_graph.f90 Normal file
View File

@ -0,0 +1,83 @@
subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
! Compute the graphical solution of the QP equation
implicit none
include 'parameters.h'
! Input variables
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) :: eta
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS,nspin)
double precision,intent(in) :: eGWlin(nBas)
! Local variables
integer :: p
integer :: nIt
integer,parameter :: maxIt = 10
double precision,parameter :: thresh = 1d-6
double precision,external :: USigmaC,dUSigmaC
double precision :: sig,dsig
double precision :: f,df
double precision :: w
! Output variables
double precision,intent(out) :: eGW(nBas)
! Run Newton's algorithm to find the root
do p=nC+1,nBas-nR
write(*,*) '-----------------'
write(*,'(A10,I3)') 'Orbital ',p
write(*,*) '-----------------'
w = eGWlin(p)
nIt = 0
f = 1d0
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f
do while (abs(f) > thresh .and. nIt < maxIt)
nIt = nIt + 1
sig = USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho)
dsig = dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho)
f = w - eHF(p) - sig
df = 1d0 - dsig
w = w - f/df
write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sig
end do
if(nIt == maxIt) then
write(*,*) 'Newton root search has not converged!'
eGW(p) = eGWlin(p)
else
eGW(p) = w
write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGW(p)*HaToeV
write(*,*)
end if
end do
end subroutine unrestricted_QP_graph

63
src/QuAcK/unrestricted_excitation_density.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ab,ERI_bb,XpY_a,XpY_b,rho)
subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho)
! Compute excitation densities for unrestricted reference
@ -14,11 +14,10 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_
integer,intent(in) :: nSa
integer,intent(in) :: nSb
integer,intent(in) :: nSt
double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: XpY_a(nSa,nSa)
double precision,intent(in) :: XpY_b(nSb,nSb)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: XpY(nSt,nSt)
! Local variables
@ -32,34 +31,34 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_
rho(:,:,:,:) = 0d0
!-------------
! alpha block
!-------------
!-------------!
! alpha block !
!-------------!
do p=nC(1)+1,nBas-nR(1)
do q=nC(1)+1,nBas-nR(1)
! Same-spin contribution
do ia=1,nSa
do ia=1,nSt
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aa(p,j,q,b)*XpY_a(ia,jb)
rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aaaa(p,j,q,b)*XpY(ia,jb)
enddo
enddo
enddo
! Opposite-spin contribution
do ia=1,nSb
jb = 0
do ia=1,nSt
jb = nSa
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
rho(p,q,nSa+ia,1) = rho(p,q,nSa+ia,1) + ERI_ab(p,j,q,b)*XpY_b(ia,jb)
rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aabb(p,j,q,b)*XpY(ia,jb)
enddo
enddo
@ -68,34 +67,34 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_
enddo
enddo
!------------
! Beta block
!------------
!------------!
! Beta block !
!------------!
do p=nC(2)+1,nBas-nR(2)
do q=nC(2)+1,nBas-nR(2)
! Same-spin contribution
do ia=1,nSb
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_bb(p,j,q,b)*XpY_b(ia,jb)
enddo
enddo
enddo
! Opposite-spin contribution
do ia=1,nSa
do ia=1,nSt
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
rho(p,q,nSb+ia,2) = rho(p,q,nSb+ia,2) + ERI_ab(j,p,b,q)*XpY_a(ia,jb)
rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_aabb(j,p,b,q)*XpY(ia,jb)
enddo
enddo
enddo
! Same-spin contribution
do ia=1,nSt
jb = nSa
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_bbbb(p,j,q,b)*XpY(ia,jb)
enddo
enddo

123
src/QuAcK/unrestricted_linear_response.f90 Normal file
View File

@ -0,0 +1,123 @@
subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, &
e,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,Omega,XpY,XmY)
! Compute linear response for unrestricted formalism
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: ispin
logical,intent(in) :: dRPA
logical,intent(in) :: TDA
logical,intent(in) :: BSE
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nSa
integer,intent(in) :: nSb
integer,intent(in) :: nSt
integer,intent(in) :: nS_sc
double precision,intent(in) :: lambda
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: OmRPA(nS_sc)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS_sc,nspin)
! Local variables
integer :: ia
double precision :: trace_matrix
double precision,allocatable :: A(:,:)
double precision,allocatable :: B(:,:)
double precision,allocatable :: ApB(:,:)
double precision,allocatable :: AmB(:,:)
double precision,allocatable :: AmBSq(:,:)
double precision,allocatable :: AmBIv(:,:)
double precision,allocatable :: Z(:,:)
! Output variables
double precision,intent(out) :: EcRPA
double precision,intent(out) :: Omega(nSt)
double precision,intent(out) :: XpY(nSt,nSt)
double precision,intent(out) :: XmY(nSt,nSt)
! Memory allocation
allocate(A(nSt,nSt),B(nSt,nSt),ApB(nSt,nSt),AmB(nSt,nSt),AmBSq(nSt,nSt),AmBIv(nSt,nSt),Z(nSt,nSt))
! Build A and B matrices
call unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A)
if(BSE) &
call unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,A)
! Tamm-Dancoff approximation
B = 0d0
if(.not. TDA) then
call unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,B)
if(BSE) &
call unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,B)
end if
! Build A + B and A - B matrices
ApB = A + B
AmB = A - B
! Diagonalize linear response matrix
call diagonalize_matrix(nSt,AmB,Omega)
if(minval(Omega) < 0d0) &
call print_warning('You may have instabilities in linear response: A-B is not positive definite!!')
do ia=1,nSt
if(Omega(ia) < 0d0) Omega(ia) = 0d0
end do
call ADAt(nSt,AmB,1d0*sqrt(Omega),AmBSq)
call ADAt(nSt,AmB,1d0/sqrt(Omega),AmBIv)
Z = matmul(AmBSq,matmul(ApB,AmBSq))
call diagonalize_matrix(nSt,Z,Omega)
if(minval(Omega) < 0d0) &
call print_warning('You may have instabilities in linear response: negative excitations!!')
do ia=1,nSt
if(Omega(ia) < 0d0) Omega(ia) = 0d0
end do
Omega = sqrt(Omega)
XpY = matmul(transpose(Z),AmBSq)
call DA(nSt,1d0/sqrt(Omega),XpY)
XmY = matmul(transpose(Z),AmBIv)
call DA(nSt,1d0*sqrt(Omega),XmY)
! Compute the RPA correlation energy
EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nSt,A))
end subroutine unrestricted_linear_response

175
src/QuAcK/unrestricted_linear_response_A_matrix.f90 Normal file
View File

@ -0,0 +1,175 @@
subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, &
e,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_lr)
! Compute linear response
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dRPA
integer,intent(in) :: ispin
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nSa
integer,intent(in) :: nSb
integer,intent(in) :: nSt
double precision,intent(in) :: lambda
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas)
! Local variables
double precision :: delta_dRPA
double precision,external :: Kronecker_delta
integer :: i,j,a,b,ia,jb
! Output variables
double precision,intent(out) :: A_lr(nSt,nSt)
! Direct RPA
delta_dRPA = 0d0
if(dRPA) delta_dRPA = 1d0
!-----------------------------------------------
! Build A matrix for spin-conserving transitions
!-----------------------------------------------
if(ispin == 1) then
! aaaa block
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
A_lr(ia,jb) = (e(a,1) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+ lambda*ERI_aaaa(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_aaaa(i,b,j,a)
end do
end do
end do
end do
! aabb block
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
A_lr(ia,nSa+jb) = lambda*ERI_aabb(i,b,a,j)
end do
end do
end do
end do
! bbaa block
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
A_lr(nSa+ia,jb) = lambda*ERI_aabb(b,i,j,a)
end do
end do
end do
end do
! bbbb block
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
A_lr(nSa+ia,nSa+jb) = (e(a,2) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+ lambda*ERI_bbbb(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_bbbb(i,b,j,a)
end do
end do
end do
end do
end if
!-----------------------------------------------
! Build A matrix for spin-flip transitions
!-----------------------------------------------
if(ispin == 2) then
A_lr(:,:) = 0d0
! abab block
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
A_lr(ia,jb) = (e(a,2) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
- (1d0 - delta_dRPA)*lambda*ERI_abab(i,b,j,a)
end do
end do
end do
end do
! baba block
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
A_lr(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
- (1d0 - delta_dRPA)*lambda*ERI_abab(b,j,i,a)
end do
end do
end do
end do
end if
end subroutine unrestricted_linear_response_A_matrix

170
src/QuAcK/unrestricted_linear_response_B_matrix.f90 Normal file
View File

@ -0,0 +1,170 @@
subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,B_lr)
! Compute linear response
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dRPA
integer,intent(in) :: ispin
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nSa
integer,intent(in) :: nSb
integer,intent(in) :: nSt
double precision,intent(in) :: lambda
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas)
! Local variables
double precision :: delta_dRPA
double precision,external :: Kronecker_delta
integer :: i,j,a,b,ia,jb
! Output variables
double precision,intent(out) :: B_lr(nSt,nSt)
! Direct RPA
delta_dRPA = 0d0
if(dRPA) delta_dRPA = 1d0
!-----------------------------------------------
! Build B matrix for spin-conserving transitions
!-----------------------------------------------
if(ispin == 1) then
! aaaa block
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
B_lr(ia,jb) = lambda*ERI_aaaa(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_aaaa(i,j,b,a)
end do
end do
end do
end do
! aabb block
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
B_lr(ia,nSa+jb) = lambda*ERI_aabb(i,j,a,b)
end do
end do
end do
end do
! bbaa block
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
B_lr(nSa+ia,jb) = lambda*ERI_aabb(j,i,b,a)
end do
end do
end do
end do
! bbbb block
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
B_lr(nSa+ia,nSa+jb) = lambda*ERI_bbbb(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_bbbb(i,j,b,a)
end do
end do
end do
end do
end if
!-----------------------------------------------
! Build B matrix for spin-flip transitions
!-----------------------------------------------
if(ispin == 2) then
B_lr(:,:) = 0d0
! abba block
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
B_lr(ia,nSa+jb) = - (1d0 - delta_dRPA)*lambda*ERI_abab(i,a,b,j)
end do
end do
end do
end do
! baab block
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
B_lr(nSa+ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI_abab(b,j,i,a)
end do
end do
end do
end do
end if
end subroutine unrestricted_linear_response_B_matrix

90
src/QuAcK/unrestricted_renormalization_factor.f90 Normal file
View File

@ -0,0 +1,90 @@
subroutine unrestricted_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,j,a,b,p,q,jb
double precision :: eps
! Output variables
double precision,intent(out) :: Z(nBas,nspin)
! Initialize
Z(:,:) = 0d0
!--------------!
! 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)
Z(p,1) = Z(p,1) + rho(p,i,jb,1)**2*(eps/(eps**2 + eta**2))**2
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)
Z(p,1) = Z(p,1) + rho(p,a,jb,1)**2*(eps/(eps**2 + eta**2))**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 i=nC(2)+1,nO(2)
do jb=1,nSt
eps = e(p,2) - e(i,2) + Omega(jb)
Z(p,2) = Z(p,2) + rho(p,i,jb,2)**2*(eps/(eps**2 + eta**2))**2
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)
Z(p,2) = Z(p,2) + rho(p,a,jb,2)**2*(eps/(eps**2 + eta**2))**2
end do
end do
end do
! Final rescaling
Z(:,:) = 1d0/(1d0 + Z(:,:))
end subroutine unrestricted_renormalization_factor

16
src/QuAcK/unrestricted_self_energy_correlation_diag.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,e,Omega,rho,SigC)
subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC)
! Compute diagonal of the correlation part of the self-energy
@ -13,8 +13,6 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nS
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nSa
integer,intent(in) :: nSb
integer,intent(in) :: nSt
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: Omega(nSt)
@ -33,9 +31,9 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nS
SigC(:,:) = 0d0
!--------------
! Spin-up part
!--------------
!--------------!
! Spin-up part !
!--------------!
! Occupied part of the correlation self-energy
@ -59,9 +57,9 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nS
end do
end do
!----------------
! Spin-down part
!----------------
!----------------!
! Spin-down part !
!----------------!
! Occupied part of the correlation self-energy

30
src/QuAcK/unrestricted_spatial_to_spin_MO_energy.f90 Normal file
View File

@ -0,0 +1,30 @@
subroutine unrestricted_spatial_to_spin_MO_energy(nBas,e,nBas2,se)
! Convert MO energies from unrestricted spatial to spin orbitals
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nBas2
double precision,intent(in) :: e(nBas,nspin)
! Local variables
integer :: p
! Output variables
double precision,intent(out) :: se(nBas2)
do p=1,nBas2,2
se(p) = e(p,1)
enddo
do p=2,nBas2,2
se(p) = e(p,2)
enddo
end subroutine unrestricted_spatial_to_spin_MO_energy

4
src/utils/read_integrals.f90
View File

@ -37,6 +37,10 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
open(unit=10,file='int/Nuc.dat')
open(unit=11,file='int/ERI.dat')
open(unit=21,file='int/x.dat')
open(unit=22,file='int/y.dat')
open(unit=23,file='int/z.dat')
! Read overlap integrals
S(:,:) = 0d0