LCPQ
/
quack
mirror of https://github.com/pfloos/quack
4
1
Fork 0
Code Releases Activity

scaling of N-centered

This commit is contained in:
Clotilde Marut 2020-09-29 09:43:01 +02:00
parent 1d3e156838
commit a96719d0c2
17 changed files with 236 additions and 198 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
input/basis
View File

@ -1,30 +1,18 @@
1 6
S 8
1 1469.0000000 0.0007660
2 220.5000000 0.0058920
3 50.2600000 0.0296710
4 14.2400000 0.1091800
5 4.5810000 0.2827890
6 1.5800000 0.4531230
7 0.5640000 0.2747740
8 0.0734500 0.0097510
S 8
1 1469.0000000 -0.0001200
2 220.5000000 -0.0009230
3 50.2600000 -0.0046890
4 14.2400000 -0.0176820
5 4.5810000 -0.0489020
6 1.5800000 -0.0960090
7 0.5640000 -0.1363800
8 0.0734500 0.5751020
1 3
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
S 1
1 0.0280500 1.0000000
P 3
1 1.5340000 0.0227840
2 0.2749000 0.1391070
3 0.0736200 0.5003750
1 0.1220000 1.0000000
P 1
1 0.0240300 1.0000000
D 1
1 0.1239000 1.0000000
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

26
input/dft
View File

@ -18,22 +18,22 @@
1
# Number of states in ensemble (nEns)
3
# occupation numbers of orbitals
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# occupation numbers of orbitals nO and nO+1
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# Ensemble weights: wEns(1),...,wEns(nEns-1)
0.00 1.00
0.00 1.00
# Ncentered ? 0 for NO
1
0
# Parameters for CC weight-dependent exchange functional
0.420431 0.069097 -0.295049
0.135075 -0.00770826 -0.028057
0.445525 0.0901503 -0.286898
0.191734 -0.0364788 -0.017035
# choice of UCC exchange coefficient : 1 for Cx1, 2 for Cx2, 3 for Cx1*Cx2
2
# GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type

5
input/molecule
View File

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

5
input/molecule.xyz
View File

@ -1,3 +1,4 @@
1
2
Li 0.0000000000 0.0000000000 0.0000000000
H 0.0000000000 0.0000000000 0.0000000000
H 0.0000000000 0.0000000000 0.7408481486

44
input/weight
View File

@ -1,30 +1,18 @@
1 6
S 8
1 1469.0000000 0.0007660
2 220.5000000 0.0058920
3 50.2600000 0.0296710
4 14.2400000 0.1091800
5 4.5810000 0.2827890
6 1.5800000 0.4531230
7 0.5640000 0.2747740
8 0.0734500 0.0097510
S 8
1 1469.0000000 -0.0001200
2 220.5000000 -0.0009230
3 50.2600000 -0.0046890
4 14.2400000 -0.0176820
5 4.5810000 -0.0489020
6 1.5800000 -0.0960090
7 0.5640000 -0.1363800
8 0.0734500 0.5751020
1 3
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
S 1
1 0.0280500 1.0000000
P 3
1 1.5340000 0.0227840
2 0.2749000 0.1391070
3 0.0736200 0.5003750
1 0.1220000 1.0000000
P 1
1 0.0240300 1.0000000
D 1
1 0.1239000 1.0000000
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

50
scan_w.sh
View File

