PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-04 05:03:49 +01:00
Code Releases Activity

revive eDFT code

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-03-14 23:00:44 +01:00
parent 6b7ebf8158
commit 0be0d5d8c3
31 changed files with 243 additions and 291 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
GoXC
View File

@ -11,8 +11,8 @@ if [ $# = 2 ]
then
cp examples/molecule."$1" input/molecule
cp examples/basis."$1"."$2" input/basis
cp basis/"$2" input/basis.qcaml
cp examples/basis."$1"."$2" input/weight
./bin/IntPak
./bin/eDFT
fi

2
examples/molecule.F2
View File

@ -2,4 +2,4 @@
2 9 9 0 0
# Znuc x y z
F 0. 0. 0.
F 0. 0. 2.640
F 0. 0. 2.1

2
examples/molecule.LiH
View File

@ -2,4 +2,4 @@
2 2 2 0 0
# Znuc x y z
Li 0. 0. 0.
H 0. 0. 3.018
H 0. 0. 3.5

34
input/basis
View File

@ -1,18 +1,26 @@
1 3
1 8
S 6
1 1264.5857000 0.0019448
2 189.9368100 0.0148351
3 43.1590890 0.0720906
4 12.0986630 0.2371542
5 3.8063232 0.4691987
6 1.2728903 0.3565202
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 3.1964631 -0.1126487
2 0.7478133 -0.2295064
3 0.2199663 1.1869167
P 3
1 3.1964631 0.0559802
2 0.7478133 0.2615506
3 0.2199663 0.7939723
S 1
1 0.1220000 1.0000000
1 0.0823099 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
1 0.0823099 1.0000000
S 1
1 0.1220000 1.0000000
1 0.0207000 1.0000000
P 1
1 0.7270000 1.0000000
1 0.0207000 1.0000000
D 1
1 0.4000000 1.0000000

8
input/dft
View File

@ -1,12 +1,12 @@
# exchange rung: Hartree = 0, LDA = 1 (S51), GGA = 2(G96,B88), hybrid = 4, Hartree-Fock = 666
666 HF
1 S51
# correlation rung: Hartree = 0, LDA = 1 (W38,VWN5,C16,LF19), GGA = 2(LYP), hybrid = 4(B3LYP), Hartree-Fock = 666
0 HF
0 LF19
# quadrature grid SG-n
1
# Number of states in ensemble (nEns)
1
3
# Ensemble weights: wEns(1),...,wEns(nEns-1)
0.00000 0.00000
0.00000 0.50000
# eKS: maxSCF thresh DIIS n_diis guess_type ortho_type
64 0.0000001 T 15 1 1

10
input/methods
View File

@ -1,16 +1,16 @@
# RHF UHF MOM
T F F
# MP2 MP3 MP2-F12
T F F
F F F
# CCD CCSD CCSD(T) ringCCD ladderCCD
F F F F F
# CIS RPA RPAx ppRPA ADC
T T T F F
F F F F F
# GF2 GF3
F F
T F
# G0W0 evGW qsGW
T F F
# G0T0 evGT qsGT
F F F
# G0T0 evGT qsGT
T F F
# MCMP2
F

5
input/molecule
View File

@ -1,5 +1,4 @@
# nAt nEla nElb nCore nRyd
2 1 1 0 0
1 2 2 0 0
# Znuc x y z
H 0. 0. 0.
H 0. 0. 2.3
Be 0.0 0.0 0.0

5
input/molecule.xyz
View File

@ -1,4 +1,3 @@
2
1
H 0.0000000000 0.0000000000 0.0000000000
H 0.0000000000 0.0000000000 1.2171076727
Be 0.0000000000 0.0000000000 0.0000000000

2
input/options
View File

@ -9,7 +9,7 @@
# GF: maxSCF thresh DIIS n_diis renormalization
256 0.00001 T 5 3
# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta
256 0.00001 T 5 F F T F F F F 0.000
1 0.00001 T 5 F F T F F F T 0.000
# ACFDT: AC Kx XBS
T F T
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift

34
input/weight
View File

@ -1,18 +1,26 @@
1 3
1 8
S 6
1 1264.5857000 0.0019448
2 189.9368100 0.0148351
3 43.1590890 0.0720906
4 12.0986630 0.2371542
5 3.8063232 0.4691987
6 1.2728903 0.3565202
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 3.1964631 -0.1126487
2 0.7478133 -0.2295064
3 0.2199663 1.1869167
P 3
1 3.1964631 0.0559802
2 0.7478133 0.2615506
3 0.2199663 0.7939723
S 1
1 0.1220000 1.0000000
1 0.0823099 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
1 0.0823099 1.0000000
S 1
1 0.1220000 1.0000000
1 0.0207000 1.0000000
P 1
1 0.7270000 1.0000000
1 0.0207000 1.0000000
D 1
1 0.4000000 1.0000000

2
scan_H2.sh
View File

@ -1,7 +1,7 @@
#! /bin/bash
MOL="H2"
BASIS="cc-pvdz"
BASIS="cc-pvqz"
R_START=1.0
R_END=2.4
DR=0.1

4
scan_LiH.sh
View File

@ -2,9 +2,9 @@
MOL="LiH"
BASIS="cc-pvqz"
R_START=2.5
R_START=3.5
R_END=3.5
DR=0.001
DR=0.1
for R in $(seq $R_START $DR $R_END)
do

55
src/QuAcK/G0T0.f90
View File

@ -1,6 +1,6 @@
subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
! Perform G0W0 calculation with a T-matrix self-energy (G0T0)
! Perform one-shot calculation with a T-matrix self-energy (G0T0)
implicit none
include 'parameters.h'
@ -57,10 +57,10 @@ subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), &
Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), &
rho1s(nBas,nO-nC,nVVs),rho2s(nBas,nV-nR,nOOs), &
rho1s(nBas,nO,nVVs),rho2s(nBas,nV,nOOs), &
Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), &
Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
rho1t(nBas,nO-nC,nVVt),rho2t(nBas,nV-nR,nOOt), &
rho1t(nBas,nO,nVVt),rho2t(nBas,nV,nOOt), &
SigT(nBas),Z(nBas),eG0T0(nBas))
!----------------------------------------------
@ -137,53 +137,4 @@ subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
call print_G0T0(nBas,nO,eHF(:),ENuc,ERHF,SigT(:),Z(:),eG0T0(:),EcRPA(:))
! Perform BSE calculation
! if(BSE) then
! ! Singlet manifold
! if(singlet_manifold) then
! ispin = 1
! EcBSE(ispin) = 0d0
! call linear_response(ispin,.false.,.false.,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI, &
! rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
! call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
! call linear_response(ispin,.false.,.false.,BSE,nBas,nC,nO,nV,nR,nS,eG0T0,ERI, &
! rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
! call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
! endif
! ! Triplet manifold
! if(triplet_manifold) then
! ispin = 2
! EcBSE(ispin) = 0d0
! call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI, &
! rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
! call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
! call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0T0,ERI, &
! rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
! call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
! endif
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A40,F15.6)') 'BSE@G0T0 correlation energy (singlet) =',EcBSE(1)
! write(*,'(2X,A40,F15.6)') 'BSE@G0T0 correlation energy (triplet) =',EcBSE(2)
! write(*,'(2X,A40,F15.6)') 'BSE@G0T0 correlation energy =',EcBSE(1) + EcBSE(2)
! write(*,'(2X,A40,F15.6)') 'BSE@G0T0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
! endif
end subroutine G0T0

16
src/QuAcK/G0W0.f90
View File

@ -1,4 +1,5 @@
subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, &
singlet_manifold,triplet_manifold,linearize,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eGW)
! Perform G0W0 calculation
@ -18,6 +19,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif
logical,intent(in) :: TDA
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: linearize
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
@ -101,15 +103,21 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif
eGWlin(:) = eHF(:) + Z(:)*SigC(:)
if(linearize) then
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)
! Find graphical solution of the QP equation
else
! call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
! Find graphical solution of the QP equation
eGW(:) = eGWlin(:)
call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
end if
! Dump results

