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

Merge branch 'master' of github.com:scemama/quantum_package

This commit is contained in:
Anthony Scemama 2018-12-19 14:34:20 +01:00
commit efc9c449e5
57 changed files with 6602 additions and 50 deletions

2
configure vendored
View File

@ -27,7 +27,7 @@ EOF
exit
}
export QP_ROOT=${PWD}
PACKAGES=""
while : ; do

View File

@ -0,0 +1,368 @@
subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center)
implicit none
BEGIN_DOC
! subroutine that returs all integrals over r of type erf(mu_in * |r-C_center|)/|r-C_center|
END_DOC
double precision, intent(in) :: mu_in,C_center(3)
double precision, intent(out) :: integrals_ao(ao_num,ao_num)
double precision :: NAI_pol_mult_erf_ao
integer :: i,j,l,k,m
do k = 1, ao_num
do m = 1, ao_num
integrals_ao(m,k) = NAI_pol_mult_erf_ao(m,k,mu_in,C_center)
enddo
enddo
end
double precision function NAI_pol_mult_erf_ao(i_ao,j_ao,mu_in,C_center)
implicit none
BEGIN_DOC
! computes the following integral :
! int[-infty;+infty] dr AO_i_ao (r) AO_j_ao(r) erf(mu_in * |r-C_center|)/|r-C_center|
END_DOC
integer, intent(in) :: i_ao,j_ao
double precision, intent(in) :: mu_in, C_center(3)
integer :: i,j,num_A,num_B, power_A(3), power_B(3), n_pt_in
double precision :: A_center(3), B_center(3),integral, alpha,beta, NAI_pol_mult_erf
num_A = ao_nucl(i_ao)
power_A(1:3)= ao_power(i_ao,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
num_B = ao_nucl(j_ao)
power_B(1:3)= ao_power(j_ao,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
n_pt_in = n_pt_max_integrals
NAI_pol_mult_erf_ao = 0.d0
do i = 1, ao_prim_num(i_ao)
alpha = ao_expo_ordered_transp(i,i_ao)
do j = 1, ao_prim_num(j_ao)
beta = ao_expo_ordered_transp(j,j_ao)
integral = NAI_pol_mult_erf(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in,mu_in)
NAI_pol_mult_erf_ao += integral * ao_coef_normalized_ordered_transp(j,j_ao)*ao_coef_normalized_ordered_transp(i,i_ao)
enddo
enddo
end
double precision function NAI_pol_mult_erf(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in,mu_in)
! function that computes the folowing integral :
! int{dr} of (x-A_x)^ax (x-B_X)^bx exp(-alpha (x-A_x)^2 - beta (x-B_x)^2 ) erf(mu_in*(r-R_c))/(r-R_c)
implicit none
integer, intent(in) :: n_pt_in
double precision,intent(in) :: C_center(3),A_center(3),B_center(3),alpha,beta,mu_in
integer, intent(in) :: power_A(3),power_B(3)
integer :: i,j,k,l,n_pt
double precision :: P_center(3)
double precision :: d(0:n_pt_in),pouet,coeff,dist,const,pouet_2,factor
double precision :: I_n_special_exact,integrate_bourrin,I_n_bibi
double precision :: V_e_n,const_factor,dist_integral,tmp
double precision :: accu,rint,p_inv,p,rho,p_inv_2
integer :: n_pt_out,lmax
include 'utils/constants.include.F'
p = alpha + beta
p_inv = 1.d0/p
p_inv_2 = 0.5d0 * p_inv
rho = alpha * beta * p_inv
dist = 0.d0
dist_integral = 0.d0
do i = 1, 3
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
dist += (A_center(i) - B_center(i))*(A_center(i) - B_center(i))
dist_integral += (P_center(i) - C_center(i))*(P_center(i) - C_center(i))
enddo
const_factor = dist*rho
if(const_factor > 80.d0)then
NAI_pol_mult_erf = 0.d0
return
endif
double precision :: p_new
p_new = mu_in/dsqrt(p+ mu_in * mu_in)
factor = dexp(-const_factor)
coeff = dtwo_pi * factor * p_inv * p_new
lmax = 20
! print*, "b"
do i = 0, n_pt_in
d(i) = 0.d0
enddo
n_pt = 2 * ( (power_A(1) + power_B(1)) +(power_A(2) + power_B(2)) +(power_A(3) + power_B(3)) )
const = p * dist_integral * p_new * p_new
if (n_pt == 0) then
pouet = rint(0,const)
NAI_pol_mult_erf = coeff * pouet
return
endif
! call give_polynom_mult_center_mono_elec_erf(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in)
p_new = p_new * p_new
call give_polynom_mult_center_mono_elec_erf_opt(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in,p,p_inv,p_inv_2,p_new,P_center)
if(n_pt_out<0)then
NAI_pol_mult_erf = 0.d0
return
endif
accu = 0.d0
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
do i =0 ,n_pt_out,2
accu += d(i) * rint(i/2,const)
enddo
NAI_pol_mult_erf = accu * coeff
end
subroutine give_polynom_mult_center_mono_elec_erf_opt(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in,p,p_inv,p_inv_2,p_new,P_center)
!!!! subroutine that returns the explicit polynom in term of the "t" variable of the following polynomw ::
!!!! I_x1(a_x, d_x,p,q) * I_x1(a_y, d_y,p,q) * I_x1(a_z, d_z,p,q)
!!!! it is for the nuclear electron atraction
implicit none
integer, intent(in) :: n_pt_in
integer,intent(out) :: n_pt_out
double precision, intent(in) :: A_center(3), B_center(3),C_center(3),p,p_inv,p_inv_2,p_new,P_center(3)
double precision, intent(in) :: alpha,beta,mu_in
integer, intent(in) :: power_A(3), power_B(3)
integer :: a_x,b_x,a_y,b_y,a_z,b_z
double precision :: d(0:n_pt_in)
double precision :: d1(0:n_pt_in)
double precision :: d2(0:n_pt_in)
double precision :: d3(0:n_pt_in)
double precision :: accu
accu = 0.d0
!COMPTEUR irp_rdtsc1 = irp_rdtsc()
ASSERT (n_pt_in > 1)
double precision :: R1x(0:2), B01(0:2), R1xp(0:2),R2x(0:2)
R1x(0) = (P_center(1) - A_center(1))
R1x(1) = 0.d0
R1x(2) = -(P_center(1) - C_center(1))* p_new
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R1xp(0) = (P_center(1) - B_center(1))
R1xp(1) = 0.d0
R1xp(2) =-(P_center(1) - C_center(1))* p_new
!R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R2x(0) = p_inv_2
R2x(1) = 0.d0
R2x(2) = -p_inv_2* p_new
!R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2
do i = 0,n_pt_in
d(i) = 0.d0
enddo
do i = 0,n_pt_in
d1(i) = 0.d0
enddo
do i = 0,n_pt_in
d2(i) = 0.d0
enddo
do i = 0,n_pt_in
d3(i) = 0.d0
enddo
integer :: n_pt1,n_pt2,n_pt3,dim,i
n_pt1 = n_pt_in
n_pt2 = n_pt_in
n_pt3 = n_pt_in
a_x = power_A(1)
b_x = power_B(1)
call I_x1_pol_mult_mono_elec(a_x,b_x,R1x,R1xp,R2x,d1,n_pt1,n_pt_in)
! print*,'passed the first I_x1'
if(n_pt1<0)then
n_pt_out = -1
do i = 0,n_pt_in
d(i) = 0.d0
enddo
return
endif
R1x(0) = (P_center(2) - A_center(2))
R1x(1) = 0.d0
R1x(2) = -(P_center(2) - C_center(2))* p_new
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R1xp(0) = (P_center(2) - B_center(2))
R1xp(1) = 0.d0
R1xp(2) =-(P_center(2) - C_center(2))* p_new
!R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
a_y = power_A(2)
b_y = power_B(2)
call I_x1_pol_mult_mono_elec(a_y,b_y,R1x,R1xp,R2x,d2,n_pt2,n_pt_in)
! print*,'passed the second I_x1'
if(n_pt2<0)then
n_pt_out = -1
do i = 0,n_pt_in
d(i) = 0.d0
enddo
return
endif
R1x(0) = (P_center(3) - A_center(3))
R1x(1) = 0.d0
R1x(2) = -(P_center(3) - C_center(3))* p_new
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R1xp(0) = (P_center(3) - B_center(3))
R1xp(1) = 0.d0
R1xp(2) =-(P_center(3) - C_center(3))* p_new
!R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2
a_z = power_A(3)
b_z = power_B(3)
! print*,'a_z = ',a_z
! print*,'b_z = ',b_z
call I_x1_pol_mult_mono_elec(a_z,b_z,R1x,R1xp,R2x,d3,n_pt3,n_pt_in)
! print*,'passed the third I_x1'
if(n_pt3<0)then
n_pt_out = -1
do i = 0,n_pt_in
d(i) = 0.d0
enddo
return
endif
integer :: n_pt_tmp
n_pt_tmp = 0
call multiply_poly(d1,n_pt1,d2,n_pt2,d,n_pt_tmp)
do i = 0,n_pt_tmp
d1(i) = 0.d0
enddo
n_pt_out = 0
call multiply_poly(d ,n_pt_tmp ,d3,n_pt3,d1,n_pt_out)
do i = 0, n_pt_out
d(i) = d1(i)
enddo
end
subroutine give_polynom_mult_center_mono_elec_erf(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in)
!!!! subroutine that returns the explicit polynom in term of the "t" variable of the following polynomw ::
!!!! I_x1(a_x, d_x,p,q) * I_x1(a_y, d_y,p,q) * I_x1(a_z, d_z,p,q)
!!!! it is for the nuclear electron atraction
implicit none
integer, intent(in) :: n_pt_in
integer,intent(out) :: n_pt_out
double precision, intent(in) :: A_center(3), B_center(3),C_center(3)
double precision, intent(in) :: alpha,beta,mu_in
integer, intent(in) :: power_A(3), power_B(3)
integer :: a_x,b_x,a_y,b_y,a_z,b_z
double precision :: d(0:n_pt_in)
double precision :: d1(0:n_pt_in)
double precision :: d2(0:n_pt_in)
double precision :: d3(0:n_pt_in)
double precision :: accu, pq_inv, p10_1, p10_2, p01_1, p01_2
double precision :: p,P_center(3),rho,p_inv,p_inv_2
accu = 0.d0
!COMPTEUR irp_rdtsc1 = irp_rdtsc()
ASSERT (n_pt_in > 1)
p = alpha+beta
p_inv = 1.d0/p
p_inv_2 = 0.5d0/p
do i =1, 3
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
enddo
double precision :: R1x(0:2), B01(0:2), R1xp(0:2),R2x(0:2)
R1x(0) = (P_center(1) - A_center(1))
R1x(1) = 0.d0
R1x(2) = -(P_center(1) - C_center(1))* mu_in**2 / (p+mu_in*mu_in)
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R1xp(0) = (P_center(1) - B_center(1))
R1xp(1) = 0.d0
R1xp(2) =-(P_center(1) - C_center(1))* mu_in**2 / (p+mu_in*mu_in)
!R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R2x(0) = p_inv_2
R2x(1) = 0.d0
R2x(2) = -p_inv_2* mu_in**2 / (p+mu_in*mu_in)
!R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2
do i = 0,n_pt_in
d(i) = 0.d0
enddo
do i = 0,n_pt_in
d1(i) = 0.d0
enddo
do i = 0,n_pt_in
d2(i) = 0.d0
enddo
do i = 0,n_pt_in
d3(i) = 0.d0
enddo
integer :: n_pt1,n_pt2,n_pt3,dim,i
n_pt1 = n_pt_in
n_pt2 = n_pt_in
n_pt3 = n_pt_in
a_x = power_A(1)
b_x = power_B(1)
call I_x1_pol_mult_mono_elec(a_x,b_x,R1x,R1xp,R2x,d1,n_pt1,n_pt_in)
! print*,'passed the first I_x1'
if(n_pt1<0)then
n_pt_out = -1
do i = 0,n_pt_in
d(i) = 0.d0
enddo
return
endif
R1x(0) = (P_center(2) - A_center(2))
R1x(1) = 0.d0
R1x(2) = -(P_center(2) - C_center(2))* mu_in**2 / (p+mu_in*mu_in)
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R1xp(0) = (P_center(2) - B_center(2))
R1xp(1) = 0.d0
R1xp(2) =-(P_center(2) - C_center(2))* mu_in**2 / (p+mu_in*mu_in)
!R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
a_y = power_A(2)
b_y = power_B(2)
call I_x1_pol_mult_mono_elec(a_y,b_y,R1x,R1xp,R2x,d2,n_pt2,n_pt_in)
! print*,'passed the second I_x1'
if(n_pt2<0)then
n_pt_out = -1
do i = 0,n_pt_in
d(i) = 0.d0
enddo
return
endif
R1x(0) = (P_center(3) - A_center(3))
R1x(1) = 0.d0
R1x(2) = -(P_center(3) - C_center(3))* mu_in**2 / (p+mu_in*mu_in)
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R1xp(0) = (P_center(3) - B_center(3))
R1xp(1) = 0.d0
R1xp(2) =-(P_center(3) - C_center(3))* mu_in**2 / (p+mu_in*mu_in)
!R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2
a_z = power_A(3)
b_z = power_B(3)
! print*,'a_z = ',a_z
! print*,'b_z = ',b_z
call I_x1_pol_mult_mono_elec(a_z,b_z,R1x,R1xp,R2x,d3,n_pt3,n_pt_in)
! print*,'passed the third I_x1'
if(n_pt3<0)then
n_pt_out = -1
do i = 0,n_pt_in
d(i) = 0.d0
enddo
return
endif
integer :: n_pt_tmp
n_pt_tmp = 0
call multiply_poly(d1,n_pt1,d2,n_pt2,d,n_pt_tmp)
do i = 0,n_pt_tmp
d1(i) = 0.d0
enddo
n_pt_out = 0
call multiply_poly(d ,n_pt_tmp ,d3,n_pt3,d1,n_pt_out)
do i = 0, n_pt_out
d(i) = d1(i)
enddo
end

View File

@ -2,7 +2,14 @@
Becke Numerical Grid
====================
Lebedev-Laikov grids, using the code distributed through CCL (http://www.ccl.net/).
This module contains all quantities needed to build the Becke's grid used in general for DFT integration. Note that it can be used for whatever integration in R^3 as long as the functions to be integrated are mostly concentrated near the atomic regions.
This grid is built as the reunion of a spherical grid around each atom. Each spherical grid contains a certain number of radial and angular points.
For a simple example of how to use the grid, see example.irp.f.
The spherical integration uses Lebedev-Laikov grids, which was used from the code distributed through CCL (http://www.ccl.net/).
See next section for explanations and citation policies.
.. code-block:: text

View File

@ -0,0 +1,71 @@
subroutine example_becke_numerical_grid
implicit none
include 'constants.include.F'
BEGIN_DOC
! subroutine that illustrates the main features available in becke_numerical_grid
END_DOC
integer :: i,j,k,ipoint
double precision :: integral_1, integral_2,alpha,center(3)
print*,''
print*,'**************'
print*,'**************'
print*,'routine that illustrates the use of the grid'
print*,'**************'
print*,'This grid is built as the reunion of a spherical grid around each atom'
print*,'Each spherical grid contains a certain number of radial and angular points'
print*,''
print*,'n_points_integration_angular = ',n_points_integration_angular
print*,'n_points_radial_grid = ',n_points_radial_grid
print*,''
print*,'As an example of the use of the grid, we will compute the integral of a 3D gaussian'
! parameter of the gaussian: center of the gaussian is set to the first nucleus
center(1:3)=nucl_coord(1,1:3)
! alpha = exponent of the gaussian
alpha = 1.d0
print*,''
print*,'The first example uses the grid points as one-dimensional array'
print*,'This is the mostly used representation of the grid'
print*,'It is the easyest way to use it with no drawback in terms of accuracy'
integral_1 = 0.d0
! you browse all the grid points as a one-dimensional array
do i = 1, n_points_final_grid
double precision :: weight, r(3)
! you get x, y and z of the ith grid point
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
weight = final_weight_functions_at_final_grid_points(i)
double precision :: distance, f_r
! you compute the function to be integrated
distance = dsqrt( (r(1) - center(1))**2 + (r(2) - center(2))**2 + (r(3) - center(3))**2 )
f_r = dexp(-alpha * distance)
! you add the contribution of the grid point to the integral
integral_1 += f_r * weight
enddo
print*,'integral_1 =',integral_1
print*,'(pi/alpha)**1.5 =',(pi / alpha)**1.5
print*,''
print*,''
print*,'The second example uses the grid points as a collection of spherical grids centered on each atom'
print*,'This is mostly useful if one needs to split contributions between radial/angular/atomic of an integral'
! you browse the nuclei
do i = 1, nucl_num
! you browse the radial points attached to each nucleus
do j = 1, n_points_radial_grid
! you browse the angular points attached to each radial point of each nucleus
do k = 1, n_points_integration_angular
r(1) = grid_points_per_atom(1,k,j,i)
r(2) = grid_points_per_atom(2,k,j,i)
r(3) = grid_points_per_atom(3,k,j,i)
weight = final_weight_functions_at_grid_points(k,j,i)
distance = dsqrt( (r(1) - center(1))**2 + (r(2) - center(2))**2 + (r(3) - center(3))**2 )
f_r = dexp(-alpha * distance)
integral_2 += f_r * weight
enddo
enddo
enddo
print*,'integral_2 =',integral_2
print*,'(pi/alpha)**1.5 =',(pi / alpha)**1.5
print*,''
end

View File

@ -123,8 +123,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, &
(n_points_integration_angular,n_points_radial_grid,nucl_num) ]
BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_integration_angular,n_points_radial_grid,nucl_num) ]
BEGIN_DOC
! Weight function at grid points : w_n(r) according to the equation (22)
! of Becke original paper (JCP, 88, 1988)
@ -158,7 +157,6 @@ BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, &
enddo
accu = 1.d0/accu
weight_functions_at_grid_points(l,k,j) = tmp_array(j) * accu
! print*,weight_functions_at_grid_points(l,k,j)
enddo
enddo
enddo

View File

@ -1,3 +1,34 @@
subroutine set_bit_to_integer(i_physical,key,Nint)
use bitmasks
BEGIN_DOC
! set to 1 the bit number i_physical in the bitstring key
END_DOC
implicit none
integer, intent(in) :: i_physical,Nint
integer(bit_kind), intent(inout) :: key(Nint)
integer :: k,j,i
k = ishft(i_physical-1,-bit_kind_shift)+1
j = i_physical-ishft(k-1,bit_kind_shift)-1
key(k) = ibset(key(k),j)
end
subroutine clear_bit_to_integer(i_physical,key,Nint)
use bitmasks
BEGIN_DOC
! set to 0 the bit number i_physical in the bitstring key
END_DOC
implicit none
integer, intent(in) :: i_physical,Nint
integer(bit_kind), intent(inout) :: key(Nint)
integer :: k,j,i
k = ishft(i_physical-1,-bit_kind_shift)+1
j = i_physical-ishft(k-1,bit_kind_shift)-1
key(k) = ibclr(key(k),j)
end
subroutine bitstring_to_list( string, list, n_elements, Nint)
use bitmasks
implicit none

81
src/bitmask/example.irp.f Normal file
View File

@ -0,0 +1,81 @@
subroutine example_bitmask
use bitmasks ! you need to include the bitmasks_module.f90 features
implicit none
BEGIN_DOC
! subroutine that illustrates the main features available in bitmask
END_DOC
integer :: i,j
print*,''
print*,'**************'
print*,'**************'
print*,'MO class: to set the various type of MO class, see the following exectuable'
print*,'qp_set_mo_class'
print*,'**************'
print*,'number of core orbitals = ',n_core_orb
print*,'list of the core orbitals '
do i = 1, n_core_orb
write(*,'(2(I3,X))')i,list_core(i)
enddo
print*,'number of inact orbitals = ',n_inact_orb
print*,'list of the inact orbitals '
do i = 1, n_inact_orb
write(*,'(2(I3,X))')i,list_inact(i)
enddo
print*,'number of act orbitals = ',n_act_orb
print*,'list of the act orbitals '
do i = 1, n_act_orb
write(*,'(2(I3,X))')i,list_act(i)
enddo
print*,'number of virt orbitals = ',n_virt_orb
print*,'list of the virt orbitals '
do i = 1, n_virt_orb
write(*,'(2(I3,X))')i,list_virt(i)
enddo
print*,''
print*,'**************'
print*,'**************'
print*,'manipulating bitstrings (usefull for determinant representation)'
print*,'**************'
integer(bit_kind), allocatable :: key(:)
print*,'Size of the integers used to represent all the orbitals '
print*,'bit_kind = ',bit_kind
print*,'Number of bits in the integers ',bit_kind_size
print*,'Number of integers to represent all the orbitals on integer'
print*,'N_int = ',N_int
allocate(key(N_int))
print*,'**** '
print*,' initialize a bistring to zero '
do i = 1, N_int
key(i) = 0_bit_kind
enddo
print*,'print a human readable representation of the bitstring'
call bitstring_to_str( output, key, N_int )
print *, trim(output)
integer :: i_orb
character*(2048) :: output
do i_orb = 1, min(4,mo_tot_num) ! you set the first four bits to 1 in key
call set_bit_to_integer(i_orb,key,N_int)
enddo
print*,'print a human readable representation of the bitstring'
call bitstring_to_str( output, key, N_int )
print *, trim(output)
print*,''
integer :: n_elements
integer, allocatable :: list_occ(:)
allocate(list_occ(N_int*bit_kind_size))
call bitstring_to_list( key, list_occ, n_elements, N_int)
print*,'number of bits set to 1 = ',n_elements
print*,'list of bits set to 1 '
do i = 1, n_elements
write(*,'(2(I3,X))')i,list_occ(i)
enddo
call clear_bit_to_integer(2,key,N_int) ! you set to 0 the second bit
print*,'print a human readable representation of the bitstring'
call bitstring_to_str( output, key, N_int )
print *, trim(output)
end

View File

@ -0,0 +1,26 @@
[data_energy_var]
type: double precision
doc: Variational energy computed with the wave function
interface: ezfio, provider
size: (determinants.n_states)
[data_energy_proj]
type: double precision
doc: Projected energy computed with the wave function
interface: ezfio, provider
size: (determinants.n_states)
[data_one_body_alpha_dm_mo]
interface: ezfio, provider
doc: Alpha one body density matrix on the MO basis computed with the wave function
type: double precision
size: (mo_basis.mo_tot_num,mo_basis.mo_tot_num,determinants.n_states)
[data_one_body_beta_dm_mo]
interface: ezfio, provider
doc: Beta one body density matrix on the MO basis computed with the wave function
type: double precision
size: (mo_basis.mo_tot_num,mo_basis.mo_tot_num,determinants.n_states)

View File

@ -0,0 +1 @@
determinants

View File

@ -0,0 +1,6 @@
============
Data energy and density
============
This module contains some global variables (such as densities and energies) which are stored in the EZFIO folder in a different place than determinants. This is used in practice to store density matrices which can be obtained from any methods, as long as they are stored in the same MO basis which is used for the calculations. In |RS-DFT| calculations, this can be done to perform damping on the density in order to speed up convergence.

View File

@ -0,0 +1,20 @@
program save_one_body_dm
implicit none
BEGIN_DOC
! programs that computes the one body density on the mo basis for alpha and beta electrons from the wave function stored in the EZFIO folder, and then save it into the EZFIO folder data_energy_and_density.
!
! Then, the global variable data_one_body_alpha_dm_mo and data_one_body_beta_dm_mo will automatically read the density in a further calculation.
!
! This can be used to perform dampin on the density in RS-DFT calculation (see the density_for_dft module).
END_DOC
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
call ezfio_set_data_energy_and_density_data_one_body_alpha_dm_mo(one_body_dm_mo_alpha)
call ezfio_set_data_energy_and_density_data_one_body_beta_dm_mo(one_body_dm_mo_beta)
end

View File

@ -0,0 +1,494 @@
BEGIN_PROVIDER [ double precision, psi_energy_bielec, (N_states) ]
implicit none
BEGIN_DOC
! Energy of the current wave function
END_DOC
integer :: i,j
call u_0_H_u_0_bielec(psi_energy_bielec,psi_coef,N_det,psi_det,N_int,N_states,psi_det_size)
do i=N_det+1,N_states
psi_energy(i) = 0.d0
enddo
END_PROVIDER
subroutine H_S2_u_0_bielec_nstates_openmp(v_0,s_0,u_0,N_st,sze)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_0 = H|u_0> and s_0 = S^2 |u_0>
!
! Assumes that the determinants are in psi_det
!
! istart, iend, ishift, istep are used in ZMQ parallelization.
END_DOC
integer, intent(in) :: N_st,sze
double precision, intent(inout) :: v_0(sze,N_st), s_0(sze,N_st), u_0(sze,N_st)
integer :: k
double precision, allocatable :: u_t(:,:), v_t(:,:), s_t(:,:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
allocate(u_t(N_st,N_det),v_t(N_st,N_det),s_t(N_st,N_det))
do k=1,N_st
call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)
enddo
v_t = 0.d0
s_t = 0.d0
call dtranspose( &
u_0, &
size(u_0, 1), &
u_t, &
size(u_t, 1), &
N_det, N_st)
call H_S2_u_0_bielec_nstates_openmp_work(v_t,s_t,u_t,N_st,sze,1,N_det,0,1)
deallocate(u_t)
call dtranspose( &
v_t, &
size(v_t, 1), &
v_0, &
size(v_0, 1), &
N_st, N_det)
call dtranspose( &
s_t, &
size(s_t, 1), &
s_0, &
size(s_0, 1), &
N_st, N_det)
deallocate(v_t,s_t)
do k=1,N_st
call dset_order(v_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
call dset_order(s_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
enddo
end
subroutine H_S2_u_0_bielec_nstates_openmp_work(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_t = H|u_t> and s_t = S^2 |u_t>
!
! Default should be 1,N_det,0,1
END_DOC
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
double precision, intent(in) :: u_t(N_st,N_det)
double precision, intent(out) :: v_t(N_st,sze), s_t(N_st,sze)
PROVIDE ref_bitmask_energy N_int
select case (N_int)
case (1)
call H_S2_u_0_bielec_nstates_openmp_work_1(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
case (2)
call H_S2_u_0_bielec_nstates_openmp_work_2(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
case (3)
call H_S2_u_0_bielec_nstates_openmp_work_3(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
case (4)
call H_S2_u_0_bielec_nstates_openmp_work_4(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
case default
call H_S2_u_0_bielec_nstates_openmp_work_N_int(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
end select
end
BEGIN_TEMPLATE
subroutine H_S2_u_0_bielec_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_t = H|u_t> and s_t = S^2 |u_t>
!
! Default should be 1,N_det,0,1
END_DOC
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
double precision, intent(in) :: u_t(N_st,N_det)
double precision, intent(out) :: v_t(N_st,sze), s_t(N_st,sze)
double precision :: hij, sij
integer :: i,j,k,l
integer :: k_a, k_b, l_a, l_b, m_a, m_b
integer :: istate
integer :: krow, kcol, krow_b, kcol_b
integer :: lrow, lcol
integer :: mrow, mcol
integer(bit_kind) :: spindet($N_int)
integer(bit_kind) :: tmp_det($N_int,2)
integer(bit_kind) :: tmp_det2($N_int,2)
integer(bit_kind) :: tmp_det3($N_int,2)
integer(bit_kind), allocatable :: buffer(:,:)
integer :: n_doubles
integer, allocatable :: doubles(:)
integer, allocatable :: singles_a(:)
integer, allocatable :: singles_b(:)
integer, allocatable :: idx(:), idx0(:)
integer :: maxab, n_singles_a, n_singles_b, kcol_prev
integer*8 :: k8
maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
allocate(idx0(maxab))
do i=1,maxab
idx0(i) = i
enddo
! Prepare the array of all alpha single excitations
! -------------------------------------------------
PROVIDE N_int nthreads_davidson
!$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) &
!$OMP SHARED(psi_bilinear_matrix_rows, N_det, &
!$OMP psi_bilinear_matrix_columns, &
!$OMP psi_det_alpha_unique, psi_det_beta_unique, &
!$OMP n_det_alpha_unique, n_det_beta_unique, N_int, &
!$OMP psi_bilinear_matrix_transp_rows, &
!$OMP psi_bilinear_matrix_transp_columns, &
!$OMP psi_bilinear_matrix_transp_order, N_st, &
!$OMP psi_bilinear_matrix_order_transp_reverse, &
!$OMP psi_bilinear_matrix_columns_loc, &
!$OMP psi_bilinear_matrix_transp_rows_loc, &
!$OMP istart, iend, istep, irp_here, v_t, s_t, &
!$OMP ishift, idx0, u_t, maxab) &
!$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i, &
!$OMP lcol, lrow, l_a, l_b, &
!$OMP buffer, doubles, n_doubles, &
!$OMP tmp_det2, hij, sij, idx, l, kcol_prev, &
!$OMP singles_a, n_singles_a, singles_b, &
!$OMP n_singles_b, k8)
! Alpha/Beta double excitations
! =============================
allocate( buffer($N_int,maxab), &
singles_a(maxab), &
singles_b(maxab), &
doubles(maxab), &
idx(maxab))
kcol_prev=-1
ASSERT (iend <= N_det)
ASSERT (istart > 0)
ASSERT (istep > 0)
!$OMP DO SCHEDULE(dynamic,64)
do k_a=istart+ishift,iend,istep
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
if (kcol /= kcol_prev) then
call get_all_spin_singles_$N_int( &
psi_det_beta_unique, idx0, &
tmp_det(1,2), N_det_beta_unique, &
singles_b, n_singles_b)
endif
kcol_prev = kcol
! Loop over singly excited beta columns
! -------------------------------------
do i=1,n_singles_b
lcol = singles_b(i)
tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol)
l_a = psi_bilinear_matrix_columns_loc(lcol)
ASSERT (l_a <= N_det)
do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow)
ASSERT (l_a <= N_det)
idx(j) = l_a
l_a = l_a+1
enddo
j = j-1
call get_all_spin_singles_$N_int( &
buffer, idx, tmp_det(1,1), j, &
singles_a, n_singles_a )
! Loop over alpha singles
! -----------------------
do k = 1,n_singles_a
l_a = singles_a(k)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
call i_H_j_double_alpha_beta(tmp_det,tmp_det2,$N_int,hij)
call get_s2(tmp_det,tmp_det2,$N_int,sij)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
s_t(l,k_a) = s_t(l,k_a) + sij * u_t(l,l_a)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP DO SCHEDULE(dynamic,64)
do k_a=istart+ishift,iend,istep
! Single and double alpha excitations
! ===================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
! Initial determinant is at k_b in beta-major representation
! ----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
ASSERT (k_b <= N_det)
spindet(1:$N_int) = tmp_det(1:$N_int,1)
! Loop inside the beta column to gather all the connected alphas
lcol = psi_bilinear_matrix_columns(k_a)
l_a = psi_bilinear_matrix_columns_loc(lcol)
do i=1,N_det_alpha_unique
if (l_a > N_det) exit
lcol = psi_bilinear_matrix_columns(l_a)
if (lcol /= kcol) exit
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow)
idx(i) = l_a
l_a = l_a+1
enddo
i = i-1
call get_all_spin_singles_and_doubles_$N_int( &
buffer, idx, spindet, i, &
singles_a, doubles, n_singles_a, n_doubles )
! Compute Hij for all alpha singles
! ----------------------------------
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
do i=1,n_singles_a
l_a = singles_a(i)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
call i_H_j_mono_spin_bielec( tmp_det, tmp_det2, $N_int, 1, hij)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
! single => sij = 0
enddo
enddo
! Compute Hij for all alpha doubles
! ----------------------------------
do i=1,n_doubles
l_a = doubles(i)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
call i_H_j_double_spin( tmp_det(1,1), psi_det_alpha_unique(1, lrow), $N_int, hij)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
! same spin => sij = 0
enddo
enddo
! Single and double beta excitations
! ==================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
kcol = psi_bilinear_matrix_columns(k_a)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
spindet(1:$N_int) = tmp_det(1:$N_int,2)
! Initial determinant is at k_b in beta-major representation
! -----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
ASSERT (k_b <= N_det)
! Loop inside the alpha row to gather all the connected betas
lrow = psi_bilinear_matrix_transp_rows(k_b)
l_b = psi_bilinear_matrix_transp_rows_loc(lrow)
do i=1,N_det_beta_unique
if (l_b > N_det) exit
lrow = psi_bilinear_matrix_transp_rows(l_b)
if (lrow /= krow) exit
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol)
idx(i) = l_b
l_b = l_b+1
enddo
i = i-1
call get_all_spin_singles_and_doubles_$N_int( &
buffer, idx, spindet, i, &
singles_b, doubles, n_singles_b, n_doubles )
! Compute Hij for all beta singles
! ----------------------------------
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
do i=1,n_singles_b
l_b = singles_b(i)
ASSERT (l_b <= N_det)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol)
call i_H_j_mono_spin_bielec( tmp_det, tmp_det2, $N_int, 2, hij)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_a <= N_det)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
! single => sij = 0
enddo
enddo
! Compute Hij for all beta doubles
! ----------------------------------
do i=1,n_doubles
l_b = doubles(i)
ASSERT (l_b <= N_det)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
call i_H_j_double_spin( tmp_det(1,2), psi_det_beta_unique(1, lcol), $N_int, hij)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_a <= N_det)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
! same spin => sij = 0
enddo
enddo
! Diagonal contribution
! =====================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
double precision, external :: diag_H_mat_elem_bielec, diag_S_mat_elem
hij = diag_H_mat_elem_bielec(tmp_det,$N_int)
sij = diag_S_mat_elem(tmp_det,$N_int)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,k_a)
s_t(l,k_a) = s_t(l,k_a) + sij * u_t(l,k_a)
enddo
end do
!$OMP END DO
deallocate(buffer, singles_a, singles_b, doubles, idx)
!$OMP END PARALLEL
end
SUBST [ N_int ]
1;;
2;;
3;;
4;;
N_int;;
END_TEMPLATE
subroutine u_0_H_u_0_bielec(e_0,u_0,n,keys_tmp,Nint,N_st,sze)
use bitmasks
implicit none
BEGIN_DOC
! Computes e_0 = <u_0|H|u_0>/<u_0|u_0>
!
! n : number of determinants
!
END_DOC
integer, intent(in) :: n,Nint, N_st, sze
double precision, intent(out) :: e_0(N_st)
double precision, intent(inout) :: u_0(sze,N_st)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
double precision, allocatable :: v_0(:,:), s_0(:,:), u_1(:,:)
double precision :: u_dot_u,u_dot_v,diag_H_mat_elem
integer :: i,j
allocate (v_0(n,N_st),s_0(n,N_st),u_1(n,N_st))
u_1(1:n,:) = u_0(1:n,:)
call H_S2_u_0_bielec_nstates_openmp(v_0,s_0,u_1,N_st,n)
u_0(1:n,:) = u_1(1:n,:)
deallocate(u_1)
double precision :: norm
do i=1,N_st
norm = u_dot_u(u_0(1,i),n)
if (norm /= 0.d0) then
e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/u_dot_u(u_0(1,i),n)
else
e_0(i) = 0.d0
endif
enddo
deallocate (s_0, v_0)
end

