LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-08-25 05:51:46 +02:00
Code Issues Releases Wiki Activity

NOFT (#73)

Browse Source 
* Fixed travis

* NOFT by T2

* Bug with n_tasks_max
This commit is contained in:
Pierre-Francois Loos 2018-09-19 10:08:07 +02:00 committed by Anthony Scemama
parent 9a77f8d342
commit e89d623be5
10 changed files with 1056 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
plugins/NOFT/.gitignore vendored Normal file
View File

@ -0,0 +1,5 @@
IRPF90_temp/
IRPF90_man/
irpf90.make
irpf90_entities
tags

19
plugins/NOFT/EZFIO.cfg Normal file
View File

@ -0,0 +1,19 @@
[do_JK_functionals]
type: logical
doc: Compute energies for JK-only functionals
interface: ezfio,provider,ocaml
default: True
[do_JKL_functionals]
type: logical
doc: Compute energies for JKL-only functionals (PNOFs)
interface: ezfio,provider,ocaml
default: True
[do_PT2_NOFT]
type: logical
doc: Compute PT2 correction for NOFT
interface: ezfio,provider,ocaml
default: False

1
plugins/NOFT/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@
Hartree_Fock Determinants

72
plugins/NOFT/NOFT.irp.f Normal file
View File

@ -0,0 +1,72 @@
program NOFT
implicit none
BEGIN_DOC
! Natural orbital functional theory module
END_DOC
PROVIDE mo_bielec_integrals_in_map
integer :: i,j
integer :: nMO,FL
double precision :: ET,EV
double precision :: integral,get_mo_bielec_integral
double precision,allocatable :: n(:)
print*, ''
print*, '*******************************'
print*, '*** NOFT functionals ***'
print*, '*******************************'
print*, ''
print*, 'SD = single determinant'
print*, 'MBB = Muller, Buijse and Baerends'
print*, 'POWER = Cioslowski and Pernal'
print*, 'BCC2 = Gritsenko and coworkers'
print*, 'CA = Csanyi and Arias'
print*, 'CGA = Csanyi, Goedecker and Arias'
print*, 'GU = Goedecker and Umrigar'
print*, 'ML = Marques and Lathiotakis'
print*, 'MLSIC = ML with self-interaction correction'
print*, 'PNOF2 = Piris natural orbital functional 2 (bug)'
print*, 'PNOF3 = Piris natural orbital functional 3'
print*, 'PNOF4 = Piris natural orbital functional 4'
print*, 'PNOF5 = Piris natural orbital functional 5 (NYI)'
print*, 'PNOF6x = Piris natural orbital functional 6 (x = d, u, h)'
print*, 'PNOF7 = Piris natural orbital functional 7 (NYI)'
print*, ''
print*, '*******************************'
print*, ''
print*, '*******************************'
print*, '*** NOFT energies ***'
print*, '*******************************'
print*, ''
! Occupation numbers
nMO = mo_tot_num
FL = elec_num/2
allocate(n(nMO))
n(1:nMO) = 0.5d0*mo_occ(1:mo_tot_num)
! Compute core energies
call NOFT_core(nMO,ET,EV,n)
! JK-only functionals
if(do_JK_functionals) call NOFT_JKfunc(nMO,FL,ET,EV,n)
! JKL-only functionals
if(do_JKL_functionals) call NOFT_JKLfunc(nMO,FL,ET,EV,n)
! PT2-NOFT correction
if(do_PT2_NOFT) call NOFT_JKLfunc(nMO,FL,n)
! End
print*, '*******************************'
print*, ''
end

484
plugins/NOFT/NOFT_JKLfunc.irp.f Normal file
View File

@ -0,0 +1,484 @@
subroutine NOFT_JKLfunc(nMO,FL,ET,EV,n)
! JKL-only functionals for NOFT
END_DOC
PROVIDE mo_bielec_integrals_in_map
! Input variables
integer,intent(in) :: nMO,FL
double precision,intent(in) :: ET,EV
double precision,intent(in) :: n(nMO)
! Local variables
integer :: i,j
double precision :: EJ_SD,EK_SD
double precision :: EJ_PNOF2,EJ_PNOF3,EJ_PNOF4,EJ_PNOF5,EJ_PNOF6d,EJ_PNOF6u,EJ_PNOF6h,EJ_PNOF7
double precision :: EK_PNOF2,EK_PNOF3,EK_PNOF4,EK_PNOF5,EK_PNOF6d,EK_PNOF6u,EK_PNOF6h,EK_PNOF7
double precision :: EL_PNOF2,EL_PNOF3,EL_PNOF4,EL_PNOF5,EL_PNOF6d,EL_PNOF6u,EL_PNOF6h,EL_PNOF7
double precision :: E_PNOF2,E_PNOF3,E_PNOF4,E_PNOF5,E_PNOF6d,E_PNOF6u,E_PNOF6h,E_PNOF7
double precision :: get_mo_bielec_integral
double precision :: SF,Sd,Su,Sh,Delta_ij,T_ij,Pi_ij
double precision,allocatable :: h(:),kappa(:),gam(:),Jint(:,:),Kint(:,:),Lint(:,:)
! memory allocation
allocate(h(nMO),kappa(nMO),gam(nMO),Jint(nMO,nMO),Kint(nMO,nMO),Lint(nMO,nMO))
! Useful quantities
h(1:nMO) = 1d0 - n(1:nMO)
SF = 0d0
do i=1,FL
SF = SF + h(i)
enddo
! Useful quantities for PNOF6
do i=1,FL
kappa(i) = h(i)*exp(-SF)
enddo
do i=FL+1,nMO
kappa(i) = n(i)*exp(-SF)
enddo
gam(:) = 0d0
do i=1,nMO
do j=i,FL
gam(i) = gam(i) + kappa(j)
enddo
gam(i) = n(i)*h(i) + kappa(i)*kappa(i) - kappa(i)*gam(i)
enddo
Sd = 0d0
do i=1,FL
Sd = Sd + gam(i)
enddo
Su = 0d0
do i=FL+1,nMO
Su = Su + gam(i)
enddo
Sh = 0.5d0*(Sd + Su)
! Coulomb, exchange and time-inversion integrals
do i=1,nMO
do j=1,nMO
Jint(i,j) = get_mo_bielec_integral(i,j,i,j,mo_integrals_map)
Kint(i,j) = get_mo_bielec_integral(i,j,j,i,mo_integrals_map)
Lint(i,j) = get_mo_bielec_integral(i,i,j,j,mo_integrals_map)
enddo
enddo
!****************************************
!*** Coulomb and exchange parts of SD ***
!****************************************
EJ_SD = +2d0*dot_product(n,matmul(Jint,n))
EK_SD = -1d0*dot_product(n,matmul(Kint,n))
! *************
! *** PNOF2 ***
! *************
EJ_PNOF2 = 0d0
EK_PNOF2 = 0d0
EL_PNOF2 = 0d0
do i=1,nMO
do j=1,nMO
if(i == j) then
Delta_ij = n(i)*n(j)
T_ij = 0d0
Pi_ij = sqrt(n(i)*n(j))
elseif(i <= FL .and. j <= FL) then
Delta_ij = h(i)*h(j)
T_ij = n(i)*n(j) - Delta_ij
Pi_ij = sqrt(n(i)*n(j)) + sqrt(h(i)*h(j)) + T_ij
elseif(i <= FL .and. j > FL) then
Delta_ij = n(j)*h(i)*(1d0 - SF)/SF
T_ij = n(i)*n(j) - Delta_ij
Pi_ij = sqrt(n(i)*n(j)) - sqrt(n(j)*h(i)) + T_ij
elseif(i > FL .and. j <= FL) then
Delta_ij = n(i)*h(j)*(1d0 - SF)/SF
T_ij = n(i)*n(j) - Delta_ij
Pi_ij = sqrt(n(i)*n(j)) - sqrt(n(i)*h(j)) + T_ij
elseif(i > FL .and. j > FL) then
Delta_ij = n(i)*n(j)
T_ij = n(i)*n(j) - Delta_ij
Pi_ij = T_ij
else
Delta_ij = 0d0
T_ij = 0d0
Pi_ij = 0d0
endif
EJ_PNOF2 = EJ_PNOF2 - 2d0*Delta_ij*Jint(i,j)
EK_PNOF2 = EK_PNOF2 + Delta_ij*Kint(i,j)
EL_PNOF2 = EL_PNOF2 + Pi_ij*Lint(i,j)
enddo
enddo
! *************
! *** PNOF3 ***
! *************
EJ_PNOF3 = 0d0
EK_PNOF3 = 0d0
EL_PNOF3 = 0d0
do i=1,nMO
do j=1,nMO
if(i == j) then
Delta_ij = n(i)*n(j)
Pi_ij = sqrt(n(i)*n(j))
elseif(i <= FL .and. j <= FL) then
Delta_ij = h(i)*h(j)
Pi_ij = n(i)*n(j) - sqrt(n(i)*n(j))
elseif(i <= FL .and. j > FL) then
Delta_ij = n(j)*h(i)*(1d0 - SF)/SF
Pi_ij = n(i)*n(j) - sqrt(n(i)*n(j)) - sqrt(n(j)*h(i))
elseif(i > FL .and. j <= FL) then
Delta_ij = n(i)*h(j)*(1d0 - SF)/SF
Pi_ij = n(i)*n(j) - sqrt(n(i)*n(j)) - sqrt(n(i)*h(j))
elseif(i > FL .and. j > FL) then
Delta_ij = n(i)*n(j)
Pi_ij = n(i)*n(j) + sqrt(n(i)*n(j))
else
Delta_ij = 0d0
Pi_ij = 0d0
endif
EJ_PNOF3 = EJ_PNOF3 - Delta_ij*Jint(i,j)
EK_PNOF3 = EK_PNOF3
EL_PNOF3 = EL_PNOF3 + Pi_ij*Lint(i,j)
enddo
enddo
! *************
! *** PNOF4 ***
! *************
EJ_PNOF4 = 0d0
EK_PNOF4 = 0d0
EL_PNOF4 = 0d0
do i=1,nMO
do j=1,nMO
if(i == j) then
Delta_ij = n(i)*n(j)
Pi_ij = sqrt(n(i)*n(j))
elseif(i <= FL .and. j <= FL) then
Delta_ij = h(i)*h(j)
Pi_ij = - sqrt(h(i)*h(j))
elseif(i <= FL .and. j > FL) then
Delta_ij = n(j)*h(i)*(1d0 - SF)/SF
Pi_ij = - sqrt( (h(i)*n(j)/SF) * (n(i)-n(j)+h(i)*n(j)/SF))
elseif(i > FL .and. j <= FL) then
Delta_ij = n(i)*h(j)*(1d0 - SF)/SF
Pi_ij = - sqrt( (h(j)*n(i)/SF) * (n(j)-n(i)+h(j)*n(i)/SF))
elseif(i > FL .and. j >= FL) then
Delta_ij = n(i)*n(j)
Pi_ij = sqrt(n(i)*n(j))
else
Delta_ij = 0d0
Pi_ij = 0d0
endif
EJ_PNOF4 = EJ_PNOF4 - 2d0*Delta_ij*Jint(i,j)
EK_PNOF4 = EK_PNOF4 + Delta_ij*Kint(i,j)
EL_PNOF4 = EL_PNOF4 + Pi_ij*Lint(i,j)
enddo
enddo
! **************
! *** PNOF6d ***
! **************
EJ_PNOF6d = 0d0
EK_PNOF6d = 0d0
EL_PNOF6d = 0d0
do i=1,nMO
do j=1,nMO
if(i == j) then
Delta_ij = n(i)*n(j)
Pi_ij = sqrt(n(i)*n(j))
elseif(i <= FL .and. j <= FL) then
Delta_ij = h(i)*h(j)*exp(-2d0*SF)
Pi_ij = - sqrt(h(i)*h(j))*exp(-SF)
elseif(i <= FL .and. j > FL) then
Delta_ij = gam(i)*gam(j)/Sd
Pi_ij = - sqrt( (n(i)*h(j) + gam(i)*gam(j)/Sd) * (n(j)*h(i) + gam(i)*gam(j)/Sd))
elseif(i > FL .and. j <= FL) then
Delta_ij = gam(i)*gam(j)/Sd
Pi_ij = - sqrt( (n(i)*h(j) + gam(i)*gam(j)/Sd) * (n(j)*h(i) + gam(i)*gam(j)/Sd))
elseif(i > FL .and. j >= FL) then
Delta_ij = n(i)*n(j)*exp(-2d0*SF)
Pi_ij = sqrt(n(i)*n(j))*exp(-SF)
else
Delta_ij = 0d0
Pi_ij = 0d0
endif
EJ_PNOF6d = EJ_PNOF6d - 2d0*Delta_ij*Jint(i,j)
EK_PNOF6d = EK_PNOF6d + Delta_ij*Kint(i,j)
EL_PNOF6d = EL_PNOF6d + Pi_ij*Lint(i,j)
enddo
enddo
! **************
! *** PNOF6u ***
! **************
EJ_PNOF6u = 0d0
EK_PNOF6u = 0d0
EL_PNOF6u = 0d0
do i=1,nMO
do j=1,nMO
if(i == j) then
Delta_ij = n(i)*n(j)
Pi_ij = sqrt(n(i)*n(j))
elseif(i <= FL .and. j <= FL) then
Delta_ij = h(i)*h(j)*exp(-2d0*SF)
Pi_ij = - sqrt(h(i)*h(j))*exp(-SF)
elseif(i <= FL .and. j > FL) then
Delta_ij = gam(i)*gam(j)/Su
Pi_ij = - sqrt( (n(i)*h(j) + gam(i)*gam(j)/Su) * (n(j)*h(i) + gam(i)*gam(j)/Su))
elseif(i > FL .and. j <= FL) then
Delta_ij = gam(i)*gam(j)/Su
Pi_ij = - sqrt( (n(i)*h(j) + gam(i)*gam(j)/Su) * (n(j)*h(i) + gam(i)*gam(j)/Su))
elseif(i > FL .and. j >= FL) then
Delta_ij = n(i)*n(j)*exp(-2d0*SF)
Pi_ij = sqrt(n(i)*n(j))*exp(-SF)
else
Delta_ij = 0d0
Pi_ij = 0d0
endif
EJ_PNOF6u = EJ_PNOF6u - 2d0*Delta_ij*Jint(i,j)
EK_PNOF6u = EK_PNOF6u + Delta_ij*Kint(i,j)
EL_PNOF6u = EL_PNOF6u + Pi_ij*Lint(i,j)
enddo
enddo
! **************
! *** PNOF6h ***
! **************
EJ_PNOF6h = 0d0
EK_PNOF6h = 0d0
EL_PNOF6h = 0d0
do i=1,nMO
do j=1,nMO
if(i == j) then
Delta_ij = n(i)*n(j)
Pi_ij = sqrt(n(i)*n(j))
elseif(i <= FL .and. j <= FL) then
Delta_ij = h(i)*h(j)*exp(-2d0*SF)
Pi_ij = - sqrt(h(i)*h(j))*exp(-SF)
elseif(i <= FL .and. j > FL) then
Delta_ij = gam(i)*gam(j)/Sh
Pi_ij = - sqrt( (n(i)*h(j) + gam(i)*gam(j)/Sh) * (n(j)*h(i) + gam(i)*gam(j)/Sh))
elseif(i > FL .and. j <= FL) then
Delta_ij = gam(i)*gam(j)/Sh
Pi_ij = - sqrt( (n(i)*h(j) + gam(i)*gam(j)/Sh) * (n(j)*h(i) + gam(i)*gam(j)/Sh))
elseif(i > FL .and. j >= FL) then
Delta_ij = n(i)*n(j)*exp(-2d0*SF)
Pi_ij = sqrt(n(i)*n(j))*exp(-SF)
else
Delta_ij = 0d0
Pi_ij = 0d0
endif
EJ_PNOF6h = EJ_PNOF6h - 2d0*Delta_ij*Jint(i,j)
EK_PNOF6h = EK_PNOF6h + Delta_ij*Kint(i,j)
EL_PNOF6h = EL_PNOF6h + Pi_ij*Lint(i,j)
enddo
enddo
! Add the SD part
EJ_PNOF2 = EJ_SD + EJ_PNOF2
EJ_PNOF3 = EJ_SD + EJ_PNOF3
EJ_PNOF4 = EJ_SD + EJ_PNOF4
! EJ_PNOF5 = EJ_SD + EJ_PNOF5
EJ_PNOF6d = EJ_SD + EJ_PNOF6d
EJ_PNOF6u = EJ_SD + EJ_PNOF6u
EJ_PNOF6h = EJ_SD + EJ_PNOF6h
! EJ_PNOF7 = EJ_SD + EJ_PNOF7
EK_PNOF2 = EK_SD + EK_PNOF2
EK_PNOF3 = EK_SD + EK_PNOF3
EK_PNOF4 = EK_SD + EK_PNOF4
! EK_PNOF5 = EK_SD + EK_PNOF5
EK_PNOF6d = EK_SD + EK_PNOF6d
EK_PNOF6u = EK_SD + EK_PNOF6u
EK_PNOF6h = EK_SD + EK_PNOF6h
! EK_PNOF7 = EK_SD + EK_PNOF7
! Compute total energies
E_PNOF2 = ET + EV + EJ_PNOF2 + EK_PNOF2 + EL_PNOF2
E_PNOF3 = ET + EV + EJ_PNOF3 + EK_PNOF3 + EL_PNOF3
E_PNOF4 = ET + EV + EJ_PNOF4 + EK_PNOF4 + EL_PNOF4
! E_PNOF5 = ET + EV + EJ_PNOF5 + EK_PNOF5 + EL_PNOF5
E_PNOF6d = ET + EV + EJ_PNOF6d + EK_PNOF6d + EL_PNOF6d
E_PNOF6u = ET + EV + EJ_PNOF6u + EK_PNOF6u + EL_PNOF6u
E_PNOF6h = ET + EV + EJ_PNOF6h + EK_PNOF6h + EL_PNOF6h
! E_PNOF7 = ET + EV + EJ_PNOF7 + EK_PNOF7 + EL_PNOF7
! Dump energies
print*, '*******************************'
print*, '*** JKL NOFT functionals ***'
print*, '*******************************'
print*, ''
print*, '*** Coulomb energies ***'
print*, 'Coulomb PNOF2 energy = ',EJ_PNOF2
print*, 'Coulomb PNOF3 energy = ',EJ_PNOF3
print*, 'Coulomb PNOF4 energy = ',EJ_PNOF4
! print*, 'Coulomb PNOF5 energy = ',EJ_PNOF5
print*, 'Coulomb PNOF6d energy = ',EJ_PNOF6d
print*, 'Coulomb PNOF6u energy = ',EJ_PNOF6u
print*, 'Coulomb PNOF6h energy = ',EJ_PNOF6h
! print*, 'Coulomb PNOF7 energy = ',EJ_PNOF7
print*, ''
print*, '*** Exchange energies ***'
print*, 'Exchange PNOF2 energy = ',EK_PNOF2
print*, 'Exchange PNOF3 energy = ',EK_PNOF3
print*, 'Exchange PNOF4 energy = ',EK_PNOF4
! print*, 'Exchange PNOF5 energy = ',EK_PNOF5
print*, 'Exchange PNOF6d energy = ',EK_PNOF6d
print*, 'Exchange PNOF6u energy = ',EK_PNOF6u
print*, 'Exchange PNOF6h energy = ',EK_PNOF6h
! print*, 'Exchange PNOF7 energy = ',EK_PNOF7
print*, ''
print*, '*** Time-inversion energies ***'
print*, 'Time-inversion PNOF2 energy = ',EL_PNOF2
print*, 'Time-inversion PNOF3 energy = ',EL_PNOF3
print*, 'Time-inversion PNOF4 energy = ',EL_PNOF4
! print*, 'Time-inversion PNOF5 energy = ',EL_PNOF5
print*, 'Time-inversion PNOF6d energy = ',EL_PNOF6d
print*, 'Time-inversion PNOF6u energy = ',EL_PNOF6u
print*, 'Time-inversion PNOF6h energy = ',EL_PNOF6h
! print*, 'Time-inversion PNOF7 energy = ',EL_PNOF7
print*, ''
print*, '*** Two-electron energies ***'
print*, 'J+K+L PNOF2 energy = ',EJ_PNOF2 + EK_PNOF2 + EL_PNOF2
print*, 'J+K+L PNOF3 energy = ',EJ_PNOF3 + EK_PNOF3 + EL_PNOF3
print*, 'J+K+L PNOF4 energy = ',EJ_PNOF4 + EK_PNOF4 + EL_PNOF4
! print*, 'J+K+L PNOF5 energy = ',EJ_PNOF5 + EK_PNOF5 + EL_PNOF5
print*, 'J+K+L PNOF6d energy = ',EJ_PNOF6d + EK_PNOF6d + EL_PNOF6d
print*, 'J+K+L PNOF6u energy = ',EJ_PNOF6u + EK_PNOF6u + EL_PNOF6u
print*, 'J+K+L PNOF6h energy = ',EJ_PNOF6h + EK_PNOF6h + EL_PNOF6h
! print*, 'J+K+L PNOF7 energy = ',EJ_PNOF7 + EK_PNOF7 + EL_PNOF7
print*, ''
print*, '*** Total energies ***'
print*, 'Total PNOF2 energy = ',E_PNOF2
print*, 'Total PNOF3 energy = ',E_PNOF3
print*, 'Total PNOF4 energy = ',E_PNOF4
! print*, 'Total PNOF5 energy = ',E_PNOF5
print*, 'Total PNOF6d energy = ',E_PNOF6d
print*, 'Total PNOF6u energy = ',E_PNOF6u
print*, 'Total PNOF6h energy = ',E_PNOF6h
! print*, 'Total PNOF7 energy = ',E_PNOF7
print*, ''
end subroutine NOFT_JKLfunc

260
plugins/NOFT/NOFT_JKfunc.irp.f Normal file
View File

@ -0,0 +1,260 @@
subroutine NOFT_JKfunc(nMO,FL,ET,EV,n)
implicit none
BEGIN_DOC
! JK-only functionals for NOFT
END_DOC
PROVIDE mo_bielec_integrals_in_map
! Input variables
integer,intent(in) :: nMO,FL
double precision,intent(in) :: ET,EV
double precision,intent(in) :: n(nMO)
! Local variables
integer :: i,j
double precision :: EJ_SD
double precision :: EK_SD,EK_MBB,EK_POWER,EK_BBC1,EK_BBC2,EK_CA,EK_CGA,EK_GU,EK_ML,EK_MLSIC
double precision :: E_SD,E_MBB,E_POWER,E_BBC1,E_BBC2,E_CA,E_CGA,E_GU,E_ML,E_MLSIC
double precision :: alpha,a0,a1,b1
double precision :: f_ij,get_mo_bielec_integral
double precision,allocatable :: Jint(:,:),Kint(:,:)
! Memory allocation
allocate(Jint(nMO,nMO),Kint(nMO,nMO))
! Coulomb, exchange and time-inversion integrals
do i=1,nMO
do j=1,nMO
Jint(i,j) = get_mo_bielec_integral(i,j,i,j,mo_integrals_map)
Kint(i,j) = get_mo_bielec_integral(i,j,j,i,mo_integrals_map)
enddo
enddo
! Compute SD Coulomb energy
EJ_SD = 2d0*dot_product(n,matmul(Jint,n))
! Compute SD exchange energy
EK_SD = dot_product(n,matmul(Kint,n))
! Compute MBB exchange energy
EK_MBB = 0d0
do i=1,nMO
do j=1,nMO
EK_MBB = EK_MBB + sqrt(n(i)*n(j))*Kint(i,j)
enddo
enddo
! Compute BBC1 exchange energy
EK_BBC1 = 0d0
do i=1,nMO
do j=1,nMO
if(i /= j .and. i > FL .and. j > FL) then
f_ij = - sqrt(n(i)*n(j))
else
f_ij = sqrt(n(i)*n(j))
endif
EK_BBC1 = EK_BBC1 + f_ij*Kint(i,j)
enddo
enddo
! Compute BBC2 exchange energy
EK_BBC2 = 0d0
do i=1,nMO
do j=1,i-1
if(i > FL .and. j > FL) then
f_ij = - sqrt(n(i)*n(j))
elseif(i <= FL .and. j <= FL) then
f_ij = n(i)*n(j)
else
f_ij = sqrt(n(i)*n(j))
endif
EK_BBC2 = EK_BBC2 + f_ij*Kint(i,j)
enddo
EK_BBC2 = EK_BBC2 + n(i)*Kint(i,i)
do j=i+1,nMO
if(i > FL .and. j > FL) then
f_ij = - sqrt(n(i)*n(j))
elseif(i <= FL .and. j <= FL) then
f_ij = n(i)*n(j)
else
f_ij = sqrt(n(i)*n(j))
endif
EK_BBC2 = EK_BBC2 + f_ij*Kint(i,j)
enddo
enddo
! Compute CA exchange energy
EK_CA = 0d0
do i=1,nMO
do j=1,nMO
EK_CA = EK_CA + (sqrt(n(i)*n(j)*(1d0 - n(i))*(1d0 - n(j))) + n(i)*n(j))*Kint(i,j)
enddo
enddo
! Compute CGA exchange energy
EK_CGA = 0d0
do i=1,nMO
do j=1,nMO
EK_CGA = EK_CGA + 0.5d0*(sqrt(n(i)*n(j)*(2d0 - n(i))*(2d0 - n(j))) + n(i)*n(j))*Kint(i,j)
enddo
enddo
! Compute ML exchange energy
EK_ML = 0d0
a0 = 126.3101d0
a1 = 2213.33d0
b1 = 2338.64d0
do i=1,nMO
do j=1,nMO
EK_ML = EK_ML + n(i)*n(j)*(a0 + a1*n(i)*n(j))/(1d0 + b1*n(i)*n(j))*Kint(i,j)
enddo
enddo
! Compute MLSIC exchange energy
EK_MLSIC = 0d0
a0 = 1298.78d0
a1 = 35114.4d0
b1 = 36412.2d0
do i=1,nMO
do j=1,i-1
EK_MLSIC = EK_MLSIC + n(i)*n(j)*(a0 + a1*n(i)*n(j))/(1d0 + b1*n(i)*n(j))*Kint(i,j)
enddo
EK_MLSIC = EK_MLSIC + n(i)*n(i)*Kint(i,i)
do j=i+1,nMO
EK_MLSIC = EK_MLSIC + n(i)*n(j)*(a0 + a1*n(i)*n(j))/(1d0 + b1*n(i)*n(j))*Kint(i,j)
enddo
enddo
! Compute GU exchange energy
EK_GU = 0d0
do i=1,nMO
do j=1,i-1
EK_GU = EK_GU + sqrt(n(i)*n(j))*Kint(i,j)
enddo
EK_GU = EK_GU + n(i)*n(i)*Kint(i,j)
do j=i+1,nMO
EK_GU = EK_GU + sqrt(n(i)*n(j))*Kint(i,j)
enddo
enddo
! Compute POWER exchange energy
EK_POWER = 0d0
alpha = 1d0/3d0
do i=1,nMO
do j=1,nMO
EK_POWER = EK_POWER + (n(i)*n(j))**alpha*Kint(i,j)
enddo
enddo
! Compute total energies
E_SD = ET + EV + EJ_SD - EK_SD
E_MBB = ET + EV + EJ_SD - EK_MBB
E_BBC1 = ET + EV + EJ_SD - EK_BBC1
E_BBC2 = ET + EV + EJ_SD - EK_BBC2
E_CA = ET + EV + EJ_SD - EK_CA
E_CGA = ET + EV + EJ_SD - EK_CGA
E_ML = ET + EV + EJ_SD - EK_ML
E_MLSIC = ET + EV + EJ_SD - EK_MLSIC
E_GU = ET + EV + EJ_SD - EK_GU
E_POWER = ET + EV + EJ_SD - EK_POWER
! Dump energies
print*, '*******************************'
print*, '*** JK NOFT functionals ***'
print*, '*******************************'
print*, ''
print*, '*** Coulomb energies ***'
print*, 'Coulomb SD energy = ',EJ_SD
print*, ''
print*, '*** Exchange energies ***'
print*, 'Exchange SD energy = ',-EK_SD
print*, 'Exchange MBB energy = ',-EK_MBB
print*, 'Exchange BBC1 energy = ',-EK_BBC1
print*, 'Exchange BBC2 energy = ',-EK_BBC2
print*, 'Exchange CA energy = ',-EK_CA
print*, 'Exchange CGA energy = ',-EK_CGA
print*, 'Exchange ML energy = ',-EK_ML
print*, 'Exchange MLSIC energy = ',-EK_MLSIC
print*, 'Exchange GU energy = ',-EK_GU
print*, 'Exchange POWER energy = ',-EK_POWER
print*, ''
print*, ''
print*, '*** Two-electron energies ***'
print*, 'J+K SD energy = ',EJ_SD - EK_SD
print*, 'J+K MBB energy = ',EJ_SD - EK_MBB
print*, 'J+K BBC1 energy = ',EJ_SD - EK_BBC1
print*, 'J+K BBC2 energy = ',EJ_SD - EK_BBC2
print*, 'J+K CA energy = ',EJ_SD - EK_CA
print*, 'J+K CGA energy = ',EJ_SD - EK_CGA
print*, 'J+K ML energy = ',EJ_SD - EK_ML
print*, 'J+K MLSIC energy = ',EJ_SD - EK_MLSIC
print*, 'J+K GU energy = ',EJ_SD - EK_GU
print*, 'J+K POWER energy = ',EJ_SD - EK_POWER
print*, ''
print*, '*** Total energies ***'
print*, 'Total SD energy = ',E_SD
print*, 'Total MBB energy = ',E_MBB
print*, 'Total BBC1 energy = ',E_BBC1
print*, 'Total BBC2 energy = ',E_BBC2
print*, 'Total CA energy = ',E_CA
print*, 'Total CGA energy = ',E_CGA
print*, 'Total ML energy = ',E_ML
print*, 'Total MLSIC energy = ',E_MLSIC
print*, 'Total GU energy = ',E_GU
print*, 'Total POWER energy = ',E_POWER
print*, ''
end subroutine NOFT_JKfunc

46
plugins/NOFT/NOFT_PT2.irp.f Normal file
View File

@ -0,0 +1,46 @@
subroutine NOFT_PT2(nMO,FL,n)
! Compute the PT2 correction from NOFT
END_DOC
PROVIDE mo_bielec_integrals_in_map
! Input variables
integer,intent(in) :: nMO,FL
double precision,intent(in) :: n(nMO)
! Local variables
integer :: i,j,a,b
double precision :: EPT1,EPT2
double precision :: get_mo_bielec_integral
! memory allocation
! Useful quantities
EPT2 = 0d0
! do i=1,FL
! do j=1,FL
! do a=FL+1,nMO
! do b=FL+1,nMO
! enddo
! enddo
! enddo
! enddo
! Dump energies
print*, '*******************************'
print*, '*** PT2 NOFT corrections ***'
print*, '*******************************'
print*, ''
print*, 'Total PT1 energy = ',E_PT1
print*, 'Total PT2 energy = ',E_PT2
print*, 'Total PT1+PT2 energy = ',E_PT1 + E_PT2
print*, ''
end subroutine NOFT_PT2

51
plugins/NOFT/NOFT_core.irp.f Normal file
View File

@ -0,0 +1,51 @@
subroutine NOFT_core(nMO,ET,EV,n)
implicit none
BEGIN_DOC
! Core energy for NOFT
END_DOC
! Input variables
integer,intent(in) :: nMO
double precision,intent(in) :: n(nMO)
! Local variables
integer :: i
! Output variables
double precision,intent(out) :: ET,EV
! Compute kinetic energy
ET = 0d0
do i=1,nMO
ET = ET + n(i)*mo_kinetic_integral(i,i)
enddo
ET = 2d0*ET
! Compute nuclear attraction energy
EV = 0d0
do i=1,nMO
EV = EV + n(i)*mo_nucl_elec_integral(i,i)
enddo
EV = 2d0*EV
! Dump energies
print*, '*******************************'
print*, '*** Core energies ***'
print*, '*******************************'
print*, ''
print*, 'Kinetic energy = ',ET
print*, 'Nuclear attraction energy = ',EV
print*, ''
end subroutine NOFT_core

12
plugins/NOFT/README.rst Normal file
View File

@ -0,0 +1,12 @@
====
NOFT
====
Needed Modules
==============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.
Documentation
=============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.

106
plugins/NOFT/ezfio_interface.irp.f Normal file
View File

@ -0,0 +1,106 @@
! DO NOT MODIFY BY HAND
! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py
! from file /home/loos/quantum_package/src/NOFT/EZFIO.cfg
BEGIN_PROVIDER [ logical, do_jk_functionals ]
implicit none
BEGIN_DOC
! Compute energies for JK-only functionals
END_DOC
logical :: has
PROVIDE ezfio_filename
if (mpi_master) then
call ezfio_has_noft_do_jk_functionals(has)
if (has) then
call ezfio_get_noft_do_jk_functionals(do_jk_functionals)
else
print *, 'noft/do_jk_functionals not found in EZFIO file'
stop 1
endif
endif
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST( do_jk_functionals, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read do_jk_functionals with MPI'
endif
IRP_ENDIF
call write_time(6)
if (mpi_master) then
write(6, *) 'Read do_jk_functionals'
endif
END_PROVIDER
BEGIN_PROVIDER [ logical, do_jkl_functionals ]
implicit none
BEGIN_DOC
! Compute energies for JKL-only functionals (PNOFs)
END_DOC
logical :: has
PROVIDE ezfio_filename
if (mpi_master) then
call ezfio_has_noft_do_jkl_functionals(has)
if (has) then
call ezfio_get_noft_do_jkl_functionals(do_jkl_functionals)
else
print *, 'noft/do_jkl_functionals not found in EZFIO file'
stop 1
endif
endif
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST( do_jkl_functionals, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read do_jkl_functionals with MPI'
endif
IRP_ENDIF
call write_time(6)
if (mpi_master) then
write(6, *) 'Read do_jkl_functionals'
endif
END_PROVIDER
BEGIN_PROVIDER [ logical, do_pt2_noft ]
implicit none
BEGIN_DOC
! Compute PT2 correction for NOFT
END_DOC
logical :: has
PROVIDE ezfio_filename
if (mpi_master) then
call ezfio_has_noft_do_pt2_noft(has)
if (has) then
call ezfio_get_noft_do_pt2_noft(do_pt2_noft)
else
print *, 'noft/do_pt2_noft not found in EZFIO file'
stop 1
endif
endif
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST( do_pt2_noft, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read do_pt2_noft with MPI'
endif
IRP_ENDIF
call write_time(6)
if (mpi_master) then
write(6, *) 'Read do_pt2_noft'
endif
END_PROVIDER