13
src/QuAcK/QuAcK.f90
View File

@ -155,13 +155,13 @@ program QuAcK
! = nO*nV
call read_molecule(nNuc,nEl(:),nO(:),nC(:),nR(:))
allocate(ZNuc(nNuc),rNuc(nNuc,3))
allocate(ZNuc(nNuc),rNuc(nNuc,ncart))
! Read geometry
call read_geometry(nNuc,ZNuc,rNuc,ENuc)
allocate(CenterShell(maxShell,3),TotAngMomShell(maxShell),KShell(maxShell), &
allocate(CenterShell(maxShell,ncart),TotAngMomShell(maxShell),KShell(maxShell), &
DShell(maxShell,maxK),ExpShell(maxShell,maxK))
!------------------------------------------------------------------------
@ -547,8 +547,9 @@ program QuAcK
if(doG0W0) then
call cpu_time(start_G0W0)
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0W0)
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, &
singlet_manifold,triplet_manifold,linearize,eta, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0W0)
call cpu_time(end_G0W0)
t_G0W0 = end_G0W0 - start_G0W0
@ -602,8 +603,8 @@ program QuAcK
if(doG0T0) then
call cpu_time(start_G0T0)
! call G0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
call G0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
! call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
call cpu_time(end_G0T0)
t_G0T0 = end_G0T0 - start_G0T0

