LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:35:36 +01:00
Code Releases Activity

occnum

Browse Source
This commit is contained in:
Clotilde Marut 2020-09-11 11:55:04 +02:00
parent dbb560f77a
commit 65ca5214b6
19 changed files with 243 additions and 297 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

108
examples/basis.C2-cc-pvdz
View File

@ -1,64 +1,58 @@
1 6
S 9
1 6.665000E+03 6.920000E-04
2 1.000000E+03 5.329000E-03
3 2.280000E+02 2.707700E-02
4 6.471000E+01 1.017180E-01
5 2.106000E+01 2.747400E-01
6 7.495000E+00 4.485640E-01
7 2.797000E+00 2.850740E-01
8 5.215000E-01 1.520400E-02
9 1.596000E-01 -3.191000E-03
S 9
1 6.665000E+03 -1.460000E-04
2 1.000000E+03 -1.154000E-03
3 2.280000E+02 -5.725000E-03
4 6.471000E+01 -2.331200E-02
5 2.106000E+01 -6.395500E-02
6 7.495000E+00 -1.499810E-01
7 2.797000E+00 -1.272620E-01
8 5.215000E-01 5.445290E-01
9 1.596000E-01 5.804960E-01
1 6
S 8
1 6665.0000000 0.0006920
2 1000.000000 0.0053290
3 228.0000000 0.0270770
4 64.7100000 0.1017180
5 21.060000 0.2747400
6 7.4950000 0.4485640
7 2.7970000 0.2850740
8 0.5215000 0.0152040
S 8
1 6665.0000000 -0.0001460
2 1000.000000 -0.0011540
3 228.0000000 -0.0057250
4 64.7100000 -0.0233120
5 21.060000 -0.0639550
6 7.4950000 -0.1499810
7 2.7970000 -0.1272620
8 0.5215000 0.5445290
S 1
1 1.596000E-01 1.000000E+00
P 4
1 9.439000E+00 3.810900E-02
2 2.002000E+00 2.094800E-01
3 5.456000E-01 5.085570E-01
4 1.517000E-01 4.688420E-01
1 0.1596000 1.0000000
P 3
1 9.4390000 0.0381090
2 2.0020000 0.2094800
3 0.5456000 0.5085570
P 1
1 1.517000E-01 1.000000E+00
1 0.1517000 1.0000000
D 1
1 5.500000E-01 1.0000000
1 0.5500000 1.0000000
2 6
S 9
1 6.665000E+03 6.920000E-04
2 1.000000E+03 5.329000E-03
3 2.280000E+02 2.707700E-02
4 6.471000E+01 1.017180E-01
5 2.106000E+01 2.747400E-01
6 7.495000E+00 4.485640E-01
7 2.797000E+00 2.850740E-01
8 5.215000E-01 1.520400E-02
9 1.596000E-01 -3.191000E-03
S 9
1 6.665000E+03 -1.460000E-04
2 1.000000E+03 -1.154000E-03
3 2.280000E+02 -5.725000E-03
4 6.471000E+01 -2.331200E-02
5 2.106000E+01 -6.395500E-02
6 7.495000E+00 -1.499810E-01
7 2.797000E+00 -1.272620E-01
8 5.215000E-01 5.445290E-01
9 1.596000E-01 5.804960E-01
S 8
1 6665.0000000 0.0006920
2 1000.000000 0.0053290
3 228.0000000 0.0270770
4 64.7100000 0.1017180
5 21.060000 0.2747400
6 7.4950000 0.4485640
7 2.7970000 0.2850740
8 0.5215000 0.0152040
S 8
1 6665.0000000 -0.0001460
2 1000.000000 -0.0011540
3 228.0000000 -0.0057250
4 64.7100000 -0.0233120
5 21.060000 -0.0639550
6 7.4950000 -0.1499810
7 2.7970000 -0.1272620
8 0.5215000 0.5445290
S 1
1 1.596000E-01 1.000000E+00
P 4
1 9.439000E+00 3.810900E-02
2 2.002000E+00 2.094800E-01
3 5.456000E-01 5.085570E-01
4 1.517000E-01 4.688420E-01
1 0.1596000 1.0000000
P 3
1 9.4390000 0.0381090
2 2.0020000 0.2094800
3 0.5456000 0.5085570
P 1
1 1.517000E-01 1.000000E+00
1 0.1517000 1.0000000
D 1
1 5.500000E-01 1.0000000
1 0.5500000 1.0000000

