LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-09-27 03:51:01 +02:00
Code Issues Releases Wiki Activity

Added @eginer's becke_numerical_grid

Browse Source
This commit is contained in:
Anthony Scemama 2018-12-18 18:14:17 +01:00
parent 4b21b8cd0e
commit 3b8f0c3338
11 changed files with 7409 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
TODO
View File

@ -1,3 +1,10 @@
* Manu : README
* data_energy_and_density
* dft_keywords
* dft_quantities_on_grid
* dft_utils_one_body
* dm_for_dft
* integrals_erf
* Molden format
* Virer tous les modules qui sont dans plugins
* Permettre aux utilisateurs de facilement deposer des plugins dans plugins via une commande
@ -7,7 +14,6 @@
* Creer une page web pas trop degueu et la mettre ici : http://lcpq.github.io/quantum_package
* Un module pour lire les integrales Moleculaires depuis un FCIDUMP
* Un module pour lire des integrales Atomiques (voir module de Mimi pour lire les AO Slater)
* Include la DFT si c'est propre
* Plus de tests:
@ -18,7 +24,7 @@
* Davidson
* Lapack
* Extrapolation
* DFT?
* DFT
* User doc:

1
data/module_gitignore
View File

@ -6,6 +6,7 @@ tags
Makefile
ao_basis
ao_one_e_integrals
becke_numerical_grid
bitmask
cis
cisd

16
src/becke_numerical_grid/EZFIO.cfg Normal file
View File

@ -0,0 +1,16 @@
[n_points_integration_angular]
type: integer
doc: Number of angular points per atom for 3d numerical integration, needed for DFT for example [ 50 | 74 | 266 | 590 | 1202 | 2030 | 5810 ]
interface: ezfio,provider,ocaml
default: 590
[n_points_radial_grid]
type: integer
doc: Number of radial points per atom for 3d numerical integration, needed for DFT for example
interface: ezfio,provider,ocaml
default: 60

1
src/becke_numerical_grid/NEED Normal file
View File

@ -0,0 +1 @@
Nuclei

52
src/becke_numerical_grid/README.rst Normal file
View File

