qmckl/qmckl_ao.org at 145f4fba403f10b2c906ea1c6f5006762fef493e

TREX/qmckl

Fork 0

mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-11-19 20:42:50 +01:00

Aurélien Delval 145f4fba40

Add the qmckl_probes interface (#2 )

* comment

* Update distance test code

The distance test has been updated to the latest version, with a first attempt at using vfc_probes inside it

* Functional implementation of vfc_probes in the distance tests

This commit has the first functional vfc_ci tests. Verificarlo tests
should be written over the existing tests, and they can be enabled with
the following configure command:

QMCKL_DEVEL=1 ./configure --prefix=$PWD/_install --enable-maintainer-mode --enable-vfc_ci CC="verificarlo-f -Mpreprocess -D VFC_CI" FC="verificarlo-f -Mpreprocess -D VFC_CI" --host=x86_64

The --enable-vfc_ci flag will trigger the linking of the vfc_ci
library. Moreover, as of now, the "-Mpreprocess" and "-D VFC_CI" flags
have to be specified directly here. There is probably an appropriate
macro to place those flags into but I couldn't find it yet, and could
only manage to build the tests this way.

When the VFC_CI preprocessor is defined, somme additional code to
register and export the test probes can be executed (see
qmckl_distance.org).

As of now, the tests are built as normal, even though they are expected
to fail :

make all
make check

From there, the test_qmckl_distance (and potentially the others)
executable can be called at will. This will typically be done
automatically by vfc_ci, but one could manually execute the executable
by defining the following env variables :

VFC_PROBES_OUTPUT="test.csv" VFC_BACKENDS="libinterflop_ieee.so"

depending on the export file and the Verificarlo backend to be used.

The next steps will be to define more tests such as this one, and to
integrate them into a Verificarlo CI workflow (by writing a
vfc_tests_config.json file and using the automatic CI setup
command).

* Error in FOrtran interface fixed

* Added missing Fortran interfaces

* Modify distance test and install process integration

All probes are now ignored using only the preprocessor (instead
of checking for a facultative argument) in the distance test.
Moreover,preprocessing can now be enabled correctly using FCFLAGS
(the issue seemed to come from the order of the arguments passed
to gfortran/verificarlo-f with the preprocessor arg having to come
first).

* Add vfc_probes to AO tests

vfc_probes have been added to qmckl_ao.org in the same way as
qmckl_distance.org, which means that it can be enabled or disabled at
compile time using the --enable-vfc_ci option.

qmckl_distance.org has been slightly modified with a better indentation,
and configure.ac now adds the "-D VFC_CI" flag to CFLAGS when vfc_ci is
enabled.

* Start work on the vfc tests config file and on a probes wrapper

The goal in the next few commits is to make the integration of
vfc_probes even easier by using a wrapper to vfc_probe dedicated to
qmckl. This will make it easier to create a call to vfc_probe that can be
ignored if VFC_CI is not defined in the preprocessor. Once this is done,
the integration will be completed by trying to create an actual workflow
to automatically build the library and execute CI tests.

* Moved qmckl_probes out of src

As of now, qmckl_probes have been moved to tools, and can be built via a
bash script. This approach seems to make more sense, as this should not
be a part of the library itself, but an additional tool to test it.

* Functional Makefile setup to enable qmckl_probes

The current setup builds qmck_probes by adding it to the main QMckl
libray (by adding it to the libtool sources). The Fortran interface's
module still need to be compiled separately.

TODO : Clean the build setup, improve integration in qmckl_tests and
update tests in qmckl_ao with the new syntax.

* New probes syntax in AO tests

* Clean the probes/Makefile setup

The Fortran module is now built a the same time than the main library.
The commit also adds a few fixes in the tests and probes wrapper.

Co-authored-by: Anthony Scemama <scemama@irsamc.ups-tlse.fr>

2021-07-23 12:01:14 +02:00

96 KiB

Raw Blame History

Context
Radial part
Polynomial part
- General functions for Powers of $x-X_i$
- General functions for Value, Gradient and Laplacian of a polynomial
Combining radial and polynomial parts

qmckl_probe(#+TITLE: Atomic Orbitals

The atomic basis set is defined as a list of shells. Each shell $s$ is centered on a nucleus $A$, possesses a given angular momentum $l$ and a radial function $R_s$. The radial function is a linear combination of \emph{primitive} functions that can be of type Slater ($p=1$) or Gaussian ($p=2$):

\[ R_s(\mathbf{r}) = \mathcal{N}_s |\mathbf{r}-\mathbf{R}_A|^{n_s} \sum_{k=1}^{N_{\text{prim}}} a_{ks}\, f_{ks} \exp \left( - \gamma_{ks} | \mathbf{r}-\mathbf{R}_A | ^p \right). \]

In the case of Gaussian functions, $n_s$ is always zero. The normalization factor $\mathcal{N}_s$ ensures that all the functions of the shell are normalized to unity. Usually, basis sets are given a combination of normalized primitives, so the normalization coefficients of the primitives, $f_{ks}$, need also to be provided.

Atomic orbitals (AOs) are defined as

\[ \chi_i (\mathbf{r}) = \mathcal{M}_i\, P_{\eta(i)}(\mathbf{r})\, R_{\theta(i)} (\mathbf{r}) \]

where $\theta(i)$ returns the shell on which the AO is expanded, and $\eta(i)$ denotes which angular function is chosen. Here, the parameter $\mathcal{M}_i$ is an extra parameter which allows the normalization of the different functions of the same shell to be different, as in GAMESS for example.

In this section we describe first how the basis set is stored in the context, and then we present the kernels used to compute the values, gradients and Laplacian of the atomic basis functions.

Context

The following arrays are stored in the context:

`type`		Gaussian (`'G'`) or Slater (`'S'`)
`shell_num`		Number of shells
`prim_num`		Total number of primitives
`nucleus_index`	`[nucl_num]`	Index of the first shell of each nucleus
`shell_ang_mom`	`[shell_num]`	Angular momentum of each shell
`shell_prim_num`	`[shell_num]`	Number of primitives in each shell
`shell_prim_index`	`[shell_num]`	Address of the first primitive of each shell in the `EXPONENT` array
`shell_factor`	`[shell_num]`	Normalization factor for each shell
`exponent`	`[prim_num]`	Array of exponents
`coefficient`	`[prim_num]`	Array of coefficients
`prim_factor`	`[prim_num]`	Normalization factors of the primtives

Computed data:

`nucleus_prim_index`	`[nucl_num]`	Index of the first primitive for each nucleus
`primitive_vgl`	`[prim_num][5][walk_num][elec_num]`	Value, gradients, Laplacian of the primitives at electron positions
`primitive_vgl_date`	`uint64_t`	Late modification date of Value, gradients, Laplacian of the primitives at electron positions
`shell_vgl`	`[prim_num][5][walk_num][elec_num]`	Value, gradients, Laplacian of the primitives at electron positions
`shell_vgl_date`	`uint64_t`	Late modification date of Value, gradients, Laplacian of the shells at electron positions
`nucl_shell_index`	`[nucl_num]`	Index of the first shell for each nucleus
`exponent_sorted`	`[prim_num]`	Array of exponents for sorted primitives
`coeff_norm_sorted`	`[prim_num]`	Array of normalized coefficients for sorted primitives
`prim_factor_sorted`	`[prim_num]`	Normalization factors of the sorted primtives
`nuclear_radius`	`[nucl_num]`	Distance beyond which all the AOs are zero

For H_2 with the following basis set,

HYDROGEN
S   5
1         3.387000E+01           6.068000E-03
2         5.095000E+00           4.530800E-02
3         1.159000E+00           2.028220E-01
4         3.258000E-01           5.039030E-01
5         1.027000E-01           3.834210E-01
S   1
1         3.258000E-01           1.000000E+00
S   1
1         1.027000E-01           1.000000E+00
P   1
1         1.407000E+00           1.000000E+00
P   1
1         3.880000E-01           1.000000E+00
D   1
1         1.057000E+00           1.0000000

we have:

type = 'G'
shell_num = 12
prim_num = 20
nucleus_index = [0 , 6]
shell_ang_mom = [0, 0, 0, 1, 1, 2, 0, 0, 0, 1, 1, 2]
shell_factor  = [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.]
shell_prim_num = [5, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1]
shell_prim_index = [0 , 5 , 6 , 7 , 8 , 9 , 10, 15, 16, 17, 18, 19]
exponent = [ 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027, 1.407,
0.388, 1.057, 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027, 1.407,
0.388, 1.057]
coefficient = [ 0.006068, 0.045308, 0.202822, 0.503903, 0.383421, 1.0, 1.0,
1.0, 1.0, 1.0, 0.006068, 0.045308, 0.202822, 0.503903, 0.383421, 1.0, 1.0,
1.0, 1.0, 1.0]
prim_factor = [ 1.0006253235944540e+01, 2.4169531573445120e+00, 7.9610924849766440e-01
3.0734305383061117e-01, 1.2929684417481876e-01, 3.0734305383061117e-01,
1.2929684417481876e-01, 2.1842769845268308e+00, 4.3649547399719840e-01,
1.8135965626177861e+00, 1.0006253235944540e+01, 2.4169531573445120e+00,
7.9610924849766440e-01, 3.0734305383061117e-01, 1.2929684417481876e-01,
3.0734305383061117e-01, 1.2929684417481876e-01, 2.1842769845268308e+00,
4.3649547399719840e-01, 1.8135965626177861e+00 ]

Data structure

typedef struct qmckl_ao_basis_struct {
int32_t   uninitialized;
int64_t   shell_num;
int64_t   prim_num;
int64_t * nucleus_index;
int64_t * nucleus_shell_num;
int32_t * shell_ang_mom;
int64_t * shell_prim_num;
int64_t * nucleus_prim_index;
int64_t * shell_prim_index;
double  * shell_factor;
double  * exponent    ;
double  * coefficient ;
double  * prim_factor ;
double  * primitive_vgl;
int64_t   primitive_vgl_date;
double  * shell_vgl;
int64_t   shell_vgl_date;
bool      provided;
char      type;
} qmckl_ao_basis_struct;

The uninitialized integer contains one bit set to one for each initialization function which has not been called. It becomes equal to zero after all initialization functions have been called. The struct is then initialized and provided == true. Some values are initialized by default, and are not concerned by this mechanism.

qmckl_exit_code qmckl_init_ao_basis(qmckl_context context);

qmckl_exit_code qmckl_init_ao_basis(qmckl_context context) {

if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
 return false;
}

qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);

ctx->ao_basis.uninitialized = (1 << 12) - 1;

/* Default values */
/* ctx->ao_basis.
  ,*/

return QMCKL_SUCCESS;
}