@ -12,17 +12,27 @@ w1=0.00
XF=$3
CF=$4
# for H
#aw1="1.49852 7.79815 25.1445"
#aw2="0.424545 -0.0382349 -0.32472"
# for He
aw1="0.420431 0.069097 -0.295049"
aw2="0.135075 -0.00770826 -0.028057"
#aw1="0.429447 0.053506 -0.339391"
#aw2="0.254939 -0.0893396 0.00765453"
# for H2
#aw1="0.445525 0.0901503 -0.286898"
#aw2="0.191734 -0.0364788 -0.017035"
aw1="0.445525 0.0901503 -0.286898"
aw2="0.191734 -0.0364788 -0.017035"
# for Li
#aw1="0.0560976 -0.00796904 -0.0238162"
#aw2="0.0360106 0.00979306 -0.0172286"
#aw1="0.055105 -0.00943825 -0.0267771"
#aw2="0.0359827 0.0096623 -0.0173542"
# for Li+
#aw1="0.503566, 0.137076, -0.348529"
#aw2="0.0553828, 0.00830375, -0.0234602"
# for B
#aw1="0.052676 -0.00624118 -0.000368825"
@ -36,6 +46,9 @@ aw2="0.135075 -0.00770826 -0.028057"
#aw1="-0.00201219 -0.00371002 -0.00212719"
#aw2="-0.00117715 0.00188738 -0.000414761"
# for Be
#aw1="0.0663282, -0.0117682, -0.0335909"
#aw2="0.0479262, 0.00966351, -0.0208712"
DATA=${MOL}_${BASIS}_${XF}_${CF}_${w2}.dat
@ -65,6 +78,15 @@ do
echo " 1" >> input/dft
echo "# Number of states in ensemble (nEns)" >> input/dft
echo " 3" >> input/dft
echo "# occupation numbers of orbitals nO and nO+1" >> input/dft
echo " 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
echo " 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
echo " " >> input/dft
echo " 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
echo " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
echo " " >> input/dft
echo " 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
echo " 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
echo "# Ensemble weights: wEns(1),...,wEns(nEns-1)" >> input/dft
echo " ${w1} ${w2} " >> input/dft
echo "# Ncentered ? 0 for NO " >> input/dft
@ -74,13 +96,6 @@ do
echo ${aw2} >> input/dft
echo "# choice of UCC exchange coefficient : 1 for Cx1, 2 for Cx2, 3 for Cx1*Cx2" >> input/dft
echo "2" >> input/dft
echo "# occupation numbers of orbitals nO and nO+1" >> input/dft
echo " 1.00 0.00 " >> input/dft
echo " 0.00 0.00 " >> input/dft
echo " 0.00 0.00 " >> input/dft
echo " 0.00 0.00 " >> input/dft
echo " 1.00 0.00 " >> input/dft
echo " 0.00 1.00 " >> input/dft
echo "# GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type" >> input/dft
echo " 1000 0.00001 T 5 1 1" >> input/dft
OUTPUT=${MOL}_${BASIS}_${XF}_${CF}_${w2}.out
@ -93,7 +108,12 @@ do
EA=`grep "Electronic Affinity" ${OUTPUT} | grep " au" | tail -1 | cut -d":" -f 2 | sed 's/au//'`
FG=`grep "Fundamental Gap" ${OUTPUT} | grep " au" | tail -1 | cut -d":" -f 2 | sed 's/au//'`
Ex=`grep "Exchange energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/au//'`
echo $w1 $w2 $Ew $E0 $E1 $E2 $IP $EA $FG $Ex
echo $w1 $w2 $Ew $E0 $E1 $E2 $IP $EA $FG $Ex >> ${DATA}
HOMOa=`grep "HOMO a energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/eV//'`
LUMOa=`grep "LUMO a energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/eV//'`
HOMOb=`grep "HOMO a energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/eV//'`
LUMOb=`grep "LUMO b energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/eV//'`
echo $w1 $w2 $Ew $E0 $E1 $E2 $IP $EA $FG $Ex $HOMOa $LUMOa $HOMOb $LUMOb
echo $w1 $w2 $Ew $E0 $E1 $E2 $IP $EA $FG $Ex $HOMOa $LUMOa $HOMOb $LUMOb >> ${DATA}
done

64
src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90
View File