13
src/QuAcK/evGW.f90
View File

@ -1,5 +1,5 @@
subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW0, &
singlet_manifold,triplet_manifold,linearize,eta, &
singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0W0)
! Perform self-consistent eigenvalue-only GW calculation
@ -25,7 +25,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
logical,intent(in) :: GW0
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: linearize
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eHF(nBas)
@ -144,18 +143,10 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
endif
! Solve the quasi-particle equation (linearized or not)
if(linearize) then
eGW(:) = eHF(:) + Z(:)*SigC(:)
else
! Solve the quasi-particle equation
eGW(:) = eHF(:) + SigC(:)
endif
! Convergence criteria
Conv = maxval(abs(eGW - eOld))

30
src/QuAcK/excitation_density_Tmatrix.f90
View File

@ -24,8 +24,8 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
! Output variables
double precision,intent(out) :: rho1(nBas,nBas,nVV)
double precision,intent(out) :: rho2(nBas,nBas,nOO)
double precision,intent(out) :: rho1(nBas,nO,nVV)
double precision,intent(out) :: rho2(nBas,nV,nOO)
rho1(:,:,:) = 0d0
rho2(:,:,:) = 0d0
@ -66,7 +66,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
end do
end do
do a=nO+1,nBas-nR
do a=1,nV-nR
do ij=1,nOO
cd = 0
@ -74,9 +74,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
do d=nO+1,c
cd = cd + 1
rho2(p,a,ij) = rho2(p,a,ij) &
! + ERI(p,a,c,d)*X2(cd,ij)
+ (ERI(p,a,c,d) + ERI(p,a,d,c))*X2(cd,ij) &
/sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(c,d)))
! + ERI(p,nO+a,c,d)*X2(cd,ij)
+ (ERI(p,nO+a,c,d) + ERI(p,nO+a,d,c))*X2(cd,ij) &
/sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(c,d)))
end do
end do
@ -86,9 +86,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
do l=nC+1,k
kl = kl + 1
rho2(p,a,ij) = rho2(p,a,ij) &
! + ERI(p,a,k,l)*Y2(kl,ij)
+ (ERI(p,a,k,l) + ERI(p,a,l,k))*Y2(kl,ij) &
/sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(k,l)))
! + ERI(p,nO+a,k,l)*Y2(kl,ij)
+ (ERI(p,nO+a,k,l) + ERI(p,nO+a,l,k))*Y2(kl,ij) &
/sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(k,l)))
end do
end do
@ -131,7 +131,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
end do
end do
do a=nO+1,nBas-nR
do a=1,nV-nR
do ij=1,nOO
cd = 0
@ -139,7 +139,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
do d=nO+1,c-1
cd = cd + 1
rho2(p,a,ij) = rho2(p,a,ij) &
+ (ERI(p,a,c,d) - ERI(p,a,d,c))*X2(cd,ij)
+ (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij)
end do
end do
@ -148,7 +148,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
do l=nC+1,k-1
kl = kl + 1
rho2(p,a,ij) = rho2(p,a,ij) &
+ (ERI(p,a,k,l) - ERI(p,a,l,k))*Y2(kl,ij)
+ (ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij)
end do
end do
@ -190,7 +190,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
end do
end do
do a=nO+1,nBas-nR
do a=1,nV-nR
do ij=1,nOO
cd = 0
@ -198,7 +198,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
do d=c+1,nBas-nR
cd = cd + 1
rho2(p,a,ij) = rho2(p,a,ij) &
+ (ERI(p,a,c,d) - ERI(p,a,d,c))*X2(cd,ij)
+ (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij)
end do
end do
@ -207,7 +207,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
do l=k+1,nO
kl = kl + 1
rho2(p,a,ij) = rho2(p,a,ij) &
+ (ERI(p,a,k,l) - ERI(p,a,l,k))*Y2(kl,ij)
+ (ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij)
end do
end do