Access functions

When all the data for the AOs have been provided, the following function returns true.

bool      qmckl_ao_basis_provided           (const qmckl_context context);

#+NAME:post

Initialization functions

To set the basis set, all the following functions need to be called.

qmckl_exit_code  qmckl_set_ao_basis_type             (qmckl_context context, const char      t);
qmckl_exit_code  qmckl_set_ao_basis_shell_num        (qmckl_context context, const int64_t   shell_num);
qmckl_exit_code  qmckl_set_ao_basis_prim_num         (qmckl_context context, const int64_t   prim_num);
qmckl_exit_code  qmckl_set_ao_basis_nucleus_index    (qmckl_context context, const int64_t * nucleus_index);
qmckl_exit_code  qmckl_set_ao_basis_nucleus_shell_num(qmckl_context context, const int64_t * nucleus_shell_num);
qmckl_exit_code  qmckl_set_ao_basis_shell_ang_mom    (qmckl_context context, const int32_t * shell_ang_mom);
qmckl_exit_code  qmckl_set_ao_basis_shell_prim_num   (qmckl_context context, const int64_t * shell_prim_num);
qmckl_exit_code  qmckl_set_ao_basis_shell_prim_index (qmckl_context context, const int64_t * shell_prim_index);
qmckl_exit_code  qmckl_set_ao_basis_shell_factor     (qmckl_context context, const double  * shell_factor);
qmckl_exit_code  qmckl_set_ao_basis_exponent         (qmckl_context context, const double  * exponent);
qmckl_exit_code  qmckl_set_ao_basis_coefficient      (qmckl_context context, const double  * coefficient);
qmckl_exit_code  qmckl_set_ao_basis_prim_factor      (qmckl_context context, const double  * prim_factor);