@ -0,0 +1,52 @@
====================
Becke Numerical Grid
====================
Lebedev-Laikov grids, using the code distributed through CCL (http://www.ccl.net/).
.. code-block:: text
This subroutine is part of a set of subroutines that generate
Lebedev grids [1-6] for integration on a sphere. The original
C-code [1] was kindly provided by Dr. Dmitri N. Laikov and
translated into fortran by Dr. Christoph van Wuellen.
This subroutine was translated using a C to fortran77 conversion
tool written by Dr. Christoph van Wuellen.
Users of this code are asked to include reference [1] in their
publications, and in the user- and programmers-manuals
describing their codes.
This code was distributed through CCL (http://www.ccl.net/).
[1] V.I. Lebedev, and D.N. Laikov
"A quadrature formula for the sphere of the 131st
algebraic order of accuracy"
Doklady Mathematics, Vol. 59, No. 3, 1999, pp. 477-481.
[2] V.I. Lebedev
"A quadrature formula for the sphere of 59th algebraic
order of accuracy"
Russian Acad. Sci. Dokl. Math., Vol. 50, 1995, pp. 283-286.
[3] V.I. Lebedev, and A.L. Skorokhodov
"Quadrature formulas of orders 41, 47, and 53 for the sphere"
Russian Acad. Sci. Dokl. Math., Vol. 45, 1992, pp. 587-592.
[4] V.I. Lebedev
"Spherical quadrature formulas exact to orders 25-29"
Siberian Mathematical Journal, Vol. 18, 1977, pp. 99-107.
[5] V.I. Lebedev
"Quadratures on a sphere"
Computational Mathematics and Mathematical Physics, Vol. 16,
1976, pp. 10-24.
[6] V.I. Lebedev
"Values of the nodes and weights of ninth to seventeenth
order Gauss-Markov quadrature formulae invariant under the
octahedron group with inversion"
Computational Mathematics and Mathematical Physics, Vol. 15,
1975, pp. 44-51.

6951
src/becke_numerical_grid/angular.f Normal file
View File

File diff suppressed because it is too large Load Diff

192
src/becke_numerical_grid/grid_becke.irp.f Normal file
View File

@ -0,0 +1,192 @@
BEGIN_PROVIDER [integer, n_points_grid_per_atom]
implicit none
BEGIN_DOC
! Number of grid points per atom
END_DOC
n_points_grid_per_atom = n_points_integration_angular * n_points_radial_grid
END_PROVIDER
BEGIN_PROVIDER [double precision, angular_quadrature_points, (n_points_integration_angular,3) ]
&BEGIN_PROVIDER [double precision, weights_angular_points, (n_points_integration_angular)]
implicit none
BEGIN_DOC
! weights and grid points for the integration on the angular variables on
! the unit sphere centered on (0,0,0)
! According to the LEBEDEV scheme
END_DOC
include 'constants.include.F'
integer :: i
double precision :: accu
double precision :: degre_rad
double precision :: x(n_points_integration_angular)
double precision :: y(n_points_integration_angular)
double precision :: z(n_points_integration_angular)
double precision :: w(n_points_integration_angular)
degre_rad = pi/180.d0
accu = 0.d0
select case (n_points_integration_angular)
case (5810)
call LD5810(X,Y,Z,W,n_points_integration_angular)
case (2030)
call LD2030(X,Y,Z,W,n_points_integration_angular)
case (1202)
call LD1202(X,Y,Z,W,n_points_integration_angular)
case (0590)
call LD0590(X,Y,Z,W,n_points_integration_angular)
case (266)
call LD0266(X,Y,Z,W,n_points_integration_angular)
case (74)
call LD0074(X,Y,Z,W,n_points_integration_angular)
case (50)
call LD0050(X,Y,Z,W,n_points_integration_angular)
case default
print *, irp_here//': wrong n_points_integration_angular. Expected:'
print *, '[ 50 | 74 | 266 | 590 | 1202 | 2030 | 5810 ]'
stop -1
end select
do i = 1, n_points_integration_angular
angular_quadrature_points(i,1) = x(i)
angular_quadrature_points(i,2) = y(i)
angular_quadrature_points(i,3) = z(i)
weights_angular_points(i) = w(i) * 4.d0 * pi
accu += w(i)
enddo
END_PROVIDER
BEGIN_PROVIDER [integer , m_knowles]
implicit none
BEGIN_DOC
! value of the "m" parameter in the equation (7) of the paper of Knowles (JCP, 104, 1996)
END_DOC
m_knowles = 3
END_PROVIDER
BEGIN_PROVIDER [double precision, grid_points_radial, (n_points_radial_grid)]
&BEGIN_PROVIDER [double precision, dr_radial_integral]
implicit none
BEGIN_DOC
! points in [0,1] to map the radial integral [0,\infty]
END_DOC
dr_radial_integral = 1.d0/dble(n_points_radial_grid-1)
integer :: i
do i = 1, n_points_radial_grid
grid_points_radial(i) = dble(i-1) * dr_radial_integral
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_angular,n_points_radial_grid,nucl_num)]
BEGIN_DOC
! x,y,z coordinates of grid points used for integration in 3d space
END_DOC
implicit none
integer :: i,j,k
double precision :: dr,x_ref,y_ref,z_ref
double precision :: knowles_function
do i = 1, nucl_num
x_ref = nucl_coord(i,1)
y_ref = nucl_coord(i,2)
z_ref = nucl_coord(i,3)
do j = 1, n_points_radial_grid-1
double precision :: x,r
! x value for the mapping of the [0, +\infty] to [0,1]
x = grid_points_radial(j)
! value of the radial coordinate for the integration
r = knowles_function(alpha_knowles(int(nucl_charge(i))),m_knowles,x)
! explicit values of the grid points centered around each atom
do k = 1, n_points_integration_angular
grid_points_per_atom(1,k,j,i) = &
x_ref + angular_quadrature_points(k,1) * r
grid_points_per_atom(2,k,j,i) = &
y_ref + angular_quadrature_points(k,2) * r
grid_points_per_atom(3,k,j,i) = &
z_ref + angular_quadrature_points(k,3) * r
enddo
enddo
enddo
END_PROVIDER
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)
!
! The "n" discrete variable represents the nucleis which in this array is
! represented by the last dimension and the points are labelled by the
! other dimensions.
END_DOC
implicit none
integer :: i,j,k,l,m
double precision :: r(3)
double precision :: accu,cell_function_becke
double precision :: tmp_array(nucl_num)
! run over all points in space
! that are referred to each atom
do j = 1, nucl_num
!for each radial grid attached to the "jth" atom
do k = 1, n_points_radial_grid -1
! for each angular point attached to the "jth" atom
do l = 1, n_points_integration_angular
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)
accu = 0.d0
! For each of these points in space, ou need to evaluate the P_n(r)
do i = 1, nucl_num
! function defined for each atom "i" by equation (13) and (21) with k == 3
tmp_array(i) = cell_function_becke(r,i) ! P_n(r)
! Then you compute the summ the P_n(r) function for each of the "r" points
accu += tmp_array(i)
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
END_PROVIDER
BEGIN_PROVIDER [double precision, final_weight_functions_at_grid_points, (n_points_integration_angular,n_points_radial_grid,nucl_num) ]
BEGIN_DOC
! Total weight on each grid point which takes into account all Lebedev, Voronoi and radial weights.
END_DOC
implicit none
integer :: i,j,k,l,m
double precision :: r(3)
double precision :: accu,cell_function_becke
double precision :: tmp_array(nucl_num)
double precision :: contrib_integration,x
double precision :: derivative_knowles_function,knowles_function
! run over all points in space
do j = 1, nucl_num ! that are referred to each atom
do i = 1, n_points_radial_grid -1 !for each radial grid attached to the "jth" atom
x = grid_points_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1]
do k = 1, n_points_integration_angular ! for each angular point attached to the "jth" atom
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
final_weight_functions_at_grid_points(k,i,j) = weights_angular_points(k) * weight_functions_at_grid_points(k,i,j) * contrib_integration * dr_radial_integral
enddo
enddo
enddo
END_PROVIDER