11
src/QuAcK/linear_response_pp.f90
View File

@ -1,6 +1,7 @@
subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,EcppRPA)
subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV, &
e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,EcppRPA)
! Compute the p-p channel of the linear response: see Scueria et al. JCP 139, 104113 (2013)
! Compute the p-p channel of the linear response: see Scuseria et al. JCP 139, 104113 (2013)
implicit none
include 'parameters.h'
@ -16,6 +17,7 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
! Local variables
logical :: dump_matrices = .false.
integer :: ab,cd,ij,kl
integer :: p,q,r,s
double precision :: trace_matrix
@ -65,6 +67,10 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
M( 1:nVV ,nVV+1:nOO+nVV) = - B(1:nVV,1:nOO)
M(nVV+1:nOO+nVV, 1:nVV) = + transpose(B(1:nVV,1:nOO))
! Dump ppRPA matrices
if(dump_matrices) then
open(unit=42,file='B.dat')
open(unit=43,file='C.dat')
open(unit=44,file='D.dat')
@ -106,6 +112,7 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
close(45)
close(46)
end if
! print*, 'pp-RPA matrix'
! call matout(nOO+nVV,nOO+nVV,M(:,:))

36
src/QuAcK/renormalization_factor_Tmatrix.f90
View File

@ -15,9 +15,9 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
integer,intent(in) :: nVVs,nVVt
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1s(nVVs),Omega1t(nVVt)
double precision,intent(in) :: rho1s(nBas,nBas,nVVs),rho1t(nBas,nBas,nVVt)
double precision,intent(in) :: rho1s(nBas,nO,nVVs),rho1t(nBas,nO,nVVt)
double precision,intent(in) :: Omega2s(nOOs),Omega2t(nOOt)
double precision,intent(in) :: rho2s(nBas,nBas,nOOs),rho2t(nBas,nBas,nOOt)
double precision,intent(in) :: rho2s(nBas,nV,nOOs),rho2t(nBas,nV,nOOt)
! Local variables
@ -40,31 +40,23 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
do p=nC+1,nBas-nR
do i=nC+1,nO
cd = 0
do c=nO+1,nBas-nR
do d=nO+1,c
cd = cd + 1
do cd=1,nVVs
eps = e(p) + e(i) - Omega1s(cd)
Z(p) = Z(p) - 2d0*rho1s(p,i,cd)**2/eps**2
enddo
enddo
enddo
enddo
! Virtual part of the T-matrix self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
kl = 0
do k=nC+1,nO
do l=nC+1,k
kl = kl + 1
eps = e(p) + e(a) - Omega2s(kl)
do a=1,nV-nR
do kl=1,nOOs
eps = e(p) + e(nO+a) - Omega2s(kl)
Z(p) = Z(p) - 2d0*rho2s(p,a,kl)**2/eps**2
enddo
enddo
enddo
enddo
!----------------------------------------------
! Triplet part of the T-matrix self-energy
@ -74,31 +66,23 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
do p=nC+1,nBas-nR
do i=nC+1,nO
cd = 0
do c=nO+1,nBas-nR
do d=nO+1,c-1
cd = cd + 1
do cd=1,nVVt
eps = e(p) + e(i) - Omega1t(cd)
Z(p) = Z(p) - 2d0*rho1t(p,i,cd)**2/eps**2
enddo
enddo
enddo
enddo
! Virtual part of the T-matrix self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
kl = 0
do k=nC+1,nO
do l=nC+1,k-1
kl = kl + 1
eps = e(p) + e(a) - Omega2t(kl)
do a=1,nV-nR
do kl=1,nOOt
eps = e(p) + e(nO+a) - Omega2t(kl)
Z(p) = Z(p) - 2d0*rho2t(p,a,kl)**2/eps**2
enddo
enddo
enddo
enddo
! Compute renormalization factor from derivative of SigT