43
input/basis
View File

@ -1,18 +1,29 @@
1 3
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
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
S 1
1 0.1220000 1.0000000
1 0.0280500 1.0000000
P 3
1 1.5340000 0.0227840
2 0.2749000 0.1391070
3 0.0736200 0.5003750
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
S 1
1 0.1220000 1.0000000
P 1
1 0.7270000 1.0000000
1 0.0240300 1.0000000
D 1
1 0.1239000 1.0000000

20
input/dft
View File

@ -18,21 +18,23 @@
1
# Number of states in ensemble (nEns)
3
# occupation numbers of orbitals
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
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
# Ensemble weights: wEns(1),...,wEns(nEns-1)
1.00 0.00
0.20 0.30
# Ncentered ? 0 for NO
0
1
# Parameters for CC weight-dependent exchange functional
0.420431 0.069097 -0.295049
0.135075 -0.00770826 -0.028057
# choice of UCC exchange coefficient : 1 for Cx1, 2 for Cx2, 3 for Cx1*Cx2
2
# occupation numbers of orbitals nO and nO+1
1.00 0.00
1.00 0.00
0.00 0.00
1.00 0.00
1.00 0.00
1.00 1.00
# GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type
1000 0.00001 T 5 1 1

5
input/molecule
View File

@ -1,5 +1,4 @@
# nAt nEla nElb nCore nRyd
2 1 1 0 0
1 2 1 0 0
# Znuc x y z
H 0. 0. 0.
H 0. 0. 1.4
Li 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 0.7408481486
Li 0.0000000000 0.0000000000 0.0000000000

43
input/weight
View File

@ -1,18 +1,29 @@
1 3
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
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
S 1
1 0.1220000 1.0000000
1 0.0280500 1.0000000
P 3
1 1.5340000 0.0227840
2 0.2749000 0.1391070
3 0.0736200 0.5003750
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
S 1
1 0.1220000 1.0000000
P 1
1 0.7270000 1.0000000
1 0.0240300 1.0000000
D 1
1 0.1239000 1.0000000

6
scan_w.sh
View File

@ -76,11 +76,11 @@ do
echo "2" >> input/dft
echo "# occupation numbers of orbitals nO and nO+1" >> input/dft
echo " 1.00 0.00 " >> 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 " 1.00 0.00 " >> input/dft
echo " 1.00 1.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

2
src/eDFT/GOK_RKS.f90
View File

@ -37,7 +37,7 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_
double precision,intent(in) :: ENuc
double precision,intent(inout):: c(nBas,nBas)
double precision,intent(in):: occnum(nspin,2,nEns)
double precision,intent(in) :: occnum(nBas,nspin,nEns)
integer,intent(in) :: Cx_choice
! Local variables

2
src/eDFT/LIM_RKS.f90
View File

@ -31,7 +31,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,aCC_w1,aCC_w2,nGr
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ENuc
double precision,intent(in) :: occnum(nspin,2,nEns)
double precision,intent(in) :: occnum(nBas,nspin,nEns)
integer,intent(in) :: Cx_choice

2
src/eDFT/MOM_RKS.f90
View File

@ -31,7 +31,7 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, &
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ENuc
double precision,intent(in) :: occnum(nspin,2,nEns)
double precision,intent(in) :: occnum(nBas,nspin,nEns)
integer,intent(in) :: Cx_choice

15
src/eDFT/eDFT.f90
View File

@ -85,12 +85,6 @@ program eDFT
call read_molecule(nNuc,nEl(:),nO(:),nC(:),nR(:))
allocate(ZNuc(nNuc),rNuc(nNuc,ncart))
ncent = dble(nEl(1) + nEl(2))
print*, 'ncent=',ncent
print*, 'N-1/N=',(ncent-1.d0)/ncent
print*, 'N+1/N=',(ncent+1.d0)/ncent
! Read geometry
call read_geometry(nNuc,ZNuc,rNuc,ENuc)
@ -109,12 +103,11 @@ program eDFT
! DFT options
!------------------------------------------------------------------------
! Allocate ensemble weights
allocate(occnum(nspin,2,maxEns))
allocate(wEns(maxEns))
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,doNcentered,ncent,occnum,Cx_choice)
allocate(wEns(maxEns),occnum(nBas,nspin,maxEns))
call read_options(nBas,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,doNcentered,occnum,Cx_choice)
!------------------------------------------------------------------------
! Read one- and two-electron integrals

