10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-02 11:25:26 +02:00

merge with toto

This commit is contained in:
Emmanuel Giner 2018-12-21 16:36:53 +01:00
commit fdb6147b0a
24 changed files with 125 additions and 983 deletions

3
TODO
View File

@ -1,6 +1,3 @@
* Virer tous les modules qui sont dans plugins
* Separer integrales ERF
# qp_module
* Mettre les fichiers de test dans le directory source

View File

@ -19,7 +19,7 @@ IRPF90_FLAGS : --ninja --align=32 -DMPI
#
[OPTION]
MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
CACHE : 1 ; Enable cache_compile.py
CACHE : 0 ; Enable cache_compile.py
OPENMP : 1 ; Append OpenMP flags
# Optimization flags

View File

@ -23,7 +23,7 @@ IRPF90_FLAGS : --ninja --align=32
#
[OPTION]
MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
CACHE : 1 ; Enable cache_compile.py
CACHE : 0 ; Enable cache_compile.py
OPENMP : 1 ; Append OpenMP flags
# Optimization flags

View File

@ -23,7 +23,7 @@ IRPF90_FLAGS : --ninja --align=32
#
[OPTION]
MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
CACHE : 1 ; Enable cache_compile.py
CACHE : 0 ; Enable cache_compile.py
OPENMP : 1 ; Append OpenMP flags
# Optimization flags

View File

@ -23,7 +23,7 @@ IRPF90_FLAGS : --ninja --align=32 -DMPI
#
[OPTION]
MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
CACHE : 1 ; Enable cache_compile.py
CACHE : 0 ; Enable cache_compile.py
OPENMP : 1 ; Append OpenMP flags
# Optimization flags

View File

@ -19,7 +19,7 @@ IRPF90_FLAGS : --ninja --align=32 --assert
#
[OPTION]
MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
CACHE : 1 ; Enable cache_compile.py
CACHE : 0 ; Enable cache_compile.py
OPENMP : 1 ; Append OpenMP flags
# Optimization flags

View File

@ -19,7 +19,7 @@ IRPF90_FLAGS : --ninja --align=32 -DMPI
#
[OPTION]
MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
CACHE : 1 ; Enable cache_compile.py
CACHE : 0 ; Enable cache_compile.py
OPENMP : 1 ; Append OpenMP flags
# Optimization flags

View File

@ -23,7 +23,7 @@ IRPF90_FLAGS : --ninja --align=32 --assert
#
[OPTION]
MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
CACHE : 1 ; Enable cache_compile.py
CACHE : 0 ; Enable cache_compile.py
OPENMP : 1 ; Append OpenMP flags
# Optimization flags

View File

@ -0,0 +1,12 @@
=====================
Creating a new plugin
=====================
Create a repository, for example :file:`qp_plugins_user`, hosted somewhere
(GitLab, GitHub, etc...), and clone the repository in the
:file:`$QP_ROOT/plugins` directory.

1
plugins/local/.gitignore vendored Normal file
View File

@ -0,0 +1 @@
*

View File

