4
1
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:35:36 +01:00

saving eDFT

This commit is contained in:
Pierre-Francois Loos 2020-05-05 16:45:10 +02:00
parent d91631f73f
commit 431eab634d
17 changed files with 124 additions and 354 deletions

View File

@ -1,27 +1,62 @@
1 5
1 14
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 82.6400000 0.0020060
2 12.4100000 0.0153430
3 2.8240000 0.0755790
S 1
1 0.1220000 1.0000000
1 0.7977000 1.0000000
S 1
1 0.0297400 1.0000000
1 0.2581000 1.0000000
S 1
1 0.0898900 1.0000000
S 1
1 0.0236300 1.0000000
P 1
1 0.7270000 1.0000000
1 2.2920000 1.0000000
P 1
1 0.1410000 1.0000000
2 5
1 0.8380000 1.0000000
P 1
1 0.2920000 1.0000000
P 1
1 0.0848000 1.0000000
D 1
1 2.0620000 1.0000000
D 1
1 0.6620000 1.0000000
D 1
1 0.1900000 1.0000000
F 1
1 1.3970000 1.0000000
F 1
1 0.3600000 1.0000000
2 14
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 82.6400000 0.0020060
2 12.4100000 0.0153430
3 2.8240000 0.0755790
S 1
1 0.1220000 1.0000000
1 0.7977000 1.0000000
S 1
1 0.0297400 1.0000000
1 0.2581000 1.0000000
S 1
1 0.0898900 1.0000000
S 1
1 0.0236300 1.0000000
P 1
1 0.7270000 1.0000000
1 2.2920000 1.0000000
P 1
1 0.1410000 1.0000000
1 0.8380000 1.0000000
P 1
1 0.2920000 1.0000000
P 1
1 0.0848000 1.0000000
D 1
1 2.0620000 1.0000000
D 1
1 0.6620000 1.0000000
D 1
1 0.1900000 1.0000000
F 1
1 1.3970000 1.0000000
F 1
1 0.3600000 1.0000000

View File

@ -1,12 +1,12 @@
# Restricted or unrestricted KS calculation
MOM-RKS
LIM-RKS
# exchange rung:
# Hartree = 0
# LDA = 1: RS51,RMFL20
# GGA = 2: RB88
# Hybrid = 4
# Hartree-Fock = 666
1 RGIC
1 RCC
# correlation rung:
# Hartree = 0
# LDA = 1: RVWN5,RMFL20

View File

@ -7,7 +7,7 @@
# drCCD rCCD lCCD pCCD
F F F F
# CIS CID CISD
T F F
F F F
# RPA RPAx ppRPA
F F F
# G0F2 evGF2 G0F3 evGF3
@ -15,6 +15,6 @@
# G0W0 evGW qsGW
T F F
# G0T0 evGT qsGT
F F F
T F F
# MCMP2
F

View File

@ -8,8 +8,8 @@
T T
# GF: maxSCF thresh DIIS n_diis lin renorm
256 0.00001 T 5 T 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 T 0.00367493
# GW/GT: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta
256 0.00001 T 5 F F T F F F T 0.0
# ACFDT: AC Kx XBS
F F T
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift

View File

@ -1,27 +1,62 @@
1 5
1 14
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 82.6400000 0.0020060
2 12.4100000 0.0153430
3 2.8240000 0.0755790
S 1
1 0.1220000 1.0000000
1 0.7977000 1.0000000
S 1
1 0.0297400 1.0000000
1 0.2581000 1.0000000
S 1
1 0.0898900 1.0000000
S 1
1 0.0236300 1.0000000
P 1
1 0.7270000 1.0000000
1 2.2920000 1.0000000
P 1
1 0.1410000 1.0000000
2 5
1 0.8380000 1.0000000
P 1
1 0.2920000 1.0000000
P 1
1 0.0848000 1.0000000
D 1
1 2.0620000 1.0000000
D 1
1 0.6620000 1.0000000
D 1
1 0.1900000 1.0000000
F 1
1 1.3970000 1.0000000
F 1
1 0.3600000 1.0000000
2 14
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 82.6400000 0.0020060
2 12.4100000 0.0153430
3 2.8240000 0.0755790
S 1
1 0.1220000 1.0000000
1 0.7977000 1.0000000
S 1
1 0.0297400 1.0000000
1 0.2581000 1.0000000
S 1
1 0.0898900 1.0000000
S 1
1 0.0236300 1.0000000
P 1
1 0.7270000 1.0000000
1 2.2920000 1.0000000
P 1
1 0.1410000 1.0000000
1 0.8380000 1.0000000
P 1
1 0.2920000 1.0000000
P 1
1 0.0848000 1.0000000
D 1
1 2.0620000 1.0000000
D 1
1 0.6620000 1.0000000
D 1
1 0.1900000 1.0000000
F 1
1 1.3970000 1.0000000
F 1
1 0.3600000 1.0000000