#+NAME:pre2

#+NAME:post2

When the basis set is completely entered, other data structures are computed to accelerate the calculations. The primitives within each contraction are sorted in ascending order of their exponents, such that as soon as a primitive is zero all the following functions vanish. Also, it is possible to compute a nuclear radius beyond which all the primitives are zero up to the numerical accuracy defined in the context.

Fortran interfaces

Radial part

General functions for Gaussian basis functions

qmckl_ao_gaussian_vgl computes the values, gradients and Laplacians at a given point of n Gaussian functions centered at the same point:

\[ v_i = \exp(-a_i |X-R|^2) \] \[ \nabla_x v_i = -2 a_i (X_x - R_x) v_i \] \[ \nabla_y v_i = -2 a_i (X_y - R_y) v_i \] \[ \nabla_z v_i = -2 a_i (X_z - R_z) v_i \] \[ \Delta v_i = a_i (4 |X-R|^2 a_i - 6) v_i \]

`context`	input	Global state
`X(3)`	input	Array containing the coordinates of the points
`R(3)`	input	Array containing the x,y,z coordinates of the center
`n`	input	Number of computed Gaussians
`A(n)`	input	Exponents of the Gaussians
`VGL(ldv,5)`	output	Value, gradients and Laplacian of the Gaussians
`ldv`	input	Leading dimension of array `VGL`

Requirements

context is not 0
n > 0
ldv >= 5
A(i) > 0 for all i
X is allocated with at least $3 \times 8$ bytes
R is allocated with at least $3 \times 8$ bytes
A is allocated with at least $n \times 8$ bytes
VGL is allocated with at least $n \times 5 \times 8$ bytes

qmckl_exit_code
qmckl_ao_gaussian_vgl(const qmckl_context context,
                 const double *X,
                 const double *R,
                 const int64_t *n,
                 const int64_t *A,
                 const double *VGL,
                 const int64_t ldv);

integer function qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv) result(info)
use qmckl
implicit none
integer*8 , intent(in)  :: context
real*8    , intent(in)  :: X(3), R(3)
integer*8 , intent(in)  :: n
real*8    , intent(in)  :: A(n)
real*8    , intent(out) :: VGL(ldv,5)
integer*8 , intent(in)  :: ldv

integer*8         :: i,j
real*8            :: Y(3), r2, t, u, v

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif

if (n <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif

if (ldv < n) then
info = QMCKL_INVALID_ARG_7
return
endif


do i=1,3
Y(i) = X(i) - R(i)
end do
r2 = Y(1)*Y(1) + Y(2)*Y(2) + Y(3)*Y(3)

do i=1,n
VGL(i,1) = dexp(-A(i) * r2)
end do

do i=1,n
VGL(i,5) = A(i) * VGL(i,1)
end do

t = -2.d0 * ( X(1) - R(1) )
u = -2.d0 * ( X(2) - R(2) )
v = -2.d0 * ( X(3) - R(3) )

do i=1,n
VGL(i,2) = t * VGL(i,5)
VGL(i,3) = u * VGL(i,5)
VGL(i,4) = v * VGL(i,5)
end do

t = 4.d0 * r2
do i=1,n
VGL(i,5) = (t * A(i) - 6.d0) *  VGL(i,5)
end do

end function qmckl_ao_gaussian_vgl_f

#ifdef VFC_CI
integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context, probes) bind(C)
#else
integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C)
#endif
use qmckl
use qmckl_probes_f