@ -1,19 +0,0 @@
[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

View File

@ -1 +0,0 @@
hartree_fock determinants

View File

@ -1,72 +0,0 @@
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

View File

@ -1,484 +0,0 @@
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

View File

@ -1,260 +0,0 @@
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

View File

@ -1,46 +0,0 @@
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

View File

@ -1,51 +0,0 @@
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

View File

@ -1,12 +0,0 @@
====
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.

View File

@ -1,5 +0,0 @@
Please create a new repository on github/gitlab for your plugins. For example, see http://gitlab.com/scemama/qp_plugins and put your QP programs there.
Some useless files might remain in the directories. Please check and create nice README files.

View File

@ -2,16 +2,17 @@
# -*- coding: utf-8 -*-
"""
Usage:
qp_plugins list [ -i | -u ]
qp_plugins list [ -i | -u | -q ]
qp_plugins download <url>
qp_plugins install <name>...
qp_plugins uninstall <name>
qp_plugins create -n <name> [<needed_modules>...]
qp_plugins create -n <name> [-r <repository>] [<needed_modules>...]
Options:
list List all the plugins
-i List only the installed plugins
-u List only the uninstalled plugins
-q List the external repositories
download <url> Download an external repository.
The URL points to a tar.gz file or a git repository:
@ -23,6 +24,7 @@ Options:
uninstall Uninstall a plugin
create -n <name> Create a new plugin named <name>
-r <repository> Name of the repository in which to create the plugin
"""
@ -84,29 +86,33 @@ def save_new_module(path, l_child):
def main(arguments):
if arguments["list"]:
# Search in src all directories with a NEED file
l_tmp = [ dirname for (dirname, _, filenames) in os.walk(QP_PLUGINS, followlinks=False) for f in filenames if f == 'NEED']
if arguments["-q"]:
l_result = [ f for f in listdir(QP_PLUGINS) if f not in [ ".gitignore", "local" ] ]
# Find directories which contain modules
l_tmp = [ os.path.split(f) for f in l_tmp ]
d_tmp = {}
for (x,_) in l_tmp:
d_tmp[x] = 1
l_repository = d_tmp.keys()
m_all_instances = ModuleHandler(l_repository)
m_instance = ModuleHandler(l_repository)
l_plugins = [ module for module in m_instance.l_module ]
l_result = l_plugins
else:
# Search in src all directories with a NEED file
l_tmp = [ dirname for (dirname, _, filenames) in os.walk(QP_PLUGINS, followlinks=False) for f in filenames if f == 'NEED']
if arguments["-i"] or arguments["-u"]:
# Search in src all symbolic links that are modules
l_installed = [ f for f in os.listdir(QP_SRC) if (os.path.islink(os.path.join(QP_SRC,f)) and f != ".gitignore") ]
# Find directories which contain modules
l_tmp = [ os.path.split(f) for f in l_tmp ]
d_tmp = {}
for (x,_) in l_tmp:
d_tmp[x] = 1
l_repository = d_tmp.keys()
m_all_instances = ModuleHandler(l_repository)
m_instance = ModuleHandler(l_repository)
l_plugins = [ module for module in m_instance.l_module ]
l_result = l_plugins
if arguments["-i"]:
l_result = [ f for f in l_plugins if f in l_installed ]
if arguments["-i"] or arguments["-u"]:
# Search in src all symbolic links that are modules
l_installed = [ f for f in listdir(QP_SRC) if (os.path.islink(os.path.join(QP_SRC,f)) and f != ".gitignore") ]
elif arguments["-u"]:
l_result = [ f for f in l_plugins if f not in l_installed ]
if arguments["-i"]:
l_result = [ f for f in l_plugins if f in l_installed ]
elif arguments["-u"]:
l_result = [ f for f in l_plugins if f not in l_installed ]
for module in sorted(l_result):
print "* {0}".format(module)
@ -114,11 +120,17 @@ def main(arguments):
if arguments["create"]:
m_instance = ModuleHandler([QP_SRC])
print arguments
l_children = arguments["<needed_modules>"]
name = arguments["<name>"][0]
path = os.path.join(QP_PLUGINS, name)
if arguments["-r"]:
repository = arguments["-r"]
else:
repository = "local"
path = os.path.join(QP_PLUGINS, repository, name)
print "Created plugin:"
print path, '\n'

68
scripts/qp_test Executable file
View File

@ -0,0 +1,68 @@
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Runs all the possible tests using bats.
Usage:
qp_test [-av]
Options:
-v verbose output
-a test all installed modules. Default is to test the current directory.
"""
import os
import subprocess
try:
from docopt import docopt
from qp_path import QP_SRC, QP_PLUGINS, QP_ROOT, QP_TESTS
except ImportError:
print "Please check if you have sourced the ${QP_ROOT}/quantum_package.rc"
print "(`source ${QP_ROOT}/quantum_package.rc`)"
print sys.exit(1)
def main(arguments):
# Fetch all *.bats files
l_bats = []
def append_bats(dirname, filenames):
for f in filenames:
if f.endswith(".bats"):
l_bats.append( os.path.join(dirname,f) )
if arguments["-a"]:
for (dirname, _, filenames) in os.walk(QP_SRC, followlinks=False) :
append_bats(dirname, filenames)
else:
for (dirname, _, filenames) in os.walk(os.getcwd(), followlinks=False) :
append_bats(dirname, filenames)
print l_bats
# Execute tests
os.chdir(QP_TESTS)
for bats_file in l_bats:
print ""
print "-~-~-~-~-~-~"
print ""
print "Running tests for %s"%(bats_file)
print ""
if arguments["-v"]:
p1 = subprocess.Popen(["python2", "bats_to_sh.py", bats_file], \
stdout=subprocess.PIPE)
p2 = subprocess.Popen(["bash"], stdin=p1.stdout)
_, _ = os.waitpid(p2.pid,0)
_, _ = os.waitpid(p1.pid,0)
else:
subprocess.check_call(["bats", bats_file])
if __name__ == '__main__':
arguments = docopt(__doc__)
main(arguments)

View File

@ -14,3 +14,4 @@ else:
QP_PLUGINS = os.path.join(QP_ROOT, "plugins")
QP_EZFIO = os.environ["QP_EZFIO"]
QP_OCAML = os.path.join(QP_ROOT, "ocaml")
QP_TESTS = os.path.join(QP_ROOT, "tests")

View File

@ -4,6 +4,7 @@ source $QP_ROOT/tests/bats/common.bats.sh
function run_init() {
cp "${QP_ROOT}/tests/input/$1" .
rm -rf -- $3
qp_create_ezfio_from_xyz $1 -o $3 $2
qp_edit -c $3
}
@ -14,7 +15,7 @@ function run_HF() {
test_exe scf || skip
qp_edit -c $1
ezfio set_file $1
ezfio set hartree_fock thresh_scf 1.e-10
ezfio set scf_utils thresh_scf 1.e-10
qp_run scf $1
energy="$(ezfio get hartree_fock energy)"
eq $energy $2 $thresh