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

scaling N-centered

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-09-29 11:47:18 +02:00
parent 9074a23700
commit 2ea01fd77e
14 changed files with 102 additions and 123 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
input/dft
View File

@ -29,8 +29,8 @@
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
# Ncentered ? 0 for NO
0
# N-centered?
T
# Parameters for CC weight-dependent exchange functional
0.445525 0.0901503 -0.286898
0.191734 -0.0364788 -0.017035

30
src/QuAcK/unrestricted_spatial_to_spin_MO_energy.f90
View File

@ -1,30 +0,0 @@
subroutine unrestricted_spatial_to_spin_MO_energy(nBas,e,nBas2,se)
! Convert MO energies from unrestricted spatial to spin orbitals
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nBas2
double precision,intent(in) :: e(nBas,nspin)
! Local variables
integer :: p
! Output variables
double precision,intent(out) :: se(nBas2)
do p=1,nBas2,2
se(p) = e(p,1)
enddo
do p=2,nBas2,2
se(p) = e(p,2)
enddo
end subroutine unrestricted_spatial_to_spin_MO_energy

15
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,Cx_choice,doNcentered,ExDD)
subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,Cx_choice,doNcentered,kappa,ExDD)
! Compute the unrestricted version of the curvature-corrected exchange ensemble derivative
@ -15,8 +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
logical,intent(in) :: doNcentered
double precision,intent(in) :: kappa(nEns)
! Local variables
@ -60,7 +60,6 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
b2 = aCC_w2(2)
c2 = aCC_w2(3)
nElw = electron_number(nGrid,weight,rhow)
@ -72,12 +71,15 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGr
w2 = wEns(3)
select case (Cx_choice)
case(1)
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
case(2)
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)))
case(3)
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))
@ -116,9 +118,6 @@ 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
if(doNcentered) ExDD(:) = kappa(:)*ExDD(:)
end subroutine UCC_lda_exchange_derivative_discontinuity

17
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,Cx_choice,doNcentered,Ex)
subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Cx_choice,doNcentered,kappa,Ex)
! Compute the unrestricted version of the curvature-corrected exchange functional
@ -16,7 +16,8 @@ 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
logical,intent(in) :: doNcentered
double precision,intent(in) :: kappa(nEns)
! Local variables
@ -61,12 +62,16 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
select case (Cx_choice)
case(1)
Cx = alpha*Fx1
case(2)
Cx = alpha*Fx2
case(3)
Cx = alpha*Fx2*Fx1
end select
nEli = electron_number(nGrid,weight,rho)
@ -87,15 +92,15 @@ 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)
if (doNcentered == 0) then
Ex = Ex - weight(iG)*dedr*r*r
else
if (doNcentered) then
Ex = Ex - weight(iG)*dedr*r*r*(nEli/nElw)
else
Ex = Ex - weight(iG)*dedr*r*r
end if
if(rI > threshold) then
if (doNcentered == 0) then
if (doNcentered) then
Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI)
else
Ex = Ex + weight(iG)*((nEli/nElw)*e_p*rI + dedr*r*rI)

42
src/eDFT/US51_lda_exchange_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered,Ex)
subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered,kappa,Ex)
! Compute the restricted version of Slater's LDA exchange individual energy
@ -11,34 +11,28 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: doNcentered
logical,intent(in) :: doNcentered
double precision,intent(in) :: kappa
! Local variables
integer :: iG
double precision :: r,rI,alpha
double precision :: e,dedr
double precision :: nEli,nElw
double precision :: Exrr,ExrI,ExrrI
! 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)
alpha = - (3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
Exrr = 0d0
ExrI = 0d0
ExrrI = 0d0
Ex = 0d0
do iG=1,nGrid
r = max(0d0,rhow(iG))
@ -49,19 +43,12 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered
e = alpha*r**(1d0/3d0)
dedr = 1d0/3d0*alpha*r**(-2d0/3d0)
if (doNcentered == 0) then
Ex = Ex - weight(iG)*dedr*r*r
else
Ex = Ex - weight(iG)*dedr*r*r*(nEli/nElw)
end if
Exrr = Exrr - weight(iG)*dedr*r*r
if(rI > threshold) then
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
ExrI = ExrI + weight(iG)*e*rI
ExrrI = ExrrI + weight(iG)*dedr*r*rI
endif
@ -69,4 +56,9 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered
enddo
Exrr = kappa*Exrr
ExrI = kappa*ExrI
Ex = Exrr + ExrI + ExrrI
end subroutine US51_lda_exchange_individual_energy

3
src/eDFT/eDFT.f90
View File

@ -50,7 +50,8 @@ program eDFT
double precision :: start_KS,end_KS,t_KS
double precision :: start_int,end_int,t_int
integer :: nEns,doNcentered
integer :: nEns
logical :: doNcentered
double precision,allocatable :: wEns(:)
integer :: maxSCF,max_diis

2
src/eDFT/eDFT_UKS.f90
View File

@ -31,7 +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
logical,intent(in) :: doNcentered
! Local variables

9
src/eDFT/read_options.f90
View File

@ -25,7 +25,7 @@ subroutine read_options(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_
character(len=12),intent(out) :: x_DFA, c_DFA
integer,intent(out) :: SGn
integer,intent(out) :: nEns
integer,intent(out) :: doNcentered
logical,intent(out) :: doNcentered
double precision,intent(out) :: wEns(maxEns)
double precision,intent(out) :: aCC_w1(3)
double precision,intent(out) :: aCC_w2(3)
@ -132,13 +132,16 @@ subroutine read_options(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_
print*,'nEl'
print*,nEl
doNcentered = .false.
read(1,*)
read(1,*) (wEns(iEns),iEns=2,nEns)
read(1,*)
read(1,*) doNcentered
read(1,*) answer
if (doNcentered == 0) then
if(answer == 'T') doNcentered = .true.
if (doNcentered) then
wEns(1) = 1d0 - wEns(2) - wEns(3)

2
src/eDFT/unrestricted_auxiliary_energy.f90
View File

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

7
src/eDFT/unrestricted_exchange_derivative_discontinuity.f90
View File

@ -1,5 +1,5 @@
subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,drhow,&
Cx_choice,doNcentered,ExDD)
Cx_choice,doNcentered,kappa,ExDD)
! Compute the exchange part of the derivative discontinuity
@ -19,7 +19,8 @@ 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
logical,intent(in) :: doNcentered
double precision,intent(in) :: kappa(nEns)
! Local variables
@ -41,7 +42,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(:),Cx_choice,doNcentered,ExDD(:))
rhow(:),Cx_choice,doNcentered,kappa,ExDD(:))
! GGA functionals

