PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-05 05:33:50 +01:00
Code Releases Activity

added occupation numbers of orbitals no and no+1 in input and N-centered weights in read_options

Browse Source
This commit is contained in:
Clotilde Marut 2020-08-01 11:45:17 +02:00
parent 9f34d27502
commit e068aab977
22 changed files with 342 additions and 175 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
input/basis
View File

@ -1,25 +1,18 @@
1 10
S 4
1 528.5000000 0.0009400
2 79.3100000 0.0072140
3 18.0500000 0.0359750
4 5.0850000 0.1277820
1 3
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
S 1
1 1.6090000 1.0000000
1 0.1220000 1.0000000
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.5363000 1.0000000
S 1
1 0.1833000 1.0000000
1 0.1220000 1.0000000
P 1
1 5.9940000 1.0000000
P 1
1 1.7450000 1.0000000
P 1
1 0.5600000 1.0000000
D 1
1 4.2990000 1.0000000
D 1
1 1.2230000 1.0000000
F 1
1 2.6800000 1.0000000
1 0.7270000 1.0000000

15
input/dft
View File

@ -19,9 +19,16 @@
# Number of states in ensemble (nEns)
3
# Ensemble weights: wEns(1),...,wEns(nEns-1)
0.00 0.00
0.00 1.00
# Parameters for CC weight-dependent exchange functional
0.420431 0.069097 -0.295049
0.135075 -0.00770826 -0.028057
-0.00201219 -0.00371002 -0.00212719
-0.00117715 0.00188738 -0.000414761
# 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
500 0.00001 T 5 1 1
1000 0.00001 T 5 1 1

5
input/molecule
View File

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

37
input/weight
View File

@ -1,25 +1,18 @@
1 10
S 4
1 528.5000000 0.0009400
2 79.3100000 0.0072140
3 18.0500000 0.0359750
4 5.0850000 0.1277820
1 3
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
S 1
1 1.6090000 1.0000000
1 0.1220000 1.0000000
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.5363000 1.0000000
S 1
1 0.1833000 1.0000000
1 0.1220000 1.0000000
P 1
1 5.9940000 1.0000000
P 1
1 1.7450000 1.0000000
P 1
1 0.5600000 1.0000000
D 1
1 4.2990000 1.0000000
D 1
1 1.2230000 1.0000000
F 1
1 2.6800000 1.0000000
1 0.7270000 1.0000000

7
src/eDFT/GOK_RKS.f90
View File

@ -1,5 +1,5 @@
subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,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,c)
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,c,occnum)
! Perform restricted Kohn-Sham calculation for ensembles
@ -37,6 +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(inout),dimension(2,2,3):: occnum
! Local variables
@ -189,7 +190,7 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_
! Compute density matrix
!------------------------------------------------------------------------
call restricted_density_matrix(nBas,nEns,nO,c(:,:),P(:,:,:))
call restricted_density_matrix(nBas,nEns,nO,c(:,:),P(:,:,:),occnum)
! Weight-dependent density matrix
@ -344,6 +345,6 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_
call restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:), &
nBas,nO,nV,T(:,:),V(:,:),ERI(:,:,:,:),ENuc,eps(:),Pw(:,:),rhow(:),drhow(:,:), &
J(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,E(:),Om(:))
J(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,E(:),Om(:),occnum)
end subroutine GOK_RKS

100
src/eDFT/GOK_UKS.f90
View File