View File

@ -328,7 +328,7 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGri
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
stop
! stop
end if

View File

@ -91,8 +91,8 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, &
write(*,'(A40)') '*************************************************'
write(*,*)
call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,nGrid,weight,maxSCF,thresh, &
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c)
! call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,nGrid,weight,maxSCF,thresh, &
! max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c)
!------------------------------------------------------------------------
! Equiensemble calculation
@ -122,7 +122,7 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, &
Om(:) = Ew(:) - Ew(1)
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' LINEAR INTERPOLATION METHOD EXCITATION ENERGIES '
write(*,'(A60)') ' MAXIMUM OVERLAP METHOD EXCITATION ENERGIES '
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A44,F16.10,A3)') ' Ensemble energy #1 ',Ew(1),' au'

View File

@ -1,85 +0,0 @@
subroutine RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rhow,ExDD)
! Compute the restricted version of the GIC exchange individual energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
! Local variables
integer :: iEns,jEns
integer :: iG
double precision :: r
double precision,allocatable :: dExdw(:)
double precision,external :: Kronecker_delta
double precision :: a,b,c,w
double precision :: dCxGICdw
! Output variables
double precision,intent(out) :: ExDD(nEns)
! Memory allocation
allocate(dExdw(nEns))
! Weight-dependent Cx coefficient for RMFL20 exchange functional
! Parameters for H2 at equilibrium
a = + 0.5739189000851961d0
b = - 0.0003469882157336496d0
c = - 0.2676338054343272d0
! Parameters for stretch H2
! a = + 0.01918229168254928d0
! b = - 0.01545313842512261d0
! c = - 0.012720073519142448d0
! Parameters for He
! a = 1.9015719148496788d0
! b = 2.5236598782764412d0
! c = 1.6652282199359842d0
w = 0.5d0*wEns(2) + wEns(3)
dCxGICdw = (0.5d0*b + (2d0*a + 0.5d0*c)*(w - 0.5d0) - (1d0 - w)*w*(3d0*b + 4d0*c*(w - 0.5d0)))
dCxGICdw = CxLDA*dCxGICdw
dExdw(:) = 0d0
do iG=1,nGrid
r = max(0d0,rhow(iG))
if(r > threshold) then
dExdw(1) = 0d0
dExdw(2) = dExdw(2) + 0.5d0*weight(iG)*dCxGICdw*r**(4d0/3d0)
dExdw(3) = dExdw(3) + 1.0d0*weight(iG)*dCxGICdw*r**(4d0/3d0)
end if
end do
ExDD(:) = 0d0
do iEns=1,nEns
do jEns=2,nEns
ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*dExdw(jEns)
end do
end do
end subroutine RGIC_lda_exchange_derivative_discontinuity

View File

@ -1,68 +0,0 @@
subroutine RGIC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
! Compute the restricted version of the GIC exchange functional
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
! Local variables
integer :: iG
double precision :: r
double precision :: a,b,c,w
double precision :: CxGIC
! Output variables
double precision :: Ex
! Weight-dependent Cx coefficient
! Weight-dependent Cx coefficient for RMFL20 exchange functional
! Parameters for H2 at equilibrium
a = + 0.5739189000851961d0
b = - 0.0003469882157336496d0
c = - 0.2676338054343272d0
! Parameters for stretch H2
! a = + 0.01918229168254928d0
! b = - 0.01545313842512261d0
! c = - 0.012720073519142448d0
! Parameters for He
! a = 1.9015719148496788d0
! b = 2.5236598782764412d0
! c = 1.6652282199359842d0
w = 0.5d0*wEns(2) + wEns(3)
CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2)
CxGIC = CxLDA*CxGIC
! Compute GIC-LDA exchange energy
Ex = 0d0
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) then
Ex = Ex + weight(iG)*CxGIC*r**(4d0/3d0)
endif
enddo
end subroutine RGIC_lda_exchange_energy