7
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,Cx_choice,doNcentered,Ex)
ERI,Pw,P,rhow,drhow,rho,drho,Cx_choice,doNcentered,kappa,Ex)
! Compute the exchange individual energy
@ -26,7 +26,8 @@ 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
logical,intent(in) :: doNcentered
double precision,intent(in) :: kappa
! Local variables
@ -51,7 +52,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,Cx_choice,doNcentered,ExLDA)
rhow,rho,Cx_choice,doNcentered,kappa,ExLDA)
Ex = ExLDA

30
src/eDFT/unrestricted_individual_energy.f90
View File

@ -43,7 +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
logical,intent(in) :: doNcentered
! Local variables
@ -68,6 +68,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
integer :: ispin,iEns,iBas
double precision,allocatable :: nEl(:)
double precision,allocatable :: kappa(:)
double precision,external :: electron_number
@ -78,12 +79,14 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
double precision,intent(out) :: Om(nEns)
allocate(nEl(nEns))
allocate(nEl(nEns),kappa(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
kappa(iEns) = nEl(iEns)/nEl(1)
end do
print*,'test1'
@ -94,10 +97,10 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
do ispin=1,nspin
do iEns=1,nEns
if (doNcentered == 0) then
ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:)))
if (doNcentered) then
ET(ispin,iEns) = kappa(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(:,:)))
ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:)))
end if
end do
end do
@ -109,10 +112,10 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
do iEns=1,nEns
do ispin=1,nspin
if (doNcentered == 0) then
EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:)))
if (doNcentered) then
EV(ispin,iEns) = kappa(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(:,:)))
EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:)))
end if
end do
end do
@ -138,9 +141,9 @@ 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
if(doNcentered) EJ(:,iEns) = kappa(iEns)*EJ(:,iEns)
end do
print*,'test4'
!------------------------------------------------------------------------
@ -166,7 +169,8 @@ 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),Cx_choice,doNcentered,Ex(ispin,iEns))
rho(:,ispin,iEns),drho(:,:,ispin,iEns),Cx_choice,doNcentered,kappa(iEns), &
Ex(ispin,iEns))
end do
end do
@ -222,7 +226,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),Cx_choice,doNcentered,ExDD(ispin,:))
rhow(:,ispin),drhow(:,:,ispin),Cx_choice,doNcentered,kappa,ExDD(ispin,:))
end do
call unrestricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,EcDD)

7
src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90
View File

@ -1,5 +1,5 @@
subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,&
Cx_choice,doNcentered,ExDD)
Cx_choice,doNcentered,kappa,ExDD)
! Compute the exchange LDA part of the derivative discontinuity
@ -18,7 +18,8 @@ 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
logical,intent(in) :: doNcentered
double precision,intent(in) :: kappa(nEns)
! Local variables
@ -38,7 +39,7 @@ 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(:),&
Cx_choice,doNcentered,ExDD(:))
Cx_choice,doNcentered,kappa,ExDD(:))
case default

10
src/eDFT/unrestricted_lda_exchange_individual_energy.f90
View File

@ -1,5 +1,5 @@
subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,&
Cx_choice,doNcentered,Ex)
Cx_choice,doNcentered,kappa,Ex)
! Compute LDA exchange energy for individual states
@ -19,7 +19,9 @@ 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
logical,intent(in) :: doNcentered
double precision,intent(in) :: kappa
! Output variables
@ -31,12 +33,12 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEn
case ('S51')
call US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered,Ex)
call US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered,kappa,Ex)
case ('CC')
call UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,&
Cx_choice,doNcentered,Ex)
Cx_choice,doNcentered,kappa,Ex)
case default