@ -1,5 +1,5 @@
subroutine 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)
nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,occnum)
! Perform unrestricted Kohn-Sham calculation for ensembles
@ -30,10 +30,12 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
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),dimension(2,2,3) :: occnum
! Local variables
integer :: xc_rung
logical :: LDA_centered = .false.
integer :: nSCF,nBasSq
integer :: n_diis
double precision :: conv
@ -127,7 +129,9 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
if(guess_type == 1) then
do ispin=1,nspin
F(:,:,ispin) = Hc(:,:)
cp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:)))
call diagonalize_matrix(nBas,cp(:,:,ispin),eps(:,ispin))
c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
end do
else if(guess_type == 2) then
@ -135,6 +139,10 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
do ispin=1,nspin
call random_number(F(:,:,ispin))
end do
else
print*,'Wrong guess option'
stop
end if
@ -171,41 +179,22 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
nSCF = nSCF + 1
! Transform Fock matrix in orthogonal basis
do ispin=1,nspin
Fp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(F(:,:,ispin),X(:,:)))
end do
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
cp(:,:,:) = Fp(:,:,:)
do ispin=1,nspin
call diagonalize_matrix(nBas,cp(:,:,ispin),eps(:,ispin))
end do
! Back-transform eigenvectors in non-orthogonal basis
do ispin=1,nspin
c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
end do
!------------------------------------------------------------------------
! Compute density matrix
!------------------------------------------------------------------------
! Compute density matrix
!------------------------------------------------------------------------
call unrestricted_density_matrix(nBas,nEns,nO(:),c(:,:,:),P(:,:,:,:))
! Weight-dependent density matrix
Pw(:,:,:) = 0d0
do iEns=1,nEns
call unrestricted_density_matrix(nBas,nEns,nO(:),c(:,:,:),P(:,:,:,:),occnum)
! Weight-dependent density matrix
Pw(:,:,:) = 0d0
do iEns=1,nEns
Pw(:,:,:) = Pw(:,:,:) + wEns(iEns)*P(:,:,:,iEns)
end do
end do
!------------------------------------------------------------------------
! Compute one-electron density and its gradient if necessary
!------------------------------------------------------------------------
!------------------------------------------------------------------------
! Compute one-electron density and its gradient if necessary
!------------------------------------------------------------------------
do ispin=1,nspin
do iEns=1,nEns
@ -248,9 +237,9 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
! Compute exchange potential
do ispin=1,nspin
call unrestricted_exchange_potential(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),aCC_w1,aCC_w2,nBas,Pw(:,:,ispin), &
ERI(:,:,:,:),AO(:,:),dAO(:,:,:),rhow(:,ispin),drhow(:,:,ispin),Fx(:,:,ispin), &
FxHF(:,:,ispin))
call unrestricted_exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas,&
Pw(:,:,ispin),ERI(:,:,:,:),AO(:,:),dAO(:,:,:),rhow(:,ispin),drhow(:,:,ispin), &
Fx(:,:,ispin),FxHF(:,:,ispin))
end do
! Compute correlation potential
@ -282,6 +271,25 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
if(minval(rcond(:)) < 1d-15) n_diis = 0
! Transform Fock matrix in orthogonal basis
do ispin=1,nspin
Fp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(F(:,:,ispin),X(:,:)))
end do
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
cp(:,:,:) = Fp(:,:,:)
do ispin=1,nspin
call diagonalize_matrix(nBas,cp(:,:,ispin),eps(:,ispin))
end do
! Back-transform eigenvectors in non-orthogonal basis
do ispin=1,nspin
c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
end do
!------------------------------------------------------------------------
! Compute KS energy
!------------------------------------------------------------------------
@ -301,19 +309,20 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
! Coulomb energy
EJ(1) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1)))
EJ(2) = trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2)))
EJ(2) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) &
+ 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1)))
EJ(3) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2)))
! Exchange energy
do ispin=1,nspin
call unrestricted_exchange_energy(x_rung,x_DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, &
Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin))
call unrestricted_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
! Correlation energy
call unrestricted_correlation_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),drhow(:,:,:),Ec)
call unrestricted_correlation_energy(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec)
! Total energy
@ -322,7 +331,7 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
! Check the grid accuracy by computing the number of electrons
do ispin=1,nspin
nEl(ispin) = electron_number(nGrid,weight(:),rhow(:,ispin))
nEl(ispin) = electron_number(nGrid,weight,rhow(:,ispin))
end do
! Dump results
@ -352,15 +361,14 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w
! Compute final KS energy
call print_UKS(nBas,nO(:),eps(:,:),c(:,:,:),ENuc,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Ew)
call print_UKS(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew)
!------------------------------------------------------------------------
! Compute individual energies from ensemble energy
!------------------------------------------------------------------------
call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, &
AO(:,:),dAO(:,:,:),nO(:),nV(:),T(:,:),V(:,:),ERI(:,:,:,:),ENuc, &
Pw(:,:,:),rhow(:,:),drhow(:,:,:),J(:,:,:),Fx(:,:,:),FxHF(:,:,:), &
Fc(:,:,:),P(:,:,:,:),rho(:,:,:),drho(:,:,:,:),E(:),Om(:))
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,nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum)
end subroutine GOK_UKS

10
src/eDFT/LIM_RKS.f90
View File

