PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 12:23:42 +01:00
Code Releases Activity

starting modifying code for parameters of CC functional

Browse Source
This commit is contained in:
Clotilde Marut 2020-07-02 15:40:30 +02:00
parent 1753f4190d
commit c081a2664a
11 changed files with 82 additions and 104 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
input/dft
View File

@ -6,7 +6,7 @@
# GGA = 2: RB88
# Hybrid = 4
# Hartree-Fock = 666
1 S51
1 US51
# correlation rung:
# Hartree = 0
# LDA = 1: RVWN5,RMFL20
@ -19,6 +19,9 @@
# Number of states in ensemble (nEns)
3
# Ensemble weights: wEns(1),...,wEns(nEns-1)
0.333 0.333
0.000 0.000
# Parameters for CC weight-dependent exchange functional
0.420243 0.0700561 -0.288301
0.135068 -0.00774769 -0.0278205
# GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type
32 0.00001 T 5 1 1

18
src/eDFT/UCC_lda_exchange_energy.f90
View File

@ -1,4 +1,4 @@
subroutine UCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
! Compute the unrestricted version of the curvature-corrected exchange functional
@ -9,6 +9,8 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
double precision,intent(in) :: aCC_w1(3)
double precision,intent(in) :: aCC_w2(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
@ -52,15 +54,15 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
! Parameters for He N -> N-1
a1 = 0.420243d0
b1 = 0.0700561d0
c1 = -0.288301d0
a1 = aCC_w1(1)
b1 = aCC_w1(2)
c1 = aCC_w1(3)
! Parameters for He N -> N+1
a2 = 0.135068d0
b2 = -0.00774769d0
c2 = -0.0278205d0
a2 = aCC_w2(1)
b2 = aCC_w2(2)
c2 = aCC_w2(3)
! Cx coefficient for unrestricted Slater LDA exchange
@ -70,7 +72,7 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
! Fx2 for states N and N+1
w1 = wEns(2)
Fx1 = 1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)
Fx1 = 1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)
w2 = wEns(3)
Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)

10
src/eDFT/eDFT.f90
View File

@ -41,6 +41,8 @@ program eDFT
integer :: nGrid
double precision,allocatable :: root(:,:)
double precision,allocatable :: weight(:)
double precision :: aCC_w1(3)
double precision :: aCC_w2(3)
double precision,allocatable :: AO(:,:)
double precision,allocatable :: dAO(:,:,:)
@ -101,8 +103,8 @@ program eDFT
! Allocate ensemble weights
allocate(wEns(maxEns))
call read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn, &
nEns,wEns,maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type)
call read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, &
maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type)
!------------------------------------------------------------------------
! Read one- and two-electron integrals
@ -233,8 +235,8 @@ program eDFT
if(method == 'eDFT-UKS') then
call cpu_time(start_KS)
call eDFT_UKS(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)
call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, &
nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew)
call cpu_time(end_KS)
t_KS = end_KS - start_KS

6
src/eDFT/eDFT_UKS.f90
View File

@ -1,4 +1,4 @@
subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thresh,max_diis,guess_type, &
subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,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
@ -12,6 +12,8 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thre
character(len=12),intent(in) :: x_DFA,c_DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
double precision,intent(in) :: aCC_w1(3)
double precision,intent(in) :: aCC_w2(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: maxSCF,max_diis,guess_type
@ -313,7 +315,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thre
! Exchange energy
do ispin=1,nspin
call exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),nBas, &
call exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin))
end do

6
src/eDFT/exchange_energy.f90
View File

@ -1,4 +1,4 @@
subroutine exchange_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,P,FxHF,rho,drho,Ex)
subroutine exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P,FxHF,rho,drho,Ex)
! Compute the exchange energy
@ -12,6 +12,8 @@ subroutine exchange_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,P,F
logical,intent(in) :: LDA_centered
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
double precision,intent(in) :: aCC_w1(3)
double precision,intent(in) :: aCC_w2(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
@ -41,7 +43,7 @@ subroutine exchange_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,P,F
case(1)
call lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rho,ExLDA)
call lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,ExLDA)
Ex = ExLDA

8
src/eDFT/lda_exchange_energy.f90
View File

@ -1,4 +1,4 @@
subroutine lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rho,Ex)
subroutine lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
! Select LDA exchange functional
@ -11,6 +11,8 @@ subroutine lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rho,Ex)
logical,intent(in) :: LDA_centered
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
double precision,intent(in) :: aCC_w1(3)
double precision,intent(in) :: aCC_w2(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
@ -37,11 +39,11 @@ subroutine lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rho,Ex)
case ('RCC')
call RCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
call RCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
case ('UCC')
call UCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
call UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
case default

2
src/eDFT/print_unrestricted_individual_energy.f90
View File