16
src/eDFT/eDFT_UKS.f90
View File

@ -30,7 +30,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ENuc
double precision,intent(in) :: occnum(nspin,2,nEns)
double precision,intent(in) :: occnum(nBas,nspin,nEns)
integer,intent(in) :: Cx_choice
! Local variables
@ -186,7 +186,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
! Compute density matrix
!------------------------------------------------------------------------
call unrestricted_density_matrix(nBas,nEns,nO(:),c(:,:,:),P(:,:,:,:),occnum)
call unrestricted_density_matrix(nBas,nEns,c(:,:,:),P(:,:,:,:),occnum(:,:,:))
! Weight-dependent density matrix
@ -337,6 +337,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
nEl(ispin) = electron_number(nGrid,weight,rhow(:,ispin))
end do
! Dump results
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
@ -362,6 +363,17 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
end if
!!!!!
do iEns=1,nEns
print*,'occnum=',occnum(1,1,iEns),occnum(2,1,iEns),occnum(1,2,iEns),occnum(2,2,iEns)
print*,'nel up and down and total=', electron_number(nGrid,weight,&
rho(:,1,iEns)),electron_number(nGrid,weight,rho(:,2,iEns)),sum(nEl(:))
end do
!!!!!
! Compute final KS energy
call print_UKS(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew)

90
src/eDFT/read_options.f90
View File

