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

moved Cx in parameters

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-03-18 16:06:29 +01:00
parent bfd5131dae
commit d0020765fc
19 changed files with 124 additions and 134 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/molecule.H2
View File

@ -2,4 +2,4 @@
2 1 1 0 0
# Znuc x y z
H 0. 0. 0.
H 0. 0. 2.3
H 0. 0. 1.399

4
include/parameters.h
View File

@ -15,3 +15,7 @@
double precision,parameter :: pmtoau = 0.0188973d0
double precision,parameter :: BoToAn = 0.529177249d0
double precision,parameter :: CxLDA = - (3d0/4d0)*(3d0/pi)**(1d0/3d0)
double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0)
double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0)

50
input/basis
View File

@ -1,26 +1,30 @@
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 3.1964631 -0.1126487
2 0.7478133 -0.2295064
3 0.2199663 1.1869167
1 6
S 8
1 2940.0000000 0.0006800
2 441.2000000 0.0052360
3 100.5000000 0.0266060
4 28.4300000 0.0999930
5 9.1690000 0.2697020
6 3.1960000 0.4514690
7 1.1590000 0.2950740
8 0.1811000 0.0125870
S 8
1 2940.0000000 -0.0001230
2 441.2000000 -0.0009660
3 100.5000000 -0.0048310
4 28.4300000 -0.0193140
5 9.1690000 -0.0532800
6 3.1960000 -0.1207230
7 1.1590000 -0.1334350
8 0.1811000 0.5307670
S 1
1 0.0589000 1.0000000
P 3
1 3.1964631 0.0559802
2 0.7478133 0.2615506
3 0.2199663 0.7939723
S 1
1 0.0823099 1.0000000
1 3.6190000 0.0291110
2 0.7110000 0.1693650
3 0.1951000 0.5134580
P 1
1 0.0823099 1.0000000
S 1
1 0.0207000 1.0000000
P 1
1 0.0207000 1.0000000
1 0.0601800 1.0000000
D 1
1 0.4000000 1.0000000
1 0.2380000 1.0000000

12
input/dft
View File

@ -2,16 +2,16 @@
LIM-RKS
# exchange rung:
# Hartree = 0
# LDA = 1: RS51,S51,RMFL20
# GGA = 2: G96,B88
# LDA = 1: RS51,RMFL20
# GGA = 2:
# Hybrid = 4
# Hartree-Fock = 666
1 RS51
# correlation rung:
# Hartree = 0
# LDA = 1: W38,VWN5,C16,RMFL20
# GGA = 2: LYP
# Hybrid = 4: B3LYP
# LDA = 1: RVWN5,RMFL20
# GGA = 2:
# Hybrid = 4:
# Hartree-Fock = 666
1 RVWN5
# quadrature grid SG-n
@ -21,4 +21,4 @@
# Ensemble weights: wEns(1),...,wEns(nEns-1)
0.50000 0.00000
# GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type
64 0.0000001 T 5 1 1
32 0.0000001 T 5 1 1

6
input/methods
View File

@ -7,10 +7,10 @@
# CIS RPA RPAx ppRPA ADC
F F F F F
# GF2 GF3
T F
F F
# G0W0 evGW qsGW
F F F
# G0T0 evGT qsGT
T F F
# G0T0 evGT qsGT
F F F
# MCMP2
F

4
input/options
View File

@ -9,8 +9,8 @@
# 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
256 0.00001 T 5 F F T F F F T 0.000
# ACFDT: AC Kx XBS
T F F
T F T
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

50
input/weight
View File