28
src/QuAcK/renormalization_factor_Tmatrix_so.f90
View File

@ -13,9 +13,9 @@ subroutine renormalization_factor_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omeg
integer,intent(in) :: nVV
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: rho1(nBas,nO,nVV)
double precision,intent(in) :: Omega2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO)
double precision,intent(in) :: rho2(nBas,nV,nOO)
! Local variables
@ -38,13 +38,13 @@ subroutine renormalization_factor_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omeg
do p=nC+1,nBas-nR
do i=nC+1,nO
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
do cd=1,nVV
! do c=nO+1,nBas-nR
! do d=c+1,nBas-nR
! cd = cd + 1
eps = e(p) + e(i) - Omega1(cd)
Z(p) = Z(p) - rho1(p,i,cd)**2/eps**2
enddo
! enddo
enddo
enddo
enddo
@ -52,14 +52,14 @@ subroutine renormalization_factor_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omeg
! Virtual part of the T-matrix self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
eps = e(p) + e(a) - Omega2(kl)
do a=1,nV-nR
do kl=1,nOO
! do k=nC+1,nO
! do l=k+1,nO
! kl = kl + 1
eps = e(p) + e(nO+a) - Omega2(kl)
Z(p) = Z(p) - rho2(p,a,kl)**2/eps**2
enddo
! enddo
enddo
enddo
enddo

4
src/QuAcK/self_energy_Tmatrix.f90
View File

@ -17,9 +17,9 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,Omega2,rho1
integer,intent(in) :: nVV
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: rho1(nBas,nO,nVV)
double precision,intent(in) :: Omega2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO)
double precision,intent(in) :: rho2(nBas,nV,nOO)
! Local variables

38
src/QuAcK/self_energy_Tmatrix_diag.f90
View File