@ -1,5 +1,5 @@
subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,aCC_w1,aCC_w2,nGrid,weight, &
maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,c)
maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,c,occnum)
! Perform restricted Kohn-Sham calculation for ensembles
@ -30,6 +30,8 @@ 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),dimension(2,2,3) :: occnum
double precision,intent(out) :: c(nBas,nBas)
@ -73,7 +75,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,aCC_w1,aCC_w2,nGr
write(*,*)
call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,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(1),c)
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c,occnum)
!------------------------------------------------------------------------
! Equiensemble calculation
@ -94,7 +96,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,aCC_w1,aCC_w2,nGr
write(*,*)
call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,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(2),c)
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c,occnum)
!------------------------------------------------------------------------
! Equiensemble calculation
@ -115,7 +117,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,aCC_w1,aCC_w2,nGr
write(*,*)
! call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,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(3),c)
! max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c,occnum)
!------------------------------------------------------------------------

10
src/eDFT/MOM_RKS.f90
View File

@ -1,5 +1,5 @@
subroutine 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,nV,S,T,V,Hc,ERI,X,ENuc,c)
aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,c,occnum)
! Perform restricted Kohn-Sham calculation for ensembles
@ -30,6 +30,8 @@ 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),dimension(2,2,3) :: occnum
double precision,intent(out) :: c(nBas,nBas)
@ -73,7 +75,7 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, &
write(*,*)
call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,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(1),c)
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c,occnum)
!------------------------------------------------------------------------
! Equiensemble calculation
@ -94,7 +96,7 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, &
write(*,*)
! call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,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(2),c)
! max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c,occnum)
!------------------------------------------------------------------------
! Equiensemble calculation
@ -115,7 +117,7 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, &
write(*,*)
call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,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(3),c)
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c,occnum)
!------------------------------------------------------------------------
! MOM excitation energies

19
src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90
View File

@ -82,16 +82,21 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
! 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))
! 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)))
! Single weight-dependency
! 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.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 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

9
src/eDFT/UCC_lda_exchange_energy.f90
View File

@ -77,7 +77,14 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
w2 = wEns(3)
Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
Cx = alpha*Fx2*Fx1
! for two-weights ensemble
! Cx = alpha*Fx2*Fx1
! for left ensemble
! Cx = alpha*Fx1
! for right ensemble
Cx = alpha*Fx2
! Compute GIC-LDA exchange energy

9
src/eDFT/UCC_lda_exchange_individual_energy.f90
View File

@ -76,7 +76,14 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
w2 = wEns(3)
Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
Cx = alpha*Fx1*Fx2
! for two-weight ensembles
! Cx = alpha*Fx1*Fx2
! for left ensembles
! Cx = alpha*Fx1
! for right ensembles
Cx = alpha*Fx2
! Compute LDA exchange matrix in the AO basis

9
src/eDFT/UCC_lda_exchange_potential.f90
View File

@ -79,7 +79,14 @@ subroutine UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,
w2 = wEns(3)
Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
Cx = alpha*Fx2*Fx1
! for two-weight ensembles
! Cx = alpha*Fx2*Fx1
! for left ensemble
! Cx = alpha*Fx1
! for right ensemble
Cx = alpha*Fx2
! Compute LDA exchange matrix in the AO basis

25
src/eDFT/eDFT.f90
View File

@ -7,7 +7,7 @@ program eDFT
integer :: nNuc,nBas
integer :: nEl(nspin),nC(nspin),nO(nspin),nV(nspin),nR(nspin)
double precision :: ENuc,Ew
double precision :: ENuc,Ew,ncent
double precision,allocatable :: ZNuc(:),rNuc(:,:)
@ -59,6 +59,8 @@ program eDFT
integer :: guess_type
integer :: ortho_type
double precision,dimension(2,2,3) :: occnum
! Hello World
write(*,*)
@ -71,7 +73,7 @@ program eDFT
! Read input information
!------------------------------------------------------------------------
! Read number of atoms, number of electroes of the system
! Read number of atoms, number of electrons of the system
! nC = number of core orbitals
! nO = number of occupied orbitals
! nV = number of virtual orbitals (see below)
@ -82,6 +84,12 @@ 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)
@ -100,11 +108,12 @@ program eDFT
! DFT options
!------------------------------------------------------------------------
! Allocate ensemble weights
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)
maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,ncent,occnum)
!------------------------------------------------------------------------
! Read one- and two-electron integrals
@ -166,7 +175,7 @@ program eDFT
call cpu_time(start_KS)
call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight, &
maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO(1),nV(1), &
S,T,V,Hc,ERI,X,ENuc,Ew,c)
S,T,V,Hc,ERI,X,ENuc,Ew,c,occnum)
call cpu_time(end_KS)
t_KS = end_KS - start_KS
@ -184,7 +193,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(:,:))
S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum)
call cpu_time(end_KS)
t_KS = end_KS - start_KS
@ -202,7 +211,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(:,:))
S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum)
call cpu_time(end_KS)
t_KS = end_KS - start_KS
@ -219,7 +228,7 @@ 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)
nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,occnum)
call cpu_time(end_KS)
t_KS = end_KS - start_KS
@ -236,7 +245,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,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew)
nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,occnum)
call cpu_time(end_KS)
t_KS = end_KS - start_KS