@ -1,5 +1,5 @@
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,doNcentered,ncent,occnum,Cx_choice)
subroutine read_options(nBas,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,doNcentered,occnum,Cx_choice)
! Read DFT options
@ -7,9 +7,16 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
! Local variables
integer :: I,J
integer :: iBas
integer :: iEns
integer :: iParam
character(len=1) :: answer
double precision,allocatable :: nEl(:)
! Output variables
@ -17,11 +24,12 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
integer,intent(out) :: x_rung,c_rung
character(len=12),intent(out) :: x_DFA, c_DFA
integer,intent(out) :: SGn
integer,intent(out) :: nEns, doNcentered
integer,intent(out) :: nEns
integer,intent(out) :: doNcentered
double precision,intent(out) :: wEns(maxEns)
double precision,intent(out) :: aCC_w1(3)
double precision,intent(out) :: aCC_w2(3)
double precision,intent(inout):: occnum(nspin,2,maxEns)
double precision,intent(out) :: occnum(nBas,nspin,maxEns)
integer,intent(out) :: maxSCF
double precision,intent(out) :: thresh
@ -29,13 +37,8 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
integer,intent(out) :: max_diis
integer,intent(out) :: guess_type
integer,intent(out) :: ortho_type
double precision,intent(in) :: ncent
integer,intent(out) :: Cx_choice
! Local variables
character(len=1) :: answer
! Open file with method specification
open(unit=1,file='input/dft')
@ -95,26 +98,60 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
write(*,*)'----------------------------------------------------------'
write(*,*)
! Read ensemble weights for unphysical (integer number of electrons) left or right N-centered ensembles:
! (alpha,0) or (0,alpha) in input
! read(1,*)
! read(1,*) (wEns(I),I=2,nEns)
! wEns(2) = (ncent/(ncent-1.d0))*wEns(2)
! wEns(3) = (ncent/(ncent+1.d0))*wEns(3)
! wEns(1) = 1d0 - ((ncent-1.d0)/ncent)*wEns(2) - ((ncent+1.d0)/ncent)*wEns(3) ! for N-centered
! Read occupation numbers for orbitals nO and nO+1
occnum(:,:,:) = 0d0
do iEns=1,nEns
read(1,*)
read(1,*) (occnum(iBas,1,iEns),iBas=1,nBas)
read(1,*) (occnum(iBas,2,iEns),iBas=1,nBas)
end do
do iEns=1,nEns
write(*,*)
write(*,*) '==============='
write(*,*) 'State n.',iEns
write(*,*) '==============='
write(*,*)
write(*,*) 'Spin-up occupation numbers'
write(*,*) (int(occnum(iBas,1,iEns)),iBas=1,nBas)
write(*,*) 'Spin-down occupation numbers'
write(*,*) (int(occnum(iBas,2,iEns)),iBas=1,nBas)
write(*,*)
end do
! Read ensemble weights for real physical (fractional number of electrons) ensemble (w1,w2)
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*,'nEl'
print*,nEl
read(1,*)
read(1,*) (wEns(I),I=2,nEns)
read(1,*) (wEns(iEns),iEns=2,nEns)
read(1,*)
read(1,*) doNcentered
if (doNcentered==0) then
wEns(1) = 1d0 - wEns(2) - wEns(3)
else
wEns(2) = (ncent/(ncent-1.d0))*wEns(2)
wEns(3) = (ncent/(ncent+1.d0))*wEns(3)
wEns(1) = 1d0 - ((ncent-1.d0)/ncent)*wEns(2) - ((ncent+1.d0)/ncent)*wEns(3) ! for N-centered
! wEns(1) = 1d0 - nEl(2)/nEl(1)*wEns(2) - nEl(3)/nEl(1)*wEns(3)
wEns(1) = 1d0 - wEns(2) - wEns(3)
wEns(2) = nEl(1)/nEl(2)*wEns(2)
wEns(3) = nEl(1)/nEl(3)*wEns(3)
end if
write(*,*)'----------------------------------------------------------'
write(*,*)' Ensemble weights '
write(*,*)'----------------------------------------------------------'
@ -123,8 +160,8 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
! Read parameters for weight-dependent functional
read(1,*)
read(1,*) (aCC_w1(I),I=1,3)
read(1,*) (aCC_w2(I),I=1,3)
read(1,*) (aCC_w1(iParam),iParam=1,3)
read(1,*) (aCC_w2(iParam),iParam=1,3)
! Read choice of exchange coefficient
read(1,*)
read(1,*) Cx_choice
@ -140,13 +177,6 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
write(*,*)'----------------------------------------------------------'
call matout(3,1,aCC_w2)
write(*,*)
! Read occupation numbers for orbitals nO and nO+1
read(1,*)
do J=1,nEns
read(1,*) (occnum(1,I,J),I=1,2)
read(1,*) (occnum(2,I,J),I=1,2)
end do
! Read KS options

2
src/eDFT/restricted_auxiliary_energy.f90
View File

@ -12,7 +12,7 @@ subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux,occnum)
integer,intent(in) :: nEns
integer,intent(in) :: nO
double precision,intent(in) :: eps(nBas)
double precision,intent(in) :: occnum(nspin,2,nEns)
double precision,intent(in) :: occnum(nBas,nspin,nEns)
! Local variables

2
src/eDFT/restricted_density_matrix.f90
View File

@ -12,7 +12,7 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P,occnum)
integer,intent(in) :: nEns
integer,intent(in) :: nO
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: occnum(nspin,2,nEns)
double precision,intent(in) :: occnum(nBas,nspin,nEns)
! Local variables

2
src/eDFT/restricted_individual_energy.f90
View File

@ -38,7 +38,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n
double precision,intent(in) :: J(nBas,nBas)
double precision :: Ew
double precision,intent(in) :: occnum(nspin,2,nEns)
double precision,intent(in) :: occnum(nBas,nspin,nEns)
integer,intent(in) :: Cx_choice

85
src/eDFT/unrestricted_auxiliary_energy.f90
View File