3
src/density_for_dft/NEED Normal file
View File

@ -0,0 +1,3 @@
determinants
dft_keywords
data_energy_and_density

View File

@ -0,0 +1,4 @@
==========
DM_for_dft
==========

View File

@ -0,0 +1,77 @@
BEGIN_PROVIDER [double precision, one_body_dm_mo_alpha_for_dft, (mo_tot_num,mo_tot_num, N_states)]
implicit none
BEGIN_DOC
! density matrix for alpha electrons in the MO basis used for all DFT calculations based on the density
END_DOC
double precision :: delta_alpha(mo_tot_num,mo_tot_num,N_states)
if(density_for_dft .EQ. "damping_rs_dft")then
delta_alpha = one_body_dm_mo_alpha - data_one_body_alpha_dm_mo
one_body_dm_mo_alpha_for_dft = data_one_body_alpha_dm_mo + damping_for_rs_dft * delta_alpha
else if (density_for_dft .EQ. "input_density")then
one_body_dm_mo_alpha_for_dft = data_one_body_alpha_dm_mo
else if (density_for_dft .EQ. "WFT")then
provide mo_coef
one_body_dm_mo_alpha_for_dft = one_body_dm_mo_alpha
endif
END_PROVIDER
BEGIN_PROVIDER [double precision, one_body_dm_mo_beta_for_dft, (mo_tot_num,mo_tot_num, N_states)]
implicit none
BEGIN_DOC
! density matrix for beta electrons in the MO basis used for all DFT calculations based on the density
END_DOC
double precision :: delta_beta(mo_tot_num,mo_tot_num,N_states)
if(density_for_dft .EQ. "damping_rs_dft")then
delta_beta = one_body_dm_mo_beta - data_one_body_beta_dm_mo
one_body_dm_mo_beta_for_dft = data_one_body_beta_dm_mo + damping_for_rs_dft * delta_beta
else if (density_for_dft .EQ. "input_density")then
one_body_dm_mo_beta_for_dft = data_one_body_beta_dm_mo
else if (density_for_dft .EQ. "WFT")then
provide mo_coef
one_body_dm_mo_beta_for_dft = one_body_dm_mo_beta
endif
END_PROVIDER
BEGIN_PROVIDER [double precision, one_body_dm_mo_for_dft, (mo_tot_num,mo_tot_num, N_states)]
implicit none
one_body_dm_mo_for_dft = one_body_dm_mo_beta_for_dft + one_body_dm_mo_alpha_for_dft
END_PROVIDER
BEGIN_PROVIDER [double precision, one_body_dm_average_mo_for_dft, (mo_tot_num,mo_tot_num)]
implicit none
integer :: i
one_body_dm_average_mo_for_dft = 0.d0
do i = 1, N_states
one_body_dm_average_mo_for_dft(:,:) += one_body_dm_mo_for_dft(:,:,i) * state_average_weight(i)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_dm_alpha_ao_for_dft, (ao_num,ao_num,N_states) ]
&BEGIN_PROVIDER [ double precision, one_body_dm_beta_ao_for_dft, (ao_num,ao_num,N_states) ]
BEGIN_DOC
! one body density matrix on the AO basis based on one_body_dm_mo_alpha_for_dft
END_DOC
implicit none
integer :: i,j,k,l,istate
double precision :: mo_alpha,mo_beta
one_body_dm_alpha_ao_for_dft = 0.d0
one_body_dm_beta_ao_for_dft = 0.d0
do k = 1, ao_num
do l = 1, ao_num
do i = 1, mo_tot_num
do j = 1, mo_tot_num
do istate = 1, N_states
mo_alpha = one_body_dm_mo_alpha_for_dft(j,i,istate)
mo_beta = one_body_dm_mo_beta_for_dft(j,i,istate)
one_body_dm_alpha_ao_for_dft(l,k,istate) += mo_coef(k,i) * mo_coef(l,j) * mo_alpha
one_body_dm_beta_ao_for_dft(l,k,istate) += mo_coef(k,i) * mo_coef(l,j) * mo_beta
enddo
enddo
enddo
enddo
enddo
END_PROVIDER

View File

@ -249,9 +249,21 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
logical :: t
double precision :: hij_elec
! output : 0 : not connected
! i : connected to determinant i of the past
! -i : is the ith determinant of the reference wf keys
BEGIN_DOC
! input : key : a given Slater determinant
!
! : keys: a list of Slater determinants
!
! : Ndet: the number of Slater determinants in keys
!
! : N_past_in the number of Slater determinants for the connectivity research
!
! output : 0 : key not connected to the N_past_in first Slater determinants in keys
!
! i : key is connected to determinant i of keys
!
! -i : key is the ith determinant of the reference wf keys
END_DOC
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
@ -349,9 +361,21 @@ integer function connected_to_ref_by_mono(key,keys,Nint,N_past_in,Ndet)
logical :: t
double precision :: hij_elec
! output : 0 : not connected
! i : connected to determinant i of the past
! -i : is the ith determinant of the refernce wf keys
BEGIN_DOC
! input : key : a given Slater determinant
!
! : keys: a list of Slater determinants
!
! : Ndet: the number of Slater determinants in keys
!
! : N_past_in the number of Slater determinants for the connectivity research
!
! output : 0 : key not connected by a MONO EXCITATION to the N_past_in first Slater determinants in keys
!
! i : key is connected by a MONO EXCITATION to determinant i of keys
!
! -i : key is the ith determinant of the reference wf keys
END_DOC
ASSERT (Nint > 0)
ASSERT (Nint == N_int)

View File

@ -65,26 +65,3 @@ logical function is_spin_flip_possible(key_in,i_flip,ispin)
endif
end
subroutine set_bit_to_integer(i_physical,key,Nint)
use bitmasks
implicit none
integer, intent(in) :: i_physical,Nint
integer(bit_kind), intent(inout) :: key(Nint)
integer :: k,j,i
k = shiftr(i_physical-1,bit_kind_shift)+1
j = i_physical-shiftl(k-1,bit_kind_shift)-1
key(k) = ibset(key(k),j)
end
subroutine clear_bit_to_integer(i_physical,key,Nint)
use bitmasks
implicit none
integer, intent(in) :: i_physical,Nint
integer(bit_kind), intent(inout) :: key(Nint)
integer :: k,j,i
k = shiftr(i_physical-1,bit_kind_shift)+1
j = i_physical-shiftl(k-1,bit_kind_shift)-1
key(k) = ibclr(key(k),j)
end

View File

@ -0,0 +1,165 @@
subroutine example_determinants
use bitmasks ! you need to include the bitmasks_module.f90 features
implicit none
BEGIN_DOC
! subroutine that illustrates the main features available in determinants
END_DOC
print*,'a determinant is stored as a binary representation of the occupancy of the spatial orbitals'
print*,'see the bitmask module for more information about that '
print*,'a spin determinant is an array of (N_int) integers of type bit_kind (see bitmask for more information)'
print*,'A determinant containing alpha and beta electrons is an array of dimension (2,N_int)'
integer(bit_kind), allocatable :: det_i(:,:)
allocate(det_i(N_int,2))
print*,'det_i(1,:) alpha spins '
print*,'det_i(2,:) beta spins '
integer :: i,j
print*,'initialize det_i to an electron occupation corresponding RHF or ROHF: ref_bitmask '
do i = 1, N_int
det_i(i,1) = ref_bitmask(i,1)
det_i(i,2) = ref_bitmask(i,2)
enddo
print*,''
print*,'print a human readable representation of the determinant '
call print_det(det_i,N_int)
print*,'doing a single excitation on top of det_i'
integer :: h1,p1,s1,i_ok
h1 = 1
p1 = elec_alpha_num + 1
s1 = 1
print*,'h1 --> p1 of spin s1'
print*,'i_ok == +1 : excitation is possible '
print*,'i_ok == -1 : excitation is NOT possible '
call do_mono_excitation(det_i,h1,p1,s1,i_ok)
print*,'h1,p1,s1,i_ok'
print*, h1,p1,s1,i_ok
if(i_ok == -1)then
print*,'excitation was not possible '
stop
endif
call debug_det(det_i,N_int)
print*,'computing the interaction between ref_determinant and det_i '
double precision :: h0i,hii,h00
call i_H_j(det_i,det_i,N_int,h0i)
print*,' < ref | H | det_i > = ',h0i
print*,'computing the diagonal Hamiltonian matrix element of det_i '
call i_H_j(ref_bitmask,det_i,N_int,hii)
print*,'< det_i | H | det_i > = ',hii
print*,'computing the first-order coefficient of det_i with H0=EN '
double precision :: c_i
call i_H_j(ref_bitmask,ref_bitmask,N_int,h00)
c_i = h0i/(h00 - hii)
print*,'c_i^{(1)} = ',c_i
print*,''
print*,'doing another single excitation on top of det_i'
h1 = elec_alpha_num
p1 = elec_alpha_num + 1
s1 = 2
call do_mono_excitation(det_i,h1,p1,s1,i_ok)
print*,'h1,p1,s1,i_ok'
print*, h1,p1,s1,i_ok
call i_H_j(det_i,det_i,N_int,h0i)
print*,' < ref | H | det_i > = ',h0i
print*,'computing the diagonal Hamiltonian matrix element of det_i '
call i_H_j(ref_bitmask,ref_bitmask,N_int,h00)
c_i = h0i/(h00 - hii)
print*,'c_i^{(1)} = ',c_i
print*,''
print*,'Finding the excitation degree between two arbitrary determinants '
integer :: exc(0:2,2,2)
double precision :: phase
integer :: h2,p2,s2,degree
call get_excitation_degree(ref_bitmask,det_i,degree,N_int)
print*,'degree = ',degree
print*,'Finding the differences in terms of holes and particles, together with the fermionic phase '
call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
print*,'Fermionic phase for the excitation from ref_bitmask to det_i'
print*,phase
print*,'put the excitation information in a human readable format'
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
print*,'s1',s1
print*,'h1,p1 = ',h1,p1
print*,'s2',s2
print*,'h2,p2 = ',h2,p2
print*,''
print*,'Finding the occupancy of det_i'
integer, allocatable :: occ(:,:)
integer :: n_occ_ab(2)
allocate(occ(N_int*bit_kind_size,2))
call bitstring_to_list_ab(det_i, occ, n_occ_ab, N_int)
print*,'alpha electrons orbital occupancy'
do i = 1, n_occ_ab(1) ! browsing the alpha electrons
print*,occ(i,1)
enddo
print*,'beta electrons orbital occupancy'
do i = 1, n_occ_ab(2) ! browsing the beta electrons
print*,occ(i,2)
enddo
end
subroutine example_determinants_psi_det
use bitmasks ! you need to include the bitmasks_module.f90 features
implicit none
BEGIN_DOC
! subroutine that illustrates the main features available in determinants using the psi_det/psi_coef
END_DOC
read_wf = .True.
touch read_wf
! you force the wave function to be set to the one in the EZFIO folder
call routine_example_psi_det
end
subroutine routine_example_psi_det
use bitmasks ! you need to include the bitmasks_module.f90 features
implicit none
BEGIN_DOC
! subroutine that illustrates the main features available in determinants using many determinants
END_DOC
integer :: i,j
integer, allocatable :: degree_list(:)
integer, allocatable :: idx(:)
allocate(degree_list(N_det),idx(0:N_det))
print*,'Number of determinants in the wave function'
print*,'N_det = ',N_det
print*,''
print*,'Printing in a human readable format all Slater determinants '
do i = 1, N_det
call debug_det(psi_det(1,1,i),N_int)
enddo
print*,''
print*,'Number of states computed '
print*,'N_states = ',N_states
print*,'Printing the coefficients for all states for all Slater determinants '
do j = 1, N_states
print*,'State = ',j
do i = 1, N_det
write(*,'(I9,X,F16.10)')i,psi_coef(i,j)
enddo
enddo
print*,''
print*,'Finding the connection through a bielectronic operator in the wave function'
print*,'You want to know the connections of the first determinant '
! wave function determinant exc degree list
call get_excitation_degree_vector( psi_det , psi_det(1,1,1),degree_list,N_int,N_det,idx)
double precision :: hij
double precision, allocatable :: i_H_psi(:)
allocate(i_H_psi(N_states))
i_H_psi = 0.d0
print*,'Computing <psi_det(1) | H | psi_det > = \sum_I c_I <psi_det(1)| H | psi_det(I)>'
do i = 1, idx(0) ! number of Slater determinants connected to the first one
print*,'Determinant connected'
call debug_det(psi_det(1,1,idx(i)),N_int)
print*,'excitation degree = ',degree_list(i)
call i_H_j(psi_det(1,1,1) , psi_det(1,1,idx(i)),hij,N_int)
do j = 1, N_states
i_H_psi(j) += hij * psi_coef(idx(i),j)
enddo
enddo
print*,'i_H_psi = ',i_H_psi
end

View File

@ -1,16 +0,0 @@
subroutine apply_mono(i_hole,i_particle,ispin_excit,key_in,Nint)
implicit none
integer, intent(in) :: i_hole,i_particle,ispin_excit,Nint
integer(bit_kind), intent(inout) :: key_in(Nint,2)
integer :: k,j
use bitmasks
! hole
k = shiftr(i_hole-1,bit_kind_shift)+1
j = i_hole-shiftl(k-1,bit_kind_shift)-1
key_in(k,ispin_excit) = ibclr(key_in(k,ispin_excit),j)
k = shiftr(i_particle-1,bit_kind_shift)+1
j = i_particle-shiftl(k-1,bit_kind_shift)-1
key_in(k,ispin_excit) = ibset(key_in(k,ispin_excit),j)
end

View File

@ -0,0 +1,26 @@
BEGIN_PROVIDER [ double precision, psi_energy_h_core, (N_states) ]
implicit none
integer :: i
integer :: j,k
double precision :: tmp(mo_tot_num,mo_tot_num),mono_ints(mo_tot_num,mo_tot_num)
BEGIN_DOC
! psi_energy_h_core = <Psi| h_{core} |Psi>
! computed using the one_body_dm_mo_alpha+one_body_dm_mo_beta and mo_mono_elec_integral
END_DOC
psi_energy_h_core = 0.d0
do i = 1, N_states
do j = 1, mo_tot_num
do k = 1, mo_tot_num
psi_energy_h_core(i) += mo_mono_elec_integral(k,j) * (one_body_dm_mo_alpha(k,j,i) + one_body_dm_mo_beta(k,j,i))
enddo
enddo
enddo
double precision :: accu
do i = 1, N_states
do j = 1, mo_tot_num
accu += one_body_dm_mo_alpha(j,j,i) + one_body_dm_mo_beta(j,j,i)
enddo
accu = (elec_alpha_num + elec_beta_num ) / accu
psi_energy_h_core(i) = psi_energy_h_core(i) * accu
enddo
END_PROVIDER

View File

@ -0,0 +1,37 @@
[exchange_functional]
type: character*(32)
doc: name of the exchange functional
interface: ezfio, provider, ocaml
default: short_range_LDA
[correlation_functional]
type: character*(32)
doc: name of the correlation functional
interface: ezfio, provider, ocaml
default: short_range_LDA
[HF_exchange]
type: double precision
doc: Percentage of HF exchange in the DFT model
interface: ezfio,provider,ocaml
default: 0.
[mu_erf]
type: double precision
doc: cutting of the interaction in the range separated model
interface: ezfio,provider,ocaml
default: 0.5
ezfio_name: mu_erf
[density_for_dft]
type: character*(32)
doc: Type of density used for DFT calculation. If WFT it uses the density of the WFT stored in terms of determinants. If input_density it uses the one-body dm stored in data_.../ . If damping_rs_dft it uses the damping density between WFT and input_density
interface: ezfio, provider, ocaml
default: WFT
[damping_for_rs_dft]
type: double precision
doc: damping factor for the density used in RSFT.
interface: ezfio,provider,ocaml
default: 0.5

1
src/dft_keywords/NEED Normal file
View File

@ -0,0 +1 @@
ezfio_files

View File

@ -0,0 +1,7 @@
============
DFT Keywords
============
This module contains all keywords which are related to a DFT calculation

View File

@ -0,0 +1,19 @@
BEGIN_PROVIDER [ character*(32), DFT_TYPE]
implicit none
BEGIN_DOC
! defines the type of DFT applied: LDA, GGA etc ...
END_DOC
logical :: is_lda
if(correlation_functional.eq."None")then
is_lda = (index(exchange_functional,"LDA") .ne. 0)
else if(exchange_functional.eq."None")then
is_lda = (index(correlation_functional,"LDA") .ne. 0)
else
is_lda = (index(correlation_functional,"LDA") .ne. 0) .and. (index(exchange_functional,"LDA") .ne. 0)
endif
if(is_lda)then
DFT_TYPE = "LDA"
else
DFT_TYPE = "GGA"
endif
END_PROVIDER

View File

@ -0,0 +1,6 @@
dft_keywords
determinants
ao_basis
mo_basis
becke_numerical_grid
density_for_dft

View File