implicit none

integer(c_int64_t), intent(in), value :: context
logical(C_BOOL) :: vfc_err

#ifdef VFC_CI
type(vfc_probes) :: probes
integer(C_INT) :: vfc_err
#endif

integer*8                     :: n, ldv, j, i
double precision              :: X(3), R(3), Y(3), r2
double precision, allocatable :: VGL(:,:), A(:)
double precision              :: epsilon

epsilon = qmckl_get_numprec_epsilon(context)

X = (/ 1.1 , 2.2 ,  3.3 /)
R = (/ 0.1 , 1.2 , -2.3 /)
Y(:) = X(:) - R(:)
r2 = Y(1)**2 + Y(2)**2 + Y(3)**2

n = 10;
ldv = 100;

allocate (A(n), VGL(ldv,5))
do i=1,n
  A(i) = 0.0013 * dble(ishft(1,i))
end do


test_qmckl_ao_gaussian_vgl = &
    qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv)

vfc_err = qmckl_probe("ao"//C_NULL_CHAR, "gaussian_vgl"//C_NULL_CHAR, &
 DBLE(test_qmckl_ao_gaussian_vgl))

if (test_qmckl_ao_gaussian_vgl /= 0) return
#endif

test_qmckl_ao_gaussian_vgl = -1

#ifndef VFC_CI
do i=1,n
  test_qmckl_ao_gaussian_vgl = -11
  if (dabs(1.d0 - VGL(i,1) / (&
       dexp(-A(i) * r2) &
       )) > epsilon ) return

  test_qmckl_ao_gaussian_vgl = -12
  if (dabs(1.d0 - VGL(i,2) / (&
       -2.d0 * A(i) * Y(1) * dexp(-A(i) * r2) &
       )) > epsilon ) return

  test_qmckl_ao_gaussian_vgl = -13
  if (dabs(1.d0 - VGL(i,3) / (&
       -2.d0 * A(i) * Y(2) * dexp(-A(i) * r2) &
       )) > epsilon ) return

  test_qmckl_ao_gaussian_vgl = -14
  if (dabs(1.d0 - VGL(i,4) / (&
       -2.d0 * A(i) * Y(3) * dexp(-A(i) * r2) &
       )) > epsilon ) return

  test_qmckl_ao_gaussian_vgl = -15
  if (dabs(1.d0 - VGL(i,5) / (&
       A(i) * (4.d0*r2*A(i) - 6.d0) * dexp(-A(i) * r2) &
       )) > epsilon ) return
end do
#endif

test_qmckl_ao_gaussian_vgl = 0

deallocate(VGL)
end function test_qmckl_ao_gaussian_vgl

TODO General functions for Slater basis functions

TODO General functions for Radial functions on a grid

DONE Computation of primitives

Get

qmckl_exit_code qmckl_get_ao_basis_primitive_vgl(qmckl_context context, double* const primitive_vgl);

Provide

Compute

qmckl_context	context	in	Global state
int64_t	prim_num	in	Number of primitives
int64_t	elec_num	in	Number of electrons
int64_t	nucl_num	in	Number of nuclei
int64_t	walk_num	in	Number of walkers
int64_t	nucleus_prim_index[nucl_num]	in	Index of the 1st primitive of each nucleus
double	elec_coord[walk_num][3][elec_num]	in	Electron coordinates
double	nucl_coord[3][elec_num]	in	Nuclear coordinates
double	expo[prim_num]	in	Exponents of the primitives
double	primitive_vgl[prim_num][5][walk_num][elec_num]	out	Value, gradients and Laplacian of the primitives

integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, &
 prim_num, elec_num, nucl_num, walk_num, &
 nucleus_prim_index, elec_coord, nucl_coord, expo, primitive_vgl) &
 result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: prim_num
integer*8             , intent(in)  :: nucl_num
integer*8             , intent(in)  :: elec_num
integer*8             , intent(in)  :: walk_num
integer*8             , intent(in)  :: nucleus_prim_index(nucl_num+1)
double precision      , intent(in)  :: elec_coord(elec_num,3,walk_num)
double precision      , intent(in)  :: nucl_coord(nucl_num,3)
double precision      , intent(in)  :: expo(prim_num)
double precision      , intent(out) :: primitive_vgl(elec_num,walk_num,5,prim_num)

integer*8 :: inucl, iprim, iwalk, ielec
double precision :: x, y, z, two_a, ar2, r2, v, cutoff

info = QMCKL_SUCCESS

! Don't compute exponentials when the result will be almost zero.
cutoff = -dlog(1.d-15)