8
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,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew)
nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,occnum)
! Perform unrestricted Kohn-Sham calculation for ensembles
@ -30,6 +30,8 @@ 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),dimension(2,2,3) :: occnum
! Local variables
@ -184,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(:,:,:,:))
call unrestricted_density_matrix(nBas,nEns,nO(:),c(:,:,:),P(:,:,:,:),occnum)
! Weight-dependent density matrix
@ -369,6 +371,6 @@ 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,nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om)
AO,dAO,nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum)
end subroutine eDFT_UKS

26
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)
maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,ncent,occnum)
! Read DFT options
@ -9,7 +9,7 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
! Local variables
integer :: I
integer :: I,J
! Output variables
@ -21,6 +21,7 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
double precision,intent(out) :: wEns(maxEns)
double precision,intent(out) :: aCC_w1(3)
double precision,intent(out) :: aCC_w2(3)
double precision,intent(out),dimension(2,2,3) :: occnum
integer,intent(out) :: maxSCF
double precision,intent(out) :: thresh
@ -28,6 +29,7 @@ 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
! Local variables
@ -92,10 +94,18 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
write(*,*)'----------------------------------------------------------'
write(*,*)
! Read ensemble weights
! 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 ensemble weights for real physical (fractional number of electrons) ensemble (w1,w2)
read(1,*)
read(1,*) (wEns(I),I=2,nEns)
wEns(1) = 1d0 - sum(wEns)
wEns(1) = 1d0 - wEns(2) - wEns(3)
write(*,*)'----------------------------------------------------------'
write(*,*)' Ensemble weights '
@ -119,6 +129,14 @@ 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(*,*)
! allocate(occnum(2,2,nEns))
! Read occupation numbers for orbitals nO and nO+1
read(1,*)
do J=1,3
read(1,*) (occnum(1,I,J),I=1,2)
read(1,*) (occnum(2,I,J),I=1,2)
end do
! Read KS options

8
src/eDFT/restricted_auxiliary_energy.f90
View File

@ -1,4 +1,4 @@
subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux,occnum)
! Compute the auxiliary KS energies
@ -12,6 +12,8 @@ subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
integer,intent(in) :: nEns
integer,intent(in) :: nO
double precision,intent(in) :: eps(nBas)
double precision,intent(in),dimension(2,2,3) :: occnum
! Local variables
@ -36,8 +38,8 @@ subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
else
Eaux(iEns) = 0d0
end if
Eaux(iEns) = Eaux(iEns) + sum(occnum(:,1,iEns))*eps(nO) + sum(occnum(:,2,iEns))*eps(nO+1)
Eaux(iEns) = Eaux(iEns) + eps(nO) + eps(nO+2)
! Doubly-excited state density matrix
@ -48,7 +50,7 @@ subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
else
Eaux(iEns) = 0d0
end if
Eaux(iEns) = Eaux(iEns) + sum(occnum(:,1,iEns))*eps(nO) + sum(occnum(:,2,iEns))*eps(nO+1)
Eaux(iEns) = Eaux(iEns) + 2d0*eps(nO+1)
end subroutine restricted_auxiliary_energy

15
src/eDFT/restricted_density_matrix.f90
View File