@ -1,26 +1,30 @@
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 3.1964631 -0.1126487
2 0.7478133 -0.2295064
3 0.2199663 1.1869167
1 6
S 8
1 2940.0000000 0.0006800
2 441.2000000 0.0052360
3 100.5000000 0.0266060
4 28.4300000 0.0999930
5 9.1690000 0.2697020
6 3.1960000 0.4514690
7 1.1590000 0.2950740
8 0.1811000 0.0125870
S 8
1 2940.0000000 -0.0001230
2 441.2000000 -0.0009660
3 100.5000000 -0.0048310
4 28.4300000 -0.0193140
5 9.1690000 -0.0532800
6 3.1960000 -0.1207230
7 1.1590000 -0.1334350
8 0.1811000 0.5307670
S 1
1 0.0589000 1.0000000
P 3
1 3.1964631 0.0559802
2 0.7478133 0.2615506
3 0.2199663 0.7939723
S 1
1 0.0823099 1.0000000
1 3.6190000 0.0291110
2 0.7110000 0.1693650
3 0.1951000 0.5134580
P 1
1 0.0823099 1.0000000
S 1
1 0.0207000 1.0000000
P 1
1 0.0207000 1.0000000
1 0.0601800 1.0000000
D 1
1 0.4000000 1.0000000
1 0.2380000 1.0000000

30
src/eDFT/GOK_RKS.f90
View File

@ -1,5 +1,5 @@
subroutine GOK_RKS(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,EwGIC)
subroutine GOK_RKS(restart,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,EwGIC,F)
! Perform restricted Kohn-Sham calculation for ensembles
@ -8,6 +8,7 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
! Input variables
logical,intent(in) :: restart
integer,intent(in) :: x_rung,c_rung
character(len=12),intent(in) :: x_DFA,c_DFA
integer,intent(in) :: nEns
@ -29,6 +30,8 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ENuc
double precision,intent(inout):: F(nBas,nBas)
! Local variables
integer :: xc_rung
@ -46,7 +49,6 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
double precision,allocatable :: c(:,:)
double precision,allocatable :: cp(:,:)
double precision,allocatable :: J(:,:)
double precision,allocatable :: F(:,:)
double precision,allocatable :: Fp(:,:)
double precision,allocatable :: Fx(:,:)
double precision,allocatable :: FxHF(:,:)
@ -118,22 +120,24 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
! Memory allocation
allocate(eps(nBas),c(nBas,nBas),cp(nBas,nBas), &
J(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
Fx(nBas,nBas),FxHF(nBas,nBas),Fc(nBas,nBas),err(nBas,nBas), &
J(nBas,nBas),Fp(nBas,nBas),Fx(nBas,nBas), &
FxHF(nBas,nBas),Fc(nBas,nBas),err(nBas,nBas), &
Pw(nBas,nBas),rhow(nGrid),drhow(ncart,nGrid), &
err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis), &
P(nBas,nBas,nEns),rho(nGrid,nEns),drho(ncart,nGrid,nEns))
! Guess coefficients and eigenvalues
if(guess_type == 1) then
F(:,:) = Hc(:,:)
else if(guess_type == 2) then
call random_number(F(:,:))
if(.not. restart) then
if(guess_type == 1) then
F(:,:) = Hc(:,:)
else if(guess_type == 2) then
call random_number(F(:,:))
end if
end if
! Initialization

12
src/eDFT/LIM_RKS.f90
View File

@ -1,5 +1,5 @@
subroutine LIM_RKS(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)
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,F)
! Perform restricted Kohn-Sham calculation for ensembles
@ -29,6 +29,8 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ENuc
double precision,intent(out) :: F(nBas,nBas)
! Local variables
integer :: iEns
@ -63,8 +65,8 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
write(*,'(A40)') '*************************************************'
write(*,*)
call GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wLIM,nGrid,weight,maxSCF,thresh, &
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwZW,EwGICZW)
call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,nEns,wLIM,nGrid,weight,maxSCF,thresh, &
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwZW,EwGICZW,F)
!------------------------------------------------------------------------
! Equiensemble calculation
@ -82,8 +84,8 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
write(*,'(A40)') '*************************************************'
write(*,*)
call GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wLIM,nGrid,weight,maxSCF,thresh, &
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwEW,EwGICEW)
call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,nEns,wLIM,nGrid,weight,maxSCF,thresh, &
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwEW,EwGICEW,F)
!------------------------------------------------------------------------
! LIM excitation energies