@ -0,0 +1,77 @@
BEGIN_PROVIDER[double precision, aos_in_r_array, (ao_num,n_points_final_grid)]
&BEGIN_PROVIDER[double precision, aos_in_r_array_transp, (n_points_final_grid,ao_num)]
implicit none
BEGIN_DOC
! aos_in_r_array(i,j) = value of the ith ao on the jth grid point
!
! aos_in_r_array_transp(i,j) = value of the jth ao on the ith grid point
END_DOC
integer :: i,j
double precision :: aos_array(ao_num), r(3)
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
call give_all_aos_at_r(r,aos_array)
do j = 1, ao_num
aos_in_r_array(j,i) = aos_array(j)
aos_in_r_array_transp(i,j) = aos_array(j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, aos_grad_in_r_array, (ao_num,n_points_final_grid,3)]
&BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp, (n_points_final_grid,ao_num,3)]
implicit none
BEGIN_DOC
! aos_grad_in_r_array(i,j,k) = value of the kth component of the gradient of ith ao on the jth grid point
!
! aos_grad_in_r_array_transp(i,j,k) = value of the kth component of the gradient of jth ao on the ith grid point
!
! k = 1 : x, k= 2, y, k 3, z
END_DOC
integer :: i,j,m
double precision :: aos_array(ao_num), r(3)
double precision :: aos_grad_array(3,ao_num)
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
call give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
do m = 1, 3
do j = 1, ao_num
aos_grad_in_r_array(j,i,m) = aos_grad_array(m,j)
aos_grad_in_r_array_transp(i,j,m) = aos_grad_array(m,j)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, aos_lapl_in_r_array, (ao_num,n_points_final_grid,3)]
&BEGIN_PROVIDER[double precision, aos_lapl_in_r_array_transp, (n_points_final_grid,ao_num,3)]
implicit none
BEGIN_DOC
! aos_lapl_in_r_array(i,j,k) = value of the kth component of the laplacian of ith ao on the jth grid point
!
! aos_lapl_in_r_array_transp(i,j,k) = value of the kth component of the laplacian of jth ao on the ith grid point
!
! k = 1 : x, k= 2, y, k 3, z
END_DOC
integer :: i,j,m
double precision :: aos_array(ao_num), r(3)
double precision :: aos_grad_array(ao_num,3)
double precision :: aos_lapl_array(ao_num,3)
do m = 1, 3
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
do j = 1, ao_num
aos_lapl_in_r_array(j,i,m) = aos_lapl_array(j,m)
aos_lapl_in_r_array_transp(i,j,m) = aos_lapl_array(j,m)
enddo
enddo
enddo
END_PROVIDER

View File