@ -19,9 +19,9 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
integer,intent(in) :: nVVs,nVVt
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1s(nVVs),Omega1t(nVVt)
double precision,intent(in) :: rho1s(nBas,nBas,nVVs),rho1t(nBas,nBas,nVVt)
double precision,intent(in) :: rho1s(nBas,nO,nVVs),rho1t(nBas,nO,nVVt)
double precision,intent(in) :: Omega2s(nOOs),Omega2t(nOOt)
double precision,intent(in) :: rho2s(nBas,nBas,nOOs),rho2t(nBas,nBas,nOOt)
double precision,intent(in) :: rho2s(nBas,nV,nOOs),rho2t(nBas,nV,nOOt)
! Local variables
@ -44,31 +44,23 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
do p=nC+1,nBas-nR
do i=nC+1,nO
cd = 0
do c=nO+1,nBas-nR
do d=nO+1,c
cd = cd + 1
do cd=1,nVVs
eps = e(p) + e(i) - Omega1s(cd)
SigT(p) = SigT(p) + 2d0*rho1s(p,i,cd)**2/eps
enddo
enddo
enddo
enddo
! Virtual part of the T-matrix self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
kl = 0
do k=nC+1,nO
do l=nC+1,k
kl = kl + 1
eps = e(p) + e(a) - Omega2s(kl)
do a=1,nV-nR
do kl=1,nOOs
eps = e(p) + e(nO+a) - Omega2s(kl)
SigT(p) = SigT(p) + 2d0*rho2s(p,a,kl)**2/eps
enddo
enddo
enddo
enddo
!----------------------------------------------
! Triplet part of the T-matrix self-energy
@ -78,30 +70,22 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
do p=nC+1,nBas-nR
do i=nC+1,nO
cd = 0
do c=nO+1,nBas-nR
do d=nO+1,c-1
cd = cd + 1
do cd=1,nVVt
eps = e(p) + e(i) - Omega1t(cd)
SigT(p) = SigT(p) + 2d0*rho1t(p,i,cd)**2/eps
enddo
enddo
enddo
enddo
! Virtual part of the T-matrix self-energy
! Virtual part of the T-matrix self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
kl = 0
do k=nC+1,nO
do l=nC+1,k-1
kl = kl + 1
eps = e(p) + e(a) - Omega2t(kl)
do a=1,nV-nR
do kl=1,nOOt
eps = e(p) + e(nO+a) - Omega2t(kl)
SigT(p) = SigT(p) + 2d0*rho2t(p,a,kl)**2/eps
enddo
enddo
enddo
enddo
end subroutine self_energy_Tmatrix_diag

28
src/QuAcK/self_energy_Tmatrix_diag_so.f90
View File

@ -17,9 +17,9 @@ subroutine self_energy_Tmatrix_diag_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho
integer,intent(in) :: nVV
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1(nVV)
double precision,intent(in) :: rho1(nBas,nBas,nVV)
double precision,intent(in) :: rho1(nBas,nO,nVV)
double precision,intent(in) :: Omega2(nOO)
double precision,intent(in) :: rho2(nBas,nBas,nOO)
double precision,intent(in) :: rho2(nBas,nV,nOO)
! Local variables
@ -42,13 +42,13 @@ subroutine self_energy_Tmatrix_diag_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho
do p=nC+1,nBas-nR
do i=nC+1,nO
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
do cd=1,nVV
! do c=nO+1,nBas-nR
! do d=c+1,nBas-nR
! cd = cd + 1
eps = e(p) + e(i) - Omega1(cd)
SigT(p) = SigT(p) + rho1(p,i,cd)**2/eps
enddo
! enddo
enddo
enddo
enddo
@ -56,14 +56,14 @@ subroutine self_energy_Tmatrix_diag_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho
! Virtual part of the T-matrix self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
eps = e(p) + e(a) - Omega2(kl)
do a=1,nV-nR
do kl=1,nOO
! do k=nC+1,nO
! do l=k+1,nO
! kl = kl + 1
eps = e(p) + e(nO+a) - Omega2(kl)
SigT(p) = SigT(p) + rho2(p,a,kl)**2/eps
enddo
! enddo
enddo
enddo
enddo

2
src/QuAcK/soG0T0.f90
View File

@ -68,7 +68,7 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
allocate(Omega1(nVV),X1(nVV,nVV),Y1(nOO,nVV), &
Omega2(nOO),X2(nVV,nOO),Y2(nOO,nOO), &
rho1(nBas2,nBas2,nVV),rho2(nBas2,nBas2,nOO), &
rho1(nBas2,nO2,nVV),rho2(nBas2,nV2,nOO), &
eG0T0(nBas2),SigT(nBas2),Z(nBas2))
!----------------------------------------------

9
src/eDFT/Kohn_Sham.f90
View File