@ -1,4 +1,4 @@
subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD,Cx_choice)
subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,Cx_choice,doNcentered,ExDD)
! Compute the unrestricted version of the curvature-corrected exchange ensemble derivative
@ -15,6 +15,8 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
integer,intent(in) :: Cx_choice
integer,intent(in) :: doNcentered
! Local variables
@ -28,51 +30,39 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
double precision :: a2,b2,c2,w2
double precision :: dCxdw1,dCxdw2
double precision :: nEli,nElw
! Output variables
double precision,intent(out) :: ExDD(nEns)
! External variable
double precision,external :: electron_number
! Memory allocation
allocate(dExdw(nEns))
! Single excitation parameters
! a1 = 0.0d0
! b1 = 0.0d0
! c1 = 0.0d0
! Parameters for H2 at equilibrium
! a2 = +0.5751782560799208d0
! b2 = -0.021108186591137282d0
! c2 = -0.36718902716347124d0
! Parameters for stretch H2
! a2 = + 0.01922622507087411d0
! b2 = - 0.01799647558018601d0
! c2 = - 0.022945430666782573d0
! Parameters for He
! a2 = 1.9125735895875828d0
! b2 = 2.715266992840757d0
! c2 = 2.1634223380633086d0
! Parameters for He N -> N-1
! Parameters for N -> N-1
a1 = aCC_w1(1)
b1 = aCC_w1(2)
c1 = aCC_w1(3)
! Parameters for He N -> N+1
! Parameters for N -> N+1
a2 = aCC_w2(1)
b2 = aCC_w2(2)
c2 = aCC_w2(3)
nElw = electron_number(nGrid,weight,rhow)
! Cx coefficient for unrestricted Slater LDA exchange
@ -96,23 +86,6 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
* (0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0)))
end select
! Double weight-dependency
! dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0))) &
! * (1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2))
! dCxdw2 = (1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)) &
! * (0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0)))
! left single-weight-dependency
! dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0)))
! dCxdw2 = 0.d0
! right single-weight-dependency
! dCxdw1 = 0.d0
! dCxdw2 =(0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0)))
dCxdw1 = alpha*dCxdw1
dCxdw2 = alpha*dCxdw2
@ -143,4 +116,9 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
end do
end do
if (doNcentered .NE. 0) then
ExDD(2) = ((nElw-1)/nElw)*ExDD(2)
ExDD(3) = ((nElw+1)/nElw)*ExDD(3)
end if
end subroutine UCC_lda_exchange_derivative_discontinuity

56
src/eDFT/UCC_lda_exchange_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex,Cx_choice)
subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Cx_choice,doNcentered,Ex)
! Compute the unrestricted version of the curvature-corrected exchange functional
@ -16,12 +16,14 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
integer,intent(in) :: doNcentered
! Local variables
integer :: iG
double precision :: r,rI,alpha
double precision :: e_p,dedr
double precision :: nEli,nElw
double precision :: a1,b1,c1,w1
double precision :: a2,b2,c2,w2
@ -31,37 +33,18 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
double precision,intent(out) :: Ex
! Single excitation parameter
! External variable
! a1 = 0.0d0
! b1 = 0.0d0
! c1 = 0.0d0
double precision,external :: electron_number
! Parameters for H2 at equilibrium
! a2 = +0.5751782560799208d0
! b2 = -0.021108186591137282d0
! c2 = -0.36718902716347124d0
! Parameters for stretch H2
! a2 = + 0.01922622507087411d0
! b2 = - 0.01799647558018601d0
! c2 = - 0.022945430666782573d0
! Parameters for He
! a2 = 1.9125735895875828d0
! b2 = 2.715266992840757d0
! c2 = 2.1634223380633086d0
! Parameters for He N -> N-1
! Parameters for N -> N-1
a1 = aCC_w1(1)
b1 = aCC_w1(2)
c1 = aCC_w1(3)
! Parameters for He N -> N+1
! Parameters for N -> N+1
a2 = aCC_w2(1)
b2 = aCC_w2(2)
@ -86,14 +69,10 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
Cx = alpha*Fx2*Fx1
end select
! for two-weight ensembles
! Cx = alpha*Fx1*Fx2
nEli = electron_number(nGrid,weight,rho)
! for left ensembles
! Cx = alpha*Fx1
! for right ensembles
! Cx = alpha*Fx2
nElw = electron_number(nGrid,weight,rhow)
! Compute LDA exchange matrix in the AO basis
@ -107,13 +86,20 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
e_p = Cx*r**(1d0/3d0)
dedr = 1d0/3d0*Cx*r**(-2d0/3d0)
Ex = Ex - weight(iG)*dedr*r*r
if (doNcentered == 0) then
Ex = Ex - weight(iG)*dedr*r*r
else
Ex = Ex - weight(iG)*dedr*r*r*(nEli/nElw)
end if
if(rI > threshold) then
Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI)
if (doNcentered == 0) then
Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI)
else
Ex = Ex + weight(iG)*((nEli/nElw)*e_p*rI + dedr*r*rI)
end if
endif
endif