@ -0,0 +1,297 @@
subroutine dm_dft_alpha_beta_at_r(r,dm_a,dm_b)
implicit none
BEGIN_DOC
! input: r(1) ==> r(1) = x, r(2) = y, r(3) = z
! output : dm_a = alpha density evaluated at r(3)
! output : dm_b = beta density evaluated at r(3)
END_DOC
double precision, intent(in) :: r(3)
double precision, intent(out) :: dm_a(N_states),dm_b(N_states)
integer :: istate
double precision :: aos_array(ao_num),aos_array_bis(ao_num),u_dot_v
call give_all_aos_at_r(r,aos_array)
do istate = 1, N_states
aos_array_bis = aos_array
! alpha density
call dgemv('N',ao_num,ao_num,1.d0,one_body_dm_alpha_ao_for_dft(1,1,istate),ao_num,aos_array,1,0.d0,aos_array_bis,1)
dm_a(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
! beta density
aos_array_bis = aos_array
call dgemv('N',ao_num,ao_num,1.d0,one_body_dm_beta_ao_for_dft(1,1,istate),ao_num,aos_array,1,0.d0,aos_array_bis,1)
dm_b(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
enddo
end
subroutine dm_dft_alpha_beta_and_all_aos_at_r(r,dm_a,dm_b,aos_array)
BEGIN_DOC
! input: r(1) ==> r(1) = x, r(2) = y, r(3) = z
! output : dm_a = alpha density evaluated at r
! output : dm_b = beta density evaluated at r
! output : aos_array(i) = ao(i) evaluated at r
END_DOC
implicit none
double precision, intent(in) :: r(3)
double precision, intent(out) :: dm_a(N_states),dm_b(N_states)
double precision, intent(out) :: aos_array(ao_num)
integer :: istate
double precision :: aos_array_bis(ao_num),u_dot_v
call give_all_aos_at_r(r,aos_array)
do istate = 1, N_states
aos_array_bis = aos_array
! alpha density
call dsymv('U',ao_num,1.d0,one_body_dm_alpha_ao_for_dft(1,1,istate),size(one_body_dm_alpha_ao_for_dft,1),aos_array,1,0.d0,aos_array_bis,1)
dm_a(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
! beta density
aos_array_bis = aos_array
call dsymv('U',ao_num,1.d0,one_body_dm_beta_ao_for_dft(1,1,istate),size(one_body_dm_beta_ao_for_dft,1),aos_array,1,0.d0,aos_array_bis,1)
dm_b(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
enddo
end
subroutine density_and_grad_alpha_beta_and_all_aos_and_grad_aos_at_r(r,dm_a,dm_b, grad_dm_a, grad_dm_b, aos_array, grad_aos_array)
implicit none
BEGIN_DOC
! input : r(1) ==> r(1) = x, r(2) = y, r(3) = z
! output : dm_a = alpha density evaluated at r
! : dm_b = beta density evaluated at r
! : aos_array(i) = ao(i) evaluated at r
! : grad_dm_a(1) = X gradient of the alpha density evaluated in r
! : grad_dm_a(1) = X gradient of the beta density evaluated in r
! : grad_aos_array(1) = X gradient of the aos(i) evaluated at r
END_DOC
double precision, intent(in) :: r(3)
double precision, intent(out) :: dm_a(N_states),dm_b(N_states)
double precision, intent(out) :: grad_dm_a(3,N_states),grad_dm_b(3,N_states)
double precision, intent(out) :: grad_aos_array(3,ao_num)
integer :: i,j,istate
double precision :: aos_array(ao_num),aos_array_bis(ao_num),u_dot_v
double precision :: aos_grad_array(ao_num,3), aos_grad_array_bis(ao_num,3)
call give_all_aos_and_grad_at_r(r,aos_array,grad_aos_array)
do i = 1, ao_num
do j = 1, 3
aos_grad_array(i,j) = grad_aos_array(j,i)
enddo
enddo
do istate = 1, N_states
! alpha density
! aos_array_bis = \rho_ao * aos_array
call dsymv('U',ao_num,1.d0,one_body_dm_alpha_ao_for_dft(1,1,istate),size(one_body_dm_alpha_ao_for_dft,1),aos_array,1,0.d0,aos_array_bis,1)
dm_a(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
! grad_dm(1) = \sum_i aos_grad_array(i,1) * aos_array_bis(i)
grad_dm_a(1,istate) = u_dot_v(aos_grad_array(1,1),aos_array_bis,ao_num)
grad_dm_a(2,istate) = u_dot_v(aos_grad_array(1,2),aos_array_bis,ao_num)
grad_dm_a(3,istate) = u_dot_v(aos_grad_array(1,3),aos_array_bis,ao_num)
grad_dm_a *= 2.d0
! aos_grad_array_bis = \rho_ao * aos_grad_array
! beta density
call dsymv('U',ao_num,1.d0,one_body_dm_beta_ao_for_dft(1,1,istate),size(one_body_dm_beta_ao_for_dft,1),aos_array,1,0.d0,aos_array_bis,1)
dm_b(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
! grad_dm(1) = \sum_i aos_grad_array(i,1) * aos_array_bis(i)
grad_dm_b(1,istate) = u_dot_v(aos_grad_array(1,1),aos_array_bis,ao_num)
grad_dm_b(2,istate) = u_dot_v(aos_grad_array(1,2),aos_array_bis,ao_num)
grad_dm_b(3,istate) = u_dot_v(aos_grad_array(1,3),aos_array_bis,ao_num)
grad_dm_b *= 2.d0
! aos_grad_array_bis = \rho_ao * aos_grad_array
enddo
end
BEGIN_PROVIDER [double precision, one_body_dm_mo_alpha_at_grid_points, (n_points_integration_angular,n_points_radial_grid,nucl_num,N_states) ]
&BEGIN_PROVIDER [double precision, one_body_dm_mo_beta_at_grid_points, (n_points_integration_angular,n_points_radial_grid,nucl_num,N_states) ]
implicit none
integer :: i,j,k,l,m,istate
double precision :: contrib
double precision :: r(3)
double precision :: aos_array(ao_num),mos_array(mo_tot_num)
do j = 1, nucl_num
do k = 1, n_points_radial_grid -1
do l = 1, n_points_integration_angular
do istate = 1, N_States
one_body_dm_mo_alpha_at_grid_points(l,k,j,istate) = 0.d0
one_body_dm_mo_beta_at_grid_points(l,k,j,istate) = 0.d0
enddo
r(1) = grid_points_per_atom(1,l,k,j)
r(2) = grid_points_per_atom(2,l,k,j)
r(3) = grid_points_per_atom(3,l,k,j)
double precision :: dm_a,dm_b
call dm_dft_alpha_beta_at_r(r,dm_a,dm_b)
one_body_dm_mo_alpha_at_grid_points(l,k,j,1) = dm_a
one_body_dm_mo_beta_at_grid_points(l,k,j,1) = dm_b
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, one_body_dm_alpha_at_r, (n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, one_body_dm_beta_at_r, (n_points_final_grid,N_states) ]
implicit none
BEGIN_DOC
! one_body_dm_alpha_at_r(i,istate) = n_alpha(r_i,istate)
! one_body_dm_beta_at_r(i,istate) = n_beta(r_i,istate)
! where r_i is the ith point of the grid and istate is the state number
END_DOC
integer :: i,istate
double precision :: r(3)
double precision, allocatable :: dm_a(:),dm_b(:)
allocate(dm_a(N_states),dm_b(N_states))
do istate = 1, N_states
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
call dm_dft_alpha_beta_at_r(r,dm_a,dm_b)
one_body_dm_alpha_at_r(i,istate) = dm_a(istate)
one_body_dm_beta_at_r(i,istate) = dm_b(istate)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, one_body_dm_alpha_and_grad_at_r, (4,n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, one_body_dm_beta_and_grad_at_r, (4,n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, one_body_grad_2_dm_alpha_at_r, (n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, one_body_grad_2_dm_beta_at_r, (n_points_final_grid,N_states) ]
BEGIN_DOC
! one_body_dm_alpha_and_grad_at_r(1,i,i_state) = d\dx n_alpha(r_i,istate)
! one_body_dm_alpha_and_grad_at_r(2,i,i_state) = d\dy n_alpha(r_i,istate)
! one_body_dm_alpha_and_grad_at_r(3,i,i_state) = d\dz n_alpha(r_i,istate)
! one_body_dm_alpha_and_grad_at_r(4,i,i_state) = n_alpha(r_i,istate)
! one_body_grad_2_dm_alpha_at_r(i,istate) = d\dx n_alpha(r_i,istate)^2 + d\dy n_alpha(r_i,istate)^2 + d\dz n_alpha(r_i,istate)^2
! where r_i is the ith point of the grid and istate is the state number
END_DOC
implicit none
integer :: i,j,k,l,m,istate
double precision :: contrib
double precision :: r(3)
double precision, allocatable :: aos_array(:),grad_aos_array(:,:)
double precision, allocatable :: dm_a(:),dm_b(:), dm_a_grad(:,:), dm_b_grad(:,:)
allocate(dm_a(N_states),dm_b(N_states), dm_a_grad(3,N_states), dm_b_grad(3,N_states))
allocate(aos_array(ao_num),grad_aos_array(3,ao_num))
do istate = 1, N_states
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
!!!! Works also with the ao basis
call density_and_grad_alpha_beta_and_all_aos_and_grad_aos_at_r(r,dm_a,dm_b, dm_a_grad, dm_b_grad, aos_array, grad_aos_array)
one_body_dm_alpha_and_grad_at_r(1,i,istate) = dm_a_grad(1,istate)
one_body_dm_alpha_and_grad_at_r(2,i,istate) = dm_a_grad(2,istate)
one_body_dm_alpha_and_grad_at_r(3,i,istate) = dm_a_grad(3,istate)
one_body_dm_alpha_and_grad_at_r(4,i,istate) = dm_a(istate)
one_body_grad_2_dm_alpha_at_r(i,istate) = dm_a_grad(1,istate) * dm_a_grad(1,istate) + dm_a_grad(2,istate) * dm_a_grad(2,istate) + dm_a_grad(3,istate) * dm_a_grad(3,istate)
one_body_dm_beta_and_grad_at_r(1,i,istate) = dm_b_grad(1,istate)
one_body_dm_beta_and_grad_at_r(2,i,istate) = dm_b_grad(2,istate)
one_body_dm_beta_and_grad_at_r(3,i,istate) = dm_b_grad(3,istate)
one_body_dm_beta_and_grad_at_r(4,i,istate) = dm_b(istate)
one_body_grad_2_dm_beta_at_r(i,istate) = dm_b_grad(1,istate) * dm_b_grad(1,istate) + dm_b_grad(2,istate) * dm_b_grad(2,istate) + dm_b_grad(3,istate) * dm_b_grad(3,istate)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, one_body_dm_mo_alpha_and_grad_at_grid_points, (4,n_points_integration_angular,n_points_radial_grid,nucl_num,N_states) ]
&BEGIN_PROVIDER [double precision, one_body_dm_mo_beta_and_grad_at_grid_points, (4,n_points_integration_angular,n_points_radial_grid,nucl_num,N_states) ]
BEGIN_DOC
! one_body_dm_mo_alpha_and_grad_at_grid_points(1,.....,i_state) = d\dx \rho_{\alpha}^{\istate}(r)
! one_body_dm_mo_alpha_and_grad_at_grid_points(2,.....,i_state) = d\dy \rho_{\alpha}^{\istate}(r)
! one_body_dm_mo_alpha_and_grad_at_grid_points(3,.....,i_state) = d\dz \rho_{\alpha}^{\istate}(r)
! one_body_dm_mo_alpha_and_grad_at_grid_points(4,.....,i_state) = \rho_{\alpha}^{\istate}(r)
END_DOC
implicit none
integer :: i,j,k,l,m,istate
double precision :: contrib
double precision :: r(3)
double precision :: aos_array(ao_num),grad_aos_array(3,ao_num)
do istate = 1, N_States
do j = 1, nucl_num
do k = 1, n_points_radial_grid -1
do l = 1, n_points_integration_angular
do m = 1, 4
one_body_dm_mo_alpha_and_grad_at_grid_points(m,l,k,j,istate) = 0.d0
one_body_dm_mo_beta_and_grad_at_grid_points(m,l,k,j,istate) = 0.d0
enddo
r(1) = grid_points_per_atom(1,l,k,j)
r(2) = grid_points_per_atom(2,l,k,j)
r(3) = grid_points_per_atom(3,l,k,j)
!!!!! Works also with the ao basis
double precision :: dm_a(N_states),dm_b(N_states), dm_a_grad(3,N_states), dm_b_grad(3,N_states)
call density_and_grad_alpha_beta_and_all_aos_and_grad_aos_at_r(r,dm_a,dm_b, dm_a_grad, dm_b_grad, aos_array, grad_aos_array)
one_body_dm_mo_alpha_and_grad_at_grid_points(1,l,k,j,istate) = dm_a_grad(1,istate)
one_body_dm_mo_alpha_and_grad_at_grid_points(2,l,k,j,istate) = dm_a_grad(2,istate)
one_body_dm_mo_alpha_and_grad_at_grid_points(3,l,k,j,istate) = dm_a_grad(3,istate)
one_body_dm_mo_alpha_and_grad_at_grid_points(4,l,k,j,istate) = dm_a(istate)
one_body_dm_mo_beta_and_grad_at_grid_points(1,l,k,j,istate) = dm_b_grad(1,istate)
one_body_dm_mo_beta_and_grad_at_grid_points(2,l,k,j,istate) = dm_b_grad(2,istate)
one_body_dm_mo_beta_and_grad_at_grid_points(3,l,k,j,istate) = dm_b_grad(3,istate)
one_body_dm_mo_beta_and_grad_at_grid_points(4,l,k,j,istate) = dm_b(istate)
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, integral_density_alpha_knowles_becke_per_atom, (nucl_num)]
&BEGIN_PROVIDER [ double precision, integral_density_beta_knowles_becke_per_atom, (nucl_num)]
implicit none
double precision :: accu
integer :: i,j,k,l,istate
double precision :: x
double precision :: integrand(n_points_integration_angular), weights(n_points_integration_angular)
double precision :: f_average_angular_alpha,f_average_angular_beta
double precision :: derivative_knowles_function,knowles_function
! Run over all nuclei in order to perform the Voronoi partition
! according ot equation (6) of the paper of Becke (JCP, (88), 1988)
! Here the m index is referred to the w_m(r) weight functions of equation (22)
! Run over all points of integrations : there are
! n_points_radial_grid (i) * n_points_integration_angular (k)
do j = 1, nucl_num
integral_density_alpha_knowles_becke_per_atom(j) = 0.d0
integral_density_beta_knowles_becke_per_atom(j) = 0.d0
do i = 1, n_points_radial_grid-1
! Angular integration over the solid angle Omega for a FIXED angular coordinate "r"
f_average_angular_alpha = 0.d0
f_average_angular_beta = 0.d0
do istate = 1, N_states
do k = 1, n_points_integration_angular
f_average_angular_alpha += weights_angular_points(k) * one_body_dm_mo_alpha_at_grid_points(k,i,j,istate) * weight_functions_at_grid_points(k,i,j)
f_average_angular_beta += weights_angular_points(k) * one_body_dm_mo_beta_at_grid_points(k,i,j,istate) * weight_functions_at_grid_points(k,i,j)
enddo
enddo
!
x = grid_points_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1]
double precision :: contrib_integration
! print*,m_knowles
contrib_integration = derivative_knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x) &
*knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x)**2
integral_density_alpha_knowles_becke_per_atom(j) += contrib_integration *f_average_angular_alpha
integral_density_beta_knowles_becke_per_atom(j) += contrib_integration *f_average_angular_beta
enddo
integral_density_alpha_knowles_becke_per_atom(j) *= dr_radial_integral
integral_density_beta_knowles_becke_per_atom(j) *= dr_radial_integral
enddo
END_PROVIDER

View File

@ -0,0 +1,56 @@
BEGIN_PROVIDER[double precision, mos_in_r_array, (mo_tot_num,n_points_final_grid)]
&BEGIN_PROVIDER[double precision, mos_in_r_array_transp,(n_points_final_grid,mo_tot_num)]
implicit none
BEGIN_DOC
! mos_in_r_array(i,j) = value of the ith mo on the jth grid point
!
! mos_in_r_array_transp(i,j) = value of the jth mo on the ith grid point
END_DOC
integer :: i,j
double precision :: mos_array(mo_tot_num), r(3)
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
call give_all_mos_at_r(r,mos_array)
do j = 1, mo_tot_num
mos_in_r_array(j,i) = mos_array(j)
mos_in_r_array_transp(i,j) = mos_array(j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, mos_grad_in_r_array,(mo_tot_num,n_points_final_grid,3)]
implicit none
BEGIN_DOC
! mos_grad_in_r_array(i,j,k) = value of the kth component of the gradient of ith mo on the jth grid point
!
! mos_grad_in_r_array_transp(i,j,k) = value of the kth component of the gradient of jth mo on the ith grid point
!
! k = 1 : x, k= 2, y, k 3, z
END_DOC
integer :: m
mos_grad_in_r_array = 0.d0
do m=1,3
call dgemm('N','N',mo_tot_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_tot_num,aos_grad_in_r_array(1,1,m),ao_num,0.d0,mos_grad_in_r_array(1,1,m),mo_tot_num)
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, mos_lapl_in_r_array,(mo_tot_num,n_points_final_grid,3)]
implicit none
BEGIN_DOC
! mos_lapl_in_r_array(i,j,k) = value of the kth component of the laplacian of ith mo on the jth grid point
!
! mos_lapl_in_r_array_transp(i,j,k) = value of the kth component of the laplacian of jth mo on the ith grid point
!
! k = 1 : x, k= 2, y, k 3, z
END_DOC
integer :: m
mos_lapl_in_r_array = 0.d0
do m=1,3
call dgemm('N','N',mo_tot_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_tot_num,aos_lapl_in_r_array(1,1,m),ao_num,0.d0,mos_lapl_in_r_array(1,1,m),mo_tot_num)
enddo
END_PROVIDER

View File

@ -0,0 +1,4 @@
density_for_dft
dft_utils_on_grid
integrals_bielec
integrals_bielec_erf

View File

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

@ -0,0 +1,125 @@
BEGIN_PROVIDER [double precision, potential_x_alpha_ao,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, potential_x_beta_ao,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, potential_c_alpha_ao,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, potential_c_beta_ao,(ao_num,ao_num,N_states)]
implicit none
BEGIN_DOC
! alpha/beta exchange/correlation potentials on the AO basis
END_DOC
if(trim(exchange_functional)=="short_range_LDA")then
potential_x_alpha_ao = potential_x_alpha_ao_LDA
potential_x_beta_ao = potential_x_beta_ao_LDA
else if(exchange_functional.EQ."short_range_PBE")then
potential_x_alpha_ao = potential_x_alpha_ao_PBE
potential_x_beta_ao = potential_x_beta_ao_PBE
else if(exchange_functional.EQ."None")then
potential_x_alpha_ao = 0.d0
potential_x_beta_ao = 0.d0
else
print*, 'Exchange functional required does not exist ...'
print*,'exchange_functional',exchange_functional
stop
endif
if(trim(correlation_functional)=="short_range_LDA")then
potential_c_alpha_ao = potential_c_alpha_ao_LDA
potential_c_beta_ao = potential_c_beta_ao_LDA
else if(correlation_functional.EQ."short_range_PBE")then
potential_c_alpha_ao = potential_c_alpha_ao_PBE
potential_c_beta_ao = potential_c_beta_ao_PBE
else if(correlation_functional.EQ."None")then
potential_c_alpha_ao = 0.d0
potential_c_beta_ao = 0.d0
else
print*, 'Correlation functional required does not ecist ...'
print*,'correlation_functional',correlation_functional
stop
endif
END_PROVIDER
BEGIN_PROVIDER [double precision, potential_x_alpha_mo,(mo_tot_num,mo_tot_num,N_states)]
&BEGIN_PROVIDER [double precision, potential_x_beta_mo,(mo_tot_num,mo_tot_num,N_states)]
&BEGIN_PROVIDER [double precision, potential_c_alpha_mo,(mo_tot_num,mo_tot_num,N_states)]
&BEGIN_PROVIDER [double precision, potential_c_beta_mo,(mo_tot_num,mo_tot_num,N_states)]
implicit none
BEGIN_DOC
! alpha/beta exchange/correlation potentials on the MO basis
END_DOC
integer :: istate
do istate = 1, N_states
call ao_to_mo( &
potential_x_alpha_ao(1,1,istate), &
size(potential_x_alpha_ao,1), &
potential_x_alpha_mo(1,1,istate), &
size(potential_x_alpha_mo,1) &
)
call ao_to_mo( &
potential_x_beta_ao(1,1,istate), &
size(potential_x_beta_ao,1), &
potential_x_beta_mo(1,1,istate), &
size(potential_x_beta_mo,1) &
)
call ao_to_mo( &
potential_c_alpha_ao(1,1,istate), &
size(potential_c_alpha_ao,1), &
potential_c_alpha_mo(1,1,istate), &
size(potential_c_alpha_mo,1) &
)
call ao_to_mo( &
potential_c_beta_ao(1,1,istate), &
size(potential_c_beta_ao,1), &
potential_c_beta_mo(1,1,istate), &
size(potential_c_beta_mo,1) &
)
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, energy_x, (N_states)]
&BEGIN_PROVIDER [double precision, energy_c, (N_states)]
implicit none
if(trim(exchange_functional)=="short_range_LDA")then
energy_x = energy_x_LDA
energy_x = energy_x_LDA
else if(exchange_functional.EQ."short_range_PBE")then
energy_x = energy_x_PBE
energy_x = energy_x_PBE
else if(exchange_functional.EQ."None")then
energy_x = 0.d0
energy_x = 0.d0
else
print*, 'Exchange functional required does not exist ...'
print*,'exchange_functional',exchange_functional
stop
endif
if(trim(correlation_functional)=="short_range_LDA")then
energy_c = energy_c_LDA
energy_c = energy_c_LDA
else if(correlation_functional.EQ."short_range_PBE")then
energy_c = energy_c_PBE
energy_c = energy_c_PBE
else if(correlation_functional.EQ."None")then
energy_c = 0.d0
energy_c = 0.d0
else
print*, 'Correlation functional required does not ecist ...'
print*,'correlation_functional',correlation_functional
stop
endif
END_PROVIDER

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,477 @@
subroutine ec_pbe_sr(mu,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,ec,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo)
!************************************************************************
! Short-range PBE correlation energy functional for erf interaction
!
!
!************************************************************************
include 'constants.include.F'
implicit none
! input
double precision, intent(in) :: rhoc,rhoo,mu
double precision, intent(in) :: sigmacc,sigmaco,sigmaoo
! output
double precision, intent(out) :: ec
double precision, intent(out) :: vrhoc,vrhoo
double precision, intent(out) :: vsigmacc,vsigmaco,vsigmaoo
! local
double precision tol
parameter(tol=1d-12)
character(len=30) namedummy
double precision eccerflda
double precision vrhoccerflda
double precision vrhoocerflda
double precision ecclda
double precision vrhocclda
double precision vrhooclda
integer i,igrad
double precision rho,drho2,rhoa,rhob
double precision ecerflda,decerfldadrho
double precision eclda,decldadrho
double precision ecerfpbe,decerfpbedrho,decerfpbedrhoo
double precision decerfpbeddrho2
double precision arglog,arglogs,arglogss,alpha,beta,betas,gamma
double precision Aa,Ab,Ac,Aas,tq,tqs,tqss,decerfpur,decpur
double precision t1,t2,t3,t4,t5,t6,t7,t8,t9,t10
double precision t11,t12,t13,t14,t15,t16,t17,t18,t19
double precision zeta,phi,phi2,phi3,phi4,phis,arglogsc
double precision dlogarglog
double precision, parameter :: f13=0.333333333333333d0
! Parameter of the modified interaction
ec = 0.d0
vrhoc = 0.d0
vrhoo = 0.d0
vsigmacc = 0.d0
vsigmaco = 0.d0
vsigmaoo = 0.d0
! First-type gradient functional
igrad=1
alpha=2.78d0
gamma=3.1091d-2
! test on density
if (dabs(rhoc).lt.tol) return
double precision :: vc_a,vc_b
! Spin polarisation
rhoa=max((rhoc+rhoo)*.5d0,1.0d-15)
rhob=max((rhoc-rhoo)*.5d0,1.0d-15)
call ec_lda_sr(mu,rhoa,rhob,eccerflda,vc_a,vc_b)
ecerflda = eccerflda
vrhoccerflda = 0.5d0 * (vc_a + vc_b)
vrhoocerflda = 0.5d0 * (vc_a - vc_b)
! Density
rho = rhoc
! Square of density gradient
drho2 = sigmacc
zeta = (rhoa-rhob)/(rhoa+rhob)
! LDA energy density
double precision :: vc_a_lda,vc_b_lda
call ec_lda(rhoa,rhob,ecclda,vc_a_lda,vc_b_lda)
eclda = ecclda
!decldadrho = 0.5d0 * (vc_a_lda+vc_b_lda)
!decldadrho = 0.5d0 * (vc_a_lda-vc_b_lda)
if ((ecerflda/eclda).le.0d0) then
beta=0d0
else
beta=6.6725d-2*(ecerflda/eclda)**alpha
endif
phi=((1d0+zeta)**(2d0/3d0)+(1d0-zeta)**(2d0/3d0))/2d0
phi2=phi*phi
phi3=phi2*phi
phi4=phi3*phi
tq=drho2*6.346820607d-2*rho**(-7d0/3d0)/phi2
Ab=dexp(-ecerflda/(rho*gamma*phi3))-1d0
if (dabs(Ab).le.dabs(beta*tol)) then
ecerfpbe=ecerflda
else
Aa=beta/(gamma*Ab)
Ac=1d0+Aa*tq+Aa**2*tq**2
if (Aa.lt.tol) Aa=tol
arglog=1d0+beta*(1d0-1d0/Ac)/(gamma*Aa)
ecerfpbe=ecerflda+rho*phi3*gamma*dlog(arglog)
end if
ec = ecerfpbe
! if(ldebug) write(*,*)"ecerfpbe=",ecerfpbe
! Derive
! LDA energy density derivative
decerfldadrho = vrhoccerflda
decldadrho = 0.5d0 * (vc_a_lda+vc_b_lda)
decerfpur=(decerfldadrho-ecerflda/rho)/rho
decpur=(decldadrho-eclda/rho)/rho
betas=alpha*beta*(decerfpur*rho/ecerflda-decpur*rho/eclda)
phis=((rhoa - rhob)*((rhoa/(rhoa + rhob))**f13 - (rhob/(rhoa + rhob))**f13))/(3d0*2d0**f13*(rhoa/(rhoa + rhob))**f13*(rhob/(rhoa + rhob))**f13*(rhoa + rhob)**2)
if (dabs(Ab).le.dabs(beta*tol)) then
decerfpbedrho=decerfldadrho
else
Aas=betas/(gamma*Ab)+Aa*(1d0+1d0/Ab)*(decerfpur/phi3-3d0*phis*ecerflda/(rho*phi4))/gamma
tqs=-7d0*tq/(3d0*rho)-2d0*tq*phis/phi
arglogs=betas*tq*(1d0+Aa*tq)/(Ac*gamma)+beta*tqs*(1d0+Aa*tq)/(Ac*gamma)-beta*tq*Aa*tq*(Aas*tq+Aa*tqs)*(2d0+Aa*tq)/(Ac**2*gamma)
dlogarglog=dlog(arglog)
decerfpbedrho=decerfldadrho+gamma*(phi3*dlogarglog+3d0*rho*phis*phi2*dlogarglog+rho*phi3*arglogs/arglog)
end if
if (dabs(Ab).le.dabs(beta*tol)) then
decerfpbeddrho2=0.0d0
else
arglogsc=Ab*(Aa+2d0*Aa*Aa*tq)/(Ac*Ac)
tqss=6.346820607d-2*rho**(-7d0/3d0)/phi2
arglogss=tqss*arglogsc
decerfpbeddrho2=rho*gamma*phi3*arglogss/arglog
end if
! LDA energy density derivative
decerfldadrho = vrhoocerflda
decldadrho = 0.5d0 * (vc_a_lda-vc_b_lda)
decerfpur=decerfldadrho/rho
decpur=decldadrho/rho
betas=alpha*beta*(decerfpur*rho/ecerflda-decpur*rho/eclda)
phis=(rhob*(rhoa/(rhoa + rhob))**(2d0*f13)-rhoa*(rhob/(rhoa + rhob))**(2d0*f13))/(3d0*2d0**f13*rhoa*rhob)
if (dabs(Ab).le.dabs(beta*tol)) then
decerfpbedrhoo=decerfldadrho
else
Aas=betas/(gamma*Ab)+Aa*(1d0+1d0/Ab)*(decerfpur/phi3-3d0*phis*ecerflda/(rho*phi4))/gamma
tqs=-2d0*tq*phis/phi
arglogs=betas*tq*(1d0+Aa*tq)/(Ac*gamma)+beta*tqs*(1d0+Aa*tq)/(Ac*gamma)-beta*tq*Aa*tq*(Aas*tq+Aa*tqs)*(2d0+Aa*tq)/(Ac**2*gamma)
decerfpbedrhoo=decerfldadrho+gamma*(3d0*rho*phis*phi2*dlog(arglog)+rho*phi3*arglogs/arglog)
end if
! derivatives
vrhoc = vrhoc + decerfpbedrho
vrhoo = vrhoo + decerfpbedrhoo
vsigmacc = vsigmacc + decerfpbeddrho2
end
subroutine ex_pbe_sr(mu,rho_a,rho_b,grd_rho_a_2,grd_rho_b_2,grd_rho_a_b,ex,vx_rho_a,vx_rho_b,vx_grd_rho_a_2,vx_grd_rho_b_2,vx_grd_rho_a_b)
BEGIN_DOC
!rho_a = density alpha
!rho_b = density beta
!grd_rho_a_2 = (gradient rho_a)^2
!grd_rho_b_2 = (gradient rho_b)^2
!grd_rho_a_b = (gradient rho_a).(gradient rho_b)
!ex = exchange energy density at point r
!vx_rho_a = d ex / d rho_a
!vx_rho_b = d ex / d rho_b
!vx_grd_rho_a_2 = d ex / d grd_rho_a_2
!vx_grd_rho_b_2 = d ex / d grd_rho_b_2
!vx_grd_rho_a_b = d ex / d grd_rho_a_b
END_DOC
implicit none
! input
double precision, intent(in) :: mu,rho_a, rho_b
double precision, intent(in) :: grd_rho_a_2, grd_rho_b_2, grd_rho_a_b
! output
double precision, intent(out) :: ex
double precision, intent(out) :: vx_rho_a, vx_rho_b
double precision, intent(out) :: vx_grd_rho_a_2, vx_grd_rho_b_2, vx_grd_rho_a_b
! function
double precision berf
double precision dberfda
! local
double precision, parameter :: tol=1d-12
double precision, parameter :: f13=0.333333333333333d0
double precision exerflda,vxerflda_a,vxerflda_b
double precision dexerfldadrho
double precision exerfpbe_a, exerfpbe_b
double precision dexerfpbedrho_a, dexerfpbedrho_b
double precision dexerfpbeddrho2_a, dexerfpbeddrho2_b
double precision rho,drho2
double precision rho_a_2, rho_b_2
double precision t1,t2,t3,t4
double precision kappa,sq,sqs,sqss,fx,fxs,ksig
! Parameter of the modified interaction
! initialization
ex=0.d0
vx_rho_a=0.d0
vx_rho_b=0.d0
vx_grd_rho_a_2=0.d0
vx_grd_rho_b_2=0.d0
vx_grd_rho_a_b=0.d0
! spin scaling relation Ex[rho_a,rho_b] = (1/2) (Ex[2rho_a,2rho_a] + Ex[2rho_b,2rho_b])
! two times spin alpha density
rho = max(rho_a,tol)*2.d0
! test on density
if (rho >= tol) then
! call srLDA Ex[2*rho_a,2*rho_a]
call ex_lda_sr(mu,rho_a,rho_a,exerflda,vxerflda_a,vxerflda_b)
dexerfldadrho = (vxerflda_a + vxerflda_b)*0.5d0
! square of two times spin alpha density gradient
drho2=max(grd_rho_a_2,0d0)*4.0d0
kappa=0.804d0
sq=drho2*2.6121172985233599567768d-2*rho**(-8d0/3d0)
fx=1d0+kappa-kappa/(1d0+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq/kappa)
exerfpbe_a=exerflda*fx
! Derivatives
sqs=-8d0*sq/(3d0*rho)
fxs=kappa**2*(-1.616204596739954813d-1*mu*rho**(-4d0*f13)/3d0*dberfda(1.616204596739954813d-1*mu*rho**(-f13))*sq+berf(1.616204596739954813d-1*mu*rho**(-f13))*sqs)/(kappa+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq)**2
dexerfpbedrho_a=dexerfldadrho*fx+exerflda*fxs
sqss=2.6121172985233599567768d-2*rho**(-8d0/3d0)
dexerfpbeddrho2_a=exerflda*berf(1.616204596739954813d-1*mu*rho**(-1.d0/3.d0))*sqss*kappa**2/(kappa+berf(1.616204596739954813d-1*mu*rho**(-1.d0/3.d0))*sq)**2
endif
! two times spin beta density
rho = max(rho_b,tol)*2.d0
! test on density
if (rho >= tol) then
! call srLDA Ex[2*rho_b,2*rho_b]
call ex_lda_sr(mu,rho_b,rho_b,exerflda,vxerflda_a,vxerflda_b)
dexerfldadrho = (vxerflda_a + vxerflda_b)*0.5d0
! square of two times spin beta density gradient
drho2=max(grd_rho_b_2,0d0)*4.0d0
kappa=0.804d0
sq=drho2*2.6121172985233599567768d-2*rho**(-8d0/3d0)
fx=1d0+kappa-kappa/(1d0+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq/kappa)
exerfpbe_b=exerflda*fx
! Derivatives
sqs=-8d0*sq/(3d0*rho)
fxs=kappa**2*(-1.616204596739954813d-1*mu*rho**(-4d0*f13)/3d0*dberfda(1.616204596739954813d-1*mu*rho**(-f13))*sq+berf(1.616204596739954813d-1*mu*rho**(-f13))*sqs)/(kappa+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq)**2
dexerfpbedrho_b=dexerfldadrho*fx+exerflda*fxs
sqss=2.6121172985233599567768d-2*rho**(-8d0/3d0)
dexerfpbeddrho2_b=exerflda*berf(1.616204596739954813d-1*mu*rho**(-1.d0/3.d0))*sqss*kappa**2/(kappa+berf(1.616204596739954813d-1*mu*rho**(-1.d0/3.d0))*sq)**2
endif
ex = (exerfpbe_a+exerfpbe_b)*0.5d0
vx_rho_a = dexerfpbedrho_a
vx_rho_b = dexerfpbedrho_a
vx_grd_rho_a_2 = 2.d0*dexerfpbeddrho2_a
vx_grd_rho_b_2 = 2.d0*dexerfpbeddrho2_b
vx_grd_rho_a_b = 0.d0
end
subroutine ex_pbe_sr_only(mu,rho_a,rho_b,grd_rho_a_2,grd_rho_b_2,grd_rho_a_b,ex)
BEGIN_DOC
!rho_a = density alpha
!rho_b = density beta
!grd_rho_a_2 = (gradient rho_a)^2
!grd_rho_b_2 = (gradient rho_b)^2
!grd_rho_a_b = (gradient rho_a).(gradient rho_b)
!ex = exchange energy density at point r
END_DOC
implicit none
! input
double precision, intent(in) :: mu,rho_a, rho_b
double precision, intent(in) :: grd_rho_a_2, grd_rho_b_2, grd_rho_a_b
! output
double precision, intent(out) :: ex
! function
double precision berf
! local
double precision, parameter :: tol=1d-12
double precision, parameter :: f13=0.333333333333333d0
double precision exerflda,vxerflda_a,vxerflda_b
double precision exerfpbe_a, exerfpbe_b
double precision rho,drho2
double precision kappa,sq,fx
! initialization
ex=0.d0
! spin scaling relation Ex[rho_a,rho_b] = (1/2) (Ex[2rho_a,2rho_a] + Ex[2rho_b,2rho_b])
! two times spin alpha density
rho = max(rho_a,tol)*2.d0
! test on density
if (rho >= tol) then
! call srLDA Ex[2*rho_a,2*rho_a]
call ex_lda_sr(mu,rho_a,rho_a,exerflda,vxerflda_a,vxerflda_b)
! square of two times spin alpha density gradient
drho2=max(grd_rho_a_2,0d0)*4.0d0
kappa=0.804d0
sq=drho2*2.6121172985233599567768d-2*rho**(-8d0/3d0)
fx=1d0+kappa-kappa/(1d0+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq/kappa)
exerfpbe_a=exerflda*fx
endif
! two times spin beta density
rho = max(rho_b,tol)*2.d0
! test on density
if (rho >= tol) then
! call srLDA Ex[2*rho_b,2*rho_b]
call ex_lda_sr(mu,rho_b,rho_b,exerflda,vxerflda_a,vxerflda_b)
! square of two times spin beta density gradient
drho2=max(grd_rho_b_2,0d0)*4.0d0
kappa=0.804d0
sq=drho2*2.6121172985233599567768d-2*rho**(-8d0/3d0)
fx=1d0+kappa-kappa/(1d0+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq/kappa)
exerfpbe_b=exerflda*fx
endif
ex = (exerfpbe_a+exerfpbe_b)*0.5d0
end
subroutine ec_pbe_only(mu,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,ec)
!************************************************************************
! Short-range PBE correlation energy functional for erf interaction
!
!
!************************************************************************
include 'constants.include.F'
implicit none
! input
double precision, intent(in) :: rhoc,rhoo,mu
double precision, intent(in) :: sigmacc,sigmaco,sigmaoo
! output
double precision, intent(out) :: ec
! local
double precision tol
parameter(tol=1d-12)
character(len=30) namedummy
double precision eccerflda
double precision vrhoccerflda
double precision vrhoocerflda
double precision ecclda
double precision vrhocclda
double precision vrhooclda
integer i,igrad
double precision rho,drho2,rhoa,rhob
double precision ecerflda
double precision eclda,decldadrho
double precision ecerfpbe
double precision arglog,alpha,beta,gamma
double precision Aa,Ab,Ac,tq
double precision zeta,phi,phi2,phi3,phi4
double precision, parameter :: f13=0.333333333333333d0
! Parameter of the modified interaction
ec = 0.d0
! First-type gradient functional
igrad=1
alpha=2.78d0
gamma=3.1091d-2
! test on density
if (dabs(rhoc).lt.tol) return
double precision :: vc_a,vc_b
! Spin polarisation
rhoa=max((rhoc+rhoo)*.5d0,1.0d-15)
rhob=max((rhoc-rhoo)*.5d0,1.0d-15)
call ec_lda_sr(mu,rhoa,rhob,eccerflda,vc_a,vc_b)
ecerflda = eccerflda
vrhoccerflda = 0.5d0 * (vc_a + vc_b)
vrhoocerflda = 0.5d0 * (vc_a - vc_b)
! Density
rho = rhoc
rho = max(rho,1.d-10)
! Square of density gradient
drho2 = sigmacc
zeta = (rhoa-rhob)/(rhoa+rhob)
zeta = max(zeta,1.d-10)
! LDA energy density
double precision :: vc_a_lda,vc_b_lda
call ec_lda(rhoa,rhob,ecclda,vc_a_lda,vc_b_lda)
eclda = ecclda
decldadrho = 0.5d0 * (vc_a_lda+vc_b_lda)
decldadrho = 0.5d0 * (vc_a_lda-vc_b_lda)
if ((ecerflda/eclda).le.0d0) then
beta=0d0
else
beta=6.6725d-2*(ecerflda/eclda)**alpha
endif
phi=((1d0+zeta)**(2d0/3d0)+(1d0-zeta)**(2d0/3d0))/2d0
phi2=phi*phi
phi3=phi2*phi
phi4=phi3*phi
tq=drho2*6.346820607d-2*rho**(-7d0/3d0)/phi2
Ab=dexp(-ecerflda/(rho*gamma*phi3))-1d0
if (dabs(Ab).le.dabs(beta*tol)) then
ecerfpbe=ecerflda
else
Aa=beta/(gamma*Ab)
Ac=1d0+Aa*tq+Aa**2*tq**2
if (Aa.lt.tol) Aa=tol
arglog=1d0+beta*(1d0-1d0/Ac)/(gamma*Aa)
arglog=max(arglog,1.d-10)
ecerfpbe=ecerflda+rho*phi3*gamma*dlog(arglog)
end if
ec = ecerfpbe
end

View File

@ -0,0 +1,29 @@
BEGIN_PROVIDER [double precision, psi_dft_energy_kinetic, (N_states) ]
&BEGIN_PROVIDER [double precision, psi_dft_energy_nuclear_elec, (N_states) ]
&BEGIN_PROVIDER [double precision, psi_dft_energy_h_core, (N_states) ]
implicit none
BEGIN_DOC
! kinetic, electron-nuclear and total h_core energy computed with the density matrix one_body_dm_mo_beta_for_dft+one_body_dm_mo_alpha_for_dft
END_DOC
integer :: i,j,istate
double precision :: accu
psi_dft_energy_kinetic = 0.d0
psi_dft_energy_nuclear_elec = 0.d0
do istate = 1, N_states
do i = 1, mo_tot_num
do j = 1, mo_tot_num
psi_dft_energy_kinetic(istate) += ( one_body_dm_mo_alpha_for_dft(j,i,istate)+one_body_dm_mo_beta_for_dft(j,i,istate)) * mo_kinetic_integral(j,i)
psi_dft_energy_nuclear_elec(istate) += ( one_body_dm_mo_alpha_for_dft(j,i,istate)+one_body_dm_mo_beta_for_dft(j,i,istate)) * mo_nucl_elec_integral(j,i)
enddo
enddo
enddo
do i = 1, N_states
do j = 1, mo_tot_num
accu += one_body_dm_mo_alpha_for_dft(j,j,i) + one_body_dm_mo_beta_for_dft(j,j,i)
enddo
accu = (elec_alpha_num + elec_beta_num ) / accu
psi_energy_h_core(i) = psi_dft_energy_h_core(i) * accu
enddo
END_PROVIDER

View File

@ -0,0 +1,272 @@
BEGIN_PROVIDER[double precision, aos_vc_alpha_LDA_w, (n_points_final_grid,ao_num,N_states)]
&BEGIN_PROVIDER[double precision, aos_vc_beta_LDA_w, (n_points_final_grid,ao_num,N_states)]
&BEGIN_PROVIDER[double precision, aos_vx_alpha_LDA_w, (n_points_final_grid,ao_num,N_states)]
&BEGIN_PROVIDER[double precision, aos_vx_beta_LDA_w, (n_points_final_grid,ao_num,N_states)]
implicit none
BEGIN_DOC
! aos_vxc_alpha_LDA_w(j,i) = ao_i(r_j) * (v^x_alpha(r_j) + v^c_alpha(r_j)) * W(r_j)
END_DOC
integer :: istate,i,j
double precision :: r(3)
double precision :: mu,weight
double precision :: e_c,vc_a,vc_b,e_x,vx_a,vx_b
double precision, allocatable :: rhoa(:),rhob(:)
allocate(rhoa(N_states), rhob(N_states))
do istate = 1, N_states
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
weight=final_weight_functions_at_final_grid_points(i)
rhoa(istate) = one_body_dm_alpha_at_r(i,istate)
rhob(istate) = one_body_dm_beta_at_r(i,istate)
call ec_LDA_sr(mu_erf,rhoa(istate),rhob(istate),e_c,vc_a,vc_b)
call ex_LDA_sr(mu_erf,rhoa(istate),rhob(istate),e_x,vx_a,vx_b)
do j =1, ao_num
aos_vc_alpha_LDA_w(i,j,istate) = vc_a * aos_in_r_array(j,i)*weight
aos_vc_beta_LDA_w(i,j,istate) = vc_b * aos_in_r_array(j,i)*weight
aos_vx_alpha_LDA_w(i,j,istate) = vx_a * aos_in_r_array(j,i)*weight
aos_vx_beta_LDA_w(i,j,istate) = vx_b * aos_in_r_array(j,i)*weight
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, energy_x_LDA, (N_states) ]
&BEGIN_PROVIDER[double precision, energy_c_LDA, (N_states) ]
implicit none
BEGIN_DOC
! exchange/correlation energy with the short range LDA functional
END_DOC
integer :: istate,i,j
double precision :: r(3)
double precision :: mu,weight
double precision :: e_c,vc_a,vc_b,e_x,vx_a,vx_b
double precision, allocatable :: rhoa(:),rhob(:)
allocate(rhoa(N_states), rhob(N_states))
energy_x_LDA = 0.d0
energy_c_LDA = 0.d0
do istate = 1, N_states
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
weight=final_weight_functions_at_final_grid_points(i)
rhoa(istate) = one_body_dm_alpha_at_r(i,istate)
rhob(istate) = one_body_dm_beta_at_r(i,istate)
call ec_LDA_sr(mu_erf,rhoa(istate),rhob(istate),e_c,vc_a,vc_b)
call ex_LDA_sr(mu_erf,rhoa(istate),rhob(istate),e_x,vx_a,vx_b)
energy_x_LDA(istate) += weight * e_x
energy_c_LDA(istate) += weight * e_c
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, potential_x_alpha_ao_LDA,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, potential_x_beta_ao_LDA,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, potential_c_alpha_ao_LDA,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, potential_c_beta_ao_LDA,(ao_num,ao_num,N_states)]
implicit none
BEGIN_DOC
! short range exchange/correlation alpha/beta potentials with LDA functional on the AO basis
END_DOC
integer :: istate
double precision :: wall_1,wall_2
call wall_time(wall_1)
do istate = 1, N_states
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0,aos_in_r_array,ao_num,aos_vc_alpha_LDA_w(1,1,istate),n_points_final_grid,0.d0,potential_c_alpha_ao_LDA(1,1,istate),ao_num)
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0,aos_in_r_array,ao_num,aos_vc_beta_LDA_w(1,1,istate) ,n_points_final_grid,0.d0,potential_c_beta_ao_LDA(1,1,istate),ao_num)
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0,aos_in_r_array,ao_num,aos_vx_alpha_LDA_w(1,1,istate),n_points_final_grid,0.d0,potential_x_alpha_ao_LDA(1,1,istate),ao_num)
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0,aos_in_r_array,ao_num,aos_vx_beta_LDA_w(1,1,istate) ,n_points_final_grid,0.d0,potential_x_beta_ao_LDA(1,1,istate),ao_num)
enddo
call wall_time(wall_2)
print*,'time to provide potential_x/c_alpha/beta_ao_LDA = ',wall_2 - wall_1
END_PROVIDER
BEGIN_PROVIDER[double precision, aos_vc_alpha_PBE_w , (ao_num,n_points_final_grid,N_states)] !(n_points_final_grid,ao_num,N_states)]
&BEGIN_PROVIDER[double precision, aos_vc_beta_PBE_w , (ao_num,n_points_final_grid,N_states)]!(n_points_final_grid,ao_num,N_states)]
&BEGIN_PROVIDER[double precision, aos_vx_alpha_PBE_w , (ao_num,n_points_final_grid,N_states)] !(n_points_final_grid,ao_num,N_states)]
&BEGIN_PROVIDER[double precision, aos_vx_beta_PBE_w , (ao_num,n_points_final_grid,N_states)]!(n_points_final_grid,ao_num,N_states)]
&BEGIN_PROVIDER[double precision, aos_dvc_alpha_PBE_w , (ao_num,n_points_final_grid,3,N_states)]
&BEGIN_PROVIDER[double precision, aos_dvc_beta_PBE_w , (ao_num,n_points_final_grid,3,N_states)]
&BEGIN_PROVIDER[double precision, aos_dvx_alpha_PBE_w , (ao_num,n_points_final_grid,3,N_states)]
&BEGIN_PROVIDER[double precision, aos_dvx_beta_PBE_w , (ao_num,n_points_final_grid,3,N_states)]
&BEGIN_PROVIDER[double precision, grad_aos_dvc_alpha_PBE_w , (ao_num,n_points_final_grid,3,N_states)]
&BEGIN_PROVIDER[double precision, grad_aos_dvc_beta_PBE_w , (ao_num,n_points_final_grid,3,N_states)]
&BEGIN_PROVIDER[double precision, grad_aos_dvx_alpha_PBE_w , (ao_num,n_points_final_grid,3,N_states)]
&BEGIN_PROVIDER[double precision, grad_aos_dvx_beta_PBE_w , (ao_num,n_points_final_grid,3,N_states)]
implicit none
BEGIN_DOC
! aos_vxc_alpha_PBE_w(j,i) = ao_i(r_j) * (v^x_alpha(r_j) + v^c_alpha(r_j)) * W(r_j)
END_DOC
integer :: istate,i,j,m
double precision :: r(3)
double precision :: mu,weight
double precision, allocatable :: ex(:), ec(:)
double precision, allocatable :: rho_a(:),rho_b(:),grad_rho_a(:,:),grad_rho_b(:,:),grad_rho_a_2(:),grad_rho_b_2(:),grad_rho_a_b(:)
double precision, allocatable :: contrib_grad_xa(:,:),contrib_grad_xb(:,:),contrib_grad_ca(:,:),contrib_grad_cb(:,:)
double precision, allocatable :: vc_rho_a(:), vc_rho_b(:), vx_rho_a(:), vx_rho_b(:)
double precision, allocatable :: vx_grad_rho_a_2(:), vx_grad_rho_b_2(:), vx_grad_rho_a_b(:), vc_grad_rho_a_2(:), vc_grad_rho_b_2(:), vc_grad_rho_a_b(:)
allocate(vc_rho_a(N_states), vc_rho_b(N_states), vx_rho_a(N_states), vx_rho_b(N_states))
allocate(vx_grad_rho_a_2(N_states), vx_grad_rho_b_2(N_states), vx_grad_rho_a_b(N_states), vc_grad_rho_a_2(N_states), vc_grad_rho_b_2(N_states), vc_grad_rho_a_b(N_states))
allocate(rho_a(N_states), rho_b(N_states),grad_rho_a(3,N_states),grad_rho_b(3,N_states))
allocate(grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states), ex(N_states), ec(N_states))
allocate(contrib_grad_xa(3,N_states),contrib_grad_xb(3,N_states),contrib_grad_ca(3,N_states),contrib_grad_cb(3,N_states))
do istate = 1, N_states
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
weight=final_weight_functions_at_final_grid_points(i)
rho_a(istate) = one_body_dm_alpha_and_grad_at_r(4,i,istate)
rho_b(istate) = one_body_dm_beta_and_grad_at_r(4,i,istate)
grad_rho_a(1:3,istate) = one_body_dm_alpha_and_grad_at_r(1:3,i,istate)
grad_rho_b(1:3,istate) = one_body_dm_beta_and_grad_at_r(1:3,i,istate)
grad_rho_a_2 = 0.d0
grad_rho_b_2 = 0.d0
grad_rho_a_b = 0.d0
do m = 1, 3
grad_rho_a_2(istate) += grad_rho_a(m,istate) * grad_rho_a(m,istate)
grad_rho_b_2(istate) += grad_rho_b(m,istate) * grad_rho_b(m,istate)
grad_rho_a_b(istate) += grad_rho_a(m,istate) * grad_rho_b(m,istate)
enddo
! inputs
call GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange
ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ! outputs correlation
ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b )
vx_rho_a(istate) *= weight
vc_rho_a(istate) *= weight
vx_rho_b(istate) *= weight
vc_rho_b(istate) *= weight
do m= 1,3
contrib_grad_ca(m,istate) = weight * (2.d0 * vc_grad_rho_a_2(istate) * grad_rho_a(m,istate) + vc_grad_rho_a_b(istate) * grad_rho_b(m,istate))
contrib_grad_xa(m,istate) = weight * (2.d0 * vx_grad_rho_a_2(istate) * grad_rho_a(m,istate) + vx_grad_rho_a_b(istate) * grad_rho_b(m,istate))
contrib_grad_cb(m,istate) = weight * (2.d0 * vc_grad_rho_b_2(istate) * grad_rho_b(m,istate) + vc_grad_rho_a_b(istate) * grad_rho_a(m,istate))
contrib_grad_xb(m,istate) = weight * (2.d0 * vx_grad_rho_b_2(istate) * grad_rho_b(m,istate) + vx_grad_rho_a_b(istate) * grad_rho_a(m,istate))
enddo
do j = 1, ao_num
aos_vc_alpha_PBE_w(j,i,istate) = vc_rho_a(istate) * aos_in_r_array(j,i)
aos_vc_beta_PBE_w (j,i,istate) = vc_rho_b(istate) * aos_in_r_array(j,i)
aos_vx_alpha_PBE_w(j,i,istate) = vx_rho_a(istate) * aos_in_r_array(j,i)
aos_vx_beta_PBE_w (j,i,istate) = vx_rho_b(istate) * aos_in_r_array(j,i)
do m = 1,3
aos_dvc_alpha_PBE_w(j,i,m,istate) = contrib_grad_ca(m,istate) * aos_in_r_array(j,i)
aos_dvc_beta_PBE_w (j,i,m,istate) = contrib_grad_cb(m,istate) * aos_in_r_array(j,i)
aos_dvx_alpha_PBE_w(j,i,m,istate) = contrib_grad_xa(m,istate) * aos_in_r_array(j,i)
aos_dvx_beta_PBE_w (j,i,m,istate) = contrib_grad_xb(m,istate) * aos_in_r_array(j,i)
grad_aos_dvc_alpha_PBE_w (j,i,m,istate) = contrib_grad_ca(m,istate) * aos_grad_in_r_array(j,i,m)
grad_aos_dvc_beta_PBE_w (j,i,m,istate) = contrib_grad_cb(m,istate) * aos_grad_in_r_array(j,i,m)
grad_aos_dvx_alpha_PBE_w (j,i,m,istate) = contrib_grad_xa(m,istate) * aos_grad_in_r_array(j,i,m)
grad_aos_dvx_beta_PBE_w (j,i,m,istate) = contrib_grad_xb(m,istate) * aos_grad_in_r_array(j,i,m)
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, energy_x_PBE, (N_states) ]
&BEGIN_PROVIDER[double precision, energy_c_PBE, (N_states) ]
implicit none
BEGIN_DOC
! exchange/correlation energy with the short range PBE functional
END_DOC
integer :: istate,i,j,m
double precision :: r(3)
double precision :: mu,weight
double precision, allocatable :: ex(:), ec(:)
double precision, allocatable :: rho_a(:),rho_b(:),grad_rho_a(:,:),grad_rho_b(:,:),grad_rho_a_2(:),grad_rho_b_2(:),grad_rho_a_b(:)
double precision, allocatable :: contrib_grad_xa(:,:),contrib_grad_xb(:,:),contrib_grad_ca(:,:),contrib_grad_cb(:,:)
double precision, allocatable :: vc_rho_a(:), vc_rho_b(:), vx_rho_a(:), vx_rho_b(:)
double precision, allocatable :: vx_grad_rho_a_2(:), vx_grad_rho_b_2(:), vx_grad_rho_a_b(:), vc_grad_rho_a_2(:), vc_grad_rho_b_2(:), vc_grad_rho_a_b(:)
allocate(vc_rho_a(N_states), vc_rho_b(N_states), vx_rho_a(N_states), vx_rho_b(N_states))
allocate(vx_grad_rho_a_2(N_states), vx_grad_rho_b_2(N_states), vx_grad_rho_a_b(N_states), vc_grad_rho_a_2(N_states), vc_grad_rho_b_2(N_states), vc_grad_rho_a_b(N_states))
allocate(rho_a(N_states), rho_b(N_states),grad_rho_a(3,N_states),grad_rho_b(3,N_states))
allocate(grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states), ex(N_states), ec(N_states))
energy_x_PBE = 0.d0
energy_c_PBE = 0.d0
do istate = 1, N_states
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
weight=final_weight_functions_at_final_grid_points(i)
rho_a(istate) = one_body_dm_alpha_and_grad_at_r(4,i,istate)
rho_b(istate) = one_body_dm_beta_and_grad_at_r(4,i,istate)
grad_rho_a(1:3,istate) = one_body_dm_alpha_and_grad_at_r(1:3,i,istate)
grad_rho_b(1:3,istate) = one_body_dm_beta_and_grad_at_r(1:3,i,istate)
grad_rho_a_2 = 0.d0
grad_rho_b_2 = 0.d0
grad_rho_a_b = 0.d0
do m = 1, 3
grad_rho_a_2(istate) += grad_rho_a(m,istate) * grad_rho_a(m,istate)
grad_rho_b_2(istate) += grad_rho_b(m,istate) * grad_rho_b(m,istate)
grad_rho_a_b(istate) += grad_rho_a(m,istate) * grad_rho_b(m,istate)
enddo
! inputs
call GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange
ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ! outputs correlation
ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b )
energy_x_PBE += ex * weight
energy_c_PBE += ec * weight
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, potential_x_alpha_ao_PBE,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, potential_x_beta_ao_PBE,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, potential_c_alpha_ao_PBE,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, potential_c_beta_ao_PBE,(ao_num,ao_num,N_states)]
implicit none
BEGIN_DOC
! exchange/correlation alpha/beta potentials with the short range PBE functional on the AO basis
END_DOC
integer :: istate, m
double precision :: wall_1,wall_2
call wall_time(wall_1)
potential_c_alpha_ao_PBE = 0.d0
potential_x_alpha_ao_PBE = 0.d0
potential_c_beta_ao_PBE = 0.d0
potential_x_beta_ao_PBE = 0.d0
do istate = 1, N_states
! correlation alpha
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_vc_alpha_PBE_w(1,1,istate),size(aos_vc_alpha_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_c_alpha_ao_PBE(1,1,istate),size(potential_c_alpha_ao_PBE,1))
! correlation beta
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_vc_beta_PBE_w(1,1,istate),size(aos_vc_beta_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_c_beta_ao_PBE(1,1,istate),size(potential_c_beta_ao_PBE,1))
! exchange alpha
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_vx_alpha_PBE_w(1,1,istate),size(aos_vx_alpha_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_x_alpha_ao_PBE(1,1,istate),size(potential_x_alpha_ao_PBE,1))
! exchange beta
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_vx_beta_PBE_w(1,1,istate),size(aos_vx_beta_PBE_w,1), aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_x_beta_ao_PBE(1,1,istate), size(potential_x_beta_ao_PBE,1))
do m= 1,3
! correlation alpha
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_dvc_alpha_PBE_w(1,1,m,istate),size(aos_dvc_alpha_PBE_w,1),aos_grad_in_r_array(1,1,m),size(aos_grad_in_r_array,1),1.d0,potential_c_alpha_ao_PBE(1,1,istate),size(potential_c_alpha_ao_PBE,1))
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,grad_aos_dvc_alpha_PBE_w(1,1,m,istate),size(grad_aos_dvc_alpha_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_c_alpha_ao_PBE(1,1,istate),size(potential_c_alpha_ao_PBE,1))
! correlation beta
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_dvc_beta_PBE_w(1,1,m,istate),size(aos_dvc_beta_PBE_w,1),aos_grad_in_r_array(1,1,m),size(aos_grad_in_r_array,1),1.d0,potential_c_beta_ao_PBE(1,1,istate),size(potential_c_beta_ao_PBE,1))
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,grad_aos_dvc_beta_PBE_w(1,1,m,istate),size(grad_aos_dvc_beta_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_c_beta_ao_PBE(1,1,istate),size(potential_c_beta_ao_PBE,1))
! exchange alpha
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_dvx_alpha_PBE_w(1,1,m,istate),size(aos_dvx_alpha_PBE_w,1),aos_grad_in_r_array(1,1,m),size(aos_grad_in_r_array,1),1.d0,potential_x_alpha_ao_PBE(1,1,istate),size(potential_x_alpha_ao_PBE,1))
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,grad_aos_dvx_alpha_PBE_w(1,1,m,istate),size(grad_aos_dvx_alpha_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_x_alpha_ao_PBE(1,1,istate),size(potential_x_alpha_ao_PBE,1))
! exchange beta
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_dvx_beta_PBE_w(1,1,m,istate),size(aos_dvx_beta_PBE_w,1),aos_grad_in_r_array(1,1,m),size(aos_grad_in_r_array,1),1.d0,potential_x_beta_ao_PBE(1,1,istate),size(potential_x_beta_ao_PBE,1))
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,grad_aos_dvx_beta_PBE_w(1,1,m,istate),size(grad_aos_dvx_beta_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_x_beta_ao_PBE(1,1,istate),size(potential_x_beta_ao_PBE,1))
enddo
enddo
call wall_time(wall_2)
END_PROVIDER