4
src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90
View File

@ -26,8 +26,8 @@ subroutine RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,r
! Weight-dependent Cx coefficient for RMFL20 exchange functional
Cx(1) = -(4d0/3d0)*(1d0/pi)**(1d0/3d0)
Cx(2) = -(176d0/105d0)*(1d0/pi)**(1d0/3d0)
Cx(1) = Cx0
Cx(2) = Cx1
! Compute correlation energy for ground, singly-excited and doubly-excited states

9
src/eDFT/RMFL20_lda_exchange_energy.f90
View File

@ -17,9 +17,6 @@ subroutine RMFL20_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
! Local variables
integer :: iG
double precision :: Cx0
double precision :: Cx1
double precision :: CxLDA
double precision :: Cxw
double precision :: r
@ -27,11 +24,7 @@ subroutine RMFL20_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
double precision :: Ex
! Cx coefficient for Slater LDA exchange
Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0)
Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0)
CxLDA = - (3d0/4d0)*(3d0/pi)**(1d0/3d0)
! Weight-denepdent Cx coefficient
Cxw = CxLDA + wEns(2)*(Cx1 - Cx0)

9
src/eDFT/RMFL20_lda_exchange_individual_energy.f90
View File

@ -17,9 +17,6 @@ subroutine RMFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho
! Local variables
integer :: iG
double precision :: Cx0
double precision :: Cx1
double precision :: CxLDA
double precision :: Cxw
double precision :: r,rI
double precision :: e,dedr
@ -30,11 +27,7 @@ subroutine RMFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho
! Weight-dependent Cx coefficient for RMFL20 exchange functional
Cx0 = -(4d0/3d0)*(1d0/pi)**(1d0/3d0)
Cx1 = -(176d0/105d0)*(1d0/pi)**(1d0/3d0)
CxLDA = -(3d0/4d0)*(3d0/pi)**(1d0/3d0)
Cxw = CxLDA + wEns(2)*(Cx1 - Cx0)
Cxw = CxLDA + wEns(2)*(Cx1 - Cx0)
! Compute LDA exchange matrix in the AO basis

9
src/eDFT/RMFL20_lda_exchange_potential.f90
View File

@ -18,9 +18,6 @@ subroutine RMFL20_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
! Local variables
integer :: mu,nu,iG
double precision :: Cx0
double precision :: Cx1
double precision :: CxLDA
double precision :: Cxw
double precision :: r,vAO
@ -30,11 +27,7 @@ subroutine RMFL20_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
! Weight-dependent Cx coefficient for RMFL20 exchange functional
Cx0 = -(4d0/3d0)*(1d0/pi)**(1d0/3d0)
Cx1 = -(176d0/105d0)*(1d0/pi)**(1d0/3d0)
CxLDA = -(3d0/4d0)*(3d0/pi)**(1d0/3d0)
Cxw = CxLDA + wEns(2)*(Cx1 - Cx0)
Cxw = CxLDA + wEns(2)*(Cx1 - Cx0)
! Compute LDA exchange matrix in the AO basis

7
src/eDFT/RS51_lda_exchange_energy.f90
View File

@ -15,16 +15,11 @@ subroutine RS51_lda_exchange_energy(nGrid,weight,rho,Ex)
integer :: iG
double precision :: r
double precision :: Cx
! Output variables
double precision :: Ex
! Cx coefficient for Slater LDA exchange
Cx = -(3d0/4d0)*(3d0/pi)**(1d0/3d0)
! Compute LDA exchange energy
Ex = 0d0
@ -34,7 +29,7 @@ subroutine RS51_lda_exchange_energy(nGrid,weight,rho,Ex)
if(r > threshold) then
Ex = Ex + weight(iG)*Cx*r**(4d0/3d0)
Ex = Ex + weight(iG)*CxLDA*r**(4d0/3d0)
endif

9
src/eDFT/RS51_lda_exchange_individual_energy.f90
View File