@ -1,5 +1,5 @@
subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,Eaux,ExDD,EcDD,ExcDD,E, &
Om,Omx,Omc,Omxc,Omaux,OmxDD,OmcDD,OmxcDD)
Om,Omx,Omc,Omxc,Omaux,OmxDD,OmcDD,OmxcDD)
! Print individual energies for eDFT calculation

18
src/eDFT/read_options.f90
View File

@ -1,4 +1,5 @@
subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type)
subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, &
maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type)
! Read DFT options
@ -14,10 +15,12 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,maxSCF,th
character(len=8),intent(out) :: method
integer,intent(out) :: x_rung,c_rung
character(len=12),intent(out) :: x_DFA, c_DFA
character(len=12),intent(out) :: x_DFA, c_DFA
integer,intent(out) :: SGn
integer,intent(out) :: nEns
double precision,intent(out) :: wEns(maxEns)
double precision,intent(out) :: aCC_w1(3)
double precision,intent(out) :: aCC_w2(3)
integer,intent(out) :: maxSCF
double precision,intent(out) :: thresh
@ -99,6 +102,17 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,maxSCF,th
write(*,*)'----------------------------------------------------------'
call matout(nEns,1,wEns)
write(*,*)
! Read parameters for weight-dependent functional
read(1,*)
read(1,*) (aCC_w1(I),I=1,3)
read(1,*) (aCC_w2(I),I=1,3)
write(*,*)'----------------------------------------------------------'
write(*,*)' Ensemble weights '
write(*,*)'----------------------------------------------------------'
call matout(nEns,1,wEns)
write(*,*)
! Read KS options

10
src/eDFT/unrestricted_auxiliary_energy.f90
View File

@ -25,10 +25,10 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
iEns = 1
do ispin=1,nspin
Eaux(ispin,iEns) = sum(eps(1:nO(ispin),ispin))
Eaux(ispin,iEns) = sum(eps(1:nO(ispin),ispin))
end do
! N-1-electron ground state
! (N-1)-electron ground state
iEns = 2
@ -40,13 +40,11 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
Eaux(2,iEns) = 0d0
end if
! N+1 ground state
! (N+1)-electron ground state
iEns = 3
Eaux(1,iEns) = sum(eps(1:nO(1)+1,1))
Eaux(2,iEns) = sum(eps(1:nO(2),2))
Eaux(2,iEns) = sum(eps(1:nO(2) ,2))
end subroutine unrestricted_auxiliary_energy

97
src/eDFT/unrestricted_individual_energy.f90
View File