View File

@ -0,0 +1,74 @@
subroutine rho_ab_to_rho_oc(rho_a,rho_b,rho_o,rho_c)
implicit none
double precision, intent(in) :: rho_a,rho_b
double precision, intent(out) :: rho_o,rho_c
rho_c=rho_a+rho_b
rho_o=rho_a-rho_b
end
subroutine rho_oc_to_rho_ab(rho_o,rho_c,rho_a,rho_b)
implicit none
double precision, intent(in) :: rho_o,rho_c
double precision, intent(out) :: rho_a,rho_b
rho_a= 0.5d0*(rho_c+rho_o)
rho_b= 0.5d0*(rho_c-rho_o)
end
subroutine grad_rho_ab_to_grad_rho_oc(grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,grad_rho_o_2,grad_rho_c_2,grad_rho_o_c)
implicit none
double precision, intent(in) :: grad_rho_a_2,grad_rho_b_2,grad_rho_a_b
double precision, intent(out) :: grad_rho_o_2,grad_rho_c_2,grad_rho_o_c
grad_rho_c_2 = grad_rho_a_2 + grad_rho_b_2 + 2d0*grad_rho_a_b
grad_rho_o_2 = grad_rho_a_2 + grad_rho_b_2 - 2d0*grad_rho_a_b
grad_rho_o_c = grad_rho_a_2 - grad_rho_b_2
end
subroutine v_rho_ab_to_v_rho_oc(v_rho_a,v_rho_b,v_rho_o,v_rho_c)
implicit none
double precision, intent(in) :: v_rho_a,v_rho_b
double precision, intent(out) :: v_rho_o,v_rho_c
v_rho_c = 0.5d0*(v_rho_a + v_rho_b)
v_rho_o = 0.5d0*(v_rho_a - v_rho_b)
end
subroutine v_rho_oc_to_v_rho_ab(v_rho_o,v_rho_c,v_rho_a,v_rho_b)
implicit none
double precision, intent(in) :: v_rho_o,v_rho_c
double precision, intent(out) :: v_rho_a,v_rho_b
v_rho_a = v_rho_c + v_rho_o
v_rho_b = v_rho_c - v_rho_o
end
subroutine v_grad_rho_oc_to_v_grad_rho_ab(v_grad_rho_o_2,v_grad_rho_c_2,v_grad_rho_o_c,v_grad_rho_a_2,v_grad_rho_b_2,v_grad_rho_a_b)
implicit none
double precision, intent(in) :: v_grad_rho_o_2,v_grad_rho_c_2,v_grad_rho_o_c
double precision, intent(out) :: v_grad_rho_a_2,v_grad_rho_b_2,v_grad_rho_a_b
v_grad_rho_a_2 = v_grad_rho_o_2 + v_grad_rho_c_2 + v_grad_rho_o_c
v_grad_rho_b_2 = v_grad_rho_o_2 + v_grad_rho_c_2 - v_grad_rho_o_c
v_grad_rho_a_b = -2d0 * v_grad_rho_o_2 + 2d0 * v_grad_rho_c_2
end

View File

@ -0,0 +1,16 @@
BEGIN_PROVIDER [double precision, shifting_constant, (N_states)]
implicit none
BEGIN_DOC
! shifting_constant = (E_{Hxc} - <\Psi | V_{Hxc} | \Psi>) / N_elec
! constant to add to the potential in order to obtain the variational energy as
! the eigenvalue of the effective long-range Hamiltonian
! (see original paper of Levy PRL 113, 113002 (2014), equation (17) )
END_DOC
integer :: istate
do istate = 1, N_states
shifting_constant(istate) = energy_x(istate) + energy_c(istate) + short_range_Hartree(istate) - Trace_v_Hxc(istate)
enddo
shifting_constant = shifting_constant / dble(elec_num)
END_PROVIDER

View File

@ -0,0 +1,118 @@
BEGIN_PROVIDER [double precision, short_range_Hartree_operator, (mo_tot_num,mo_tot_num,N_states)]
&BEGIN_PROVIDER [double precision, short_range_Hartree, (N_states)]
implicit none
BEGIN_DOC
! short_range_Hartree_operator(i,j) = \int dr i(r)j(r) \int r' \rho(r') W_{ee}^{sr}
! short_range_Hartree = 0.5 * \sum_{i,j} \rho_{ij} short_range_Hartree_operator(i,j)
! = 0.5 * \int dr \int r' \rho(r) \rho(r') W_{ee}^{sr}
END_DOC
integer :: i,j,k,l,m,n,istate
double precision :: get_mo_bielec_integral,get_mo_bielec_integral_erf
double precision :: integral, integral_erf, contrib
double precision :: integrals_array(mo_tot_num,mo_tot_num),integrals_erf_array(mo_tot_num,mo_tot_num)
short_range_Hartree_operator = 0.d0
short_range_Hartree = 0.d0
do i = 1, mo_tot_num
do j = 1, mo_tot_num
if(dabs(one_body_dm_average_mo_for_dft(j,i)).le.1.d-12)cycle
call get_mo_bielec_integrals_i1j1(i,j,mo_tot_num,integrals_array,mo_integrals_map)
call get_mo_bielec_integrals_erf_i1j1(i,j,mo_tot_num,integrals_erf_array,mo_integrals_erf_map)
do istate = 1, N_states
do k = 1, mo_tot_num
do l = 1, mo_tot_num
integral = integrals_array(l,k)
integral_erf = integrals_erf_array(l,k)
contrib = one_body_dm_mo_for_dft(i,j,istate) * (integral - integral_erf)
short_range_Hartree_operator(l,k,istate) += contrib
short_range_Hartree(istate) += contrib * one_body_dm_mo_for_dft(k,l,istate)
enddo
enddo
enddo
enddo
enddo
short_range_Hartree = short_range_Hartree * 0.5d0
print*, 'short_range_Hartree',short_range_Hartree
END_PROVIDER
BEGIN_PROVIDER [double precision, effective_one_e_potential, (mo_tot_num, mo_tot_num,N_states)]
&BEGIN_PROVIDER [double precision, effective_one_e_potential_without_kin, (mo_tot_num, mo_tot_num,N_states)]
&BEGIN_PROVIDER [double precision, shifted_effective_one_e_potential_without_kin, (mo_tot_num, mo_tot_num,N_states)]
implicit none
integer :: i,j,istate
effective_one_e_potential = 0.d0
BEGIN_DOC
! effective_one_e_potential(i,j) = <i| v_{H}^{sr} |j> + <i| h_{core} |j> + <i|v_{xc} |j>
! Taking the expectation value does not provide any energy
! but effective_one_e_potential(i,j) is the potential coupling DFT and WFT part to be used in any WFT calculation
! shifted_effective_one_e_potential_without_kin = effective_one_e_potential_without_kin + shifting_constant on the diagonal
END_DOC
do istate = 1, N_states
do i = 1, mo_tot_num
do j = 1, mo_tot_num
effective_one_e_potential(i,j,istate) = short_range_Hartree_operator(i,j,istate) + mo_nucl_elec_integral(i,j) + mo_kinetic_integral(i,j) &
+ 0.5d0 * (potential_x_alpha_mo(i,j,istate) + potential_c_alpha_mo(i,j,istate) &
+ potential_x_beta_mo(i,j,istate) + potential_c_beta_mo(i,j,istate) )
effective_one_e_potential_without_kin(i,j,istate) = short_range_Hartree_operator(i,j,istate) + mo_nucl_elec_integral(i,j) &
+ 0.5d0 * (potential_x_alpha_mo(i,j,istate) + potential_c_alpha_mo(i,j,istate) &
+ potential_x_beta_mo(i,j,istate) + potential_c_beta_mo(i,j,istate) )
shifted_effective_one_e_potential_without_kin(j,i,istate) = effective_one_e_potential_without_kin(j,i,istate)
enddo
enddo
do i = 1, mo_tot_num
shifted_effective_one_e_potential_without_kin(i,i,istate) += shifting_constant(istate)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, Fock_matrix_expectation_value]
implicit none
call get_average(effective_one_e_potential,one_body_dm_average_mo_for_dft,Fock_matrix_expectation_value)
END_PROVIDER
BEGIN_PROVIDER [double precision, Trace_v_xc, (N_states)]
&BEGIN_PROVIDER [double precision, Trace_v_H, (N_states)]
&BEGIN_PROVIDER [double precision, Trace_v_Hxc, (N_states)]
implicit none
integer :: i,j,istate
double precision :: dm
BEGIN_DOC
! Trace_v_xc = \sum_{i,j} (rho_{ij}_\alpha v^{xc}_{ij}^\alpha + rho_{ij}_\beta v^{xc}_{ij}^\beta)
! Trace_v_Hxc = \sum_{i,j} v^{H}_{ij} (rho_{ij}_\alpha + rho_{ij}_\beta)
! Trace_v_Hxc = \sum_{i,j} rho_{ij} v^{Hxc}_{ij}
END_DOC
do istate = 1, N_states
Trace_v_xc(istate) = 0.d0
Trace_v_H(istate) = 0.d0
do i = 1, mo_tot_num
do j = 1, mo_tot_num
Trace_v_xc(istate) += (potential_x_alpha_mo(j,i,istate) + potential_c_alpha_mo(j,i,istate)) * one_body_dm_mo_alpha_for_dft(j,i,istate)
Trace_v_xc(istate) += (potential_x_beta_mo(j,i,istate) + potential_c_beta_mo(j,i,istate) ) * one_body_dm_mo_beta_for_dft(j,i,istate)
dm = one_body_dm_mo_alpha_for_dft(j,i,istate) + one_body_dm_mo_beta_for_dft(j,i,istate)
Trace_v_H(istate) += dm * short_range_Hartree_operator(j,i,istate)
enddo
enddo
Trace_v_Hxc(istate) = Trace_v_xc(istate) + Trace_v_H(istate)
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, DFT_one_e_energy_potential, (mo_tot_num, mo_tot_num,N_states)]
implicit none
integer :: i,j,istate
BEGIN_DOC
! one_e_energy_potential(i,j) = <i|h_{core}|j> + \int dr i(r)j(r) \int r' \rho(r') W_{ee}^{sr}
! If one take the expectation value over Psi, one gets the total one body energy
END_DOC
do istate = 1, N_states
do i = 1, mo_tot_num
do j = 1, mo_tot_num
DFT_one_e_energy_potential(j,i,istate) = mo_nucl_elec_integral(j,i) + mo_kinetic_integral(j,i) + short_range_Hartree_operator(j,i,istate) * 0.5d0
enddo
enddo
enddo
END_PROVIDER

View File

@ -0,0 +1,51 @@
subroutine GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, &
ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, &
ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b )
implicit none
double precision, intent(in) :: r(3),rho_a(N_states),rho_b(N_states),grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states)
double precision, intent(out) :: ex(N_states),vx_rho_a(N_states),vx_rho_b(N_states),vx_grad_rho_a_2(N_states),vx_grad_rho_b_2(N_states),vx_grad_rho_a_b(N_states)
double precision, intent(out) :: ec(N_states),vc_rho_a(N_states),vc_rho_b(N_states),vc_grad_rho_a_2(N_states),vc_grad_rho_b_2(N_states),vc_grad_rho_a_b(N_states)
integer :: istate
double precision :: r2(3),dr2(3), local_potential,r12,dx2,mu
do istate = 1, N_states
if(exchange_functional.EQ."short_range_PBE")then
call ex_pbe_sr(mu_erf,rho_a(istate),rho_b(istate),grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),ex(istate),vx_rho_a(istate),vx_rho_b(istate),vx_grad_rho_a_2(istate),vx_grad_rho_b_2(istate),vx_grad_rho_a_b(istate))
else if(exchange_functional.EQ."None")then
ex = 0.d0
vx_rho_a = 0.d0
vx_rho_b = 0.d0
vx_grad_rho_a_2 = 0.d0
vx_grad_rho_a_b = 0.d0
vx_grad_rho_b_2 = 0.d0
else
print*, 'Exchange functional required does not exist ...'
print*,'exchange_functional',exchange_functional
stop
endif
double precision :: rhoc,rhoo,sigmacc,sigmaco,sigmaoo,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo
if(correlation_functional.EQ."short_range_PBE")then
! convertion from (alpha,beta) formalism to (closed, open) formalism
call rho_ab_to_rho_oc(rho_a(istate),rho_b(istate),rhoo,rhoc)
call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),sigmaoo,sigmacc,sigmaco)
call ec_pbe_sr(mu_erf,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,ec(istate),vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo)
call v_rho_oc_to_v_rho_ab(vrhoo,vrhoc,vc_rho_a(istate),vc_rho_b(istate))
call v_grad_rho_oc_to_v_grad_rho_ab(vsigmaoo,vsigmacc,vsigmaco,vc_grad_rho_a_2(istate),vc_grad_rho_b_2(istate),vc_grad_rho_a_b(istate))
else if(correlation_functional.EQ."None")then
ec = 0.d0
vc_rho_a = 0.d0
vc_rho_b = 0.d0
vc_grad_rho_a_2 = 0.d0
vc_grad_rho_a_b = 0.d0
vc_grad_rho_b_2 = 0.d0
else
print*, 'Correlation functional required does not exist ...'
print*, 'correlation_functional',correlation_functional
stop
endif
enddo
end

View File

@ -8,6 +8,7 @@ davidson
davidson_dressed
davidson_undressed
determinants
dft_utils_one_body
dressing
electrons
ezfio_files

View File

@ -0,0 +1,25 @@
[disk_access_ao_integrals_erf]
type: Disk_access
doc: Read/Write AO integrals with the long range interaction from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[disk_access_mo_integrals_erf]
type: Disk_access
doc: Read/Write MO integrals with the long range interaction from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[disk_access_mo_integrals_sr]
type: Disk_access
doc: Read/Write MO integrals with the short range interaction from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[disk_access_ao_integrals_sr]
type: Disk_access
doc: Read/Write AO integrals with the short range interaction from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None

View File

@ -0,0 +1,5 @@
pseudo
bitmask
zmq
integrals_bielec
dft_keywords

View File