@ -15,7 +15,6 @@ subroutine RS51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex)
! Local variables
integer :: iG
double precision :: Cx
double precision :: r,rI
double precision :: e,dedr
@ -23,10 +22,6 @@ subroutine RS51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex)
double precision,intent(out) :: Ex
! Cx coefficient for Slater LDA exchange
Cx = -(3d0/4d0)*(3d0/pi)**(1d0/3d0)
! Compute LDA exchange matrix in the AO basis
Ex = 0d0
@ -37,8 +32,8 @@ subroutine RS51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex)
if(r > threshold .and. rI > threshold) then
e = Cx*r**(1d0/3d0)
dedr = 1d0/3d0*Cx*r**(-2d0/3d0)
e = CxLDA*r**(1d0/3d0)
dedr = 1d0/3d0*CxLDA*r**(-2d0/3d0)
Ex = Ex + weight(iG)*(e*rI + dedr*r*rI - dedr*r*r)
endif

7
src/eDFT/RS51_lda_exchange_potential.f90
View File

@ -16,17 +16,12 @@ subroutine RS51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx)
! Local variables
integer :: mu,nu,iG
double precision :: Cx
double precision :: r,vAO
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Cx coefficient for Slater LDA exchange
Cx = -(3d0/4d0)*(3d0/pi)**(1d0/3d0)
! Compute LDA exchange matrix in the AO basis
Fx(:,:) = 0d0
@ -39,7 +34,7 @@ subroutine RS51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx)
if(r > threshold) then
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*Cx*r**(1d0/3d0)
Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*CxLDA*r**(1d0/3d0)
endif

20
src/eDFT/eDFT.f90
View File

@ -18,8 +18,13 @@ program eDFT
double precision,allocatable :: DShell(:,:)
double precision,allocatable :: ExpShell(:,:)
double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),X(:,:)
double precision,allocatable :: S(:,:)
double precision,allocatable :: T(:,:)
double precision,allocatable :: V(:,:)
double precision,allocatable :: Hc(:,:)
double precision,allocatable :: X(:,:)
double precision,allocatable :: ERI(:,:,:,:)
double precision,allocatable :: F(:,:)
character(len=7) :: method
integer :: x_rung,c_rung
@ -95,7 +100,8 @@ program eDFT
! Memory allocation for one- and two-electron integrals
allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI(nBas,nBas,nBas,nBas))
allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas), &
X(nBas,nBas),ERI(nBas,nBas,nBas,nBas),F(nBas,nBas))
! Read integrals
@ -137,8 +143,9 @@ program eDFT
if(method == 'GOK-RKS') then
call cpu_time(start_KS)
call GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, &
nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,EwGIC)
call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, &
nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc, &
Ew,EwGIC,F(:,:))
call cpu_time(end_KS)
t_KS = end_KS - start_KS
@ -154,8 +161,9 @@ program eDFT
if(method == 'LIM-RKS') then
call cpu_time(start_KS)
call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, &
nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc)
call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, &
nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc, &
F(:,:))
call cpu_time(end_KS)
t_KS = end_KS - start_KS

6
src/eDFT/fock_exchange_potential.f90
View File

@ -21,9 +21,9 @@ subroutine fock_exchange_potential(nBas,P,ERI,Fx)
! Compute HF exchange matrix
Fx(:,:) = 0d0
do nu=1,nBas
do si=1,nBas
do la=1,nBas
do si=1,nBas
do la=1,nBas
do nu=1,nBas
do mu=1,nBas
Fx(mu,nu) = Fx(mu,nu) - P(la,si)*ERI(mu,la,si,nu)
enddo

8
src/eDFT/hartree_coulomb.f90
View File

@ -19,10 +19,10 @@ subroutine hartree_coulomb(nBas,P,ERI,J)
double precision,intent(out) :: J(nBas,nBas)
J = 0d0
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do si=1,nBas
do si=1,nBas
do la=1,nBas
do nu=1,nBas
do mu=1,nBas
J(mu,nu) = J(mu,nu) + P(la,si)*ERI(mu,la,nu,si)
enddo
enddo