@ -1,5 +1,5 @@
subroutine Kohn_Sham(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thresh,max_diis, &
guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew)
subroutine Kohn_Sham(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thresh,max_diis,guess_type, &
nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew)
! Perform unrestricted Kohn-Sham calculation for ensembles
@ -224,7 +224,7 @@ subroutine Kohn_Sham(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thr
do ispin=1,nspin
do iEns=1,nEns
call gradient_density(nGrid,nBas,P(:,:,ispin,iEns),AO(:,:),dAO(:,:,:),drho(:,:,ispin,iEns))
! call gradient_density(nGrid,nBas,P(:,:,ispin,iEns),AO(:,:),dAO(:,:,:),drho(:,:,ispin,iEns))
end do
end do
@ -271,7 +271,8 @@ subroutine Kohn_Sham(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thr
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))
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

10
src/eDFT/Makefile
View File

@ -14,19 +14,19 @@ ifeq ($(PROFIL),1)
FC += -p -fno-inline
endif
LIBS = ~/Dropbox/quack/lib/*.a
#LIBS = -lblas -llapack
LIBS = ../../lib/*.a
LIBS += -lblas -llapack
SRCF90 = $(wildcard *.f90)
SRC = $(wildcard *.f)
SRC = $(wildcard *.F)
OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.f,$(ODIR)/%.o,$(SRC))
OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.F,$(ODIR)/%.o,$(SRC))
$(ODIR)/%.o: %.f90
$(FC) -c -o $@ $< $(FFLAGS)
$(ODIR)/%.o: %.f
$(ODIR)/%.o: %.F
$(FC) -c -o $@ $< $(FFLAGS)
$(BDIR)/eDFT: $(OBJ)

13
src/eDFT/eDFT.f90
View File

@ -31,6 +31,7 @@ program eDFT
double precision,allocatable :: dAO(:,:,:)
double precision :: start_KS,end_KS,t_KS
double precision :: start_int,end_int,t_int
integer :: nEns
double precision,allocatable :: wEns(:)
@ -96,8 +97,18 @@ program eDFT
! Read integrals
call cpu_time(start_int)
call system('./GoQCaml')
call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI)
call cpu_time(end_int)
t_int = end_int - start_int
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for reading integrals = ',t_int,' seconds'
write(*,*)
! Orthogonalization X = S^(-1/2)
call orthogonalization_matrix(ortho_type,nBas,S,X)
@ -122,7 +133,7 @@ program eDFT
!------------------------------------------------------------------------
call cpu_time(start_KS)
call Kohn_Sham(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(1:nEns),nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, &
call Kohn_Sham(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, &
nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,EKS)
call cpu_time(end_KS)

2
src/eDFT/exchange_energy.f90
View File

@ -17,7 +17,7 @@ subroutine exchange_energy(rung,DFA,nEns,wEns,nGrid,weight,nBas,P,FxHF,rho,drho,
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: FxHF(nBas,nBas)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(3,nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
! Local variables

2
src/eDFT/fock_exchange_energy.f90
View File

@ -20,6 +20,6 @@ subroutine fock_exchange_energy(nBas,P,Fx,Ex)
! Compute HF exchange energy
Ex = trace_matrix(nBas,matmul(P,Fx))
Ex = 0.5d0*trace_matrix(nBas,matmul(P,Fx))
end subroutine fock_exchange_energy

2
src/eDFT/fock_exchange_potential.f90
View File

@ -25,7 +25,7 @@ subroutine fock_exchange_potential(nBas,P,ERI,Fx)
do si=1,nBas
do la=1,nBas
do mu=1,nBas
Fx(mu,nu) = Fx(mu,nu) - P(la,si)*ERI(mu,si,la,nu)
Fx(mu,nu) = Fx(mu,nu) - P(la,si)*ERI(mu,la,si,nu)
enddo
enddo
enddo

2
src/eDFT/hartree_coulomb.f90
View File

@ -23,7 +23,7 @@ subroutine hartree_coulomb(nBas,P,ERI,J)
do nu=1,nBas
do la=1,nBas
do si=1,nBas
J(mu,nu) = J(mu,nu) + P(la,si)*ERI(mu,nu,la,si)
J(mu,nu) = J(mu,nu) + P(la,si)*ERI(mu,la,nu,si)
enddo
enddo
enddo