@ -0,0 +1,649 @@
double precision function ao_bielec_integral_erf(i,j,k,l)
implicit none
BEGIN_DOC
! integral of the AO basis <ik|jl> or (ij|kl)
! i(r1) j(r1) 1/r12 k(r2) l(r2)
END_DOC
integer,intent(in) :: i,j,k,l
integer :: p,q,r,s
double precision :: I_center(3),J_center(3),K_center(3),L_center(3)
integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
double precision :: integral
include 'utils/constants.include.F'
double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp
double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
integer :: iorder_p(3), iorder_q(3)
double precision :: ao_bielec_integral_schwartz_accel_erf
if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
ao_bielec_integral_erf = ao_bielec_integral_schwartz_accel_erf(i,j,k,l)
return
endif
dim1 = n_pt_max_integrals
num_i = ao_nucl(i)
num_j = ao_nucl(j)
num_k = ao_nucl(k)
num_l = ao_nucl(l)
ao_bielec_integral_erf = 0.d0
if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then
do p = 1, 3
I_power(p) = ao_power(i,p)
J_power(p) = ao_power(j,p)
K_power(p) = ao_power(k,p)
L_power(p) = ao_power(l,p)
I_center(p) = nucl_coord(num_i,p)
J_center(p) = nucl_coord(num_j,p)
K_center(p) = nucl_coord(num_k,p)
L_center(p) = nucl_coord(num_l,p)
enddo
double precision :: coef1, coef2, coef3, coef4
double precision :: p_inv,q_inv
double precision :: general_primitive_integral_erf
do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p,i)
do q = 1, ao_prim_num(j)
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), &
I_power,J_power,I_center,J_center,dim1)
p_inv = 1.d0/pp
do r = 1, ao_prim_num(k)
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
do s = 1, ao_prim_num(l)
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), &
K_power,L_power,K_center,L_center,dim1)
q_inv = 1.d0/qq
integral = general_primitive_integral_erf(dim1, &
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
ao_bielec_integral_erf = ao_bielec_integral_erf + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
else
do p = 1, 3
I_power(p) = ao_power(i,p)
J_power(p) = ao_power(j,p)
K_power(p) = ao_power(k,p)
L_power(p) = ao_power(l,p)
enddo
double precision :: ERI_erf
do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p,i)
do q = 1, ao_prim_num(j)
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
do r = 1, ao_prim_num(k)
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
do s = 1, ao_prim_num(l)
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
integral = ERI_erf( &
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),&
I_power(1),J_power(1),K_power(1),L_power(1), &
I_power(2),J_power(2),K_power(2),L_power(2), &
I_power(3),J_power(3),K_power(3),L_power(3))
ao_bielec_integral_erf = ao_bielec_integral_erf + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
endif
end
double precision function ao_bielec_integral_schwartz_accel_erf(i,j,k,l)
implicit none
BEGIN_DOC
! integral of the AO basis <ik|jl> or (ij|kl)
! i(r1) j(r1) 1/r12 k(r2) l(r2)
END_DOC
integer,intent(in) :: i,j,k,l
integer :: p,q,r,s
double precision :: I_center(3),J_center(3),K_center(3),L_center(3)
integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
double precision :: integral
include 'utils/constants.include.F'
double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp
double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
integer :: iorder_p(3), iorder_q(3)
double precision, allocatable :: schwartz_kl(:,:)
double precision :: schwartz_ij
dim1 = n_pt_max_integrals
num_i = ao_nucl(i)
num_j = ao_nucl(j)
num_k = ao_nucl(k)
num_l = ao_nucl(l)
ao_bielec_integral_schwartz_accel_erf = 0.d0
double precision :: thr
thr = ao_integrals_threshold*ao_integrals_threshold
allocate(schwartz_kl(0:ao_prim_num(l),0:ao_prim_num(k)))
double precision :: coef3
double precision :: coef2
double precision :: p_inv,q_inv
double precision :: coef1
double precision :: coef4
if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then
do p = 1, 3
I_power(p) = ao_power(i,p)
J_power(p) = ao_power(j,p)
K_power(p) = ao_power(k,p)
L_power(p) = ao_power(l,p)
I_center(p) = nucl_coord(num_i,p)
J_center(p) = nucl_coord(num_j,p)
K_center(p) = nucl_coord(num_k,p)
L_center(p) = nucl_coord(num_l,p)
enddo
schwartz_kl(0,0) = 0.d0
do r = 1, ao_prim_num(k)
coef1 = ao_coef_normalized_ordered_transp(r,k)*ao_coef_normalized_ordered_transp(r,k)
schwartz_kl(0,r) = 0.d0
do s = 1, ao_prim_num(l)
coef2 = coef1 * ao_coef_normalized_ordered_transp(s,l) * ao_coef_normalized_ordered_transp(s,l)
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), &
K_power,L_power,K_center,L_center,dim1)
q_inv = 1.d0/qq
schwartz_kl(s,r) = general_primitive_integral_erf(dim1, &
Q_new,Q_center,fact_q,qq,q_inv,iorder_q, &
Q_new,Q_center,fact_q,qq,q_inv,iorder_q) &
* coef2
schwartz_kl(0,r) = max(schwartz_kl(0,r),schwartz_kl(s,r))
enddo
schwartz_kl(0,0) = max(schwartz_kl(0,r),schwartz_kl(0,0))
enddo
do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p,i)
do q = 1, ao_prim_num(j)
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), &
I_power,J_power,I_center,J_center,dim1)
p_inv = 1.d0/pp
schwartz_ij = general_primitive_integral_erf(dim1, &
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
P_new,P_center,fact_p,pp,p_inv,iorder_p) * &
coef2*coef2
if (schwartz_kl(0,0)*schwartz_ij < thr) then
cycle
endif
do r = 1, ao_prim_num(k)
if (schwartz_kl(0,r)*schwartz_ij < thr) then
cycle
endif
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
do s = 1, ao_prim_num(l)
if (schwartz_kl(s,r)*schwartz_ij < thr) then
cycle
endif
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
double precision :: general_primitive_integral_erf
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), &
K_power,L_power,K_center,L_center,dim1)
q_inv = 1.d0/qq
integral = general_primitive_integral_erf(dim1, &
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
ao_bielec_integral_schwartz_accel_erf = ao_bielec_integral_schwartz_accel_erf + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
else
do p = 1, 3
I_power(p) = ao_power(i,p)
J_power(p) = ao_power(j,p)
K_power(p) = ao_power(k,p)
L_power(p) = ao_power(l,p)
enddo
double precision :: ERI_erf
schwartz_kl(0,0) = 0.d0
do r = 1, ao_prim_num(k)
coef1 = ao_coef_normalized_ordered_transp(r,k)*ao_coef_normalized_ordered_transp(r,k)
schwartz_kl(0,r) = 0.d0
do s = 1, ao_prim_num(l)
coef2 = coef1*ao_coef_normalized_ordered_transp(s,l)*ao_coef_normalized_ordered_transp(s,l)
schwartz_kl(s,r) = ERI_erf( &
ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),&
K_power(1),L_power(1),K_power(1),L_power(1), &
K_power(2),L_power(2),K_power(2),L_power(2), &
K_power(3),L_power(3),K_power(3),L_power(3)) * &
coef2
schwartz_kl(0,r) = max(schwartz_kl(0,r),schwartz_kl(s,r))
enddo
schwartz_kl(0,0) = max(schwartz_kl(0,r),schwartz_kl(0,0))
enddo
do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p,i)
do q = 1, ao_prim_num(j)
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
schwartz_ij = ERI_erf( &
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),&
I_power(1),J_power(1),I_power(1),J_power(1), &
I_power(2),J_power(2),I_power(2),J_power(2), &
I_power(3),J_power(3),I_power(3),J_power(3))*coef2*coef2
if (schwartz_kl(0,0)*schwartz_ij < thr) then
cycle
endif
do r = 1, ao_prim_num(k)
if (schwartz_kl(0,r)*schwartz_ij < thr) then
cycle
endif
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
do s = 1, ao_prim_num(l)
if (schwartz_kl(s,r)*schwartz_ij < thr) then
cycle
endif
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
integral = ERI_erf( &
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),&
I_power(1),J_power(1),K_power(1),L_power(1), &
I_power(2),J_power(2),K_power(2),L_power(2), &
I_power(3),J_power(3),K_power(3),L_power(3))
ao_bielec_integral_schwartz_accel_erf = ao_bielec_integral_schwartz_accel_erf + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
endif
deallocate (schwartz_kl)
end
subroutine compute_ao_bielec_integrals_erf(j,k,l,sze,buffer_value)
implicit none
use map_module
BEGIN_DOC
! Compute AO 1/r12 integrals for all i and fixed j,k,l
END_DOC
include 'utils/constants.include.F'
integer, intent(in) :: j,k,l,sze
real(integral_kind), intent(out) :: buffer_value(sze)
double precision :: ao_bielec_integral_erf
integer :: i
if (ao_overlap_abs(j,l) < thresh) then
buffer_value = 0._integral_kind
return
endif
if (ao_bielec_integral_erf_schwartz(j,l) < thresh ) then
buffer_value = 0._integral_kind
return
endif
do i = 1, ao_num
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then
buffer_value(i) = 0._integral_kind
cycle
endif
if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < thresh ) then
buffer_value(i) = 0._integral_kind
cycle
endif
!DIR$ FORCEINLINE
buffer_value(i) = ao_bielec_integral_erf(i,k,j,l)
enddo
end
double precision function general_primitive_integral_erf(dim, &
P_new,P_center,fact_p,p,p_inv,iorder_p, &
Q_new,Q_center,fact_q,q,q_inv,iorder_q)
implicit none
BEGIN_DOC
! Computes the integral <pq|rs> where p,q,r,s are Gaussian primitives
END_DOC
integer,intent(in) :: dim
include 'utils/constants.include.F'
double precision, intent(in) :: P_new(0:max_dim,3),P_center(3),fact_p,p,p_inv
double precision, intent(in) :: Q_new(0:max_dim,3),Q_center(3),fact_q,q,q_inv
integer, intent(in) :: iorder_p(3)
integer, intent(in) :: iorder_q(3)
double precision :: r_cut,gama_r_cut,rho,dist
double precision :: dx(0:max_dim),Ix_pol(0:max_dim),dy(0:max_dim),Iy_pol(0:max_dim),dz(0:max_dim),Iz_pol(0:max_dim)
integer :: n_Ix,n_Iy,n_Iz,nx,ny,nz
double precision :: bla
integer :: ix,iy,iz,jx,jy,jz,i
double precision :: a,b,c,d,e,f,accu,pq,const
double precision :: pq_inv, p10_1, p10_2, p01_1, p01_2,pq_inv_2
integer :: n_pt_tmp,n_pt_out, iorder
double precision :: d1(0:max_dim),d_poly(0:max_dim),rint,d1_screened(0:max_dim)
general_primitive_integral_erf = 0.d0
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: dx,Ix_pol,dy,Iy_pol,dz,Iz_pol
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: d1, d_poly
! Gaussian Product
! ----------------
double precision :: p_plus_q
p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf)
pq = p_inv*0.5d0*q_inv
pq_inv = 0.5d0/p_plus_q
p10_1 = q*pq ! 1/(2p)
p01_1 = p*pq ! 1/(2q)
pq_inv_2 = pq_inv+pq_inv
p10_2 = pq_inv_2 * p10_1*q !0.5d0*q/(pq + p*p)
p01_2 = pq_inv_2 * p01_1*p !0.5d0*p/(q*q + pq)
accu = 0.d0
iorder = iorder_p(1)+iorder_q(1)+iorder_p(1)+iorder_q(1)
!DIR$ VECTOR ALIGNED
do ix=0,iorder
Ix_pol(ix) = 0.d0
enddo
n_Ix = 0
do ix = 0, iorder_p(1)
if (abs(P_new(ix,1)) < thresh) cycle
a = P_new(ix,1)
do jx = 0, iorder_q(1)
d = a*Q_new(jx,1)
if (abs(d) < thresh) cycle
!DEC$ FORCEINLINE
call give_polynom_mult_center_x(P_center(1),Q_center(1),ix,jx,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dx,nx)
!DEC$ FORCEINLINE
call add_poly_multiply(dx,nx,d,Ix_pol,n_Ix)
enddo
enddo
if (n_Ix == -1) then
return
endif
iorder = iorder_p(2)+iorder_q(2)+iorder_p(2)+iorder_q(2)
!DIR$ VECTOR ALIGNED
do ix=0, iorder
Iy_pol(ix) = 0.d0
enddo
n_Iy = 0
do iy = 0, iorder_p(2)
if (abs(P_new(iy,2)) > thresh) then
b = P_new(iy,2)
do jy = 0, iorder_q(2)
e = b*Q_new(jy,2)
if (abs(e) < thresh) cycle
!DEC$ FORCEINLINE
call give_polynom_mult_center_x(P_center(2),Q_center(2),iy,jy,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dy,ny)
!DEC$ FORCEINLINE
call add_poly_multiply(dy,ny,e,Iy_pol,n_Iy)
enddo
endif
enddo
if (n_Iy == -1) then
return
endif
iorder = iorder_p(3)+iorder_q(3)+iorder_p(3)+iorder_q(3)
do ix=0,iorder
Iz_pol(ix) = 0.d0
enddo
n_Iz = 0
do iz = 0, iorder_p(3)
if (abs(P_new(iz,3)) > thresh) then
c = P_new(iz,3)
do jz = 0, iorder_q(3)
f = c*Q_new(jz,3)
if (abs(f) < thresh) cycle
!DEC$ FORCEINLINE
call give_polynom_mult_center_x(P_center(3),Q_center(3),iz,jz,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dz,nz)
!DEC$ FORCEINLINE
call add_poly_multiply(dz,nz,f,Iz_pol,n_Iz)
enddo
endif
enddo
if (n_Iz == -1) then
return
endif
rho = p*q *pq_inv_2 ! le rho qui va bien
dist = (P_center(1) - Q_center(1))*(P_center(1) - Q_center(1)) + &
(P_center(2) - Q_center(2))*(P_center(2) - Q_center(2)) + &
(P_center(3) - Q_center(3))*(P_center(3) - Q_center(3))
const = dist*rho
n_pt_tmp = n_Ix+n_Iy
do i=0,n_pt_tmp
d_poly(i)=0.d0
enddo
!DEC$ FORCEINLINE
call multiply_poly(Ix_pol,n_Ix,Iy_pol,n_Iy,d_poly,n_pt_tmp)
if (n_pt_tmp == -1) then
return
endif
n_pt_out = n_pt_tmp+n_Iz
do i=0,n_pt_out
d1(i)=0.d0
enddo
!DEC$ FORCEINLINE
call multiply_poly(d_poly ,n_pt_tmp ,Iz_pol,n_Iz,d1,n_pt_out)
double precision :: rint_sum
accu = accu + rint_sum(n_pt_out,const,d1)
! change p+q in dsqrt
general_primitive_integral_erf = fact_p * fact_q * accu *pi_5_2*p_inv*q_inv/dsqrt(p_plus_q)
end
double precision function ERI_erf(alpha,beta,delta,gama,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z)
implicit none
BEGIN_DOC
! ATOMIC PRIMTIVE bielectronic integral between the 4 primitives ::
! primitive_1 = x1**(a_x) y1**(a_y) z1**(a_z) exp(-alpha * r1**2)
! primitive_2 = x1**(b_x) y1**(b_y) z1**(b_z) exp(- beta * r1**2)
! primitive_3 = x2**(c_x) y2**(c_y) z2**(c_z) exp(-delta * r2**2)
! primitive_4 = x2**(d_x) y2**(d_y) z2**(d_z) exp(- gama * r2**2)
END_DOC
double precision, intent(in) :: delta,gama,alpha,beta
integer, intent(in) :: a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z
integer :: a_x_2,b_x_2,c_x_2,d_x_2,a_y_2,b_y_2,c_y_2,d_y_2,a_z_2,b_z_2,c_z_2,d_z_2
integer :: i,j,k,l,n_pt
integer :: n_pt_sup
double precision :: p,q,denom,coeff
double precision :: I_f
integer :: nx,ny,nz
include 'utils/constants.include.F'
nx = a_x+b_x+c_x+d_x
if(iand(nx,1) == 1) then
ERI_erf = 0.d0
return
endif
ny = a_y+b_y+c_y+d_y
if(iand(ny,1) == 1) then
ERI_erf = 0.d0
return
endif
nz = a_z+b_z+c_z+d_z
if(iand(nz,1) == 1) then
ERI_erf = 0.d0
return
endif
ASSERT (alpha >= 0.d0)
ASSERT (beta >= 0.d0)
ASSERT (delta >= 0.d0)
ASSERT (gama >= 0.d0)
p = alpha + beta
q = delta + gama
double precision :: p_plus_q
p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf)
ASSERT (p+q >= 0.d0)
n_pt = ishft( nx+ny+nz,1 )
coeff = pi_5_2 / (p * q * dsqrt(p_plus_q))
if (n_pt == 0) then
ERI_erf = coeff
return
endif
call integrale_new_erf(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q,n_pt)
ERI_erf = I_f * coeff
end
subroutine integrale_new_erf(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q,n_pt)
BEGIN_DOC
! calculate the integral of the polynom ::
! I_x1(a_x+b_x, c_x+d_x,p,q) * I_x1(a_y+b_y, c_y+d_y,p,q) * I_x1(a_z+b_z, c_z+d_z,p,q)
! between ( 0 ; 1)
END_DOC
implicit none
include 'utils/constants.include.F'
double precision :: p,q
integer :: a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z
integer :: i, n_iter, n_pt, j
double precision :: I_f, pq_inv, p10_1, p10_2, p01_1, p01_2,rho,pq_inv_2
integer :: ix,iy,iz, jx,jy,jz, sx,sy,sz
j = ishft(n_pt,-1)
ASSERT (n_pt > 1)
double precision :: p_plus_q
p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf)
pq_inv = 0.5d0/(p_plus_q)
pq_inv_2 = pq_inv + pq_inv
p10_1 = 0.5d0/p
p01_1 = 0.5d0/q
p10_2 = 0.5d0 * q /(p * p_plus_q)
p01_2 = 0.5d0 * p /(q * p_plus_q)
double precision :: B00(n_pt_max_integrals)
double precision :: B10(n_pt_max_integrals), B01(n_pt_max_integrals)
double precision :: t1(n_pt_max_integrals), t2(n_pt_max_integrals)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: t1, t2, B10, B01, B00
ix = a_x+b_x
jx = c_x+d_x
iy = a_y+b_y
jy = c_y+d_y
iz = a_z+b_z
jz = c_z+d_z
sx = ix+jx
sy = iy+jy
sz = iz+jz
!DIR$ VECTOR ALIGNED
do i = 1,n_pt
B10(i) = p10_1 - gauleg_t2(i,j)* p10_2
B01(i) = p01_1 - gauleg_t2(i,j)* p01_2
B00(i) = gauleg_t2(i,j)*pq_inv
enddo
if (sx > 0) then
call I_x1_new(ix,jx,B10,B01,B00,t1,n_pt)
else
!DIR$ VECTOR ALIGNED
do i = 1,n_pt
t1(i) = 1.d0
enddo
endif
if (sy > 0) then
call I_x1_new(iy,jy,B10,B01,B00,t2,n_pt)
!DIR$ VECTOR ALIGNED
do i = 1,n_pt
t1(i) = t1(i)*t2(i)
enddo
endif
if (sz > 0) then
call I_x1_new(iz,jz,B10,B01,B00,t2,n_pt)
!DIR$ VECTOR ALIGNED
do i = 1,n_pt
t1(i) = t1(i)*t2(i)
enddo
endif
I_f= 0.d0
!DIR$ VECTOR ALIGNED
do i = 1,n_pt
I_f += gauleg_w(i,j)*t1(i)
enddo
end
subroutine compute_ao_integrals_erf_jl(j,l,n_integrals,buffer_i,buffer_value)
implicit none
use map_module
BEGIN_DOC
! Parallel client for AO integrals
END_DOC
integer, intent(in) :: j,l
integer,intent(out) :: n_integrals
integer(key_kind),intent(out) :: buffer_i(ao_num*ao_num)
real(integral_kind),intent(out) :: buffer_value(ao_num*ao_num)
integer :: i,k
double precision :: ao_bielec_integral_erf,cpu_1,cpu_2, wall_1, wall_2
double precision :: integral, wall_0
double precision :: thr
integer :: kk, m, j1, i1
thr = ao_integrals_threshold
n_integrals = 0
j1 = j+ishft(l*l-l,-1)
do k = 1, ao_num ! r1
i1 = ishft(k*k-k,-1)
if (i1 > j1) then
exit
endif
do i = 1, k
i1 += 1
if (i1 > j1) then
exit
endif
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thr) then
cycle
endif
if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < thr ) then
cycle
endif
!DIR$ FORCEINLINE
integral = ao_bielec_integral_erf(i,k,j,l) ! i,k : r1 j,l : r2
if (abs(integral) < thr) then
cycle
endif
n_integrals += 1
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,buffer_i(n_integrals))
buffer_value(n_integrals) = integral
enddo
enddo
end

View File

@ -0,0 +1,194 @@
subroutine ao_bielec_integrals_erf_in_map_slave_tcp(i)
implicit none
integer, intent(in) :: i
BEGIN_DOC
! Computes a buffer of integrals. i is the ID of the current thread.
END_DOC
call ao_bielec_integrals_erf_in_map_slave(0,i)
end
subroutine ao_bielec_integrals_erf_in_map_slave_inproc(i)
implicit none
integer, intent(in) :: i
BEGIN_DOC
! Computes a buffer of integrals. i is the ID of the current thread.
END_DOC
call ao_bielec_integrals_erf_in_map_slave(1,i)
end
subroutine ao_bielec_integrals_erf_in_map_slave(thread,iproc)
use map_module
use f77_zmq
implicit none
BEGIN_DOC
! Computes a buffer of integrals
END_DOC
integer, intent(in) :: thread, iproc
integer :: j,l,n_integrals
integer :: rc
real(integral_kind), allocatable :: buffer_value(:)
integer(key_kind), allocatable :: buffer_i(:)
integer :: worker_id, task_id
character*(512) :: task
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
integer(ZMQ_PTR), external :: new_zmq_push_socket
integer(ZMQ_PTR) :: zmq_socket_push
character*(64) :: state
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
integer, external :: connect_to_taskserver
if (connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) == -1) then
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
return
endif
zmq_socket_push = new_zmq_push_socket(thread)
allocate ( buffer_i(ao_num*ao_num), buffer_value(ao_num*ao_num) )
do
integer, external :: get_task_from_taskserver
if (get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task) == -1) then
exit
endif
if (task_id == 0) exit
read(task,*) j, l
integer, external :: task_done_to_taskserver
call compute_ao_integrals_erf_jl(j,l,n_integrals,buffer_i,buffer_value)
if (task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id) == -1) then
stop 'Unable to send task_done'
endif
call push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, task_id)
enddo
integer, external :: disconnect_from_taskserver
if (disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id) == -1) then
continue
endif
deallocate( buffer_i, buffer_value )
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call end_zmq_push_socket(zmq_socket_push,thread)
end
subroutine ao_bielec_integrals_erf_in_map_collector(zmq_socket_pull)
use map_module
use f77_zmq
implicit none
BEGIN_DOC
! Collects results from the AO integral calculation
END_DOC
integer(ZMQ_PTR), intent(in) :: zmq_socket_pull
integer :: j,l,n_integrals
integer :: rc
real(integral_kind), allocatable :: buffer_value(:)
integer(key_kind), allocatable :: buffer_i(:)
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
integer(ZMQ_PTR), external :: new_zmq_pull_socket
integer*8 :: control, accu, sze
integer :: task_id, more
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
sze = ao_num*ao_num
allocate ( buffer_i(sze), buffer_value(sze) )
accu = 0_8
more = 1
do while (more == 1)
rc = f77_zmq_recv( zmq_socket_pull, n_integrals, 4, 0)
if (rc == -1) then
n_integrals = 0
return
endif
if (rc /= 4) then
print *, irp_here, ': f77_zmq_recv( zmq_socket_pull, n_integrals, 4, 0)'
stop 'error'
endif
if (n_integrals >= 0) then
if (n_integrals > sze) then
deallocate (buffer_value, buffer_i)
sze = n_integrals
allocate (buffer_value(sze), buffer_i(sze))
endif
rc = f77_zmq_recv( zmq_socket_pull, buffer_i, key_kind*n_integrals, 0)
if (rc /= key_kind*n_integrals) then
print *, rc, key_kind, n_integrals
print *, irp_here, ': f77_zmq_recv( zmq_socket_pull, buffer_i, key_kind*n_integrals, 0)'
stop 'error'
endif
rc = f77_zmq_recv( zmq_socket_pull, buffer_value, integral_kind*n_integrals, 0)
if (rc /= integral_kind*n_integrals) then
print *, irp_here, ': f77_zmq_recv( zmq_socket_pull, buffer_value, integral_kind*n_integrals, 0)'
stop 'error'
endif
rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0)
IRP_IF ZMQ_PUSH
IRP_ELSE
rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)
if (rc /= 4) then
print *, irp_here, ' : f77_zmq_send (zmq_socket_pull,...'
stop 'error'
endif
IRP_ENDIF
call insert_into_ao_integrals_erf_map(n_integrals,buffer_i,buffer_value)
accu += n_integrals
if (task_id /= 0) then
integer, external :: zmq_delete_task
if (zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more) == -1) then
stop 'Unable to delete task'
endif
endif
endif
enddo
deallocate( buffer_i, buffer_value )
integer (map_size_kind) :: get_ao_erf_map_size
control = get_ao_erf_map_size(ao_integrals_erf_map)
if (control /= accu) then
print *, ''
print *, irp_here
print *, 'Control : ', control
print *, 'Accu : ', accu
print *, 'Some integrals were lost during the parallel computation.'
print *, 'Try to reduce the number of threads.'
stop
endif
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
end

View File

@ -0,0 +1,47 @@
BEGIN_PROVIDER [double precision, big_array_coulomb_integrals_erf, (mo_tot_num,mo_tot_num, mo_tot_num)]
&BEGIN_PROVIDER [double precision, big_array_exchange_integrals_erf,(mo_tot_num,mo_tot_num, mo_tot_num)]
implicit none
integer :: i,j,k,l
double precision :: get_mo_bielec_integral_erf
double precision :: integral
do k = 1, mo_tot_num
do i = 1, mo_tot_num
do j = 1, mo_tot_num
l = j
integral = get_mo_bielec_integral_erf(i,j,k,l,mo_integrals_erf_map)
big_array_coulomb_integrals_erf(j,i,k) = integral
l = j
integral = get_mo_bielec_integral_erf(i,j,l,k,mo_integrals_erf_map)
big_array_exchange_integrals_erf(j,i,k) = integral
enddo
enddo
enddo
END_PROVIDER
!BEGIN_PROVIDER [double precision, big_array_coulomb_integrals_sr, (mo_tot_num,mo_tot_num, mo_tot_num)]
!&BEGIN_PROVIDER [double precision, big_array_exchange_integrals_sr,(mo_tot_num,mo_tot_num, mo_tot_num)]
!implicit none
!integer :: i,j,k,l
!double precision :: get_mo_bielec_integral_sr
!double precision :: integral
!do k = 1, mo_tot_num
! do i = 1, mo_tot_num
! do j = 1, mo_tot_num
! l = j
! integral = get_mo_bielec_integral_sr(i,j,k,l,mo_integrals_sr_map)
! big_array_coulomb_integrals_sr(j,i,k) = integral
! l = j
! integral = get_mo_bielec_integral_sr(i,j,l,k,mo_integrals_sr_map)
! big_array_exchange_integrals_sr(j,i,k) = integral
! enddo
! enddo
!enddo
!ND_PROVIDER