View File

@ -1,72 +0,0 @@
subroutine RGIC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,Ex)
! Compute the restricted version of the GIC exchange functional
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: rho(nGrid)
! Local variables
integer :: iG
double precision :: CxGIC
double precision :: r,rI
double precision :: e_p,dedr
double precision :: a,b,c,w
! Output variables
double precision,intent(out) :: Ex
! Weight-dependent Cx coefficient for RMFL20 exchange functional
! Parameters for H2 at equilibrium
a = + 0.5739189000851961d0
b = - 0.0003469882157336496d0
c = - 0.2676338054343272d0
! Parameters for stretch H2
! a = + 0.01918229168254928d0
! b = - 0.01545313842512261d0
! c = - 0.012720073519142448d0
! Parameters for He
! a = 1.9015719148496788d0
! b = 2.5236598782764412d0
! c = 1.6652282199359842d0
w = 0.5d0*wEns(2) + wEns(3)
CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2)
CxGIC = CxLDA*CxGIC
! Compute LDA exchange matrix in the AO basis
Ex = 0d0
do iG=1,nGrid
r = max(0d0,rhow(iG))
rI = max(0d0,rho(iG))
if(r > threshold .and. rI > threshold) then
e_p = CxGIC*r**(1d0/3d0)
dedr = 1d0/3d0*CxGIC*r**(-2d0/3d0)
Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI - dedr*r*r)
endif
enddo
end subroutine RGIC_lda_exchange_individual_energy

View File

@ -1,75 +0,0 @@
subroutine RGIC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
! Compute the restricted version of the GIC exchange potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
! Local variables
integer :: mu,nu,iG
double precision :: r,vAO
double precision :: a,b,c,w
double precision :: CxGIC
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Weight-dependent Cx coefficient for RMFL20 exchange functional
! Parameters for H2 at equilibrium
a = + 0.5739189000851961d0
b = - 0.0003469882157336496d0
c = - 0.2676338054343272d0
! Parameters for stretch H2
! a = + 0.01918229168254928d0
! b = - 0.01545313842512261d0
! c = - 0.012720073519142448d0
! Parameters for He
! a = 1.9015719148496788d0
! b = 2.5236598782764412d0
! c = 1.6652282199359842d0
w = 0.5d0*wEns(2) + wEns(3)
CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2)
CxGIC = CxLDA*CxGIC
! Compute LDA exchange matrix in the AO basis
Fx(:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) then
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*CxGIC*r**(1d0/3d0)
endif
enddo
enddo
enddo
end subroutine RGIC_lda_exchange_potential

View File

@ -37,9 +37,9 @@ subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow
call RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:))
case ('RGIC')
case ('RCC')
call RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:))
call RCC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:))
case default

View File

@ -35,9 +35,9 @@ subroutine lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rho,Ex)
call RMFL20_lda_exchange_energy(LDA_centered,nEns,wEns,nGrid,weight,rho,Ex)
case ('RGIC')
case ('RCC')
call RGIC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
call RCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
case default

View File

@ -32,9 +32,9 @@ subroutine lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weigh
call RMFL20_lda_exchange_individual_energy(LDA_centered,nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex)
case ('RGIC')
case ('RCC')
call RGIC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex)
call RCC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex)
case default

View File

@ -38,9 +38,9 @@ subroutine lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,A
call RMFL20_lda_exchange_potential(LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
case ('RGIC')
case ('RCC')
call RGIC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
call RCC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
case default

View File

@ -37,7 +37,7 @@ subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
Eaux(iEns) = 0d0
end if
Eaux(iEns) = Eaux(iEns) + eps(nO) + eps(nO+1)
Eaux(iEns) = Eaux(iEns) + eps(nO) + eps(nO+2)
! Doubly-excited state density matrix

View File

@ -38,7 +38,7 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P)
end if
P(:,:,iEns) = P(:,:,iEns) + 1d0*matmul(c(:,nO :nO ),transpose(c(:,nO :nO ))) &
+ 1d0*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1)))
+ 1d0*matmul(c(:,nO+2:nO+2),transpose(c(:,nO+2:nO+2)))
! Doubly-excited state density matrix