26
src/eDFT/US51_lda_exchange_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex)
subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered,Ex)
! Compute the restricted version of Slater's LDA exchange individual energy
@ -11,17 +11,29 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex)
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: doNcentered
! Local variables
integer :: iG
double precision :: r,rI,alpha
double precision :: e,dedr
double precision :: nEli,nElw
! Output variables
double precision,intent(out) :: Ex
! External variable
double precision,external :: electron_number
nEli = electron_number(nGrid,weight,rho)
nElw = electron_number(nGrid,weight,rhow)
! Compute LDA exchange matrix in the AO basis
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
@ -37,11 +49,19 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex)
e = alpha*r**(1d0/3d0)
dedr = 1d0/3d0*alpha*r**(-2d0/3d0)
Ex = Ex - weight(iG)*dedr*r*r
if (doNcentered == 0) then
Ex = Ex - weight(iG)*dedr*r*r
else
Ex = Ex - weight(iG)*dedr*r*r*(nEli/nElw)
end if
if(rI > threshold) then
Ex = Ex + weight(iG)*(e*rI + dedr*r*rI)
if (doNcentered == 0) then
Ex = Ex + weight(iG)*(e*rI + dedr*r*rI)
else
Ex = Ex + weight(iG)*((nEli/nElw)*e*rI + dedr*r*rI)
end if
endif

9
src/eDFT/eDFT.f90
View File

@ -187,7 +187,7 @@ program eDFT
call cpu_time(start_KS)
call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), &
aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), &
S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice)
S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice,doNcentered)
call cpu_time(end_KS)
t_KS = end_KS - start_KS
@ -205,7 +205,7 @@ program eDFT
call cpu_time(start_KS)
call MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), &
aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), &
S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice)
S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice,doNcentered)
call cpu_time(end_KS)
t_KS = end_KS - start_KS
@ -222,7 +222,8 @@ program eDFT
call cpu_time(start_KS)
call GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type, &
nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,occnum,Cx_choice)
nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,occnum, &
Cx_choice,doNcentered)
call cpu_time(end_KS)
t_KS = end_KS - start_KS
@ -239,7 +240,7 @@ program eDFT
call cpu_time(start_KS)
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,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice)
nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice,doNcentered)
call cpu_time(end_KS)
t_KS = end_KS - start_KS

6
src/eDFT/eDFT_UKS.f90
View File

@ -1,5 +1,5 @@
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,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice)
nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice,doNcentered)
! Perform unrestricted Kohn-Sham calculation for ensembles
@ -31,6 +31,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
double precision,intent(in) :: ENuc
double precision,intent(in) :: occnum(nBas,nspin,nEns)
integer,intent(in) :: Cx_choice
integer,intent(in) :: doNcentered
! Local variables
@ -380,6 +381,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
!------------------------------------------------------------------------
call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
AO,dAO,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum,Cx_choice)
AO,dAO,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum, &
Cx_choice,doNcentered)
end subroutine eDFT_UKS