View File

@ -0,0 +1,682 @@
use map_module
!! AO Map
!! ======
BEGIN_PROVIDER [ type(map_type), ao_integrals_erf_map ]
implicit none
BEGIN_DOC
! AO integrals
END_DOC
integer(key_kind) :: key_max
integer(map_size_kind) :: sze
call bielec_integrals_index(ao_num,ao_num,ao_num,ao_num,key_max)
sze = key_max
call map_init(ao_integrals_erf_map,sze)
print*, 'AO map initialized : ', sze
END_PROVIDER
BEGIN_PROVIDER [ integer, ao_integrals_erf_cache_min ]
&BEGIN_PROVIDER [ integer, ao_integrals_erf_cache_max ]
implicit none
BEGIN_DOC
! Min and max values of the AOs for which the integrals are in the cache
END_DOC
ao_integrals_erf_cache_min = max(1,ao_num - 63)
ao_integrals_erf_cache_max = ao_num
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_integrals_erf_cache, (0:64*64*64*64) ]
use map_module
implicit none
BEGIN_DOC
! Cache of AO integrals for fast access
END_DOC
PROVIDE ao_bielec_integrals_erf_in_map
integer :: i,j,k,l,ii
integer(key_kind) :: idx
real(integral_kind) :: integral
!$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
do l=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
do k=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
do j=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
do i=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(ao_integrals_erf_map,idx,integral)
ii = l-ao_integrals_erf_cache_min
ii = ior( ishft(ii,6), k-ao_integrals_erf_cache_min)
ii = ior( ishft(ii,6), j-ao_integrals_erf_cache_min)
ii = ior( ishft(ii,6), i-ao_integrals_erf_cache_min)
ao_integrals_erf_cache(ii) = integral
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
double precision function get_ao_bielec_integral_erf(i,j,k,l,map) result(result)
use map_module
implicit none
BEGIN_DOC
! Gets one AO bi-electronic integral from the AO map
END_DOC
integer, intent(in) :: i,j,k,l
integer(key_kind) :: idx
type(map_type), intent(inout) :: map
integer :: ii
real(integral_kind) :: tmp
PROVIDE ao_bielec_integrals_erf_in_map ao_integrals_erf_cache ao_integrals_erf_cache_min
!DIR$ FORCEINLINE
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then
tmp = 0.d0
else if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < ao_integrals_threshold) then
tmp = 0.d0
else
ii = l-ao_integrals_erf_cache_min
ii = ior(ii, k-ao_integrals_erf_cache_min)
ii = ior(ii, j-ao_integrals_erf_cache_min)
ii = ior(ii, i-ao_integrals_erf_cache_min)
if (iand(ii, -64) /= 0) then
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(map,idx,tmp)
tmp = tmp
else
ii = l-ao_integrals_erf_cache_min
ii = ior( ishft(ii,6), k-ao_integrals_erf_cache_min)
ii = ior( ishft(ii,6), j-ao_integrals_erf_cache_min)
ii = ior( ishft(ii,6), i-ao_integrals_erf_cache_min)
tmp = ao_integrals_erf_cache(ii)
endif
endif
result = tmp
end
subroutine get_ao_bielec_integrals_erf(j,k,l,sze,out_val)
use map_module
BEGIN_DOC
! Gets multiple AO bi-electronic integral from the AO map .
! All i are retrieved for j,k,l fixed.
END_DOC
implicit none
integer, intent(in) :: j,k,l, sze
real(integral_kind), intent(out) :: out_val(sze)
integer :: i
integer(key_kind) :: hash
double precision :: thresh
PROVIDE ao_bielec_integrals_erf_in_map ao_integrals_erf_map
thresh = ao_integrals_threshold
if (ao_overlap_abs(j,l) < thresh) then
out_val = 0.d0
return
endif
double precision :: get_ao_bielec_integral_erf
do i=1,sze
out_val(i) = get_ao_bielec_integral_erf(i,j,k,l,ao_integrals_erf_map)
enddo
end
subroutine get_ao_bielec_integrals_erf_non_zero(j,k,l,sze,out_val,out_val_index,non_zero_int)
use map_module
implicit none
BEGIN_DOC
! Gets multiple AO bi-electronic integral from the AO map .
! All non-zero i are retrieved for j,k,l fixed.
END_DOC
integer, intent(in) :: j,k,l, sze
real(integral_kind), intent(out) :: out_val(sze)
integer, intent(out) :: out_val_index(sze),non_zero_int
integer :: i
integer(key_kind) :: hash
double precision :: thresh,tmp
PROVIDE ao_bielec_integrals_erf_in_map
thresh = ao_integrals_threshold
non_zero_int = 0
if (ao_overlap_abs(j,l) < thresh) then
out_val = 0.d0
return
endif
non_zero_int = 0
do i=1,sze
integer, external :: ao_l4
double precision, external :: ao_bielec_integral_erf
!DIR$ FORCEINLINE
if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < thresh) then
cycle
endif
call bielec_integrals_index(i,j,k,l,hash)
call map_get(ao_integrals_erf_map, hash,tmp)
if (dabs(tmp) < thresh ) cycle
non_zero_int = non_zero_int+1
out_val_index(non_zero_int) = i
out_val(non_zero_int) = tmp
enddo
end
function get_ao_erf_map_size()
implicit none
integer (map_size_kind) :: get_ao_erf_map_size
BEGIN_DOC
! Returns the number of elements in the AO map
END_DOC
get_ao_erf_map_size = ao_integrals_erf_map % n_elements
end
subroutine clear_ao_erf_map
implicit none
BEGIN_DOC
! Frees the memory of the AO map
END_DOC
call map_deinit(ao_integrals_erf_map)
FREE ao_integrals_erf_map
end
BEGIN_TEMPLATE
subroutine dump_$ao_integrals(filename)
use map_module
implicit none
BEGIN_DOC
! Save to disk the $ao integrals
END_DOC
character*(*), intent(in) :: filename
integer(cache_key_kind), pointer :: key(:)
real(integral_kind), pointer :: val(:)
integer*8 :: i,j, n
call ezfio_set_work_empty(.False.)
open(unit=66,file=filename,FORM='unformatted')
write(66) integral_kind, key_kind
write(66) $ao_integrals_map%sorted, $ao_integrals_map%map_size, &
$ao_integrals_map%n_elements
do i=0_8,$ao_integrals_map%map_size
write(66) $ao_integrals_map%map(i)%sorted, $ao_integrals_map%map(i)%map_size,&
$ao_integrals_map%map(i)%n_elements
enddo
do i=0_8,$ao_integrals_map%map_size
key => $ao_integrals_map%map(i)%key
val => $ao_integrals_map%map(i)%value
n = $ao_integrals_map%map(i)%n_elements
write(66) (key(j), j=1,n), (val(j), j=1,n)
enddo
close(66)
end
IRP_IF COARRAY
subroutine communicate_$ao_integrals()
use map_module
implicit none
BEGIN_DOC
! Communicate the $ao integrals with co-array
END_DOC
integer(cache_key_kind), pointer :: key(:)
real(integral_kind), pointer :: val(:)
integer*8 :: i,j, k, nmax
integer*8, save :: n[*]
integer :: copy_n
real(integral_kind), allocatable :: buffer_val(:)[:]
integer(cache_key_kind), allocatable :: buffer_key(:)[:]
real(integral_kind), allocatable :: copy_val(:)
integer(key_kind), allocatable :: copy_key(:)
n = 0_8
do i=0_8,$ao_integrals_map%map_size
n = max(n,$ao_integrals_map%map(i)%n_elements)
enddo
sync all
nmax = 0_8
do j=1,num_images()
nmax = max(nmax,n[j])
enddo
allocate( buffer_key(nmax)[*], buffer_val(nmax)[*])
allocate( copy_key(nmax), copy_val(nmax))
do i=0_8,$ao_integrals_map%map_size
key => $ao_integrals_map%map(i)%key
val => $ao_integrals_map%map(i)%value
n = $ao_integrals_map%map(i)%n_elements
do j=1,n
buffer_key(j) = key(j)
buffer_val(j) = val(j)
enddo
sync all
do j=1,num_images()
if (j /= this_image()) then
copy_n = n[j]
do k=1,copy_n
copy_val(k) = buffer_val(k)[j]
copy_key(k) = buffer_key(k)[j]
copy_key(k) = copy_key(k)+ishft(i,-map_shift)
enddo
call map_append($ao_integrals_map, copy_key, copy_val, copy_n )
endif
enddo
sync all
enddo
deallocate( buffer_key, buffer_val, copy_val, copy_key)
end
IRP_ENDIF
integer function load_$ao_integrals(filename)
implicit none
BEGIN_DOC
! Read from disk the $ao integrals
END_DOC
character*(*), intent(in) :: filename
integer*8 :: i
integer(cache_key_kind), pointer :: key(:)
real(integral_kind), pointer :: val(:)
integer :: iknd, kknd
integer*8 :: n, j
load_$ao_integrals = 1
open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN')
read(66,err=98,end=98) iknd, kknd
if (iknd /= integral_kind) then
print *, 'Wrong integrals kind in file :', iknd
stop 1
endif
if (kknd /= key_kind) then
print *, 'Wrong key kind in file :', kknd
stop 1
endif
read(66,err=98,end=98) $ao_integrals_map%sorted, $ao_integrals_map%map_size,&
$ao_integrals_map%n_elements
do i=0_8, $ao_integrals_map%map_size
read(66,err=99,end=99) $ao_integrals_map%map(i)%sorted, &
$ao_integrals_map%map(i)%map_size, $ao_integrals_map%map(i)%n_elements
call cache_map_reallocate($ao_integrals_map%map(i),$ao_integrals_map%map(i)%map_size)
enddo
do i=0_8, $ao_integrals_map%map_size
key => $ao_integrals_map%map(i)%key
val => $ao_integrals_map%map(i)%value
n = $ao_integrals_map%map(i)%n_elements
read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n)
enddo
call map_sort($ao_integrals_map)
load_$ao_integrals = 0
return
99 continue
call map_deinit($ao_integrals_map)
98 continue
stop 'Problem reading $ao_integrals_map file in work/'
end
SUBST [ ao_integrals_map, ao_integrals, ao_num ]
ao_integrals_erf_map ; ao_integrals_erf ; ao_num ;;
mo_integrals_erf_map ; mo_integrals_erf ; mo_tot_num;;
END_TEMPLATE
BEGIN_PROVIDER [ type(map_type), mo_integrals_erf_map ]
implicit none
BEGIN_DOC
! MO integrals
END_DOC
integer(key_kind) :: key_max
integer(map_size_kind) :: sze
call bielec_integrals_index(mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num,key_max)
sze = key_max
call map_init(mo_integrals_erf_map,sze)
print*, 'MO ERF map initialized'
END_PROVIDER
subroutine insert_into_ao_integrals_erf_map(n_integrals,buffer_i, buffer_values)
use map_module
implicit none
BEGIN_DOC
! Create new entry into AO map
END_DOC
integer, intent(in) :: n_integrals
integer(key_kind), intent(inout) :: buffer_i(n_integrals)
real(integral_kind), intent(inout) :: buffer_values(n_integrals)
call map_append(ao_integrals_erf_map, buffer_i, buffer_values, n_integrals)
end
subroutine insert_into_mo_integrals_erf_map(n_integrals, &
buffer_i, buffer_values, thr)
use map_module
implicit none
BEGIN_DOC
! Create new entry into MO map, or accumulate in an existing entry
END_DOC
integer, intent(in) :: n_integrals
integer(key_kind), intent(inout) :: buffer_i(n_integrals)
real(integral_kind), intent(inout) :: buffer_values(n_integrals)
real(integral_kind), intent(in) :: thr
call map_update(mo_integrals_erf_map, buffer_i, buffer_values, n_integrals, thr)
end
BEGIN_PROVIDER [ integer, mo_integrals_erf_cache_min ]
&BEGIN_PROVIDER [ integer, mo_integrals_erf_cache_max ]
implicit none
BEGIN_DOC
! Min and max values of the MOs for which the integrals are in the cache
END_DOC
mo_integrals_erf_cache_min = max(1,elec_alpha_num - 31)
mo_integrals_erf_cache_max = min(mo_tot_num,mo_integrals_erf_cache_min+63)
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_integrals_erf_cache, (0:64*64*64*64) ]
implicit none
BEGIN_DOC
! Cache of MO integrals for fast access
END_DOC
PROVIDE mo_bielec_integrals_erf_in_map
integer :: i,j,k,l
integer :: ii
integer(key_kind) :: idx
real(integral_kind) :: integral
FREE ao_integrals_erf_cache
!$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
do l=mo_integrals_erf_cache_min,mo_integrals_erf_cache_max
do k=mo_integrals_erf_cache_min,mo_integrals_erf_cache_max
do j=mo_integrals_erf_cache_min,mo_integrals_erf_cache_max
do i=mo_integrals_erf_cache_min,mo_integrals_erf_cache_max
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(mo_integrals_erf_map,idx,integral)
ii = l-mo_integrals_erf_cache_min
ii = ior( ishft(ii,6), k-mo_integrals_erf_cache_min)
ii = ior( ishft(ii,6), j-mo_integrals_erf_cache_min)
ii = ior( ishft(ii,6), i-mo_integrals_erf_cache_min)
mo_integrals_erf_cache(ii) = integral
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
double precision function get_mo_bielec_integral_erf(i,j,k,l,map)
use map_module
implicit none
BEGIN_DOC
! Returns one integral <ij|kl> in the MO basis
END_DOC
integer, intent(in) :: i,j,k,l
integer(key_kind) :: idx
integer :: ii
type(map_type), intent(inout) :: map
real(integral_kind) :: tmp
PROVIDE mo_bielec_integrals_erf_in_map mo_integrals_erf_cache
ii = l-mo_integrals_erf_cache_min
ii = ior(ii, k-mo_integrals_erf_cache_min)
ii = ior(ii, j-mo_integrals_erf_cache_min)
ii = ior(ii, i-mo_integrals_erf_cache_min)
if (iand(ii, -64) /= 0) then
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(map,idx,tmp)
get_mo_bielec_integral_erf = dble(tmp)
else
ii = l-mo_integrals_erf_cache_min
ii = ior( ishft(ii,6), k-mo_integrals_erf_cache_min)
ii = ior( ishft(ii,6), j-mo_integrals_erf_cache_min)
ii = ior( ishft(ii,6), i-mo_integrals_erf_cache_min)
get_mo_bielec_integral_erf = mo_integrals_erf_cache(ii)
endif
end
double precision function mo_bielec_integral_erf(i,j,k,l)
implicit none
BEGIN_DOC
! Returns one integral <ij|kl> in the MO basis
END_DOC
integer, intent(in) :: i,j,k,l
double precision :: get_mo_bielec_integral_erf
PROVIDE mo_bielec_integrals_erf_in_map mo_integrals_erf_cache
!DIR$ FORCEINLINE
PROVIDE mo_bielec_integrals_erf_in_map
mo_bielec_integral_erf = get_mo_bielec_integral_erf(i,j,k,l,mo_integrals_erf_map)
return
end
subroutine get_mo_bielec_integrals_erf(j,k,l,sze,out_val,map)
use map_module
implicit none
BEGIN_DOC
! Returns multiple integrals <ij|kl> in the MO basis, all
! i for j,k,l fixed.
END_DOC
integer, intent(in) :: j,k,l, sze
double precision, intent(out) :: out_val(sze)
type(map_type), intent(inout) :: map
integer :: i
integer(key_kind) :: hash(sze)
real(integral_kind) :: tmp_val(sze)
PROVIDE mo_bielec_integrals_erf_in_map
do i=1,sze
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,hash(i))
enddo
if (key_kind == 8) then
call map_get_many(map, hash, out_val, sze)
else
call map_get_many(map, hash, tmp_val, sze)
! Conversion to double precision
do i=1,sze
out_val(i) = dble(tmp_val(i))
enddo
endif
end
subroutine get_mo_bielec_integrals_erf_ij(k,l,sze,out_array,map)
use map_module
implicit none
BEGIN_DOC
! Returns multiple integrals <ij|kl> in the MO basis, all
! i(1)j(2) 1/r12 k(1)l(2)
! i, j for k,l fixed.
END_DOC
integer, intent(in) :: k,l, sze
double precision, intent(out) :: out_array(sze,sze)
type(map_type), intent(inout) :: map
integer :: i,j,kk,ll,m
integer(key_kind),allocatable :: hash(:)
integer ,allocatable :: pairs(:,:), iorder(:)
real(integral_kind), allocatable :: tmp_val(:)
PROVIDE mo_bielec_integrals_erf_in_map
allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze), &
tmp_val(sze*sze))
kk=0
out_array = 0.d0
do j=1,sze
do i=1,sze
kk += 1
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,hash(kk))
pairs(1,kk) = i
pairs(2,kk) = j
iorder(kk) = kk
enddo
enddo
logical :: integral_is_in_map
if (key_kind == 8) then
call i8radix_sort(hash,iorder,kk,-1)
else if (key_kind == 4) then
call iradix_sort(hash,iorder,kk,-1)
else if (key_kind == 2) then
call i2radix_sort(hash,iorder,kk,-1)
endif
call map_get_many(mo_integrals_erf_map, hash, tmp_val, kk)
do ll=1,kk
m = iorder(ll)
i=pairs(1,m)
j=pairs(2,m)
out_array(i,j) = tmp_val(ll)
enddo
deallocate(pairs,hash,iorder,tmp_val)
end
subroutine get_mo_bielec_integrals_erf_i1j1(k,l,sze,out_array,map)
use map_module
implicit none
BEGIN_DOC
! Returns multiple integrals <ik|jl> in the MO basis, all
! i(1)j(1) erf(mu_erf * r12) /r12 k(2)l(2)
! i, j for k,l fixed.
END_DOC
integer, intent(in) :: k,l, sze
double precision, intent(out) :: out_array(sze,sze)
type(map_type), intent(inout) :: map
integer :: i,j,kk,ll,m
integer(key_kind),allocatable :: hash(:)
integer ,allocatable :: pairs(:,:), iorder(:)
real(integral_kind), allocatable :: tmp_val(:)
PROVIDE mo_bielec_integrals_erf_in_map
allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze), &
tmp_val(sze*sze))
kk=0
out_array = 0.d0
do j=1,sze
do i=1,sze
kk += 1
!DIR$ FORCEINLINE
call bielec_integrals_index(i,k,j,l,hash(kk))
pairs(1,kk) = i
pairs(2,kk) = j
iorder(kk) = kk
enddo
enddo
logical :: integral_is_in_map
if (key_kind == 8) then
call i8radix_sort(hash,iorder,kk,-1)
else if (key_kind == 4) then
call iradix_sort(hash,iorder,kk,-1)
else if (key_kind == 2) then
call i2radix_sort(hash,iorder,kk,-1)
endif
call map_get_many(mo_integrals_erf_map, hash, tmp_val, kk)
do ll=1,kk
m = iorder(ll)
i=pairs(1,m)
j=pairs(2,m)
out_array(i,j) = tmp_val(ll)
enddo
deallocate(pairs,hash,iorder,tmp_val)
end
subroutine get_mo_bielec_integrals_erf_coulomb_ii(k,l,sze,out_val,map)
use map_module
implicit none
BEGIN_DOC
! Returns multiple integrals <ki|li>
! k(1)i(2) 1/r12 l(1)i(2) :: out_val(i1)
! for k,l fixed.
END_DOC
integer, intent(in) :: k,l, sze
double precision, intent(out) :: out_val(sze)
type(map_type), intent(inout) :: map
integer :: i
integer(key_kind) :: hash(sze)
real(integral_kind) :: tmp_val(sze)
PROVIDE mo_bielec_integrals_erf_in_map
integer :: kk
do i=1,sze
!DIR$ FORCEINLINE
call bielec_integrals_index(k,i,l,i,hash(i))
enddo
if (key_kind == 8) then
call map_get_many(map, hash, out_val, sze)
else
call map_get_many(map, hash, tmp_val, sze)
! Conversion to double precision
do i=1,sze
out_val(i) = dble(tmp_val(i))
enddo
endif
end
subroutine get_mo_bielec_integrals_erf_exch_ii(k,l,sze,out_val,map)
use map_module
implicit none
BEGIN_DOC
! Returns multiple integrals <ki|il>
! k(1)i(2) 1/r12 i(1)l(2) :: out_val(i1)
! for k,l fixed.
END_DOC
integer, intent(in) :: k,l, sze
double precision, intent(out) :: out_val(sze)
type(map_type), intent(inout) :: map
integer :: i
integer(key_kind) :: hash(sze)
real(integral_kind) :: tmp_val(sze)
PROVIDE mo_bielec_integrals_erf_in_map
integer :: kk
do i=1,sze
!DIR$ FORCEINLINE
call bielec_integrals_index(k,i,i,l,hash(i))
enddo
if (key_kind == 8) then
call map_get_many(map, hash, out_val, sze)
else
call map_get_many(map, hash, tmp_val, sze)
! Conversion to double precision
do i=1,sze
out_val(i) = dble(tmp_val(i))
enddo
endif
end
integer*8 function get_mo_erf_map_size()
implicit none
BEGIN_DOC
! Return the number of elements in the MO map
END_DOC
get_mo_erf_map_size = mo_integrals_erf_map % n_elements
end

View File