@ -12,85 +12,30 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux,occnum)
integer,intent(in) :: nEns
integer,intent(in) :: nO(nspin)
double precision,intent(in) :: eps(nBas,nspin)
double precision,intent(in) :: occnum(nspin,2,nEns)
double precision,intent(in) :: occnum(nBas,nspin,nEns)
! Local variables
integer :: iEns,ispin
integer :: iEns
integer :: ispin
integer :: p
! Output variables
double precision,intent(out) :: Eaux(nspin,nEns)
!----------------------------------------------------------
!-------------- GOK-UKS ----------------------------------
!----------------------------------------------------------
iEns = 1
do ispin=1,nspin
Eaux(ispin,iEns) = sum(eps(1:nO(ispin),ispin))
! Compute auxiliary energies for each state of the ensemble based on occupation numbers
Eaux(:,:) = 0d0
do iEns=1,nEns
do ispin=1,nspin
do p=1,nBas
Eaux(ispin,iEns) = Eaux(ispin,iEns) + occnum(p,ispin,iEns)*eps(p,ispin)
end do
end do
end do
iEns = 2
if(nO(2) > 1) then
Eaux(2,iEns) = sum(eps(1:nO(2)-1,2))
else
Eaux(2,iEns) = 0d0
end if
Eaux(2,iEns) = sum(eps(1:nO(2)-1,2)) + occnum(2,1,iEns)*eps(nO(2),2) &
+ occnum(2,2,iEns)*eps(nO(2)+1,2)
if(nO(1) > 1) then
Eaux(1,iEns) = sum(eps(1:nO(1)-1,1))
else
Eaux(1,iEns) = 0d0
end if
Eaux(1,iEns) = Eaux(1,iEns) + occnum(1,1,iEns)*eps(nO(1),1) &
+ occnum(1,2,iEns)*eps(nO(1)+1,1)
iEns = 3
if(nO(1) > 1) then
Eaux(1,iEns) = sum(eps(1:nO(1)-1,1))
else
Eaux(1,iEns) = 0d0
end if
Eaux(1,iEns) = Eaux(1,iEns) + occnum(1,1,iEns)*eps(nO(1),1) &
+ occnum(1,2,iEns)*eps(nO(1)+1,1)
if(nO(2) > 1) then
Eaux(2,iEns) = sum(eps(1:nO(2)-1,2))
else
Eaux(2,iEns) = 0d0
end if
Eaux(2,iEns) = Eaux(2,iEns) + occnum(2,1,iEns)*eps(nO(2),2) &
+ occnum(2,2,iEns)* eps(nO(2)+1,2)
!----------------------------------------------------------
!-------------- eDFT-UKS ----------------------------------
!----------------------------------------------------------
! N-electron ground state
! iEns = 1
! do ispin=1,nspin
! Eaux(ispin,iEns) = sum(eps(1:nO(ispin),ispin))
! end do
! (N-1)-electron ground state
! iEns = 2
! Eaux(2,iEns) = sum(eps(1:nO(2),2))
! if(nO(1) > 1) then
! Eaux(1,iEns) = sum(eps(1:nO(1)-1,1))
! else
! Eaux(1,iEns) = 0d0
! end if
! (N+1)-electron ground state
! iEns = 3
! Eaux(2,iEns) = sum(eps(1:nO(2)+1,2))
! Eaux(1,iEns) = sum(eps(1:nO(1),1))
end subroutine unrestricted_auxiliary_energy