71
src/becke_numerical_grid/integration_radial.irp.f Normal file
View File

@ -0,0 +1,71 @@
double precision function knowles_function(alpha,m,x)
implicit none
BEGIN_DOC
! Function proposed by Knowles (JCP, 104, 1996) for distributing the radial points :
! the Log "m" function ( equation (7) in the paper )
END_DOC
double precision, intent(in) :: alpha,x
integer, intent(in) :: m
!print*, x
knowles_function = -alpha * dlog(1.d0-x**m)
end
double precision function derivative_knowles_function(alpha,m,x)
implicit none
BEGIN_DOC
! Derivative of the function proposed by Knowles (JCP, 104, 1996) for distributing the radial points
END_DOC
double precision, intent(in) :: alpha,x
integer, intent(in) :: m
double precision :: f
f = x**(m-1)
derivative_knowles_function = alpha * dble(m) * f / (1.d0 - x*f)
end
BEGIN_PROVIDER [double precision, alpha_knowles, (100)]
implicit none
integer :: i
BEGIN_DOC
! Recommended values for the alpha parameters according to the paper of Knowles (JCP, 104, 1996)
! as a function of the nuclear charge
END_DOC
! H-He
alpha_knowles(1) = 5.d0
alpha_knowles(2) = 5.d0
! Li-Be
alpha_knowles(3) = 7.d0
alpha_knowles(4) = 7.d0
! B-Ne
do i = 5, 10
alpha_knowles(i) = 5.d0
enddo
! Na-Mg
do i = 11, 12
alpha_knowles(i) = 7.d0
enddo
! Al-Ar
do i = 13, 18
alpha_knowles(i) = 5.d0
enddo
! K-Ca
do i = 19, 20
alpha_knowles(i) = 7.d0
enddo
! Sc-Zn
do i = 21, 30
alpha_knowles(i) = 5.d0
enddo
! Ga-Kr
do i = 31, 36
alpha_knowles(i) = 7.d0
enddo
END_PROVIDER

65
src/becke_numerical_grid/selected_points.irp.f Normal file
View File