@ -1,5 +1,5 @@
subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,dAO, &
nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om)
subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,dAO, &
nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om)
! Compute unrestricted individual energies as well as excitation energies
@ -48,16 +48,16 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
double precision :: EJ(nsp,nEns)
double precision :: Ex(nspin,nEns)
double precision :: Ec(nsp,nEns)
double precision :: Exc(nEns) !!!!!
double precision :: Eaux(nspin,nEns) !!!!!
double precision :: Exc(nEns)
double precision :: Eaux(nspin,nEns)
double precision :: ExDD(nspin,nEns) !!!!!
double precision :: ExDD(nspin,nEns)
double precision :: EcDD(nsp,nEns)
double precision :: ExcDD(nspin,nEns) !!!!!
double precision :: ExcDD(nsp,nEns)
double precision :: Omx(nEns), Omc(nEns), Omxc(nEns) !!!!!
double precision :: Omaux(nEns) !!!!!
double precision :: OmxDD(nEns),OmcDD(nEns),OmxcDD(nEns) !!!!!
double precision :: Omx(nEns),Omc(nEns),Omxc(nEns)
double precision :: Omaux(nEns)
double precision :: OmxDD(nEns),OmcDD(nEns),OmxcDD(nEns)
double precision,external :: trace_matrix
@ -129,7 +129,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
do iEns=1,nEns
call unrestricted_correlation_individual_energy(c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight, &
rhow,drhow,rho(:,:,iEns),drho(:,:,:,iEns),Ec(:,iEns))
end do
end do
!------------------------------------------------------------------------
! Compute auxiliary energies
@ -140,75 +140,29 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
!------------------------------------------------------------------------
! Compute derivative discontinuities
!------------------------------------------------------------------------
do iEns=1,nEns
do ispin=1,nspin
call exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:), &
rhow(:,ispin),drhow(:,:,ispin),ExDD(ispin,iEns))
do ispin=1,nspin
call restricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:), &
rhow(:,ispin),drhow(:,:,ispin),EcDD(:,iEns))
! EcDD(ispin,:) = 0.d0 !!!!!!!!!
ExcDD(ispin,iEns) = ExDD(ispin,iEns) +sum(EcDD(:,iEns))
end do
end do
call exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns,nGrid,weight, &
rhow(:,ispin),drhow(:,:,ispin),ExDD(ispin,:))
end do
call unrestricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,EcDD)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!------------------------------------------------------------------------
! Exchange energy
!------------------------------------------------------------------------
! do iEns=1,nEns
! do ispin=1,nspin
! call exchange_potential(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,P(:,:,ispin,iEns),ERI(:,:,:,:), &
! AO(:,:),dAO(:,:,:),rho(:,ispin,iEns),drho(:,:,ispin,iEns),Fx(:,:,ispin),FxHF(:,:,ispin))
! call exchange_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,P(:,:,ispin,iEns),FxHF(:,:,ispin), &
! rho(:,ispin,iEns),drho(:,:,ispin,iEns),Ex(ispin,iEns))
! end do
! end do
!------------------------------------------------------------------------
! Correlation energy
!------------------------------------------------------------------------
! do iEns=1,nEns
! call correlation_individual_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),drhow(:,:,:), &
! rho(:,:,iEns),drho(:,:,:,iEns),Ec(:,iEns))
! end do
!------------------------------------------------------------------------
! Compute derivative discontinuities
!------------------------------------------------------------------------
! call correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),drhow(:,:,:),EcDD(:,:))
ExcDD(1,:) = ExDD(1,:) + EcDD(1,:)
ExcDD(2,:) = EcDD(2,:)
ExcDD(3,:) = ExDD(2,:) + EcDD(3,:)
!------------------------------------------------------------------------
! Total energy
!------------------------------------------------------------------------
! do iEns=1,nEns
! E(iEns) = ENuc + sum(ET(:,iEns)) + sum(EV(:,iEns)) + sum(EJ(:,iEns)) &
! + sum(Ex(:,iEns)) + sum(Ec(:,iEns)) + sum(EcDD(:,iEns))
! end do
do iEns=1,nEns
Exc(iEns) = sum(Ex(:,iEns)) + sum(Ec(:,iEns))
E(iEns) =sum( ET(:,iEns)) + sum( EV(:,iEns)) + sum(EJ(:,iEns)) &
E(iEns) = sum(ET(:,iEns)) + sum(EV(:,iEns)) + sum(EJ(:,iEns)) &
+ sum(Ex(:,iEns)) + sum(Ec(:,iEns)) + sum(ExcDD(:,iEns))
end do
!------------------------------------------------------------------------
! Excitation energies
!------------------------------------------------------------------------
@ -216,11 +170,11 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
do iEns=1,nEns
Om(iEns) = E(iEns) - E(1)
Omx(iEns) = sum(Ex(:,iEns)) -sum(Ex(:,1))
Omc(iEns) = sum(Ec(:,iEns)) -sum(Ec(:,1))
Omxc(iEns) = Exc(iEns) - Exc(1)
Omx(iEns) = sum(Ex(:,iEns)) - sum(Ex(:,1))
Omc(iEns) = sum(Ec(:,iEns)) - sum(Ec(:,1))
Omxc(iEns) = Exc(iEns) - Exc(1)
Omaux(iEns) = sum(Eaux(:,iEns)) -sum(Eaux(:,1))
Omaux(iEns) = sum(Eaux(:,iEns)) - sum(Eaux(:,1))
OmxDD(iEns) = sum(ExDD(:,iEns)) - sum(ExDD(:,1))
OmcDD(iEns) = sum(EcDD(:,iEns)) - sum(EcDD(:,1))
@ -233,8 +187,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
! Dump results
!------------------------------------------------------------------------
call print_unrestricted_individual_energy(nEns,ENuc,Ew,ET(:,:),EV(:,:),EJ(:,:),Ex(:,:),Ec(:,:),Exc(:), &
Eaux(:,:),ExDD(:,:),EcDD(:,:),ExcDD(:,:),E(:), &
Om(:),Omx(:),Omc(:),Omxc(:),Omaux(:),OmxDD(:),OmcDD(:),OmxcDD(:))
call print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,Eaux,ExDD,EcDD,ExcDD,E, &
Om,Omx,Omc,Omxc,Omaux,OmxDD,OmcDD,OmxcDD)
end subroutine unrestricted_individual_energy

4
src/eDFT/unrestricted_lda_correlation_derivative_discontinuity.f90
View File

@ -28,13 +28,13 @@ subroutine unrestricted_lda_correlation_derivative_discontinuity(DFA,nEns,wEns,n
! Wigner's LDA correlation functional: Wigner, Trans. Faraday Soc. 34 (1938) 678
case ('RW38')
case ('UW38')
Ec(:,:) = 0d0
! Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
case ('RVWN5')
case ('UVWN5')
Ec(:,:) = 0d0