@ -1,4 +1,4 @@
subroutine restricted_density_matrix(nBas,nEns,nO,c,P)
subroutine restricted_density_matrix(nBas,nEns,nO,c,P,occnum)
! Calculate density matrices
@ -12,6 +12,8 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P)
integer,intent(in) :: nEns
integer,intent(in) :: nO
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in),dimension(2,2,3) :: occnum
! Local variables
@ -25,7 +27,7 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P)
iEns = 1
P(:,:,iEns) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
P(:,:,iEns) = 2.d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! Doubly-excited state density matrix
@ -36,9 +38,8 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P)
else
P(:,:,iEns) = 0d0
end if
P(:,:,iEns) = P(:,:,iEns) + 1d0*matmul(c(:,nO :nO ),transpose(c(:,nO :nO ))) &
+ 1d0*matmul(c(:,nO+2:nO+2),transpose(c(:,nO+2:nO+2)))
P(:,:,iEns) = P(:,:,iEns) + sum(occnum(:,1,iEns))*matmul(c(:,nO:nO),transpose(c(:,nO:nO))) &
+ sum(occnum(:,2,iEns))*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1)))
! Doubly-excited state density matrix
@ -49,7 +50,7 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P)
else
P(:,:,iEns) = 0d0
end if
P(:,:,iEns) = P(:,:,iEns) + 2d0*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1)))
P(:,:,iEns) = P(:,:,iEns) + sum(occnum(:,1,iEns))*matmul(c(:,nO:nO),transpose(c(:,nO:nO))) &
+ sum(occnum(:,2,iEns))*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1)))
end subroutine restricted_density_matrix

6
src/eDFT/restricted_individual_energy.f90
View File

@ -1,5 +1,5 @@
subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,P,rho,drho,Ew,E,Om)
nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,P,rho,drho,Ew,E,Om,occnum)
! Compute individual energies as well as excitation energies
@ -38,6 +38,8 @@ 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),dimension(2,2,3) :: occnum
! Local variables
@ -108,7 +110,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n
! Compute auxiliary energies
!------------------------------------------------------------------------
call restricted_auxiliary_energy(nBas,nEns,nO,eps(:),Eaux(:))
call restricted_auxiliary_energy(nBas,nEns,nO,eps(:),Eaux(:),occnum)
!------------------------------------------------------------------------
! Compute derivative discontinuities

72
src/eDFT/unrestricted_auxiliary_energy.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
subroutine unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux,occnum)
! Compute the auxiliary KS energies
@ -12,6 +12,8 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
integer,intent(in) :: nEns
integer,intent(in) :: nO(nspin)
double precision,intent(in) :: eps(nBas,nspin)
double precision,intent(in),dimension(2,2,3) :: occnum
! Local variables
@ -21,30 +23,74 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
double precision,intent(out) :: Eaux(nspin,nEns)
! N-electron ground state
!----------------------------------------------------------
!-------------- GOK-UKS ----------------------------------
!----------------------------------------------------------
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
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
! iEns = 3
Eaux(2,iEns) = sum(eps(1:nO(2)+1,2))
Eaux(1,iEns) = sum(eps(1:nO(1),1))
! Eaux(2,iEns) = sum(eps(1:nO(2)+1,2))
! Eaux(1,iEns) = sum(eps(1:nO(1),1))
end subroutine unrestricted_auxiliary_energy

79
src/eDFT/unrestricted_density_matrix.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_density_matrix(nBas,nEns,nO,c,P)
subroutine unrestricted_density_matrix(nBas,nEns,nO,c,P,occnum)
! Calculate density matrices
@ -12,6 +12,8 @@ subroutine unrestricted_density_matrix(nBas,nEns,nO,c,P)
integer,intent(in) :: nEns
integer,intent(in) :: nO(nspin)
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in),dimension(2,2,3) :: occnum
! Local variables
@ -22,26 +24,75 @@ subroutine unrestricted_density_matrix(nBas,nEns,nO,c,P)
double precision,intent(out) :: P(nBas,nBas,nspin,nEns)
! N-electron ground state
!-------------------------------------------------------
!------------------------ GOK-UKS ----------------------
!-------------------------------------------------------
iEns = 1
do ispin=1,nspin
P(:,:,ispin,iEns) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
P(:,:,ispin,iEns) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
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
! 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)))
! 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

6
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,aCC_w1,aCC_w2,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)
nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum)
! Compute unrestricted individual energies as well as excitation energies
@ -41,6 +41,8 @@ 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),dimension(2,2,3) :: occnum
! Local variables
@ -182,7 +184,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
! Compute auxiliary energies
!------------------------------------------------------------------------
call unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
call unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux,occnum)
!------------------------------------------------------------------------
! Compute derivative discontinuities