@ -0,0 +1,549 @@
subroutine mo_bielec_integrals_erf_index(i,j,k,l,i1)
use map_module
implicit none
BEGIN_DOC
! Computes an unique index for i,j,k,l integrals
END_DOC
integer, intent(in) :: i,j,k,l
integer(key_kind), intent(out) :: i1
integer(key_kind) :: p,q,r,s,i2
p = min(i,k)
r = max(i,k)
p = p+ishft(r*r-r,-1)
q = min(j,l)
s = max(j,l)
q = q+ishft(s*s-s,-1)
i1 = min(p,q)
i2 = max(p,q)
i1 = i1+ishft(i2*i2-i2,-1)
end
BEGIN_PROVIDER [ logical, mo_bielec_integrals_erf_in_map ]
use map_module
implicit none
integer(bit_kind) :: mask_ijkl(N_int,4)
integer(bit_kind) :: mask_ijk(N_int,3)
BEGIN_DOC
! If True, the map of MO bielectronic integrals is provided
END_DOC
real :: map_mb
mo_bielec_integrals_erf_in_map = .True.
if (read_mo_integrals_erf) then
print*,'Reading the MO integrals_erf'
call map_load_from_disk(trim(ezfio_filename)//'/work/mo_ints_erf',mo_integrals_erf_map)
print*, 'MO integrals_erf provided'
return
else
PROVIDE ao_bielec_integrals_erf_in_map
endif
! call four_index_transform_block(ao_integrals_erf_map,mo_integrals_erf_map, &
! mo_coef, size(mo_coef,1), &
! 1, 1, 1, 1, ao_num, ao_num, ao_num, ao_num, &
! 1, 1, 1, 1, mo_num, mo_num, mo_num, mo_num)
call add_integrals_to_map_erf(full_ijkl_bitmask_4)
integer*8 :: get_mo_erf_map_size, mo_erf_map_size
mo_erf_map_size = get_mo_erf_map_size()
! print*,'Molecular integrals ERF provided:'
! print*,' Size of MO ERF map ', map_mb(mo_integrals_erf_map) ,'MB'
! print*,' Number of MO ERF integrals: ', mo_erf_map_size
if (write_mo_integrals_erf) then
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_erf',mo_integrals_erf_map)
call ezfio_set_integrals_bielec_erf_disk_access_mo_integrals_erf("Read")
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj_from_ao, (mo_tot_num,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj_exchange_from_ao, (mo_tot_num,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj_anti_from_ao, (mo_tot_num,mo_tot_num) ]
BEGIN_DOC
! mo_bielec_integral_jj_from_ao(i,j) = J_ij
! mo_bielec_integral_jj_exchange_from_ao(i,j) = J_ij
! mo_bielec_integral_jj_anti_from_ao(i,j) = J_ij - K_ij
END_DOC
implicit none
integer :: i,j,p,q,r,s
double precision :: c
real(integral_kind) :: integral
integer :: n, pp
real(integral_kind), allocatable :: int_value(:)
integer, allocatable :: int_idx(:)
double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:)
if (.not.do_direct_integrals) then
PROVIDE ao_bielec_integrals_erf_in_map mo_coef
endif
mo_bielec_integral_erf_jj_from_ao = 0.d0
mo_bielec_integral_erf_jj_exchange_from_ao = 0.d0
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: iqrs, iqsr
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE (i,j,p,q,r,s,integral,c,n,pp,int_value,int_idx, &
!$OMP iqrs, iqsr,iqri,iqis) &
!$OMP SHARED(mo_tot_num,mo_coef_transp,ao_num,&
!$OMP ao_integrals_threshold,do_direct_integrals) &
!$OMP REDUCTION(+:mo_bielec_integral_erf_jj_from_ao,mo_bielec_integral_erf_jj_exchange_from_ao)
allocate( int_value(ao_num), int_idx(ao_num), &
iqrs(mo_tot_num,ao_num), iqis(mo_tot_num), iqri(mo_tot_num),&
iqsr(mo_tot_num,ao_num) )
!$OMP DO SCHEDULE (guided)
do s=1,ao_num
do q=1,ao_num
do j=1,ao_num
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqrs(i,j) = 0.d0
iqsr(i,j) = 0.d0
enddo
enddo
if (do_direct_integrals) then
double precision :: ao_bielec_integral_erf
do r=1,ao_num
call compute_ao_bielec_integrals_erf(q,r,s,ao_num,int_value)
do p=1,ao_num
integral = int_value(p)
if (abs(integral) > ao_integrals_threshold) then
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqrs(i,r) += mo_coef_transp(i,p) * integral
enddo
endif
enddo
call compute_ao_bielec_integrals_erf(q,s,r,ao_num,int_value)
do p=1,ao_num
integral = int_value(p)
if (abs(integral) > ao_integrals_threshold) then
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqsr(i,r) += mo_coef_transp(i,p) * integral
enddo
endif
enddo
enddo
else
do r=1,ao_num
call get_ao_bielec_integrals_erf_non_zero(q,r,s,ao_num,int_value,int_idx,n)
do pp=1,n
p = int_idx(pp)
integral = int_value(pp)
if (abs(integral) > ao_integrals_threshold) then
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqrs(i,r) += mo_coef_transp(i,p) * integral
enddo
endif
enddo
call get_ao_bielec_integrals_erf_non_zero(q,s,r,ao_num,int_value,int_idx,n)
do pp=1,n
p = int_idx(pp)
integral = int_value(pp)
if (abs(integral) > ao_integrals_threshold) then
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqsr(i,r) += mo_coef_transp(i,p) * integral
enddo
endif
enddo
enddo
endif
iqis = 0.d0
iqri = 0.d0
do r=1,ao_num
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqis(i) += mo_coef_transp(i,r) * iqrs(i,r)
iqri(i) += mo_coef_transp(i,r) * iqsr(i,r)
enddo
enddo
do i=1,mo_tot_num
!DIR$ VECTOR ALIGNED
do j=1,mo_tot_num
c = mo_coef_transp(j,q)*mo_coef_transp(j,s)
mo_bielec_integral_erf_jj_from_ao(j,i) += c * iqis(i)
mo_bielec_integral_erf_jj_exchange_from_ao(j,i) += c * iqri(i)
enddo
enddo
enddo
enddo
!$OMP END DO NOWAIT
deallocate(iqrs,iqsr,int_value,int_idx)
!$OMP END PARALLEL
mo_bielec_integral_erf_jj_anti_from_ao = mo_bielec_integral_erf_jj_from_ao - mo_bielec_integral_erf_jj_exchange_from_ao
! end
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj, (mo_tot_num,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj_exchange, (mo_tot_num,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj_anti, (mo_tot_num,mo_tot_num) ]
implicit none
BEGIN_DOC
! mo_bielec_integral_jj(i,j) = J_ij
! mo_bielec_integral_jj_exchange(i,j) = K_ij
! mo_bielec_integral_jj_anti(i,j) = J_ij - K_ij
END_DOC
integer :: i,j
double precision :: get_mo_bielec_integral_erf
PROVIDE mo_bielec_integrals_erf_in_map
mo_bielec_integral_erf_jj = 0.d0
mo_bielec_integral_erf_jj_exchange = 0.d0
do j=1,mo_tot_num
do i=1,mo_tot_num
mo_bielec_integral_erf_jj(i,j) = get_mo_bielec_integral_erf(i,j,i,j,mo_integrals_erf_map)
mo_bielec_integral_erf_jj_exchange(i,j) = get_mo_bielec_integral_erf(i,j,j,i,mo_integrals_erf_map)
mo_bielec_integral_erf_jj_anti(i,j) = mo_bielec_integral_erf_jj(i,j) - mo_bielec_integral_erf_jj_exchange(i,j)
enddo
enddo
END_PROVIDER
subroutine clear_mo_erf_map
implicit none
BEGIN_DOC
! Frees the memory of the MO map
END_DOC
call map_deinit(mo_integrals_erf_map)
FREE mo_integrals_erf_map mo_bielec_integral_erf_jj mo_bielec_integral_erf_jj_anti
FREE mo_bielec_integral_Erf_jj_exchange mo_bielec_integrals_erf_in_map
end
subroutine provide_all_mo_integrals_erf
implicit none
provide mo_integrals_erf_map mo_bielec_integral_erf_jj mo_bielec_integral_erf_jj_anti
provide mo_bielec_integral_erf_jj_exchange mo_bielec_integrals_erf_in_map
end
subroutine add_integrals_to_map_erf(mask_ijkl)
use bitmasks
implicit none
BEGIN_DOC
! Adds integrals to tha MO map according to some bitmask
END_DOC
integer(bit_kind), intent(in) :: mask_ijkl(N_int,4)
integer :: i,j,k,l
integer :: i0,j0,k0,l0
double precision :: c, cpu_1, cpu_2, wall_1, wall_2, wall_0
integer, allocatable :: list_ijkl(:,:)
integer :: n_i, n_j, n_k, n_l
integer, allocatable :: bielec_tmp_0_idx(:)
real(integral_kind), allocatable :: bielec_tmp_0(:,:)
double precision, allocatable :: bielec_tmp_1(:)
double precision, allocatable :: bielec_tmp_2(:,:)
double precision, allocatable :: bielec_tmp_3(:,:,:)
!DIR$ ATTRIBUTES ALIGN : 64 :: bielec_tmp_1, bielec_tmp_2, bielec_tmp_3
integer :: n_integrals
integer :: size_buffer
integer(key_kind),allocatable :: buffer_i(:)
real(integral_kind),allocatable :: buffer_value(:)
double precision :: map_mb
integer :: i1,j1,k1,l1, ii1, kmax, thread_num
integer :: i2,i3,i4
double precision,parameter :: thr_coef = 1.d-10
PROVIDE ao_bielec_integrals_in_map mo_coef
!Get list of MOs for i,j,k and l
!-------------------------------
allocate(list_ijkl(mo_tot_num,4))
call bitstring_to_list( mask_ijkl(1,1), list_ijkl(1,1), n_i, N_int )
call bitstring_to_list( mask_ijkl(1,2), list_ijkl(1,2), n_j, N_int )
call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int )
call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int )
character*(2048) :: output(1)
print*, 'i'
call bitstring_to_str( output(1), mask_ijkl(1,1), N_int )
print *, trim(output(1))
j = 0
do i = 1, N_int
j += popcnt(mask_ijkl(i,1))
enddo
if(j==0)then
return
endif
print*, 'j'
call bitstring_to_str( output(1), mask_ijkl(1,2), N_int )
print *, trim(output(1))
j = 0
do i = 1, N_int
j += popcnt(mask_ijkl(i,2))
enddo
if(j==0)then
return
endif
print*, 'k'
call bitstring_to_str( output(1), mask_ijkl(1,3), N_int )
print *, trim(output(1))
j = 0
do i = 1, N_int
j += popcnt(mask_ijkl(i,3))
enddo
if(j==0)then
return
endif
print*, 'l'
call bitstring_to_str( output(1), mask_ijkl(1,4), N_int )
print *, trim(output(1))
j = 0
do i = 1, N_int
j += popcnt(mask_ijkl(i,4))
enddo
if(j==0)then
return
endif
size_buffer = min(ao_num*ao_num*ao_num,16000000)
print*, 'Providing the ERF molecular integrals '
print*, 'Buffers : ', 8.*(mo_tot_num*(n_j)*(n_k+1) + mo_tot_num+&
ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core'
call wall_time(wall_1)
call cpu_time(cpu_1)
double precision :: accu_bis
accu_bis = 0.d0
!$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, &
!$OMP bielec_tmp_0_idx, bielec_tmp_0, bielec_tmp_1,bielec_tmp_2,bielec_tmp_3,&
!$OMP buffer_i,buffer_value,n_integrals,wall_2,i0,j0,k0,l0, &
!$OMP wall_0,thread_num,accu_bis) &
!$OMP DEFAULT(NONE) &
!$OMP SHARED(size_buffer,ao_num,mo_tot_num,n_i,n_j,n_k,n_l, &
!$OMP mo_coef_transp, &
!$OMP mo_coef_transp_is_built, list_ijkl, &
!$OMP mo_coef_is_built, wall_1, &
!$OMP mo_coef,mo_integrals_threshold,mo_integrals_erf_map)
n_integrals = 0
wall_0 = wall_1
allocate(bielec_tmp_3(mo_tot_num, n_j, n_k), &
bielec_tmp_1(mo_tot_num), &
bielec_tmp_0(ao_num,ao_num), &
bielec_tmp_0_idx(ao_num), &
bielec_tmp_2(mo_tot_num, n_j), &
buffer_i(size_buffer), &
buffer_value(size_buffer) )
thread_num = 0
!$ thread_num = omp_get_thread_num()
!$OMP DO SCHEDULE(guided)
do l1 = 1,ao_num
bielec_tmp_3 = 0.d0
do k1 = 1,ao_num
bielec_tmp_2 = 0.d0
do j1 = 1,ao_num
call get_ao_bielec_integrals_erf(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) ! all integrals for a given l1, k1
! call compute_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1))
enddo
do j1 = 1,ao_num
kmax = 0
do i1 = 1,ao_num
c = bielec_tmp_0(i1,j1)
if (c == 0.d0) then
cycle
endif
kmax += 1
bielec_tmp_0(kmax,j1) = c
bielec_tmp_0_idx(kmax) = i1
enddo
if (kmax==0) then
cycle
endif
bielec_tmp_1 = 0.d0
ii1=1
! sum_m c_m^i (m)
do ii1 = 1,kmax-4,4
i1 = bielec_tmp_0_idx(ii1)
i2 = bielec_tmp_0_idx(ii1+1)
i3 = bielec_tmp_0_idx(ii1+2)
i4 = bielec_tmp_0_idx(ii1+3)
do i = list_ijkl(1,1), list_ijkl(n_i,1)
bielec_tmp_1(i) = bielec_tmp_1(i) + &
mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + &
mo_coef_transp(i,i2) * bielec_tmp_0(ii1+1,j1) + &
mo_coef_transp(i,i3) * bielec_tmp_0(ii1+2,j1) + &
mo_coef_transp(i,i4) * bielec_tmp_0(ii1+3,j1)
enddo ! i
enddo ! ii1
i2 = ii1
do ii1 = i2,kmax
i1 = bielec_tmp_0_idx(ii1)
do i = list_ijkl(1,1), list_ijkl(n_i,1)
bielec_tmp_1(i) = bielec_tmp_1(i) + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1)
enddo ! i
enddo ! ii1
c = 0.d0
do i = list_ijkl(1,1), list_ijkl(n_i,1)
c = max(c,abs(bielec_tmp_1(i)))
if (c>mo_integrals_threshold) exit
enddo
if ( c < mo_integrals_threshold ) then
cycle
endif
do j0 = 1, n_j
j = list_ijkl(j0,2)
c = mo_coef_transp(j,j1)
if (abs(c) < thr_coef) then
cycle
endif
do i = list_ijkl(1,1), list_ijkl(n_i,1)
bielec_tmp_2(i,j0) = bielec_tmp_2(i,j0) + c * bielec_tmp_1(i)
enddo ! i
enddo ! j
enddo !j1
if ( maxval(abs(bielec_tmp_2)) < mo_integrals_threshold ) then
cycle
endif
do k0 = 1, n_k
k = list_ijkl(k0,3)
c = mo_coef_transp(k,k1)
if (abs(c) < thr_coef) then
cycle
endif
do j0 = 1, n_j
j = list_ijkl(j0,2)
do i = list_ijkl(1,1), k
bielec_tmp_3(i,j0,k0) = bielec_tmp_3(i,j0,k0) + c* bielec_tmp_2(i,j0)
enddo!i
enddo !j
enddo !k
enddo !k1
do l0 = 1,n_l
l = list_ijkl(l0,4)
c = mo_coef_transp(l,l1)
if (abs(c) < thr_coef) then
cycle
endif
j1 = ishft((l*l-l),-1)
do j0 = 1, n_j
j = list_ijkl(j0,2)
if (j > l) then
exit
endif
j1 += 1
do k0 = 1, n_k
k = list_ijkl(k0,3)
i1 = ishft((k*k-k),-1)
if (i1<=j1) then
continue
else
exit
endif
bielec_tmp_1 = 0.d0
do i0 = 1, n_i
i = list_ijkl(i0,1)
if (i>k) then
exit
endif
bielec_tmp_1(i) = c*bielec_tmp_3(i,j0,k0)
! i1+=1
enddo
do i0 = 1, n_i
i = list_ijkl(i0,1)
if(i> min(k,j1-i1+list_ijkl(1,1)-1))then
exit
endif
if (abs(bielec_tmp_1(i)) < mo_integrals_threshold) then
cycle
endif
n_integrals += 1
buffer_value(n_integrals) = bielec_tmp_1(i)
!DIR$ FORCEINLINE
call mo_bielec_integrals_index(i,j,k,l,buffer_i(n_integrals))
if (n_integrals == size_buffer) then
call insert_into_mo_integrals_erf_map(n_integrals,buffer_i,buffer_value,&
real(mo_integrals_threshold,integral_kind))
n_integrals = 0
endif
enddo
enddo
enddo
enddo
call wall_time(wall_2)
if (thread_num == 0) then
if (wall_2 - wall_0 > 1.d0) then
wall_0 = wall_2
print*, 100.*float(l1)/float(ao_num), '% in ', &
wall_2-wall_1, 's', map_mb(mo_integrals_erf_map) ,'MB'
endif
endif
enddo
!$OMP END DO NOWAIT
deallocate (bielec_tmp_1,bielec_tmp_2,bielec_tmp_3)
integer :: index_needed
call insert_into_mo_integrals_erf_map(n_integrals,buffer_i,buffer_value,&
real(mo_integrals_threshold,integral_kind))
deallocate(buffer_i, buffer_value)
!$OMP END PARALLEL
call map_merge(mo_integrals_erf_map)
call wall_time(wall_2)
call cpu_time(cpu_2)
integer*8 :: get_mo_erf_map_size, mo_map_size
mo_map_size = get_mo_erf_map_size()
deallocate(list_ijkl)
print*,'Molecular ERF integrals provided:'
print*,' Size of MO ERF map ', map_mb(mo_integrals_erf_map) ,'MB'
print*,' Number of MO ERF integrals: ', mo_map_size
print*,' cpu time :',cpu_2 - cpu_1, 's'
print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')'
end

View File

@ -0,0 +1,125 @@
BEGIN_PROVIDER [ logical, ao_bielec_integrals_erf_in_map ]
implicit none
use f77_zmq
use map_module
BEGIN_DOC
! Map of Atomic integrals
! i(r1) j(r2) 1/r12 k(r1) l(r2)
END_DOC
integer :: i,j,k,l
double precision :: ao_bielec_integral_erf,cpu_1,cpu_2, wall_1, wall_2
double precision :: integral, wall_0
include 'utils/constants.include.F'
! For integrals file
integer(key_kind),allocatable :: buffer_i(:)
integer,parameter :: size_buffer = 1024*64
real(integral_kind),allocatable :: buffer_value(:)
integer :: n_integrals, rc
integer :: kk, m, j1, i1, lmax
character*(64) :: fmt
integral = ao_bielec_integral_erf(1,1,1,1)
double precision :: map_mb
PROVIDE read_ao_integrals_erf disk_access_ao_integrals_erf
if (read_ao_integrals_erf) then
print*,'Reading the AO ERF integrals'
call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map)
print*, 'AO ERF integrals provided'
ao_bielec_integrals_erf_in_map = .True.
return
endif
print*, 'Providing the AO ERF integrals'
call wall_time(wall_0)
call wall_time(wall_1)
call cpu_time(cpu_1)
integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull
call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'ao_integrals_erf')
character(len=:), allocatable :: task
allocate(character(len=ao_num*12) :: task)
write(fmt,*) '(', ao_num, '(I5,X,I5,''|''))'
do l=1,ao_num
write(task,fmt) (i,l, i=1,l)
integer, external :: add_task_to_taskserver
if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task)) == -1) then
stop 'Unable to add task to server'
endif
enddo
deallocate(task)
integer, external :: zmq_set_running
if (zmq_set_running(zmq_to_qp_run_socket) == -1) then
print *, irp_here, ': Failed in zmq_set_running'
endif
PROVIDE nproc
!$OMP PARALLEL DEFAULT(shared) private(i) num_threads(nproc+1)
i = omp_get_thread_num()
if (i==0) then
call ao_bielec_integrals_erf_in_map_collector(zmq_socket_pull)
else
call ao_bielec_integrals_erf_in_map_slave_inproc(i)
endif
!$OMP END PARALLEL
call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'ao_integrals_erf')
print*, 'Sorting the map'
call map_sort(ao_integrals_erf_map)
call cpu_time(cpu_2)
call wall_time(wall_2)
integer(map_size_kind) :: get_ao_erf_map_size, ao_erf_map_size
ao_erf_map_size = get_ao_erf_map_size()
print*, 'AO ERF integrals provided:'
print*, ' Size of AO ERF map : ', map_mb(ao_integrals_erf_map) ,'MB'
print*, ' Number of AO ERF integrals :', ao_erf_map_size
print*, ' cpu time :',cpu_2 - cpu_1, 's'
print*, ' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )'
ao_bielec_integrals_erf_in_map = .True.
if (write_ao_integrals_erf) then
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map)
call ezfio_set_integrals_bielec_erf_disk_access_ao_integrals_erf("Read")
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_bielec_integral_erf_schwartz,(ao_num,ao_num) ]
implicit none
BEGIN_DOC
! Needed to compute Schwartz inequalities
END_DOC
integer :: i,k
double precision :: ao_bielec_integral_erf,cpu_1,cpu_2, wall_1, wall_2
ao_bielec_integral_erf_schwartz(1,1) = ao_bielec_integral_erf(1,1,1,1)
!$OMP PARALLEL DO PRIVATE(i,k) &
!$OMP DEFAULT(NONE) &
!$OMP SHARED (ao_num,ao_bielec_integral_erf_schwartz) &
!$OMP SCHEDULE(dynamic)
do i=1,ao_num
do k=1,i
ao_bielec_integral_erf_schwartz(i,k) = dsqrt(ao_bielec_integral_erf(i,k,i,k))
ao_bielec_integral_erf_schwartz(k,i) = ao_bielec_integral_erf_schwartz(i,k)
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER

View File

@ -0,0 +1,32 @@
program qp_ao_ints
use omp_lib
implicit none
BEGIN_DOC
! Increments a running calculation to compute AO integral_erfs
END_DOC
integer :: i
call switch_qp_run_to_master
zmq_context = f77_zmq_ctx_new ()
! Set the state of the ZMQ
zmq_state = 'ao_integral_erfs'
! Provide everything needed
double precision :: integral_erf, ao_bielec_integral_erf
integral_erf = ao_bielec_integral_erf(1,1,1,1)
character*(64) :: state
call wait_for_state(zmq_state,state)
do while (state /= 'Stopped')
!$OMP PARALLEL DEFAULT(PRIVATE) PRIVATE(i)
i = omp_get_thread_num()
call ao_bielec_integrals_erf_in_map_slave_tcp(i)
!$OMP END PARALLEL
call wait_for_state(zmq_state,state)
enddo
print *, 'Done'
end

View File

@ -0,0 +1,47 @@
BEGIN_PROVIDER [ logical, read_ao_integrals_erf ]
&BEGIN_PROVIDER [ logical, read_mo_integrals_erf ]
&BEGIN_PROVIDER [ logical, write_ao_integrals_erf ]
&BEGIN_PROVIDER [ logical, write_mo_integrals_erf ]
BEGIN_DOC
! One level of abstraction for disk_access_ao_integrals_erf and disk_access_mo_integrals_erf
END_DOC
implicit none
if (disk_access_ao_integrals_erf.EQ.'Read') then
read_ao_integrals_erf = .True.
write_ao_integrals_erf = .False.
else if (disk_access_ao_integrals_erf.EQ.'Write') then
read_ao_integrals_erf = .False.
write_ao_integrals_erf = .True.
else if (disk_access_ao_integrals_erf.EQ.'None') then
read_ao_integrals_erf = .False.
write_ao_integrals_erf = .False.
else
print *, 'bielec_integrals_erf/disk_access_ao_integrals_erf has a wrong type'
stop 1
endif
if (disk_access_mo_integrals_erf.EQ.'Read') then
read_mo_integrals_erf = .True.
write_mo_integrals_erf = .False.
else if (disk_access_mo_integrals_erf.EQ.'Write') then
read_mo_integrals_erf = .False.
write_mo_integrals_erf = .True.
else if (disk_access_mo_integrals_erf.EQ.'None') then
read_mo_integrals_erf = .False.
write_mo_integrals_erf = .False.
else
print *, 'bielec_integrals_erf/disk_access_mo_integrals_erf has a wrong type'
stop 1
endif
END_PROVIDER

View File

@ -0,0 +1,36 @@
subroutine save_erf_bi_elec_integrals_mo
implicit none
integer :: i,j,k,l
PROVIDE mo_bielec_integrals_erf_in_map
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_erf',mo_integrals_erf_map)
call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read")
end
subroutine save_erf_bi_elec_integrals_ao
implicit none
integer :: i,j,k,l
PROVIDE ao_bielec_integrals_erf_in_map
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map)
call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read")
end
subroutine save_erf_bielec_ints_mo_into_ints_mo
implicit none
integer :: i,j,k,l
PROVIDE mo_bielec_integrals_erf_in_map
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_erf_map)
call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read")
end
subroutine save_erf_bielec_ints_mo_into_ints_ao
implicit none
integer :: i,j,k,l
PROVIDE ao_bielec_integrals_erf_in_map
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_erf_map)
call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read")
end

View File

@ -0,0 +1,20 @@
program write_integrals
implicit none
BEGIN_DOC
! This program saves the bielec erf integrals into the EZFIO folder
END_DOC
disk_access_mo_integrals = "None"
touch disk_access_mo_integrals
disk_access_ao_integrals = "None"
touch disk_access_ao_integrals
call routine
end
subroutine routine
implicit none
call save_erf_bi_elec_integrals_ao
call save_erf_bi_elec_integrals_mo
end

View File

@ -0,0 +1,21 @@
program write_integrals_erf_to_regular_integrals
implicit none
BEGIN_DOC
! This program saves the bielec erf integrals into the EZFIO folder but at the regular bielec integrals place.
! Threfore, if the user runs a future calculation with a reading of the integrals, the calculation will be performed with the erf bielec integrals instead of the regular bielec integrals
END_DOC
disk_access_mo_integrals = "None"
touch disk_access_mo_integrals
disk_access_ao_integrals = "None"
touch disk_access_ao_integrals
call routine
end
subroutine routine
implicit none
call save_erf_bielec_ints_mo_into_ints_ao
call save_erf_bielec_ints_mo_into_ints_mo
end

View File

@ -0,0 +1,55 @@
subroutine give_all_mos_at_r(r,mos_array)
implicit none
double precision, intent(in) :: r(3)
double precision, intent(out) :: mos_array(mo_tot_num)
double precision :: aos_array(ao_num)
call give_all_aos_at_r(r,aos_array)
call dgemv('N',mo_tot_num,ao_num,1.d0,mo_coef_transp,mo_tot_num,aos_array,1,0.d0,mos_array,1)
end
subroutine give_all_mos_and_grad_at_r(r,mos_array,mos_grad_array)
implicit none
double precision, intent(in) :: r(3)
double precision, intent(out) :: mos_array(mo_tot_num)
double precision, intent(out) :: mos_grad_array(mo_tot_num,3)
integer :: i,j,k
double precision :: aos_array(ao_num),aos_grad_array(ao_num,3)
call give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
mos_array=0d0
mos_grad_array=0d0
do j = 1, mo_tot_num
do k=1, ao_num
mos_array(j) += mo_coef(k,j)*aos_array(k)
mos_grad_array(j,1) += mo_coef(k,j)*aos_grad_array(k,1)
mos_grad_array(j,2) += mo_coef(k,j)*aos_grad_array(k,2)
mos_grad_array(j,3) += mo_coef(k,j)*aos_grad_array(k,3)
enddo
enddo
end
subroutine give_all_mos_and_grad_and_lapl_at_r(r,mos_array,mos_grad_array,mos_lapl_array)
implicit none
double precision, intent(in) :: r(3)
double precision, intent(out) :: mos_array(mo_tot_num)
double precision, intent(out) :: mos_grad_array(mo_tot_num,3),mos_lapl_array(mo_tot_num,3)
integer :: i,j,k
double precision :: aos_array(ao_num),aos_grad_array(ao_num,3),aos_lapl_array(ao_num,3)
call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
mos_array=0d0
mos_grad_array=0d0
mos_lapl_array=0d0
do j = 1, mo_tot_num
do k=1, ao_num
mos_array(j) += mo_coef(k,j)*aos_array(k)
mos_grad_array(j,1) += mo_coef(k,j)*aos_grad_array(k,1)
mos_grad_array(j,2) += mo_coef(k,j)*aos_grad_array(k,2)
mos_grad_array(j,3) += mo_coef(k,j)*aos_grad_array(k,3)
mos_lapl_array(j,1) += mo_coef(k,j)*aos_lapl_array(k,1)
mos_lapl_array(j,2) += mo_coef(k,j)*aos_lapl_array(k,2)
mos_lapl_array(j,3) += mo_coef(k,j)*aos_lapl_array(k,3)
enddo
enddo
end

View File

@ -90,6 +90,7 @@ END_PROVIDER
&BEGIN_PROVIDER [ double precision, nucl_dist_vec_y, (nucl_num,nucl_num) ]
&BEGIN_PROVIDER [ double precision, nucl_dist_vec_z, (nucl_num,nucl_num) ]
&BEGIN_PROVIDER [ double precision, nucl_dist, (nucl_num,nucl_num) ]
&BEGIN_PROVIDER [ double precision, nucl_dist_inv, (nucl_num,nucl_num) ]
implicit none
BEGIN_DOC
! nucl_dist : Nucleus-nucleus distances
@ -97,6 +98,8 @@ END_PROVIDER
! nucl_dist_2 : Nucleus-nucleus distances squared
! nucl_dist_vec : Nucleus-nucleus distances vectors
!
! nucl_dist_inv(I,J) = inversed of the distance between nucleus I and nucleus J
END_DOC
integer :: ie1, ie2, l
@ -124,6 +127,13 @@ END_PROVIDER
ASSERT (nucl_dist(ie1,ie2) > 0.d0)
enddo
enddo
do ie1 = 1, nucl_num
do ie2 = 1, nucl_num
if(ie1 == ie2)cycle
nucl_dist_inv(ie2,ie1) = 1.d0/nucl_dist(ie2,ie1)
enddo
enddo
END_PROVIDER

View File

@ -10,3 +10,9 @@ double precision, parameter :: dtwo_pi = 2.d0*dacos(-1.d0)
double precision, parameter :: inv_sq_pi = 1.d0/dsqrt(dacos(-1.d0))
double precision, parameter :: inv_sq_pi_2 = 0.5d0/dsqrt(dacos(-1.d0))
double precision, parameter :: thresh = 1.d-15
double precision, parameter :: cx_lda = -0.73855876638202234d0
double precision, parameter :: c_2_4_3 = 2.5198420997897464d0
double precision, parameter :: cst_lda = -0.93052573634909996d0
double precision, parameter :: c_4_3 = 1.3333333333333333d0
double precision, parameter :: c_1_3 = 0.3333333333333333d0