21
src/eDFT/unrestricted_auxiliary_energy.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_auxiliary_energy(nBas,nEns,eps,Eaux,occnum)
subroutine unrestricted_auxiliary_energy(nBas,nEns,eps,occnum,doNcentered,Eaux)
! Compute the auxiliary KS energies
@ -12,18 +12,28 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,eps,Eaux,occnum)
integer,intent(in) :: nEns
double precision,intent(in) :: eps(nBas,nspin)
double precision,intent(in) :: occnum(nBas,nspin,nEns)
integer,intent(in) :: doNcentered
! Local variables
integer :: iEns
integer :: iEns,iBas
integer :: ispin
integer :: p
double precision,allocatable :: nEl(:)
! Output variables
double precision,intent(out) :: Eaux(nspin,nEns)
allocate(nEl(nEns))
nEl(:) = 0d0
do iEns=1,nEns
do iBas=1,nBas
nEl(iEns) = nEl(iEns) + occnum(iBas,1,iEns) + occnum(iBas,2,iEns)
end do
end do
! Compute auxiliary energies for each state of the ensemble based on occupation numbers
Eaux(:,:) = 0d0
@ -35,6 +45,11 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,eps,Eaux,occnum)
end do
end do
! if (doNcentered .NE. 0) then
! Eaux(:,iEns) = Eaux(:,iEns)*(nEl(iEns)/nEl(1))
! end if
end do
end subroutine unrestricted_auxiliary_energy

6
src/eDFT/unrestricted_exchange_derivative_discontinuity.f90
View File

@ -1,4 +1,5 @@
subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,drhow,ExDD,Cx_choice)
subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,drhow,&
Cx_choice,doNcentered,ExDD)
! Compute the exchange part of the derivative discontinuity
@ -18,6 +19,7 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid)
integer,intent(in) :: Cx_choice
integer,intent(in) :: doNcentered
! Local variables
@ -39,7 +41,7 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC
case(1)
call unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),&
rhow(:),ExDD(:),Cx_choice)
rhow(:),Cx_choice,doNcentered,ExDD(:))
! GGA functionals

5
src/eDFT/unrestricted_exchange_individual_energy.f90
View File

@ -1,5 +1,5 @@
subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
ERI,Pw,P,rhow,drhow,rho,drho,Ex,Cx_choice)
ERI,Pw,P,rhow,drhow,rho,drho,Cx_choice,doNcentered,Ex)
! Compute the exchange individual energy
@ -26,6 +26,7 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
integer,intent(in) :: Cx_choice
integer,intent(in) :: doNcentered
! Local variables
@ -50,7 +51,7 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE
case(1)
call unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,&
rhow,rho,ExLDA,Cx_choice)
rhow,rho,Cx_choice,doNcentered,ExLDA)
Ex = ExLDA

49
src/eDFT/unrestricted_individual_energy.f90
View File