do inucl=1,nucl_num
 ! C is zero-based, so shift bounds by one
 do iprim = nucleus_prim_index(inucl)+1, nucleus_prim_index(inucl+1)
    do iwalk = 1, walk_num
       do ielec = 1, elec_num
          x = elec_coord(ielec,1,iwalk) - nucl_coord(inucl,1)
          y = elec_coord(ielec,2,iwalk) - nucl_coord(inucl,2)
          z = elec_coord(ielec,3,iwalk) - nucl_coord(inucl,3)

          r2 = x*x + y*y + z*z
          ar2 = expo(iprim)*r2
          if (ar2 > cutoff) cycle

          v = dexp(-ar2)
          two_a = -2.d0 * expo(iprim) * v

          primitive_vgl(ielec, iwalk, 1, iprim) = v
          primitive_vgl(ielec, iwalk, 2, iprim) = two_a * x
          primitive_vgl(ielec, iwalk, 3, iprim) = two_a * y
          primitive_vgl(ielec, iwalk, 4, iprim) = two_a * z
          primitive_vgl(ielec, iwalk, 5, iprim) = two_a * (3.d0 - 2.d0*ar2)

       end do
    end do
 end do
end do

end function qmckl_compute_ao_basis_primitive_gaussian_vgl_f

Test

Ideas for improvement

// m : walkers
// j : electrons
// l : primitives

k=0;
for (m=0 ; m<walk_num ; ++m) {
for (j=0 ; j<elec_num ; ++j) {
for (i=0 ; i<nucl_num ; ++i) {

  r2 = nucl_elec_dist[i][j];

  if (r2 < nucl_radius2[i]) {

    for (l=0 ; l<prim_num ; ++l) {
      tmp[k].i   = i;
      tmp[k].j   = j;
      tmp[k].m   = m;
      tmp[k].ar2 = -expo[l] *r2;
      ++k;
    }
  }
}
}
}

// sort(tmp) in increasing ar2;
// Identify first ar2 above numerical accuracy threshold
// Compute vectorized exponentials on significant values

Computation of shells

Get

qmckl_exit_code qmckl_get_ao_basis_shell_vgl(qmckl_context context, double* const shell_vgl);

Provide

Compute

`qmckl_context`	`context`	in	Global state
`int64_t`	`prim_num`	in	Number of primitives
`int64_t`	`shell_num`	in	Number of shells
`int64_t`	`elec_num`	in	Number of electrons
`int64_t`	`nucl_num`	in	Number of nuclei
`int64_t`	`walk_num`	in	Number of walkers
`int64_t`	`nucleus_shell_num[nucl_num]`	in	Number of shells for each nucleus
`int64_t`	`nucleus_index[nucl_num]`	in	Index of the 1st shell of each nucleus
`int64_t`	`shell_prim_index[shell_num]`	in	Index of the 1st primitive of each shell
`int64_t`	`shell_prim_num[shell_num]`	in	Number of primitives per shell
`double`	`elec_coord[walk_num][3][elec_num]`	in	Electron coordinates
`double`	`nucl_coord[3][elec_num]`	in	Nuclear coordinates
`double`	`expo[prim_num]`	in	Exponents of the primitives
`double`	`coef[prim_num]`	in	Coefficients of the primitives
`double`	`shell_vgl[shell_num][5][walk_num][elec_num]`	out	Value, gradients and Laplacian of the shells

integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, &
 prim_num, shell_num, elec_num, nucl_num, walk_num, &
 nucleus_shell_num, nucleus_index, shell_prim_index, shell_prim_num, &
 elec_coord, nucl_coord, expo, coef, shell_vgl) &
 result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: prim_num
integer*8             , intent(in)  :: shell_num
integer*8             , intent(in)  :: nucl_num
integer*8             , intent(in)  :: elec_num
integer*8             , intent(in)  :: walk_num
integer*8             , intent(in)  :: nucleus_shell_num(nucl_num)
integer*8             , intent(in)  :: nucleus_index(nucl_num)
integer*8             , intent(in)  :: shell_prim_index(shell_num)
integer*8             , intent(in)  :: shell_prim_num(shell_num)
double precision      , intent(in)  :: elec_coord(elec_num,3,walk_num)
double precision      , intent(in)  :: nucl_coord(nucl_num,3)
double precision      , intent(in)  :: expo(prim_num)
double precision      , intent(in)  :: coef(prim_num)
double precision      , intent(out) :: shell_vgl(elec_num,walk_num,5,shell_num)

integer*8 :: inucl, iprim, iwalk, ielec, ishell
double precision :: x, y, z, two_a, ar2, r2, v, cutoff

info = QMCKL_SUCCESS

! Don't compute exponentials when the result will be almost zero.
! TODO : Use numerical precision here
cutoff = -dlog(1.d-15)