@ -0,0 +1,65 @@
BEGIN_PROVIDER [integer, n_points_final_grid]
BEGIN_DOC
! Number of points which are non zero
END_DOC
integer :: i,j,k,l
double precision :: r(3)
n_points_final_grid = 0
do j = 1, nucl_num
do i = 1, n_points_radial_grid -1
do k = 1, n_points_integration_angular
r(1) = grid_points_per_atom(1,k,i,j)
r(2) = grid_points_per_atom(2,k,i,j)
r(3) = grid_points_per_atom(3,k,i,j)
! condition to maybe select points
if(.True.)then
n_points_final_grid += 1
endif
enddo
enddo
enddo
print*,'n_points_final_grid = ',n_points_final_grid
print*,'n max point = ',n_points_integration_angular*n_points_radial_grid*nucl_num
END_PROVIDER
BEGIN_PROVIDER [double precision, final_grid_points, (3,n_points_final_grid)]
&BEGIN_PROVIDER [double precision, final_weight_functions_at_final_grid_points, (n_points_final_grid) ]
&BEGIN_PROVIDER [integer, index_final_points, (3,n_points_final_grid) ]
&BEGIN_PROVIDER [integer, index_final_points_reverse, (n_points_integration_angular,n_points_radial_grid,nucl_num) ]
implicit none
BEGIN_DOC
! final_grid_points(1:3,j) = (/ x, y, z /) of the jth grid point
!
! final_weight_functions_at_final_grid_points(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions
!
! index_final_points(1:3,i) = gives the angular, radial and atomic indices associated to the ith grid point
!
! index_final_points_reverse(i,j,k) = index of the grid point having i as angular, j as radial and l as atomic indices
END_DOC
integer :: i,j,k,l,i_count
double precision :: r(3)
i_count = 0
do j = 1, nucl_num
do i = 1, n_points_radial_grid -1
do k = 1, n_points_integration_angular
r(1) = grid_points_per_atom(1,k,i,j)
r(2) = grid_points_per_atom(2,k,i,j)
r(3) = grid_points_per_atom(3,k,i,j)
! condition to maybe select points
if(.True.)then
i_count += 1
final_grid_points(1,i_count) = grid_points_per_atom(1,k,i,j)
final_grid_points(2,i_count) = grid_points_per_atom(2,k,i,j)
final_grid_points(3,i_count) = grid_points_per_atom(3,k,i,j)
final_weight_functions_at_final_grid_points(i_count) = final_weight_functions_at_grid_points(k,i,j)
index_final_points(1,i_count) = k
index_final_points(2,i_count) = i
index_final_points(3,i_count) = j
index_final_points_reverse(k,i,j) = i_count
endif
enddo
enddo
enddo
END_PROVIDER

50
src/becke_numerical_grid/step_function_becke.irp.f Normal file
View File

@ -0,0 +1,50 @@
double precision function step_function_becke(x)
implicit none
BEGIN_DOC
! Step function of the Becke paper (1988, JCP,88(4))
END_DOC
double precision, intent(in) :: x
double precision :: f_function_becke
integer :: i,n_max_becke
step_function_becke = f_function_becke(x)
do i = 1, 4
step_function_becke = f_function_becke(step_function_becke)
enddo
step_function_becke = 0.5d0*(1.d0 - step_function_becke)
end
double precision function f_function_becke(x)
implicit none
double precision, intent(in) :: x
f_function_becke = 1.5d0 * x - 0.5d0 * x*x*x
end
double precision function cell_function_becke(r,atom_number)
implicit none
double precision, intent(in) :: r(3)
integer, intent(in) :: atom_number
BEGIN_DOC
! atom_number :: atom on which the cell function of Becke (1988, JCP,88(4))
! r(1:3) :: x,y,z coordinantes of the current point
END_DOC
double precision :: mu_ij,nu_ij
double precision :: distance_i,distance_j,step_function_becke
integer :: j
distance_i = (r(1) - nucl_coord_transp(1,atom_number) ) * (r(1) - nucl_coord_transp(1,atom_number))
distance_i += (r(2) - nucl_coord_transp(2,atom_number) ) * (r(2) - nucl_coord_transp(2,atom_number))
distance_i += (r(3) - nucl_coord_transp(3,atom_number) ) * (r(3) - nucl_coord_transp(3,atom_number))
distance_i = dsqrt(distance_i)
cell_function_becke = 1.d0
do j = 1, nucl_num
if(j==atom_number)cycle
distance_j = (r(1) - nucl_coord_transp(1,j) ) * (r(1) - nucl_coord_transp(1,j))
distance_j+= (r(2) - nucl_coord_transp(2,j) ) * (r(2) - nucl_coord_transp(2,j))
distance_j+= (r(3) - nucl_coord_transp(3,j) ) * (r(3) - nucl_coord_transp(3,j))
distance_j = dsqrt(distance_j)
mu_ij = (distance_i - distance_j)*nucl_dist_inv(atom_number,j)
nu_ij = mu_ij + slater_bragg_type_inter_distance_ua(atom_number,j) * (1.d0 - mu_ij*mu_ij)
cell_function_becke *= step_function_becke(nu_ij)
enddo
end

3
src/dummy/NEED
View File

@ -1,5 +1,6 @@
ao_basis
ao_one_e_integrals
becke_numerical_grid
bitmask
cis
cisd
@ -20,6 +21,7 @@ mo_basis
mo_guess
mo_one_e_integrals
mpi
mrpt_utils
nuclei
perturbation
pseudo
@ -31,4 +33,3 @@ selectors_utils
single_ref_method
utils
zmq
mrpt_utils