90
src/eDFT/unrestricted_density_matrix.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_density_matrix(nBas,nEns,nO,c,P,occnum)
subroutine unrestricted_density_matrix(nBas,nEns,c,P,occnum)
! Calculate density matrices
@ -10,89 +10,39 @@ subroutine unrestricted_density_matrix(nBas,nEns,nO,c,P,occnum)
integer,intent(in) :: nBas
integer,intent(in) :: nEns
integer,intent(in) :: nO(nspin)
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: occnum(nspin,2,nEns)
double precision,intent(in) :: occnum(nBas,nspin,nEns)
! Local variables
integer :: ispin
integer :: iEns
integer :: q
integer :: mu,nu
! Output variables
double precision,intent(out) :: P(nBas,nBas,nspin,nEns)
!-------------------------------------------------------
!------------------------ GOK-UKS ----------------------
!-------------------------------------------------------
iEns = 1
do ispin=1,nspin
P(:,:,ispin,iEns) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
! Compute density matrix for each state of the ensemble based on occupation numbers
P(:,:,:,:) = 0d0
do iEns=1,nEns
do ispin=1,nspin
do mu=1,nBas
do nu=1,nBas
do q=1,nBas
P(mu,nu,ispin,iEns) = P(mu,nu,ispin,iEns) &
+ occnum(q,ispin,iEns)*c(mu,q,ispin)*c(nu,q,ispin)
end do
end do
end do
end do
end do
iEns = 2
if(nO(1) > 1) then
P(:,:,1,iEns) = matmul(c(:,1:nO(1)-1,1),transpose(c(:,1:nO(1)-1,1)))
else
P(:,:,1,iEns)=0.d0
end if
P(:,:,1,iEns) = P(:,:,1,iEns) + occnum(1,1,iEns)* matmul(c(:,nO(1):nO(1),1),transpose(c(:,nO(1):nO(1),1))) &
+ occnum(1,2,iEns)* matmul(c(:,nO(1)+1:nO(1)+1,1),transpose(c(:,nO(1)+1:nO(1)+1,1)))
if(nO(2) > 1) then
P(:,:,2,iEns) = matmul(c(:,1:nO(2)-1,2),transpose(c(:,1:nO(2)-1,2)))
else
P(:,:,2,iEns)=0.d0
end if
P(:,:,2,iEns) = P(:,:,2,iEns) + occnum(2,1,iEns)* matmul(c(:,1:nO(2),2),transpose(c(:,1:nO(2),2))) &
+ occnum(2,2,iEns)*matmul(c(:,1:nO(2)+1,2),transpose(c(:,1:nO(2)+1,2)))
iEns = 3
if(nO(1) > 1) then
P(:,:,1,iEns) = matmul(c(:,1:nO(1)-1,1),transpose(c(:,1:nO(1)-1,1)))
else
P(:,:,1,iEns)=0.d0
end if
P(:,:,1,iEns) = P(:,:,1,iEns) + occnum(1,1,iEns)* matmul(c(:,nO(1):nO(1),1),transpose(c(:,nO(1):nO(1),1))) &
+ occnum(1,2,iEns)* matmul(c(:,nO(1)+1:nO(1)+1,1),transpose(c(:,nO(1)+1:nO(1)+1,1)))
if(nO(2) > 1) then
P(:,:,2,iEns) = matmul(c(:,1:nO(2)-1,2),transpose(c(:,1:nO(2)-1,2)))
else
P(:,:,2,iEns)=0.d0
end if
P(:,:,2,iEns) = P(:,:,2,iEns) + occnum(2,1,iEns)* matmul(c(:,nO(2):nO(2),2),transpose(c(:,nO(2):nO(2),2))) &
+ occnum(2,2,iEns)*matmul(c(:,nO(2)+1:nO(2)+1,2),transpose(c(:,nO(2)+1:nO(2)+1,2)))
!-------------------------------------------------------------------
!--------------- For eDFT_UKS (fundamental gap)---------------------
!-------------------------------------------------------------------
! N-electron ground state
! iEns = 1
! do ispin=1,nspin
! P(:,:,ispin,iEns) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
! end do
! (N-1)-electron state: remove spin-up electrons
! iEns = 2
! P(:,:,2,iEns) = matmul(c(:,1:nO(2),2),transpose(c(:,1:nO(2),2)))
! if (nO(1) > 1) then
! P(:,:,1,iEns) = matmul(c(:,1:nO(1)-1,1),transpose(c(:,1:nO(1)-1,1)))
! else
! P(:,:,1,iEns) = 0.d0
! end if
! (N+1)-electron state: remove spin-down electrons
! iEns = 3
! P(:,:,2,iEns) = matmul(c(:,1:nO(2)+1,2),transpose(c(:,1:nO(2)+1,2)))
! P(:,:,1,iEns) = matmul(c(:,1:nO(1),1),transpose(c(:,1:nO(1),1)))
end subroutine unrestricted_density_matrix

2
src/eDFT/unrestricted_individual_energy.f90
View File

@ -41,7 +41,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
double precision,intent(in) :: FxHF(nBas,nBas,nspin)
double precision,intent(in) :: Fc(nBas,nBas,nspin)
double precision :: Ew
double precision,intent(in) :: occnum(nspin,2,nEns)
double precision,intent(in) :: occnum(nBas,nspin,nEns)
integer,intent(in) :: Cx_choice