do inucl=1,nucl_num
 do ishell=nucleus_index(inucl)+1, nucleus_index(inucl)+nucleus_shell_num(inucl)
    ! C is zero-based, so shift bounds by one

    do iwalk = 1, walk_num
       do ielec = 1, elec_num

          shell_vgl(ielec, iwalk, 1:5, ishell) = 0.d0

          x = elec_coord(ielec,1,iwalk) - nucl_coord(inucl,1)
          y = elec_coord(ielec,2,iwalk) - nucl_coord(inucl,2)
          z = elec_coord(ielec,3,iwalk) - nucl_coord(inucl,3)

          r2 = x*x + y*y + z*z

          do iprim = shell_prim_index(ishell)+1, shell_prim_index(ishell)+shell_prim_num(ishell)

             ar2 = expo(iprim)*r2
             if (ar2 > cutoff) then
                cycle
             end if

             v = coef(iprim) * dexp(-ar2)
             two_a = -2.d0 * expo(iprim) * v

             shell_vgl(ielec, iwalk, 1, ishell) = &
                  shell_vgl(ielec, iwalk, 1, ishell) + v

             shell_vgl(ielec, iwalk, 2, ishell) = &
                  shell_vgl(ielec, iwalk, 2, ishell) + two_a * x

             shell_vgl(ielec, iwalk, 3, ishell) = &
                  shell_vgl(ielec, iwalk, 3, ishell) + two_a * y

             shell_vgl(ielec, iwalk, 4, ishell) = &
                  shell_vgl(ielec, iwalk, 4, ishell) + two_a * z

             shell_vgl(ielec, iwalk, 5, ishell) = &
                  shell_vgl(ielec, iwalk, 5, ishell) + two_a * (3.d0 - 2.d0*ar2)

          end do

       end do
    end do

 end do
end do

end function qmckl_compute_ao_basis_shell_gaussian_vgl_f

Test

Polynomial part

General functions for Powers of $x-X_i$

The qmckl_ao_power function computes all the powers of the n input data up to the given maximum value given in input for each of the $n$ points:

\[ P_{ik} = X_i^k \]

qmckl_context	context	in	Global state
int64_t	n	in	Number of values
double	X[n]	in	Array containing the input values
int32_t	LMAX[n]	in	Array containing the maximum power for each value
double	P[n][ldp]	out	Array containing all the powers of `X`
int64_t	ldp	in	Leading dimension of array `P`

Requirements

context is not QMCKL_NULL_CONTEXT
n > 0
X is allocated with at least $n \times 8$ bytes
LMAX is allocated with at least $n \times 4$ bytes
P is allocated with at least $n \times \max_i \text{LMAX}_i \times 8$ bytes
LDP >= $\max_i$ LMAX[i]

C Header

qmckl_exit_code qmckl_ao_power (
      const qmckl_context context,
      const int64_t n,
      const double* X,
      const int32_t* LMAX,
      double* const P,
      const int64_t ldp );

Source

integer function qmckl_ao_power_f(context, n, X, LMAX, P, ldp) result(info)
use qmckl
implicit none
integer*8 , intent(in)  :: context
integer*8 , intent(in)  :: n
real*8    , intent(in)  :: X(n)
integer   , intent(in)  :: LMAX(n)
real*8    , intent(out) :: P(ldp,n)
integer*8 , intent(in)  :: ldp

integer*8  :: i,k

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
 info = QMCKL_INVALID_CONTEXT
 return
endif

if (n <= ldp) then
 info = QMCKL_INVALID_ARG_2
 return
endif

k = MAXVAL(LMAX)
if (LDP < k) then
 info = QMCKL_INVALID_ARG_6
 return
endif

if (k <= 0) then
 info = QMCKL_INVALID_ARG_4
 return
endif

do i=1,n
 P(1,i) = X(i)
 do k=2,LMAX(i)
    P(k,i) = P(k-1,i) * X(i)
 end do
end do

end function qmckl_ao_power_f

C interface

Fortran interface

Test

#ifdef VFC_CI
integer(c_int32_t) function test_qmckl_ao_power(context, probes) bind(C)
#else
integer(c_int32_t) function test_qmckl_ao_power(context) bind(C)

use qmckl
use qmckl_probes_f

implicit none

logical(C_BOOL) :: vfc_err


integer(qmckl_context), intent(in), value :: context

integer*8                     :: n, LDP
integer, allocatable          :: LMAX(:)
double precision, allocatable :: X(:), P(:,:)
integer*8                     :: i,j
double precision              :: epsilon

epsilon = qmckl_get_numprec_epsilon(context)

n = 100;
LDP = 10;

allocate(X(n), P(LDP,n), LMAX(n))

do j=1,n
 X(j) = -5.d0 + 0.1d0 * dble(j)
 LMAX(j) = 1 + int(mod(j, 5),4)
end do

test_qmckl_ao_power = qmckl_ao_power(context, n, X, LMAX, P, LDP)

vfc_err = qmckl_probe("ao"//C_NULL_CHAR, "power"//C_NULL_CHAR, &
DBLE(test_qmckl_ao_power))

if (test_qmckl_ao_power /= QMCKL_SUCCESS) return
#endif

test_qmckl_ao_power = QMCKL_FAILURE

do j=1,n
 do i=1,LMAX(j)
    if ( X(j)**i == 0.d0 ) then
       if ( P(i,j) /= 0.d0) return
    else
       if ( dabs(1.d0 - P(i,j) / (X(j)**i)) > epsilon ) return
    end if
 end do
end do

test_qmckl_ao_power = QMCKL_SUCCESS
deallocate(X,P,LMAX)
end function test_qmckl_ao_power

General functions for Value, Gradient and Laplacian of a polynomial

A polynomial is centered on a nucleus $\mathbf{R}_i$

\[ P_l(\mathbf{r},\mathbf{R}_i) = (x-X_i)^a (y-Y_i)^b (z-Z_i)^c \]

The gradients with respect to electron coordinates are