@ -1,5 +1,6 @@
subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,dAO, &
T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum,Cx_choice)
T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum,&
Cx_choice,doNcentered)
! Compute unrestricted individual energies as well as excitation energies
@ -42,6 +43,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
double precision :: Ew
double precision,intent(in) :: occnum(nBas,nspin,nEns)
integer,intent(in) :: Cx_choice
integer,intent(in) :: doNcentered
! Local variables
@ -64,7 +66,9 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
double precision,external :: trace_matrix
integer :: ispin,iEns
integer :: ispin,iEns,iBas
double precision,allocatable :: nEl(:)
double precision,external :: electron_number
@ -73,15 +77,31 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
double precision,intent(out) :: E(nEns)
double precision,intent(out) :: Om(nEns)
allocate(nEl(nEns))
nEl(:) = 0d0
do iEns=1,nEns
do iBas=1,nBas
nEl(iEns) = nEl(iEns) + occnum(iBas,1,iEns) + occnum(iBas,2,iEns)
end do
end do
print*,'test1'
!------------------------------------------------------------------------
! Kinetic energy
!------------------------------------------------------------------------
do ispin=1,nspin
do iEns=1,nEns
ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:)))
if (doNcentered == 0) then
ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:)))
else
ET(ispin,iEns) = (nEl(iEns)/nEl(1))*trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:)))
end if
end do
end do
print*,'test2'
!------------------------------------------------------------------------
! Potential energy
@ -89,10 +109,15 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
do iEns=1,nEns
do ispin=1,nspin
EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:)))
end do
if (doNcentered == 0) then
EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:)))
else
EV(ispin,iEns) = (nEl(iEns)/nEl(1))*trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:)))
end if
end do
end do
print*,'test3'
!------------------------------------------------------------------------
! Individual Hartree energy
!------------------------------------------------------------------------
@ -113,9 +138,11 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
EJ(3,iEns) = trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) &
- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2)))
if (doNcentered .NE. 0) then
EJ(:,iEns) = (nEl(iEns)/nEl(1))*EJ(:,iEns)
end if
end do
print*,'test4'
!------------------------------------------------------------------------
! Checking Hartree contributions for each individual states
!------------------------------------------------------------------------
@ -139,10 +166,12 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
do ispin=1,nspin
call unrestricted_exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,ERI, &
Pw(:,:,ispin),P(:,:,ispin,iEns),rhow(:,ispin),drhow(:,:,ispin), &
rho(:,ispin,iEns),drho(:,:,ispin,iEns),Ex(ispin,iEns),Cx_choice)
rho(:,ispin,iEns),drho(:,:,ispin,iEns),Cx_choice,doNcentered,Ex(ispin,iEns))
end do
end do
print*,'test5'
!------------------------------------------------------------------------
! Checking exchange contributions for each individual states
!------------------------------------------------------------------------
@ -184,7 +213,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
! Compute auxiliary energies
!------------------------------------------------------------------------
call unrestricted_auxiliary_energy(nBas,nEns,eps,Eaux,occnum)
call unrestricted_auxiliary_energy(nBas,nEns,eps,occnum,doNcentered,Eaux)
!------------------------------------------------------------------------
! Compute derivative discontinuities
@ -193,7 +222,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
do ispin=1,nspin
call unrestricted_exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight, &
rhow(:,ispin),drhow(:,:,ispin),ExDD(ispin,:),Cx_choice)
rhow(:,ispin),drhow(:,:,ispin),Cx_choice,doNcentered,ExDD(ispin,:))
end do
call unrestricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,EcDD)

9
src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90
View File

@ -1,4 +1,5 @@
subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD,Cx_choice)
subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,&
Cx_choice,doNcentered,ExDD)
! Compute the exchange LDA part of the derivative discontinuity
@ -17,6 +18,7 @@ subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
integer,intent(in) :: Cx_choice
integer,intent(in) :: doNcentered
! Local variables
@ -25,7 +27,7 @@ subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_
double precision,intent(out) :: ExDD(nEns)
! Select correlation functional
! Select exchange functional
select case (DFA)
@ -35,7 +37,8 @@ subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_
case ('CC')
call UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),ExDD(:),Cx_choice)
call UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),&
Cx_choice,doNcentered,ExDD(:))
case default

9
src/eDFT/unrestricted_lda_exchange_individual_energy.f90
View File

@ -1,4 +1,5 @@
subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex,Cx_choice)
subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,&
Cx_choice,doNcentered,Ex)
! Compute LDA exchange energy for individual states
@ -18,6 +19,7 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEn
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
integer,intent(in) ::doNcentered
! Output variables
@ -29,11 +31,12 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEn
case ('S51')
call US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex)
call US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered,Ex)
case ('CC')
call UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex,Cx_choice)
call UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,&
Cx_choice,doNcentered,Ex)
case default