\begin{eqnarray*} \frac{\partial }{\partial x} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & a (x-X_i)^{a-1} (y-Y_i)^b (z-Z_i)^c \\ \frac{\partial }{\partial y} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & b (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c \\ \frac{\partial }{\partial z} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & c (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} \\ \end{eqnarray*}

and the Laplacian is

\begin{eqnarray*} \left( \frac{\partial }{\partial x^2} + \frac{\partial }{\partial y^2} + \frac{\partial }{\partial z^2} \right) P_l \left(\mathbf{r},\mathbf{R}_i \right) & = & a(a-1) (x-X_i)^{a-2} (y-Y_i)^b (z-Z_i)^c + \\ && b(b-1) (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c + \\ && c(c-1) (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1}. \end{eqnarray*}

qmckl_ao_polynomial_vgl computes the values, gradients and Laplacians at a given point in space, of all polynomials with an angular momentum up to lmax.

qmckl_context	context	in	Global state
double	X[3]	in	Array containing the coordinates of the points
double	R[3]	in	Array containing the x,y,z coordinates of the center
int32_t	lmax	in	Maximum angular momentum
int64_t	n	inout	Number of computed polynomials
int32_t	L[n][ldl]	out	Contains a,b,c for all `n` results
int64_t	ldl	in	Leading dimension of `L`
double	VGL[n][ldv]	out	Value, gradients and Laplacian of the polynomials
int64_t	ldv	in	Leading dimension of array `VGL`

Requirements

context is not QMCKL_NULL_CONTEXT
n > 0
lmax >= 0
ldl >= 3
ldv >= 5
X is allocated with at least $3 \times 8$ bytes
R is allocated with at least $3 \times 8$ bytes
n >= (lmax+1)(lmax+2)(lmax+3)/6
L is allocated with at least $3 \times n \times 4$ bytes
VGL is allocated with at least $5 \times n \times 8$ bytes
On output, n should be equal to (lmax+1)(lmax+2)(lmax+3)/6
On output, the powers are given in the following order (l=a+b+c):
- Increasing values of l
- Within a given value of l, alphabetical order of the string made by a*"x" + b*"y" + c*"z" (in Python notation). For example, with a=0, b=2 and c=1 the string is "yyz"

C Header

qmckl_exit_code qmckl_ao_polynomial_vgl (
      const qmckl_context context,
      const double* X,
      const double* R,
      const int32_t lmax,
      int64_t* n,
      int32_t* const L,
      const int64_t ldl,
      double* const VGL,
      const int64_t ldv );

Source

integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, ldv) result(info)
use qmckl
implicit none
integer*8 , intent(in)  :: context
real*8    , intent(in)  :: X(3), R(3)
integer   , intent(in)  :: lmax
integer*8 , intent(out) :: n
integer   , intent(out) :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6)
integer*8 , intent(in)  :: ldl
real*8    , intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6)
integer*8 , intent(in)  :: ldv

integer*8         :: i,j
integer           :: a,b,c,d
real*8            :: Y(3)
integer           :: lmax_array(3)
real*8            :: pows(-2:lmax,3)
integer, external :: qmckl_ao_power_f
double precision  :: xy, yz, xz
double precision  :: da, db, dc, dd

info = 0

if (context == QMCKL_NULL_CONTEXT) then
 info = QMCKL_INVALID_CONTEXT
 return
endif

if (lmax < 0) then
 info = QMCKL_INVALID_ARG_4
 return
endif

if (ldl < 3) then
 info = QMCKL_INVALID_ARG_7
 return
endif

if (ldv < 5) then
 info = QMCKL_INVALID_ARG_9
 return
endif


do i=1,3
 Y(i) = X(i) - R(i)
end do

lmax_array(1:3) = lmax
if (lmax == 0) then
 VGL(1,1) = 1.d0
 vgL(2:5,1) = 0.d0
 l(1:3,1) = 0
 n=1
else if (lmax > 0) then
 pows(-2:0,1:3) = 1.d0
 do i=1,lmax
    pows(i,1) = pows(i-1,1) * Y(1)
    pows(i,2) = pows(i-1,2) * Y(2)
    pows(i,3) = pows(i-1,3) * Y(3)
 end do

 VGL(1:5,1:4) = 0.d0
 l  (1:3,1:4) = 0

 VGL(1  ,1  ) = 1.d0
 vgl(1:5,2:4) = 0.d0

 l  (1,2) = 1
 vgl(1,2) = pows(1,1)
 vgL(2,2) = 1.d0

 l  (2,3) = 1
 vgl(1,3) = pows(1,2)
 vgL(3,3) = 1.d0

 l  (3,4) = 1
 vgl(1,4) = pows(1,3)
 vgL(4,4) = 1.d0

 n=4
endif

! l>=2
dd = 2.d0
do d=2,lmax
 da = dd
 do a=d,0,-1
    db = dd-da
    do b=d-a,0,-1
       c  = d  - a  - b
       dc = dd - da - db
       n = n+1

       l(1,n) = a
       l(2,n) = b
       l(3,n) = c

       xy = pows(a,1) * pows(b,2)
       yz = pows(b,2) * pows(c,3)
       xz = pows(a,1) * pows(c,3)

       vgl(1,n) = xy * pows(c,3)

       xy = dc * xy
       xz = db * xz
       yz = da * yz

       vgl(2,n) = pows(a-1,1) * yz
       vgl(3,n) = pows(b-1,2) * xz
       vgl(4,n) = pows(c-1,3) * xy

       vgl(5,n) = &
            (da-1.d0) * pows(a-2,1) * yz + &
            (db-1.d0) * pows(b-2,2) * xz + &
            (dc-1.d0) * pows(c-2,3) * xy

       db = db - 1.d0
    end do
    da = da - 1.d0
 end do
 dd = dd + 1.d0
end do

info = QMCKL_SUCCESS

end function qmckl_ao_polynomial_vgl_f

C interface

Fortran interface

Test

#ifdef VFC_CI
integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context, probes) bind(C)
#else
integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
#endif
use qmckl
use qmckl_probes_f

implicit none

integer(c_int64_t), intent(in), value :: context
logical(C_BOOL) :: vfc_err

#ifdef VFC_CI
type(vfc_probes) :: probes
integer(C_INT) :: vfc_err
#endif

integer                       :: lmax, d, i
integer, allocatable          :: L(:,:)
integer*8                     :: n, ldl, ldv, j
double precision              :: X(3), R(3), Y(3)
double precision, allocatable :: VGL(:,:)
double precision              :: w
double precision              :: epsilon

epsilon = qmckl_get_numprec_epsilon(context)

X = (/ 1.1 , 2.2 ,  3.3 /)
R = (/ 0.1 , 1.2 , -2.3 /)
Y(:) = X(:) - R(:)

lmax = 4;
ldl = 3;
ldv = 100;

d = (lmax+1)*(lmax+2)*(lmax+3)/6

allocate (L(ldl,d), VGL(ldv,d))

test_qmckl_ao_polynomial_vgl = &
   qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv)

vfc_err = qmckl_probe("ao"//C_NULL_CHAR, "polynomial_vgl"//C_NULL_CHAR, &
DBLE(test_qmckl_ao_polynomial_vgl))

if (test_qmckl_ao_polynomial_vgl /= QMCKL_SUCCESS) return
if (n /= d) return
#endif

#ifdef VFC_CI
do j=1,n
 test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
 do i=1,3
    if (L(i,j) < 0) return
 end do
 test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
 if (dabs(1.d0 - VGL(1,j) / (&
      Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**L(3,j)  &
      )) > epsilon ) return

 test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
 if (L(1,j) < 1) then
    if (VGL(2,j) /= 0.d0) return
 else
    if (dabs(1.d0 - VGL(2,j) / (&
         L(1,j) * Y(1)**(L(1,j)-1) * Y(2)**L(2,j) * Y(3)**L(3,j) &
         )) > epsilon ) return
 end if

 test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
 if (L(2,j) < 1) then
    if (VGL(3,j) /= 0.d0) return
 else
    if (dabs(1.d0 - VGL(3,j) / (&
         L(2,j) * Y(1)**L(1,j) * Y(2)**(L(2,j)-1) * Y(3)**L(3,j) &
         )) > epsilon ) return
 end if

 test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
 if (L(3,j) < 1) then
    if (VGL(4,j) /= 0.d0) return
 else
    if (dabs(1.d0 - VGL(4,j) / (&
         L(3,j) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-1) &
         )) > epsilon ) return
 end if

 test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
 w = 0.d0
 if (L(1,j) > 1) then
    w = w + L(1,j) * (L(1,j)-1) * Y(1)**(L(1,j)-2) * Y(2)**L(2,j) * Y(3)**L(3,j)
 end if
 if (L(2,j) > 1) then
    w = w + L(2,j) * (L(2,j)-1) * Y(1)**L(1,j) * Y(2)**(L(2,j)-2) * Y(3)**L(3,j)
 end if
 if (L(3,j) > 1) then
    w = w + L(3,j) * (L(3,j)-1) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-2)
 end if
 if (dabs(1.d0 - VGL(5,j) / w) > epsilon ) return
end do
#endif

test_qmckl_ao_polynomial_vgl = QMCKL_SUCCESS

deallocate(L,VGL)
end function test_qmckl_ao_polynomial_vgl

#ifdef VFC_CI
int  test_qmckl_ao_polynomial_vgl(qmckl_context context, vfc_probes * probes);
assert(0 == test_qmckl_ao_polynomial_vgl(context, &probes));
#else
int  test_qmckl_ao_polynomial_vgl(qmckl_context context);
assert(0 == test_qmckl_ao_polynomial_vgl(context));
#endif

96 KiB Raw Blame History

Context

Data structure

Access functions

Initialization functions

Fortran interfaces

Radial part

General functions for Gaussian basis functions

TODO General functions for Slater basis functions

TODO General functions for Radial functions on a grid

DONE Computation of primitives

Get

Provide

Compute

Test

Ideas for improvement

Computation of shells

Get

Provide

Compute

Test

Polynomial part

General functions for Powers of $x-X_i$

Requirements

C Header

Source

C interface

Fortran interface

Test

General functions for Value, Gradient and Laplacian of a polynomial

Requirements

C Header

Source

C interface

Fortran interface

Test

Combining radial and polynomial parts

96 KiB

Raw Blame History