mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-08-14 16:58:38 +02:00
78fdbfca77
The previous value of epsilon, which was pulled form the qmckl_context and reused fot the qmckl probes, was arout 4.10^-16. It's a good value to use for the tests with IEEE arithmetic, however, MCA introduces some noise which is at this order of magnitude, and caused one probe in particular to fail, while its accuracy was actually satisfying. This value is now multiplied by 16 if VFC_CI is defined (to ask for 4 less significant digits in base 2), which makes more sense for MCA.
3835 lines
130 KiB
Org Mode
3835 lines
130 KiB
Org Mode
#+TITLE: Atomic Orbitals
|
|
#+SETUPFILE: ../tools/theme.setup
|
|
#+INCLUDE: ../tools/lib.org
|
|
|
|
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.
|
|
|
|
* Headers :noexport:
|
|
#+begin_src elisp :noexport :results none
|
|
(org-babel-lob-ingest "../tools/lib.org")
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :tangle (eval h_private_type)
|
|
#ifndef QMCKL_AO_HPT
|
|
#define QMCKL_AO_HPT
|
|
|
|
#include <stdbool.h>
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval c_test) :noweb yes
|
|
#include "qmckl.h"
|
|
#include "assert.h"
|
|
#ifdef HAVE_CONFIG_H
|
|
#include "config.h"
|
|
#endif
|
|
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
#include "chbrclf.h"
|
|
#include "qmckl_ao_private_func.h"
|
|
|
|
|
|
int main() {
|
|
qmckl_context context;
|
|
context = qmckl_context_create();
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval c)
|
|
#ifdef HAVE_CONFIG_H
|
|
#include "config.h"
|
|
#endif
|
|
|
|
#ifdef HAVE_STDINT_H
|
|
#include <stdint.h>
|
|
#elif HAVE_INTTYPES_H
|
|
#include <inttypes.h>
|
|
#endif
|
|
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
#include <stdbool.h>
|
|
#include <assert.h>
|
|
|
|
#include "qmckl.h"
|
|
#include "qmckl_context_private_type.h"
|
|
#include "qmckl_memory_private_type.h"
|
|
#include "qmckl_memory_private_func.h"
|
|
#include "qmckl_ao_private_type.h"
|
|
#include "qmckl_ao_private_func.h"
|
|
#+end_src
|
|
|
|
* 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 |
|
|
| ~nucleus_shell_num~ | ~[nucl_num]~ | Number of shells per 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 |
|
|
| ~ao_num~ | | Number of AOs |
|
|
| ~ao_cartesian~ | | If true, use polynomials. Otherwise, use spherical harmonics |
|
|
| ~ao_factor~ | ~[ao_num]~ | Normalization factor of the AO |
|
|
| ~ao_shell~ | ~[ao_num]~ | For each AO, specify to which shell it belongs |
|
|
|
|
Computed data:
|
|
|
|
|--------------------------+--------------------------------------+-----------------------------------------------------------------------------------------------|
|
|
| ~coefficient_normalized~ | ~[prim_num]~ | Normalized primitive coefficients |
|
|
| ~nucleus_prim_index~ | ~[nucl_num]~ | Index of the first primitive for each nucleus |
|
|
| ~nucleus_max_ang_mom~ | ~[nucl_num]~ | Maximum angular momentum for each nucleus |
|
|
| ~nucleus_range~ | ~[nucl_num]~ | Distance beyond which all the AOs are zero |
|
|
|--------------------------+--------------------------------------+-----------------------------------------------------------------------------------------------|
|
|
| ~primitive_vgl~ | ~[5][walk_num][elec_num][prim_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~ | ~[5][walk_num][elec_num][shell_num]~ | Value, gradients, Laplacian of the primitives at electron positions |
|
|
| ~shell_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the AOs at electron positions |
|
|
| ~ao_vgl~ | ~[5][walk_num][elec_num][ao_num]~ | Value, gradients, Laplacian of the primitives at electron positions |
|
|
| ~ao_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the AOs 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 |
|
|
|
|
For H_2 with the following basis set,
|
|
|
|
#+BEGIN_EXAMPLE
|
|
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
|
|
#+END_EXAMPLE
|
|
|
|
we have:
|
|
|
|
#+BEGIN_EXAMPLE
|
|
type = 'G'
|
|
shell_num = 12
|
|
prim_num = 20
|
|
ao_num = 38
|
|
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 ]
|
|
#+END_EXAMPLE
|
|
|
|
** Data structure
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_type)
|
|
typedef struct qmckl_ao_basis_struct {
|
|
int64_t shell_num;
|
|
int64_t prim_num;
|
|
int64_t ao_num;
|
|
int64_t * nucleus_index;
|
|
int64_t * nucleus_shell_num;
|
|
int32_t * shell_ang_mom;
|
|
int64_t * shell_prim_num;
|
|
int64_t * shell_prim_index;
|
|
double * shell_factor;
|
|
double * exponent;
|
|
double * coefficient;
|
|
double * prim_factor;
|
|
double * ao_factor;
|
|
|
|
int64_t * nucleus_prim_index;
|
|
double * coefficient_normalized;
|
|
int32_t * nucleus_max_ang_mom;
|
|
double * nucleus_range;
|
|
double * primitive_vgl;
|
|
int64_t primitive_vgl_date;
|
|
double * shell_vgl;
|
|
int64_t shell_vgl_date;
|
|
double * ao_vgl;
|
|
int64_t ao_vgl_date;
|
|
|
|
int32_t uninitialized;
|
|
bool provided;
|
|
bool ao_cartesian;
|
|
char type;
|
|
} qmckl_ao_basis_struct;
|
|
#+end_src
|
|
|
|
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.
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func)
|
|
qmckl_exit_code qmckl_init_ao_basis(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c)
|
|
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 << 14) - 1;
|
|
|
|
/* Default values */
|
|
ctx->ao_basis.ao_cartesian = true;
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
** Access functions
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :exports none
|
|
char qmckl_get_ao_basis_type (const qmckl_context context);
|
|
int64_t qmckl_get_ao_basis_shell_num (const qmckl_context context);
|
|
int64_t qmckl_get_ao_basis_prim_num (const qmckl_context context);
|
|
int64_t qmckl_get_ao_basis_ao_num (const qmckl_context context);
|
|
int64_t* qmckl_get_ao_basis_nucleus_index (const qmckl_context context);
|
|
int64_t* qmckl_get_ao_basis_nucleus_shell_num(const qmckl_context context);
|
|
int32_t* qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context);
|
|
int64_t* qmckl_get_ao_basis_shell_prim_num (const qmckl_context context);
|
|
int64_t* qmckl_get_ao_basis_shell_prim_index (const qmckl_context context);
|
|
double* qmckl_get_ao_basis_shell_factor (const qmckl_context context);
|
|
double* qmckl_get_ao_basis_exponent (const qmckl_context context);
|
|
double* qmckl_get_ao_basis_coefficient (const qmckl_context context);
|
|
double* qmckl_get_ao_basis_prim_factor (const qmckl_context context);
|
|
double* qmckl_get_ao_basis_ao_factor (const qmckl_context context);
|
|
#+end_src
|
|
|
|
When all the data for the AOs have been provided, the following
|
|
function returns ~true~.
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
bool qmckl_ao_basis_provided (const qmckl_context context);
|
|
#+end_src
|
|
|
|
#+NAME:post
|
|
#+begin_src c :exports none
|
|
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
|
|
return NULL;
|
|
}
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
char qmckl_get_ao_basis_type (const qmckl_context context) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return (char) 0;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1;
|
|
|
|
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
|
|
return (char) 0;
|
|
}
|
|
|
|
assert (ctx->ao_basis.type != (char) 0);
|
|
return ctx->ao_basis.type;
|
|
}
|
|
|
|
|
|
int64_t qmckl_get_ao_basis_shell_num (const qmckl_context context) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return (int64_t) 0;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 1;
|
|
|
|
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
|
|
return (int64_t) 0;
|
|
}
|
|
|
|
assert (ctx->ao_basis.shell_num > (int64_t) 0);
|
|
return ctx->ao_basis.shell_num;
|
|
}
|
|
|
|
|
|
int64_t qmckl_get_ao_basis_prim_num (const qmckl_context context) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return (int64_t) 0;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 2;
|
|
|
|
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
|
|
return (int64_t) 0;
|
|
}
|
|
|
|
assert (ctx->ao_basis.prim_num > (int64_t) 0);
|
|
return ctx->ao_basis.prim_num;
|
|
}
|
|
|
|
|
|
int64_t* qmckl_get_ao_basis_nucleus_shell_num (const qmckl_context context) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return NULL;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 3;
|
|
|
|
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
|
|
return NULL;
|
|
}
|
|
|
|
assert (ctx->ao_basis.nucleus_shell_num != NULL);
|
|
return ctx->ao_basis.nucleus_shell_num ;
|
|
}
|
|
|
|
int64_t* qmckl_get_ao_basis_nucleus_index (const qmckl_context context) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return NULL;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 4;
|
|
|
|
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
|
|
return NULL;
|
|
}
|
|
|
|
assert (ctx->ao_basis.nucleus_index != NULL);
|
|
return ctx->ao_basis.nucleus_index ;
|
|
}
|
|
|
|
|
|
int32_t* qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return NULL;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 5;
|
|
|
|
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
|
|
return NULL;
|
|
}
|
|
|
|
assert (ctx->ao_basis.shell_ang_mom != NULL);
|
|
return ctx->ao_basis.shell_ang_mom;
|
|
}
|
|
|
|
|
|
int64_t* qmckl_get_ao_basis_shell_prim_num (const qmckl_context context) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return NULL;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 6;
|
|
|
|
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
|
|
return NULL;
|
|
}
|
|
|
|
assert (ctx->ao_basis.shell_prim_num != NULL);
|
|
return ctx->ao_basis.shell_prim_num;
|
|
}
|
|
|
|
|
|
int64_t* qmckl_get_ao_basis_shell_prim_index (const qmckl_context context) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return NULL;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 7;
|
|
|
|
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
|
|
return NULL;
|
|
}
|
|
|
|
assert (ctx->ao_basis.shell_prim_index != NULL);
|
|
return ctx->ao_basis.shell_prim_index;
|
|
}
|
|
|
|
|
|
double* qmckl_get_ao_basis_shell_factor (const qmckl_context context) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return NULL;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 8;
|
|
|
|
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
|
|
return NULL;
|
|
}
|
|
|
|
assert (ctx->ao_basis.shell_factor != NULL);
|
|
return ctx->ao_basis.shell_factor;
|
|
}
|
|
|
|
|
|
double* qmckl_get_ao_basis_exponent (const qmckl_context context) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return NULL;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
|
|
int32_t mask = 1 << 9;
|
|
|
|
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
|
|
return NULL;
|
|
}
|
|
|
|
assert (ctx->ao_basis.exponent != NULL);
|
|
return ctx->ao_basis.exponent;
|
|
}
|
|
|
|
|
|
double* qmckl_get_ao_basis_coefficient (const qmckl_context context) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return NULL;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 10;
|
|
|
|
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
|
|
return NULL;
|
|
}
|
|
|
|
assert (ctx->ao_basis.coefficient != NULL);
|
|
return ctx->ao_basis.coefficient;
|
|
}
|
|
|
|
|
|
double* qmckl_get_ao_basis_prim_factor (const qmckl_context context) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return NULL;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 11;
|
|
|
|
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
|
|
return NULL;
|
|
}
|
|
|
|
assert (ctx->ao_basis.prim_factor != NULL);
|
|
return ctx->ao_basis.prim_factor;
|
|
}
|
|
|
|
|
|
int64_t qmckl_get_ao_basis_ao_num (const qmckl_context context) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return (int64_t) 0;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 12;
|
|
|
|
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
|
|
return (int64_t) 0;
|
|
}
|
|
|
|
assert (ctx->ao_basis.ao_num > (int64_t) 0);
|
|
return ctx->ao_basis.ao_num;
|
|
}
|
|
|
|
|
|
double* qmckl_get_ao_basis_ao_factor (const qmckl_context context) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return NULL;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 13;
|
|
|
|
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
|
|
return NULL;
|
|
}
|
|
|
|
assert (ctx->ao_basis.ao_factor != NULL);
|
|
return ctx->ao_basis.ao_factor;
|
|
}
|
|
|
|
|
|
bool qmckl_ao_basis_provided(const 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);
|
|
|
|
return ctx->ao_basis.provided;
|
|
}
|
|
|
|
|
|
|
|
#+end_src
|
|
|
|
** Initialization functions
|
|
|
|
To set the basis set, all the following functions need to be
|
|
called.
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
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_ao_num (qmckl_context context, const int64_t ao_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);
|
|
qmckl_exit_code qmckl_set_ao_basis_ao_factor (qmckl_context context, const double * ao_factor);
|
|
qmckl_exit_code qmckl_set_ao_basis_cartesian (qmckl_context context, const bool cartesian);
|
|
#+end_src
|
|
|
|
#+NAME:pre2
|
|
#+begin_src c :exports none
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
#+end_src
|
|
|
|
#+NAME:post2
|
|
#+begin_src c :exports none
|
|
ctx->ao_basis.uninitialized &= ~mask;
|
|
ctx->ao_basis.provided = (ctx->ao_basis.uninitialized == 0);
|
|
if (ctx->ao_basis.provided) {
|
|
qmckl_exit_code rc_ = qmckl_finalize_basis(context);
|
|
if (rc_ != QMCKL_SUCCESS) return rc_;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_set_ao_basis_type(qmckl_context context, const char t) {
|
|
<<pre2>>
|
|
|
|
if (t != 'G' && t != 'S') {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_ao_basis_type",
|
|
NULL);
|
|
}
|
|
|
|
int32_t mask = 1;
|
|
ctx->ao_basis.type = t;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
|
|
qmckl_exit_code qmckl_set_ao_basis_shell_num(qmckl_context context, const int64_t shell_num) {
|
|
<<pre2>>
|
|
|
|
if (shell_num <= 0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_ao_basis_shell_num",
|
|
"shell_num <= 0");
|
|
}
|
|
|
|
int64_t prim_num = qmckl_get_ao_basis_prim_num(context);
|
|
|
|
if (0L < prim_num && prim_num < shell_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_ao_basis_shell_num",
|
|
"shell_num > prim_num");
|
|
}
|
|
|
|
int32_t mask = 1 << 1;
|
|
ctx->ao_basis.shell_num = shell_num;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
|
|
qmckl_exit_code qmckl_set_ao_basis_prim_num(qmckl_context context, const int64_t prim_num) {
|
|
<<pre2>>
|
|
|
|
if (prim_num <= 0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_ao_basis_shell_num",
|
|
"prim_num must be positive");
|
|
}
|
|
|
|
int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
|
|
|
|
if (prim_num < shell_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_ao_basis_shell_num",
|
|
"prim_num < shell_num");
|
|
}
|
|
|
|
int32_t mask = 1 << 2;
|
|
ctx->ao_basis.prim_num = prim_num;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
|
|
qmckl_exit_code qmckl_set_ao_basis_nucleus_shell_num(qmckl_context context, const int64_t* nucleus_shell_num) {
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 3;
|
|
|
|
const int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
|
|
if (shell_num == 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_nucleus_shell_num",
|
|
"shell_num is not set");
|
|
}
|
|
|
|
if (ctx->ao_basis.nucleus_shell_num != NULL) {
|
|
qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.nucleus_shell_num);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc,
|
|
"qmckl_set_ao_basis_nucleus_shell_num",
|
|
NULL);
|
|
}
|
|
}
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = shell_num * sizeof(int64_t);
|
|
int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);
|
|
|
|
if (new_array == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_set_ao_basis_nucleus_shell_num",
|
|
NULL);
|
|
}
|
|
|
|
memcpy(new_array, nucleus_shell_num, mem_info.size);
|
|
|
|
ctx->ao_basis.nucleus_shell_num = new_array;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
qmckl_exit_code qmckl_set_ao_basis_nucleus_index(qmckl_context context, const int64_t* nucleus_index) {
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 4;
|
|
|
|
const int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
|
|
if (shell_num == 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_nucleus_index",
|
|
"shell_num is not set");
|
|
}
|
|
|
|
if (ctx->ao_basis.nucleus_index != NULL) {
|
|
qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.nucleus_index);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc,
|
|
"qmckl_set_ao_basis_nucleus_index",
|
|
NULL);
|
|
}
|
|
}
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = shell_num * sizeof(int64_t);
|
|
int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);
|
|
|
|
if (new_array == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_set_ao_basis_nucleus_index",
|
|
NULL);
|
|
}
|
|
|
|
memcpy(new_array, nucleus_index, mem_info.size);
|
|
|
|
ctx->ao_basis.nucleus_index = new_array;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
|
|
qmckl_exit_code qmckl_set_ao_basis_shell_ang_mom(qmckl_context context, const int32_t* shell_ang_mom) {
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 5;
|
|
|
|
const int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
|
|
if (shell_num == 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_shell_ang_mom",
|
|
"shell_num is not set");
|
|
}
|
|
|
|
if (ctx->ao_basis.shell_ang_mom != NULL) {
|
|
qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_ang_mom);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc,
|
|
"qmckl_set_ao_basis_shell_ang_mom",
|
|
NULL);
|
|
}
|
|
}
|
|
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = shell_num * sizeof(int32_t);
|
|
int32_t * new_array = (int32_t*) qmckl_malloc(context, mem_info);
|
|
|
|
if (new_array == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_set_ao_basis_shell_ang_mom",
|
|
NULL);
|
|
}
|
|
|
|
memcpy(new_array, shell_ang_mom, mem_info.size);
|
|
|
|
ctx->ao_basis.shell_ang_mom = new_array;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
|
|
qmckl_exit_code qmckl_set_ao_basis_shell_prim_num(qmckl_context context, const int64_t* shell_prim_num) {
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 6;
|
|
|
|
const int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
|
|
if (shell_num == 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_shell_prim_num",
|
|
"shell_num is not set");
|
|
}
|
|
|
|
if (ctx->ao_basis.shell_prim_num != NULL) {
|
|
qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_prim_num);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc,
|
|
"qmckl_set_ao_basis_shell_prim_num",
|
|
NULL);
|
|
}
|
|
}
|
|
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = shell_num * sizeof(int64_t);
|
|
int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);
|
|
|
|
if (new_array == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_set_ao_basis_shell_prim_num",
|
|
NULL);
|
|
}
|
|
|
|
memcpy(new_array, shell_prim_num, mem_info.size);
|
|
|
|
ctx->ao_basis.shell_prim_num = new_array;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
|
|
qmckl_exit_code qmckl_set_ao_basis_shell_prim_index(qmckl_context context, const int64_t* shell_prim_index) {
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 7;
|
|
|
|
const int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
|
|
if (shell_num == 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_shell_prim_index",
|
|
"shell_num is not set");
|
|
}
|
|
|
|
if (ctx->ao_basis.shell_prim_index != NULL) {
|
|
qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_prim_index);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc,
|
|
"qmckl_set_ao_basis_shell_prim_index",
|
|
NULL);
|
|
}
|
|
}
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = shell_num * sizeof(int64_t);
|
|
int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);
|
|
|
|
if (new_array == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_set_ao_basis_shell_prim_index",
|
|
NULL);
|
|
}
|
|
|
|
memcpy(new_array, shell_prim_index, mem_info.size);
|
|
|
|
ctx->ao_basis.shell_prim_index = new_array;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
|
|
qmckl_exit_code qmckl_set_ao_basis_shell_factor(qmckl_context context, const double* shell_factor) {
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 8;
|
|
|
|
const int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
|
|
if (shell_num == 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_shell_factor",
|
|
"shell_num is not set");
|
|
}
|
|
|
|
|
|
if (ctx->ao_basis.shell_factor != NULL) {
|
|
qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_factor);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc,
|
|
"qmckl_set_ao_basis_shell_factor",
|
|
NULL);
|
|
}
|
|
}
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = shell_num * sizeof(double);
|
|
double* new_array = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (new_array == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_set_ao_basis_shell_factor",
|
|
NULL);
|
|
}
|
|
|
|
memcpy(new_array, shell_factor, mem_info.size);
|
|
|
|
ctx->ao_basis.shell_factor = new_array;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
qmckl_exit_code qmckl_set_ao_basis_exponent(qmckl_context context, const double* exponent) {
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 9;
|
|
|
|
const int64_t prim_num = qmckl_get_ao_basis_prim_num(context);
|
|
if (prim_num == 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_exponent",
|
|
"prim_num is not set");
|
|
}
|
|
|
|
if (ctx->ao_basis.exponent != NULL) {
|
|
qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.exponent);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc,
|
|
"qmckl_set_ao_basis_exponent",
|
|
NULL);
|
|
}
|
|
}
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = prim_num * sizeof(double);
|
|
double* new_array = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (new_array == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_set_ao_basis_exponent",
|
|
NULL);
|
|
}
|
|
|
|
memcpy(new_array, exponent, mem_info.size);
|
|
|
|
ctx->ao_basis.exponent = new_array;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
qmckl_exit_code qmckl_set_ao_basis_coefficient(qmckl_context context, const double* coefficient) {
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 10;
|
|
|
|
const int64_t prim_num = qmckl_get_ao_basis_prim_num(context);
|
|
if (prim_num == 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_coefficient",
|
|
"prim_num is not set");
|
|
}
|
|
|
|
if (ctx->ao_basis.coefficient != NULL) {
|
|
qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.coefficient);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc,
|
|
"qmckl_set_ao_basis_coefficient",
|
|
NULL);
|
|
}
|
|
}
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = prim_num * sizeof(double);
|
|
double* new_array = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (new_array == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_set_ao_basis_coefficient",
|
|
NULL);
|
|
}
|
|
|
|
memcpy(new_array, coefficient, mem_info.size);
|
|
|
|
ctx->ao_basis.coefficient = new_array;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
|
|
qmckl_exit_code qmckl_set_ao_basis_prim_factor(qmckl_context context, const double* prim_factor) {
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 11;
|
|
|
|
const int64_t prim_num = qmckl_get_ao_basis_prim_num(context);
|
|
if (prim_num == 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_prim_factor",
|
|
"prim_num is not set");
|
|
}
|
|
|
|
|
|
if (ctx->ao_basis.prim_factor != NULL) {
|
|
qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.prim_factor);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc,
|
|
"qmckl_set_ao_basis_prim_factor",
|
|
NULL);
|
|
}
|
|
}
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = prim_num * sizeof(double);
|
|
double* new_array = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (new_array == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_set_ao_basis_prim_factor",
|
|
NULL);
|
|
}
|
|
|
|
memcpy(new_array, prim_factor, mem_info.size);
|
|
|
|
ctx->ao_basis.prim_factor = new_array;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
|
|
qmckl_exit_code qmckl_set_ao_basis_ao_num(qmckl_context context, const int64_t ao_num) {
|
|
<<pre2>>
|
|
|
|
if (ao_num <= 0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_ao_basis_shell_num",
|
|
"ao_num must be positive");
|
|
}
|
|
|
|
int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
|
|
|
|
if (ao_num < shell_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_ao_basis_shell_num",
|
|
"ao_num < shell_num");
|
|
}
|
|
|
|
int32_t mask = 1 << 12;
|
|
ctx->ao_basis.ao_num = ao_num;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
|
|
qmckl_exit_code qmckl_set_ao_basis_ao_factor(qmckl_context context, const double* ao_factor) {
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 13;
|
|
|
|
const int64_t ao_num = qmckl_get_ao_basis_ao_num(context);
|
|
if (ao_num == 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_ao_factor",
|
|
"ao_num is not set");
|
|
}
|
|
|
|
|
|
if (ctx->ao_basis.ao_factor != NULL) {
|
|
qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.ao_factor);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc,
|
|
"qmckl_set_ao_basis_ao_factor",
|
|
NULL);
|
|
}
|
|
}
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = ao_num * sizeof(double);
|
|
double* new_array = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (new_array == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_set_ao_basis_ao_factor",
|
|
NULL);
|
|
}
|
|
|
|
memcpy(new_array, ao_factor, mem_info.size);
|
|
|
|
ctx->ao_basis.ao_factor = new_array;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
|
|
qmckl_exit_code qmckl_set_ao_basis_cartesian(qmckl_context context, const bool t) {
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1;
|
|
ctx->ao_basis.ao_cartesian = t;
|
|
|
|
<<post2>>
|
|
}
|
|
#+end_src
|
|
|
|
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.
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_finalize_basis(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_finalize_basis(qmckl_context context) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_INVALID_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int64_t nucl_num = 0;
|
|
qmckl_exit_code rc = QMCKL_FAILURE;
|
|
|
|
rc = qmckl_get_nucleus_num(context, &nucl_num);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* nucleus_prim_index */
|
|
{
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = (ctx->nucleus.num + (int64_t) 1) * sizeof(int64_t);
|
|
|
|
ctx->ao_basis.nucleus_prim_index = (int64_t*) qmckl_malloc(context, mem_info);
|
|
|
|
if (ctx->ao_basis.nucleus_prim_index == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"ao_basis.nucleus_prim_index",
|
|
NULL);
|
|
}
|
|
|
|
for (int64_t i=0 ; i<nucl_num ; ++i) {
|
|
int64_t shell_idx = ctx->ao_basis.nucleus_index[i];
|
|
ctx->ao_basis.nucleus_prim_index[i] = ctx->ao_basis.shell_prim_index[shell_idx];
|
|
}
|
|
ctx->ao_basis.nucleus_prim_index[nucl_num] = ctx->ao_basis.prim_num;
|
|
}
|
|
|
|
|
|
/* Normalize coefficients */
|
|
{
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = ctx->ao_basis.prim_num * sizeof(double);
|
|
|
|
ctx->ao_basis.coefficient_normalized = (double *) qmckl_malloc(context, mem_info);
|
|
|
|
if (ctx->ao_basis.coefficient_normalized == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"ao_basis.coefficient_normalized",
|
|
NULL);
|
|
}
|
|
|
|
for (int64_t ishell=0 ; ishell < ctx->ao_basis.shell_num ; ++ishell) {
|
|
for (int64_t iprim=ctx->ao_basis.shell_prim_index[ishell] ;
|
|
iprim < ctx->ao_basis.shell_prim_index[ishell]+ctx->ao_basis.shell_prim_num[ishell] ;
|
|
++iprim) {
|
|
ctx->ao_basis.coefficient_normalized[iprim] =
|
|
ctx->ao_basis.coefficient[iprim] * ctx->ao_basis.prim_factor[iprim] *
|
|
ctx->ao_basis.shell_factor[ishell];
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
/* Find max angular momentum on each nucleus */
|
|
{
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = ctx->nucleus.num * sizeof(int32_t);
|
|
|
|
ctx->ao_basis.nucleus_max_ang_mom = (int32_t *) qmckl_malloc(context, mem_info);
|
|
|
|
if (ctx->ao_basis.nucleus_max_ang_mom == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"ao_basis.nucleus_max_ang_mom",
|
|
NULL);
|
|
}
|
|
|
|
for (int64_t inucl=0 ; inucl < nucl_num ; ++inucl) {
|
|
ctx->ao_basis.nucleus_max_ang_mom[inucl] = 0;
|
|
for (int64_t ishell=ctx->ao_basis.nucleus_index[inucl] ;
|
|
ishell < ctx->ao_basis.nucleus_index[inucl] + ctx->ao_basis.nucleus_shell_num[inucl] ;
|
|
++ishell) {
|
|
ctx->ao_basis.nucleus_max_ang_mom[inucl] =
|
|
ctx->ao_basis.nucleus_max_ang_mom[inucl] > ctx->ao_basis.shell_ang_mom[ishell] ?
|
|
ctx->ao_basis.nucleus_max_ang_mom[inucl] : ctx->ao_basis.shell_ang_mom[ishell] ;
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Find distance beyond which all AOs are zero.
|
|
The distance is obtained by sqrt(log(cutoff)*range) */
|
|
{
|
|
if (ctx->ao_basis.type == 'G') {
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = ctx->nucleus.num * sizeof(double);
|
|
|
|
ctx->ao_basis.nucleus_range = (double *) qmckl_malloc(context, mem_info);
|
|
|
|
if (ctx->ao_basis.nucleus_range == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"ao_basis.nucleus_range",
|
|
NULL);
|
|
}
|
|
|
|
for (int64_t inucl=0 ; inucl < ctx->nucleus.num ; ++inucl) {
|
|
ctx->ao_basis.nucleus_range[inucl] = 0.;
|
|
for (int64_t ishell=ctx->ao_basis.nucleus_index[inucl] ;
|
|
ishell < ctx->ao_basis.nucleus_index[inucl] + ctx->ao_basis.nucleus_shell_num[inucl] ;
|
|
++ishell) {
|
|
for (int64_t iprim=ctx->ao_basis.shell_prim_index[ishell] ;
|
|
iprim < ctx->ao_basis.shell_prim_index[ishell] + ctx->ao_basis.shell_prim_num[ishell] ;
|
|
++iprim) {
|
|
double range = 1./ctx->ao_basis.exponent[iprim];
|
|
ctx->ao_basis.nucleus_range[inucl] =
|
|
ctx->ao_basis.nucleus_range[inucl] > range ?
|
|
ctx->ao_basis.nucleus_range[inucl] : range;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
/* TODO : sort the basis set here */
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
** Fortran interfaces
|
|
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_type (context, t) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
character(c_char) , intent(in) , value :: t
|
|
end function
|
|
end interface
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_shell_num(context, num) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: num
|
|
end function
|
|
end interface
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_prim_num(context, num) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: num
|
|
end function
|
|
end interface
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_nucleus_index(context, idx) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) :: idx(*)
|
|
end function
|
|
end interface
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_nucleus_shell_num(context,shell_num) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) :: shell_num(*)
|
|
end function
|
|
end interface
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_shell_ang_mom(context,shell_ang_mom) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int32_t) , intent(in) :: shell_ang_mom(*)
|
|
end function
|
|
end interface
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_shell_prim_num(context,shell_prim_num) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) :: shell_prim_num(*)
|
|
end function
|
|
end interface
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_shell_prim_index(context,shell_prim_index) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) :: shell_prim_index(*)
|
|
end function
|
|
end interface
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_shell_factor(context,shell_factor) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double) , intent(in) :: shell_factor(*)
|
|
end function
|
|
end interface
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_exponent(context,exponent) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double) , intent(in) :: exponent(*)
|
|
end function
|
|
end interface
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_coefficient(context,coefficient) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double) , intent(in) :: coefficient(*)
|
|
end function
|
|
end interface
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_prim_factor(context,prim_factor) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double) , intent(in) :: prim_factor(*)
|
|
end function
|
|
end interface
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_ao_num(context, num) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: num
|
|
end function
|
|
end interface
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_ao_factor(context,ao_factor) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double) , intent(in) :: ao_factor(*)
|
|
end function
|
|
end interface
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_cartesian(context,cartesian) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
logical (c_bool) , intent(in) , value :: cartesian
|
|
end function
|
|
end interface
|
|
#+end_src
|
|
|
|
** Test :noexport:
|
|
|
|
#+begin_src c :tangle (eval c_test) :exports none :exports none
|
|
const int64_t nucl_num = chbrclf_nucl_num;
|
|
const double* nucl_charge = chbrclf_charge;
|
|
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
|
|
|
|
qmckl_exit_code rc;
|
|
rc = qmckl_set_nucleus_num (context, nucl_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_nucleus_charge(context, nucl_charge);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
assert(qmckl_nucleus_provided(context));
|
|
|
|
|
|
const int64_t shell_num = chbrclf_shell_num;
|
|
const int64_t prim_num = chbrclf_prim_num;
|
|
const int64_t ao_num = chbrclf_ao_num;
|
|
const int64_t * nucleus_index = &(chbrclf_basis_nucleus_index[0]);
|
|
const int64_t * nucleus_shell_num = &(chbrclf_basis_nucleus_shell_num[0]);
|
|
const int32_t * shell_ang_mom = &(chbrclf_basis_shell_ang_mom[0]);
|
|
const int64_t * shell_prim_num = &(chbrclf_basis_shell_prim_num[0]);
|
|
const int64_t * shell_prim_index = &(chbrclf_basis_shell_prim_index[0]);
|
|
const double * shell_factor = &(chbrclf_basis_shell_factor[0]);
|
|
const double * exponent = &(chbrclf_basis_exponent[0]);
|
|
const double * coefficient = &(chbrclf_basis_coefficient[0]);
|
|
const double * prim_factor = &(chbrclf_basis_prim_factor[0]);
|
|
const double * ao_factor = &(chbrclf_basis_ao_factor[0]);
|
|
|
|
const char typ = 'G';
|
|
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_type (context, typ);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_shell_num (context, shell_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_prim_num (context, prim_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_shell_factor (context, shell_factor);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_exponent (context, exponent);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_coefficient (context, coefficient);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_prim_factor (context, prim_factor);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_ao_basis_ao_num(context, ao_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_ao_basis_ao_factor (context, ao_factor);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
assert(qmckl_ao_basis_provided(context));
|
|
|
|
int64_t shell_num_test ;
|
|
int64_t prim_num_test ;
|
|
int64_t ao_num_test ;
|
|
int64_t * nucleus_index_test ;
|
|
int64_t * nucleus_shell_num_test;
|
|
int32_t * shell_ang_mom_test ;
|
|
int64_t * shell_prim_num_test ;
|
|
int64_t * shell_prim_index_test ;
|
|
double * shell_factor_test ;
|
|
double * exponent_test ;
|
|
double * coefficient_test ;
|
|
double * prim_factor_test ;
|
|
double * ao_factor_test ;
|
|
char typ_test ;
|
|
|
|
|
|
typ_test = qmckl_get_ao_basis_type (context);
|
|
assert(typ == typ_test);
|
|
|
|
shell_num_test = qmckl_get_ao_basis_shell_num (context);
|
|
assert(shell_num == shell_num_test);
|
|
|
|
prim_num_test = qmckl_get_ao_basis_prim_num (context);
|
|
assert(prim_num == prim_num_test);
|
|
|
|
nucleus_index_test = qmckl_get_ao_basis_nucleus_index (context);
|
|
for (int64_t i=0 ; i < nucl_num ; ++i) {
|
|
assert(nucleus_index_test[i] == nucleus_index[i]);
|
|
}
|
|
|
|
nucleus_shell_num_test = qmckl_get_ao_basis_nucleus_shell_num (context);
|
|
for (int64_t i=0 ; i < nucl_num ; ++i) {
|
|
assert(nucleus_shell_num_test[i] == nucleus_shell_num[i]);
|
|
}
|
|
|
|
shell_ang_mom_test = qmckl_get_ao_basis_shell_ang_mom (context);
|
|
for (int64_t i=0 ; i < shell_num ; ++i) {
|
|
assert(shell_ang_mom_test[i] == shell_ang_mom[i]);
|
|
}
|
|
|
|
shell_factor_test = qmckl_get_ao_basis_shell_factor (context);
|
|
for (int64_t i=0 ; i < shell_num ; ++i) {
|
|
assert(shell_factor_test[i] == shell_factor[i]);
|
|
}
|
|
|
|
shell_prim_num_test = qmckl_get_ao_basis_shell_prim_num (context);
|
|
for (int64_t i=0 ; i < shell_num ; ++i) {
|
|
assert(shell_prim_num_test[i] == shell_prim_num[i]);
|
|
}
|
|
|
|
shell_prim_index_test = qmckl_get_ao_basis_shell_prim_index (context);
|
|
for (int64_t i=0 ; i < shell_num ; ++i) {
|
|
assert(shell_prim_index_test[i] == shell_prim_index[i]);
|
|
}
|
|
|
|
exponent_test = qmckl_get_ao_basis_exponent(context);
|
|
for (int64_t i=0 ; i < prim_num ; ++i) {
|
|
assert(exponent_test[i] == exponent[i]);
|
|
}
|
|
|
|
coefficient_test = qmckl_get_ao_basis_coefficient(context);
|
|
for (int64_t i=0 ; i < prim_num ; ++i) {
|
|
assert(coefficient_test[i] == coefficient[i]);
|
|
}
|
|
|
|
prim_factor_test = qmckl_get_ao_basis_prim_factor (context);
|
|
for (int64_t i=0 ; i < prim_num ; ++i) {
|
|
assert(prim_factor_test[i] == prim_factor[i]);
|
|
}
|
|
|
|
ao_num_test = qmckl_get_ao_basis_ao_num(context);
|
|
assert(ao_num == ao_num_test);
|
|
|
|
ao_factor_test = qmckl_get_ao_basis_ao_factor (context);
|
|
for (int64_t i=0 ; i < ao_num ; ++i) {
|
|
assert(ao_factor_test[i] == ao_factor[i]);
|
|
}
|
|
|
|
#+end_src
|
|
|
|
* Radial part
|
|
|
|
** TODO Helper functions to accelerate calculations
|
|
** 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
|
|
|
|
#+begin_src c :tangle (eval h_func)
|
|
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);
|
|
#+end_src
|
|
|
|
#+begin_src f90 :tangle (eval f)
|
|
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
|
|
#+end_src
|
|
|
|
#+begin_src f90 :tangle (eval f) :exports none
|
|
integer(c_int32_t) function qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) &
|
|
bind(C) result(info)
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double) , intent(in) :: X(3), R(3)
|
|
integer (c_int64_t) , intent(in) , value :: n
|
|
real (c_double) , intent(in) :: A(n)
|
|
real (c_double) , intent(out) :: VGL(ldv,5)
|
|
integer (c_int64_t) , intent(in) , value :: ldv
|
|
|
|
integer, external :: qmckl_ao_gaussian_vgl_f
|
|
info = qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv)
|
|
end function qmckl_ao_gaussian_vgl
|
|
#+end_src
|
|
|
|
#+begin_src f90 :tangle (eval fh_func) :exports none
|
|
interface
|
|
integer(c_int32_t) function qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: ldv
|
|
integer (c_int64_t) , intent(in) , value :: n
|
|
real (c_double) , intent(in) :: X(3), R(3), A(n)
|
|
real (c_double) , intent(out) :: VGL(ldv,5)
|
|
end function qmckl_ao_gaussian_vgl
|
|
end interface
|
|
#+end_src
|
|
|
|
# Test
|
|
#+begin_src f90 :tangle (eval f_test)
|
|
integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C)
|
|
|
|
use qmckl
|
|
use qmckl_probes_f
|
|
|
|
implicit none
|
|
|
|
integer(c_int64_t), intent(in), value :: context
|
|
logical(C_BOOL) :: vfc_err
|
|
|
|
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)
|
|
|
|
#ifdef VFC_CI
|
|
! Multplying epsilon by 16 = 2^4 is equivalent to asking 4 significant digits
|
|
! less. This makes sense because we are adding noise with MCA so we can't be
|
|
! as strict on the accuracy target.
|
|
epsilon = epsilon * 16
|
|
#endif
|
|
|
|
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_check("ao"//C_NULL_CHAR, "gaussian_vgl_2_1"//C_NULL_CHAR, &
|
|
DBLE(VGL(2,1)), DBLE(0), DBLE(epsilon))
|
|
|
|
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "gaussian_vgl_2_2"//C_NULL_CHAR, &
|
|
DBLE(VGL(2,2)), DBLE(0), DBLE(epsilon))
|
|
|
|
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "gaussian_vgl_2_3"//C_NULL_CHAR, &
|
|
DBLE(VGL(2,3)), DBLE(0), DBLE(epsilon))
|
|
|
|
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "gaussian_vgl_2_4"//C_NULL_CHAR, &
|
|
DBLE(VGL(2,4)), DBLE(0), DBLE(epsilon))
|
|
|
|
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "gaussian_vgl_2_5"//C_NULL_CHAR, &
|
|
DBLE(VGL(2,5)), DBLE(0), DBLE(epsilon))
|
|
|
|
#ifndef VFC_CI
|
|
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
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval c_test) :exports none
|
|
int test_qmckl_ao_gaussian_vgl(qmckl_context context);
|
|
assert(0 == test_qmckl_ao_gaussian_vgl(context));
|
|
#+end_src
|
|
|
|
** TODO General functions for Slater basis functions
|
|
** TODO General functions for Radial functions on a grid
|
|
** Computation of primitives
|
|
|
|
*** Get
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_ao_basis_primitive_vgl(qmckl_context context, double* const primitive_vgl);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_ao_basis_primitive_vgl(qmckl_context context, double* const primitive_vgl) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_ao_basis_primitive_vgl(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = ctx->ao_basis.prim_num * 5 * ctx->electron.num * ctx->electron.walk_num;
|
|
memcpy(primitive_vgl, ctx->ao_basis.primitive_vgl, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
if (!ctx->ao_basis.provided) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_ao_basis_primitive_vgl",
|
|
NULL);
|
|
}
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->electron.coord_new_date > ctx->ao_basis.primitive_vgl_date) {
|
|
|
|
/* Allocate array */
|
|
if (ctx->ao_basis.primitive_vgl == NULL) {
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = ctx->ao_basis.prim_num * 5 * ctx->electron.num *
|
|
ctx->electron.walk_num * sizeof(double);
|
|
double* primitive_vgl = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (primitive_vgl == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_ao_basis_primitive_vgl",
|
|
NULL);
|
|
}
|
|
ctx->ao_basis.primitive_vgl = primitive_vgl;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
if (ctx->ao_basis.type == 'G') {
|
|
rc = qmckl_compute_ao_basis_primitive_gaussian_vgl(context,
|
|
ctx->ao_basis.prim_num,
|
|
ctx->electron.num,
|
|
ctx->nucleus.num,
|
|
ctx->electron.walk_num,
|
|
ctx->ao_basis.nucleus_prim_index,
|
|
ctx->electron.coord_new,
|
|
ctx->nucleus.coord,
|
|
ctx->ao_basis.exponent,
|
|
ctx->ao_basis.primitive_vgl);
|
|
} else {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"compute_ao_basis_primitive_vgl",
|
|
"Not yet implemented");
|
|
}
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->ao_basis.primitive_vgl_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_ao_basis_primitive_gaussian_vgl
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_ao_basis_primitive_gaussian_vgl_args
|
|
| 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[5][walk_num][elec_num][prim_num] | out | Value, gradients and Laplacian of the primitives |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
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(prim_num,elec_num,walk_num,5)
|
|
|
|
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(iprim, ielec, iwalk, 1) = v
|
|
primitive_vgl(iprim, ielec, iwalk, 2) = two_a * x
|
|
primitive_vgl(iprim, ielec, iwalk, 3) = two_a * y
|
|
primitive_vgl(iprim, ielec, iwalk, 4) = two_a * z
|
|
primitive_vgl(iprim, ielec, iwalk, 5) = 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
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval h_private_func) :comments org :exports none
|
|
qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl(
|
|
const qmckl_context context,
|
|
const int64_t prim_num,
|
|
const int64_t elec_num,
|
|
const int64_t nucl_num,
|
|
const int64_t walk_num,
|
|
const int64_t* nucleus_prim_index,
|
|
const double* elec_coord,
|
|
const double* nucl_coord,
|
|
const double* expo,
|
|
double* const primitive_vgl);
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_ao_basis_primitive_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_basis_primitive_gaussian_vgl"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_ao_basis_primitive_gaussian_vgl &
|
|
(context, prim_num, elec_num, nucl_num, walk_num, nucleus_prim_index, elec_coord, nucl_coord, expo, primitive_vgl) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: prim_num
|
|
integer (c_int64_t) , intent(in) , value :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
integer (c_int64_t) , intent(in) :: nucleus_prim_index(nucl_num)
|
|
real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num)
|
|
real (c_double ) , intent(in) :: nucl_coord(elec_num,3)
|
|
real (c_double ) , intent(in) :: expo(prim_num)
|
|
real (c_double ) , intent(out) :: primitive_vgl(elec_num,walk_num,5,prim_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_ao_basis_primitive_gaussian_vgl_f
|
|
info = 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)
|
|
|
|
end function qmckl_compute_ao_basis_primitive_gaussian_vgl
|
|
#+end_src
|
|
|
|
#+begin_src python :results output :exports none
|
|
import numpy as np
|
|
|
|
def f(a,x,y):
|
|
return np.exp( -a*(np.linalg.norm(x-y))**2 )
|
|
|
|
def df(a,x,y,n):
|
|
h0 = 1.e-6
|
|
if n == 1: h = np.array([h0,0.,0.])
|
|
elif n == 2: h = np.array([0.,h0,0.])
|
|
elif n == 3: h = np.array([0.,0.,h0])
|
|
return ( f(a,x+h,y) - f(a,x-h,y) ) / (2.*h0)
|
|
|
|
def d2f(a,x,y,n):
|
|
h0 = 1.e-6
|
|
if n == 1: h = np.array([h0,0.,0.])
|
|
elif n == 2: h = np.array([0.,h0,0.])
|
|
elif n == 3: h = np.array([0.,0.,h0])
|
|
return ( f(a,x+h,y) - 2.*f(a,x,y) + f(a,x-h,y) ) / h0**2
|
|
|
|
def lf(a,x,y):
|
|
return d2f(a,x,y,1) + d2f(a,x,y,2) + d2f(a,x,y,3)
|
|
|
|
elec_26_w1 = np.array( [ 1.49050402641, 2.90106987953, -1.05920815468 ] )
|
|
elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] )
|
|
nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] )
|
|
nucl_2 = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] )
|
|
|
|
#double prim_vgl[prim_num][5][walk_num][elec_num];
|
|
a = 0.9059; x = elec_26_w1 ; y = nucl_1
|
|
print ( "[7][0][0][26] : %e"% f(a,x,y))
|
|
print ( "[7][1][0][26] : %e"% df(a,x,y,1))
|
|
print ( "[7][2][0][26] : %e"% df(a,x,y,2))
|
|
print ( "[7][3][0][26] : %e"% df(a,x,y,3))
|
|
print ( "[7][4][0][26] : %e"% lf(a,x,y))
|
|
|
|
a = 0.32578; x = elec_15_w2 ; y = nucl_2
|
|
print ( "[39][0][1][15] : %e"% f(a,x,y))
|
|
print ( "[39][1][1][15] : %e"% df(a,x,y,1))
|
|
print ( "[39][2][1][15] : %e"% df(a,x,y,2))
|
|
print ( "[39][3][1][15] : %e"% df(a,x,y,3))
|
|
print ( "[39][4][1][15] : %e"% lf(a,x,y))
|
|
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
#+begin_example
|
|
[7][0][0][26] : 1.050157e-03
|
|
[7][1][0][26] : -7.501497e-04
|
|
[7][2][0][26] : -3.825069e-03
|
|
[7][3][0][26] : 3.495056e-03
|
|
[7][4][0][26] : 2.040013e-02
|
|
[39][0][1][15] : 1.083038e-03
|
|
[39][1][1][15] : 2.378275e-03
|
|
[39][2][1][15] : 2.143086e-03
|
|
[39][3][1][15] : 4.332750e-04
|
|
[39][4][1][15] : 7.514605e-03
|
|
#+end_example
|
|
|
|
*** Test
|
|
|
|
#+begin_src c :tangle (eval c_test) :exports none
|
|
{
|
|
#define walk_num chbrclf_walk_num
|
|
#define elec_num chbrclf_elec_num
|
|
#define prim_num chbrclf_prim_num
|
|
|
|
int64_t elec_up_num = chbrclf_elec_up_num;
|
|
int64_t elec_dn_num = chbrclf_elec_dn_num;
|
|
double* elec_coord = &(chbrclf_elec_coord[0][0][0]);
|
|
|
|
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_electron_walk_num (context, walk_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
assert(qmckl_electron_provided(context));
|
|
|
|
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
|
|
|
|
double prim_vgl[5][walk_num][elec_num][prim_num];
|
|
|
|
rc = qmckl_get_ao_basis_primitive_vgl(context, &(prim_vgl[0][0][0][0]));
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
assert( fabs(prim_vgl[0][0][26][7] - ( 1.0501570432064878E-003)) < 1.e-14 );
|
|
assert( fabs(prim_vgl[1][0][26][7] - (-7.5014974095310560E-004)) < 1.e-14 );
|
|
assert( fabs(prim_vgl[2][0][26][7] - (-3.8250692897610380E-003)) < 1.e-14 );
|
|
assert( fabs(prim_vgl[3][0][26][7] - ( 3.4950559194080275E-003)) < 1.e-14 );
|
|
assert( fabs(prim_vgl[4][0][26][7] - ( 2.0392163767356572E-002)) < 1.e-14 );
|
|
|
|
assert( fabs(prim_vgl[0][1][15][39] - ( 1.0825844173157661E-003)) < 1.e-14 );
|
|
assert( fabs(prim_vgl[1][1][15][39] - ( 2.3774237611651531E-003)) < 1.e-14 );
|
|
assert( fabs(prim_vgl[2][1][15][39] - ( 2.1423191526963063E-003)) < 1.e-14 );
|
|
assert( fabs(prim_vgl[3][1][15][39] - ( 4.3312003523048492E-004)) < 1.e-14 );
|
|
assert( fabs(prim_vgl[4][1][15][39] - ( 7.5174404780004771E-003)) < 1.e-14 );
|
|
|
|
}
|
|
|
|
|
|
#+end_src
|
|
|
|
*** Ideas for improvement
|
|
|
|
#+begin_src c
|
|
// 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
|
|
#+end_src
|
|
|
|
** Computation of shells
|
|
|
|
*** Get
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_ao_basis_shell_vgl(qmckl_context context, double* const shell_vgl);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_ao_basis_shell_vgl(qmckl_context context, double* const shell_vgl) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_ao_basis_shell_vgl(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = ctx->ao_basis.shell_num * 5 * ctx->electron.num * ctx->electron.walk_num;
|
|
memcpy(shell_vgl, ctx->ao_basis.shell_vgl, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_ao_basis_shell_vgl (context, shell_vgl) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
double precision, intent(out) :: shell_vgl(*)
|
|
end function
|
|
end interface
|
|
#+end_src
|
|
|
|
*** Provide
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
if (!ctx->ao_basis.provided) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_ao_basis_shell_vgl",
|
|
NULL);
|
|
}
|
|
|
|
if(!(ctx->electron.provided)) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_electron",
|
|
NULL);
|
|
}
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->electron.coord_new_date > ctx->ao_basis.shell_vgl_date) {
|
|
|
|
/* Allocate array */
|
|
if (ctx->ao_basis.shell_vgl == NULL) {
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = ctx->ao_basis.shell_num * 5 * ctx->electron.num *
|
|
ctx->electron.walk_num * sizeof(double);
|
|
double* shell_vgl = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (shell_vgl == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_ao_basis_shell_vgl",
|
|
NULL);
|
|
}
|
|
ctx->ao_basis.shell_vgl = shell_vgl;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
if (ctx->ao_basis.type == 'G') {
|
|
rc = qmckl_compute_ao_basis_shell_gaussian_vgl(context,
|
|
ctx->ao_basis.prim_num,
|
|
ctx->ao_basis.shell_num,
|
|
ctx->electron.num,
|
|
ctx->nucleus.num,
|
|
ctx->electron.walk_num,
|
|
ctx->ao_basis.nucleus_shell_num,
|
|
ctx->ao_basis.nucleus_index,
|
|
ctx->ao_basis.shell_prim_index,
|
|
ctx->ao_basis.shell_prim_num,
|
|
ctx->electron.coord_new,
|
|
ctx->nucleus.coord,
|
|
ctx->ao_basis.exponent,
|
|
ctx->ao_basis.coefficient_normalized,
|
|
ctx->ao_basis.shell_vgl);
|
|
} else {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"compute_ao_basis_shell_vgl",
|
|
"Not yet implemented");
|
|
}
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->ao_basis.shell_vgl_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_ao_basis_shell_gaussian_vgl
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_ao_basis_shell_gaussian_vgl_args
|
|
| ~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_normalized[prim_num]~ | in | Coefficients of the primitives |
|
|
| ~double~ | ~shell_vgl[5][walk_num][elec_num][shell_num]~ | out | Value, gradients and Laplacian of the shells |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
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_normalized, 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_normalized(prim_num)
|
|
double precision , intent(out) :: shell_vgl(shell_num,elec_num,walk_num,5)
|
|
|
|
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 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
|
|
|
|
do ishell=nucleus_index(inucl)+1, nucleus_index(inucl)+nucleus_shell_num(inucl)
|
|
! C is zero-based, so shift bounds by one
|
|
|
|
shell_vgl(ishell, ielec, iwalk, 1) = 0.d0
|
|
shell_vgl(ishell, ielec, iwalk, 2) = 0.d0
|
|
shell_vgl(ishell, ielec, iwalk, 3) = 0.d0
|
|
shell_vgl(ishell, ielec, iwalk, 4) = 0.d0
|
|
shell_vgl(ishell, ielec, iwalk, 5) = 0.d0
|
|
|
|
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_normalized(iprim) * dexp(-ar2)
|
|
two_a = -2.d0 * expo(iprim) * v
|
|
|
|
shell_vgl(ishell, ielec, iwalk, 1) = &
|
|
shell_vgl(ishell, ielec, iwalk, 1) + v
|
|
|
|
shell_vgl(ishell, ielec, iwalk, 2) = &
|
|
shell_vgl(ishell, ielec, iwalk, 2) + two_a * x
|
|
|
|
shell_vgl(ishell, ielec, iwalk, 3) = &
|
|
shell_vgl(ishell, ielec, iwalk, 3) + two_a * y
|
|
|
|
shell_vgl(ishell, ielec, iwalk, 4) = &
|
|
shell_vgl(ishell, ielec, iwalk, 4) + two_a * z
|
|
|
|
shell_vgl(ishell, ielec, iwalk, 5) = &
|
|
shell_vgl(ishell, ielec, iwalk, 5) + 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
|
|
#+end_src
|
|
|
|
# #+CALL: generate_c_header(table=qmckl_ao_basis_shell_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_basis_shell_gaussian_vgl"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_private_func) :comments org
|
|
qmckl_exit_code qmckl_compute_ao_basis_shell_gaussian_vgl (
|
|
const qmckl_context context,
|
|
const int64_t prim_num,
|
|
const int64_t shell_num,
|
|
const int64_t elec_num,
|
|
const int64_t nucl_num,
|
|
const int64_t walk_num,
|
|
const int64_t* nucleus_shell_num,
|
|
const int64_t* nucleus_index,
|
|
const int64_t* shell_prim_index,
|
|
const int64_t* shell_prim_num,
|
|
const double* elec_coord,
|
|
const double* nucl_coord,
|
|
const double* expo,
|
|
const double* coef_normalized,
|
|
double* const shell_vgl );
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_ao_basis_shell_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_basis_shell_gaussian_vgl"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_ao_basis_shell_gaussian_vgl &
|
|
(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_normalized, &
|
|
shell_vgl) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: prim_num
|
|
integer (c_int64_t) , intent(in) , value :: shell_num
|
|
integer (c_int64_t) , intent(in) , value :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num)
|
|
integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num)
|
|
integer (c_int64_t) , intent(in) :: shell_prim_index(shell_num)
|
|
integer (c_int64_t) , intent(in) :: shell_prim_num(shell_num)
|
|
real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num)
|
|
real (c_double ) , intent(in) :: nucl_coord(elec_num,3)
|
|
real (c_double ) , intent(in) :: expo(prim_num)
|
|
real (c_double ) , intent(in) :: coef_normalized(prim_num)
|
|
real (c_double ) , intent(out) :: shell_vgl(shell_num,elec_num,walk_num,5)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_ao_basis_shell_gaussian_vgl_f
|
|
info = 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_normalized, &
|
|
shell_vgl)
|
|
|
|
end function qmckl_compute_ao_basis_shell_gaussian_vgl
|
|
#+end_src
|
|
|
|
#+begin_src python :results output :exports none
|
|
import numpy as np
|
|
|
|
def f(a,x,y):
|
|
return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )
|
|
|
|
def df(a,x,y,n):
|
|
h0 = 1.e-6
|
|
if n == 1: h = np.array([h0,0.,0.])
|
|
elif n == 2: h = np.array([0.,h0,0.])
|
|
elif n == 3: h = np.array([0.,0.,h0])
|
|
return ( f(a,x+h,y) - f(a,x-h,y) ) / (2.*h0)
|
|
|
|
def d2f(a,x,y,n):
|
|
h0 = 1.e-6
|
|
if n == 1: h = np.array([h0,0.,0.])
|
|
elif n == 2: h = np.array([0.,h0,0.])
|
|
elif n == 3: h = np.array([0.,0.,h0])
|
|
return ( f(a,x+h,y) - 2.*f(a,x,y) + f(a,x-h,y) ) / h0**2
|
|
|
|
def lf(a,x,y):
|
|
return d2f(a,x,y,1) + d2f(a,x,y,2) + d2f(a,x,y,3)
|
|
|
|
elec_26_w1 = np.array( [ 1.49050402641, 2.90106987953, -1.05920815468 ] )
|
|
elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] )
|
|
nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] )
|
|
nucl_2 = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] )
|
|
|
|
#double prim_vgl[prim_num][5][walk_num][elec_num];
|
|
x = elec_26_w1 ; y = nucl_1
|
|
a = [( 8.236000E+03, -1.130000E-04 * 6.1616545431994848e+02 ),
|
|
( 1.235000E+03, -8.780000E-04 * 1.4847738511079908e+02 ),
|
|
( 2.808000E+02, -4.540000E-03 * 4.8888635917437597e+01 ),
|
|
( 7.927000E+01, -1.813300E-02 * 1.8933972232608955e+01 ),
|
|
( 2.559000E+01, -5.576000E-02 * 8.1089160941724145e+00 ),
|
|
( 8.997000E+00, -1.268950E-01 * 3.7024003863155635e+00 ),
|
|
( 3.319000E+00, -1.703520E-01 * 1.7525302846177560e+00 ),
|
|
( 9.059000E-01, 1.403820E-01 * 6.6179013183966806e-01 ),
|
|
( 3.643000E-01, 5.986840E-01 * 3.3419848027174592e-01 ),
|
|
( 1.285000E-01, 3.953890E-01 * 1.5296336817449557e-01 )]
|
|
|
|
print ( "[1][0][0][26] : %25.15e"% f(a,x,y))
|
|
print ( "[1][1][0][26] : %25.15e"% df(a,x,y,1))
|
|
print ( "[1][2][0][26] : %25.15e"% df(a,x,y,2))
|
|
print ( "[1][3][0][26] : %25.15e"% df(a,x,y,3))
|
|
print ( "[1][4][0][26] : %25.15e"% lf(a,x,y))
|
|
|
|
x = elec_15_w2 ; y = nucl_2
|
|
a = [(3.387000E+01, 6.068000E-03 *1.0006253235944540e+01),
|
|
(5.095000E+00, 4.530800E-02 *2.4169531573445120e+00),
|
|
(1.159000E+00, 2.028220E-01 *7.9610924849766440e-01),
|
|
(3.258000E-01, 5.039030E-01 *3.0734305383061117e-01),
|
|
(1.027000E-01, 3.834210E-01 *1.2929684417481876e-01)]
|
|
|
|
print ( "[0][1][15][14] : %25.15e"% f(a,x,y))
|
|
print ( "[1][1][15][14] : %25.15e"% df(a,x,y,1))
|
|
print ( "[2][1][15][14] : %25.15e"% df(a,x,y,2))
|
|
print ( "[3][1][15][14] : %25.15e"% df(a,x,y,3))
|
|
print ( "[4][1][15][14] : %25.15e"% lf(a,x,y))
|
|
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
#+begin_example
|
|
[1][0][0][26] : 3.564393437193867e-02
|
|
[1][1][0][26] : -6.030177988891605e-03
|
|
[1][2][0][26] : -3.074832579871845e-02
|
|
[1][3][0][26] : 2.809546963133958e-02
|
|
[1][4][0][26] : 1.903338597841753e-02
|
|
[0][1][15][14] : 5.928089771361000e-03
|
|
[1][1][15][14] : 4.355862298893037e-03
|
|
[2][1][15][14] : 3.925108924950765e-03
|
|
[3][1][15][14] : 7.935527764416084e-04
|
|
[4][1][15][14] : 2.697495005143935e-03
|
|
#+end_example
|
|
|
|
*** Test
|
|
|
|
#+begin_src c :tangle (eval c_test) :exports none
|
|
{
|
|
#define walk_num chbrclf_walk_num
|
|
#define elec_num chbrclf_elec_num
|
|
#define shell_num chbrclf_shell_num
|
|
|
|
int64_t elec_up_num = chbrclf_elec_up_num;
|
|
int64_t elec_dn_num = chbrclf_elec_dn_num;
|
|
double* elec_coord = &(chbrclf_elec_coord[0][0][0]);
|
|
|
|
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_electron_walk_num (context, walk_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
assert(qmckl_electron_provided(context));
|
|
|
|
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
|
|
double shell_vgl[5][walk_num][elec_num][shell_num];
|
|
|
|
rc = qmckl_get_ao_basis_shell_vgl(context, &(shell_vgl[0][0][0][0]));
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
printf(" shell_vgl[1][0][0][26] %25.15e\n", shell_vgl[0][0][26][1]);
|
|
printf(" shell_vgl[1][1][0][26] %25.15e\n", shell_vgl[1][0][26][1]);
|
|
printf(" shell_vgl[1][2][0][26] %25.15e\n", shell_vgl[2][0][26][1]);
|
|
printf(" shell_vgl[1][3][0][26] %25.15e\n", shell_vgl[3][0][26][1]);
|
|
printf(" shell_vgl[1][4][0][26] %25.15e\n", shell_vgl[4][0][26][1]);
|
|
|
|
printf(" shell_vgl[14][0][1][15] %25.15e\n", shell_vgl[0][1][15][14]);
|
|
printf(" shell_vgl[14][1][1][15] %25.15e\n", shell_vgl[1][1][15][14]);
|
|
printf(" shell_vgl[14][2][1][15] %25.15e\n", shell_vgl[2][1][15][14]);
|
|
printf(" shell_vgl[14][3][1][15] %25.15e\n", shell_vgl[3][1][15][14]);
|
|
printf(" shell_vgl[14][4][1][15] %25.15e\n", shell_vgl[4][1][15][14]);
|
|
|
|
assert( fabs(shell_vgl[0][0][26][1] - ( 3.564393437193868e-02)) < 1.e-14 );
|
|
assert( fabs(shell_vgl[1][0][26][1] - (-6.030177987072189e-03)) < 1.e-14 );
|
|
assert( fabs(shell_vgl[2][0][26][1] - (-3.074832579537582e-02)) < 1.e-14 );
|
|
assert( fabs(shell_vgl[3][0][26][1] - ( 2.809546963519935e-02)) < 1.e-14 );
|
|
assert( fabs(shell_vgl[4][0][26][1] - ( 1.896046117183968e-02)) < 1.e-14 );
|
|
|
|
assert( fabs(shell_vgl[0][1][15][14] - ( 5.928089771361000e-03)) < 1.e-14 );
|
|
assert( fabs(shell_vgl[1][1][15][14] - ( 4.355862296021654e-03)) < 1.e-14 );
|
|
assert( fabs(shell_vgl[2][1][15][14] - ( 3.925108924923650e-03)) < 1.e-14 );
|
|
assert( fabs(shell_vgl[3][1][15][14] - ( 7.935527784022099e-04)) < 1.e-14 );
|
|
assert( fabs(shell_vgl[4][1][15][14] - ( 2.708246573703548e-03)) < 1.e-14 );
|
|
}
|
|
|
|
#+end_src
|
|
|
|
* Polynomial part
|
|
|
|
Going from the atomic basis set to AOs implies a systematic
|
|
construction of all the angular functions of each shell. We
|
|
consider two cases for the angular functions: the real-valued
|
|
spherical harmonics, and the polynomials in Cartesian coordinates.
|
|
In the case of spherical harmonics, the AOs are ordered in
|
|
increasing magnetic quantum number ($-l \le m \le l$), and in the
|
|
case of polynomials we choose the canonical ordering, i.e
|
|
|
|
\begin{eqnarray}
|
|
p & : & p_x, p_y, p_z \nonumber \\
|
|
d & : & d_{xx}, d_{xy}, d_{xz}, d_{yy}, d_{yz}, d_{zz} \nonumber \\
|
|
f & : & f_{xxx}, f_{xxy}, f_{xxz}, f_{xyy}, f_{xyz}, f_{xzz}, f_{yyy}, f_{yyz}, f_{yzz}, f_{zzz} \nonumber \\
|
|
{\rm etc.} \nonumber
|
|
\end{eqnarray}
|
|
|
|
** General functions for Powers of $x-X_i$
|
|
:PROPERTIES:
|
|
:Name: qmckl_ao_power
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
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 \]
|
|
|
|
#+NAME: qmckl_ao_power_args
|
|
| 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
|
|
|
|
#+CALL: generate_c_header(table=qmckl_ao_power_args,rettyp=get_value("CRetType"),fname="qmckl_ao_power")
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
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 );
|
|
#+end_src
|
|
|
|
*** Source
|
|
|
|
#+begin_src f90 :tangle (eval f)
|
|
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
|
|
#+end_src
|
|
|
|
*** C interface
|
|
#+CALL: generate_c_interface(table=qmckl_ao_power_args,rettyp=get_value("CRetType"),fname="qmckl_ao_power")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_ao_power &
|
|
(context, n, X, LMAX, P, ldp) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: n
|
|
real (c_double ) , intent(in) :: X(n)
|
|
integer (c_int32_t) , intent(in) :: LMAX(n)
|
|
real (c_double ) , intent(out) :: P(ldp,n)
|
|
integer (c_int64_t) , intent(in) , value :: ldp
|
|
|
|
integer(c_int32_t), external :: qmckl_ao_power_f
|
|
info = qmckl_ao_power_f &
|
|
(context, n, X, LMAX, P, ldp)
|
|
|
|
end function qmckl_ao_power
|
|
#+end_src
|
|
|
|
*** Fortran interface
|
|
|
|
#+CALL: generate_f_interface(table=qmckl_ao_power_args,rettyp=get_value("CRetType"),fname="qmckl_ao_power")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(c_int32_t) function qmckl_ao_power &
|
|
(context, n, X, LMAX, P, ldp) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: n
|
|
real (c_double ) , intent(in) :: X(n)
|
|
integer (c_int32_t) , intent(in) :: LMAX(n)
|
|
real (c_double ) , intent(out) :: P(ldp,n)
|
|
integer (c_int64_t) , intent(in) , value :: ldp
|
|
|
|
end function qmckl_ao_power
|
|
end interface
|
|
#+end_src
|
|
|
|
*** Test
|
|
|
|
#+begin_src f90 :tangle (eval f_test)
|
|
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)
|
|
|
|
#ifdef VFC_CI
|
|
! Multplying epsilon by 16 = 2^4 is equivalent to asking 4 significant digits
|
|
! less. This makes sense because we are adding noise with MCA so we can't be
|
|
! as strict on the accuracy target.
|
|
epsilon = epsilon * 16
|
|
#endif
|
|
|
|
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_check("ao"//C_NULL_CHAR, "power_2_2"//C_NULL_CHAR, &
|
|
DBLE(P(2,2)), DBLE(0), DBLE(epsilon))
|
|
|
|
if (test_qmckl_ao_power /= QMCKL_SUCCESS) return
|
|
|
|
test_qmckl_ao_power = QMCKL_FAILURE
|
|
|
|
#ifndef VFC_CI
|
|
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
|
|
#endif
|
|
|
|
test_qmckl_ao_power = QMCKL_SUCCESS
|
|
deallocate(X,P,LMAX)
|
|
end function test_qmckl_ao_power
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval c_test) :exports none
|
|
int test_qmckl_ao_power(qmckl_context context);
|
|
assert(0 == test_qmckl_ao_power(context));
|
|
#+end_src
|
|
|
|
** General functions for Value, Gradient and Laplacian of a polynomial
|
|
:PROPERTIES:
|
|
:Name: qmckl_ao_polynomial_vgl
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
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~.
|
|
|
|
#+NAME: qmckl_ao_polynomial_vgl_args
|
|
| 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
|
|
|
|
#+CALL: generate_c_header(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
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 );
|
|
#+end_src
|
|
|
|
*** Source
|
|
#+begin_src f90 :tangle (eval f)
|
|
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)
|
|
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
|
|
#+end_src
|
|
|
|
*** C interface
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_ao_polynomial_vgl &
|
|
(context, X, R, lmax, n, L, ldl, VGL, ldv) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double ) , intent(in) :: X(3)
|
|
real (c_double ) , intent(in) :: R(3)
|
|
integer (c_int32_t) , intent(in) , value :: lmax
|
|
integer (c_int64_t) , intent(inout) :: n
|
|
integer (c_int32_t) , intent(out) :: L(ldl,n)
|
|
integer (c_int64_t) , intent(in) , value :: ldl
|
|
real (c_double ) , intent(out) :: VGL(ldv,n)
|
|
integer (c_int64_t) , intent(in) , value :: ldv
|
|
|
|
integer(c_int32_t), external :: qmckl_ao_polynomial_vgl_f
|
|
info = qmckl_ao_polynomial_vgl_f &
|
|
(context, X, R, lmax, n, L, ldl, VGL, ldv)
|
|
|
|
end function qmckl_ao_polynomial_vgl
|
|
#+end_src
|
|
|
|
*** Fortran interface
|
|
|
|
#+CALL: generate_f_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(c_int32_t) function qmckl_ao_polynomial_vgl &
|
|
(context, X, R, lmax, n, L, ldl, VGL, ldv) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double ) , intent(in) :: X(3)
|
|
real (c_double ) , intent(in) :: R(3)
|
|
integer (c_int32_t) , intent(in) , value :: lmax
|
|
integer (c_int64_t) , intent(inout) :: n
|
|
integer (c_int32_t) , intent(out) :: L(ldl,n)
|
|
integer (c_int64_t) , intent(in) , value :: ldl
|
|
real (c_double ) , intent(out) :: VGL(ldv,n)
|
|
integer (c_int64_t) , intent(in) , value :: ldv
|
|
|
|
end function qmckl_ao_polynomial_vgl
|
|
end interface
|
|
#+end_src
|
|
|
|
*** Test
|
|
|
|
#+begin_src f90 :tangle (eval f_test)
|
|
integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
|
|
use qmckl
|
|
use qmckl_probes_f
|
|
|
|
implicit none
|
|
|
|
integer(c_int64_t), intent(in), value :: context
|
|
logical(C_BOOL) :: vfc_err
|
|
|
|
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_check("ao"//C_NULL_CHAR, "polynomial_vgl_1_2"//C_NULL_CHAR, &
|
|
DBLE(VGL(1,2)), DBLE(0), DBLE(epsilon))
|
|
|
|
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "polynomial_vgl_2_2"//C_NULL_CHAR, &
|
|
DBLE(VGL(2,2)), DBLE(0), DBLE(epsilon))
|
|
|
|
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "polynomial_vgl_3_2"//C_NULL_CHAR, &
|
|
DBLE(VGL(3,2)), DBLE(0), DBLE(epsilon))
|
|
|
|
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "polynomial_vgl_4_2"//C_NULL_CHAR, &
|
|
DBLE(VGL(4,2)), DBLE(0), DBLE(epsilon))
|
|
|
|
vfc_err = qmckl_probe_check("ao"//C_NULL_CHAR, "polynomial_vgl_5_2"//C_NULL_CHAR, &
|
|
DBLE(VGL(5,2)), DBLE(0), DBLE(epsilon))
|
|
|
|
if (test_qmckl_ao_polynomial_vgl /= QMCKL_SUCCESS) return
|
|
if (n /= d) return
|
|
|
|
#ifndef 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
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
int test_qmckl_ao_polynomial_vgl(qmckl_context context);
|
|
assert(0 == test_qmckl_ao_polynomial_vgl(context));
|
|
#+end_src
|
|
|
|
* Combining radial and polynomial parts
|
|
|
|
*** Get
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_ao_vgl(qmckl_context context, double* const ao_vgl);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_ao_vgl(qmckl_context context, double* const ao_vgl) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_ao_vgl(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = ctx->ao_basis.ao_num * 5 * ctx->electron.num * ctx->electron.walk_num;
|
|
memcpy(ao_vgl, ctx->ao_basis.ao_vgl, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_ao_vgl (context, ao_vgl) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
double precision, intent(out) :: ao_vgl(*)
|
|
end function
|
|
end interface
|
|
#+end_src
|
|
|
|
*** Provide
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
if (!ctx->ao_basis.provided) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_ao_vgl",
|
|
NULL);
|
|
}
|
|
|
|
if(!(ctx->electron.provided)) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_electron",
|
|
NULL);
|
|
}
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->electron.coord_new_date > ctx->ao_basis.ao_vgl_date) {
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
/* Provide required data */
|
|
rc = qmckl_provide_ao_basis_shell_vgl(context);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_shell_vgl", NULL);
|
|
}
|
|
|
|
/* Allocate array */
|
|
if (ctx->ao_basis.ao_vgl == NULL) {
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = ctx->ao_basis.ao_num * 5 * ctx->electron.num *
|
|
ctx->electron.walk_num * sizeof(double);
|
|
double* ao_vgl = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (ao_vgl == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_ao_basis_ao_vgl",
|
|
NULL);
|
|
}
|
|
ctx->ao_basis.ao_vgl = ao_vgl;
|
|
}
|
|
|
|
rc = qmckl_compute_ao_vgl(context,
|
|
ctx->ao_basis.ao_num,
|
|
ctx->ao_basis.shell_num,
|
|
ctx->electron.num,
|
|
ctx->nucleus.num,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.coord_new,
|
|
ctx->nucleus.coord,
|
|
ctx->ao_basis.nucleus_index,
|
|
ctx->ao_basis.nucleus_shell_num,
|
|
ctx->ao_basis.nucleus_range,
|
|
ctx->ao_basis.nucleus_max_ang_mom,
|
|
ctx->ao_basis.shell_ang_mom,
|
|
ctx->ao_basis.ao_factor,
|
|
ctx->ao_basis.shell_vgl,
|
|
ctx->ao_basis.ao_vgl);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->ao_basis.ao_vgl_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_ao_vgl
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_ao_vgl_args
|
|
| ~qmckl_context~ | ~context~ | in | Global state |
|
|
| ~int64_t~ | ~ao_num~ | in | Number of AOs |
|
|
| ~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 |
|
|
| ~double~ | ~elec_coord[walk_num][3][elec_num]~ | in | Electron coordinates |
|
|
| ~double~ | ~nucl_coord[3][nucl_num]~ | in | Nuclear coordinates |
|
|
| ~int64_t~ | ~nucleus_index[nucl_num]~ | in | Index of the 1st shell of each nucleus |
|
|
| ~int64_t~ | ~nucleus_shell_num[nucl_num]~ | in | Number of shells per nucleus |
|
|
| ~double~ | ~nucleus_range[nucl_num]~ | in | Range beyond which all is zero |
|
|
| ~int32_t~ | ~nucleus_max_ang_mom[nucl_num]~ | in | Maximum angular momentum per nucleus |
|
|
| ~int32_t~ | ~shell_ang_mom[shell_num]~ | in | Angular momentum of each shell |
|
|
| ~double~ | ~ao_factor[ao_num]~ | in | Normalization factor of the AOs |
|
|
| ~double~ | ~shell_vgl[5][walk_num][elec_num][shell_num]~ | in | Value, gradients and Laplacian of the shells |
|
|
| ~double~ | ~ao_vgl[5][walk_num][elec_num][ao_num]~ | out | Value, gradients and Laplacian of the AOs |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_ao_vgl_f(context, &
|
|
ao_num, shell_num, elec_num, nucl_num, walk_num, &
|
|
elec_coord, nucl_coord, nucleus_index, nucleus_shell_num, &
|
|
nucleus_range, nucleus_max_ang_mom, shell_ang_mom, &
|
|
ao_factor, shell_vgl, ao_vgl) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: ao_num
|
|
integer*8 , intent(in) :: shell_num
|
|
integer*8 , intent(in) :: elec_num
|
|
integer*8 , intent(in) :: nucl_num
|
|
integer*8 , intent(in) :: walk_num
|
|
double precision , intent(in) :: elec_coord(elec_num,3,walk_num)
|
|
double precision , intent(in) :: nucl_coord(nucl_num,3)
|
|
integer*8 , intent(in) :: nucleus_index(nucl_num)
|
|
integer*8 , intent(in) :: nucleus_shell_num(nucl_num)
|
|
double precision , intent(in) :: nucleus_range(nucl_num)
|
|
integer , intent(in) :: nucleus_max_ang_mom(nucl_num)
|
|
integer , intent(in) :: shell_ang_mom(shell_num)
|
|
double precision , intent(in) :: ao_factor(ao_num)
|
|
double precision , intent(in) :: shell_vgl(shell_num,elec_num,walk_num,5)
|
|
double precision , intent(out) :: ao_vgl(ao_num,elec_num,walk_num,5)
|
|
|
|
double precision :: e_coord(3), n_coord(3)
|
|
integer*8 :: n_poly
|
|
integer :: l, il, k
|
|
integer*8 :: ielec, inucl, ishell, iwalk
|
|
integer :: lstart(0:20)
|
|
double precision :: x, y, z, r2
|
|
double precision :: cutoff
|
|
integer, external :: qmckl_ao_polynomial_vgl_f
|
|
|
|
double precision, allocatable :: poly_vgl(:,:)
|
|
integer , allocatable :: powers(:,:)
|
|
|
|
allocate(poly_vgl(5,ao_num), powers(3,ao_num))
|
|
|
|
! Pre-computed data
|
|
do l=0,20
|
|
lstart(l) = l*(l+1)*(l+2)/6 +1
|
|
end do
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
! Don't compute polynomials when the radial part is zero.
|
|
! TODO : Use numerical precision here
|
|
cutoff = -dlog(1.d-15)
|
|
|
|
do iwalk = 1,walk_num
|
|
do ielec = 1, elec_num
|
|
e_coord(1) = elec_coord(ielec,1,iwalk)
|
|
e_coord(2) = elec_coord(ielec,2,iwalk)
|
|
e_coord(3) = elec_coord(ielec,3,iwalk)
|
|
k=1
|
|
do inucl=1,nucl_num
|
|
n_coord(1) = nucl_coord(inucl,1)
|
|
n_coord(2) = nucl_coord(inucl,2)
|
|
n_coord(3) = nucl_coord(inucl,3)
|
|
|
|
! Test if the electron is in the range of the nucleus
|
|
x = e_coord(1) - n_coord(1)
|
|
y = e_coord(2) - n_coord(2)
|
|
z = e_coord(3) - n_coord(3)
|
|
|
|
r2 = x*x + z*z + z*z
|
|
|
|
if (r2 > cutoff*nucleus_range(inucl)) then
|
|
cycle
|
|
end if
|
|
|
|
! Compute polynomials
|
|
info = qmckl_ao_polynomial_vgl_f(context, e_coord, n_coord, &
|
|
nucleus_max_ang_mom(inucl), n_poly, powers, 3_8, &
|
|
poly_vgl, 5_8)
|
|
|
|
! Loop over shells
|
|
do ishell = nucleus_index(inucl)+1, nucleus_index(inucl)+nucleus_shell_num(inucl)
|
|
l = shell_ang_mom(ishell)
|
|
do il = lstart(l), lstart(l+1)-1
|
|
! Value
|
|
ao_vgl(k,ielec,iwalk,1) = &
|
|
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,1) * ao_factor(k)
|
|
|
|
! Grad_x
|
|
ao_vgl(k,ielec,iwalk,2) = ( &
|
|
poly_vgl(2,il) * shell_vgl(ishell,ielec,iwalk,1) + &
|
|
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,2) &
|
|
) * ao_factor(k)
|
|
|
|
! Grad_y
|
|
ao_vgl(k,ielec,iwalk,3) = ( &
|
|
poly_vgl(3,il) * shell_vgl(ishell,ielec,iwalk,1) + &
|
|
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,3) &
|
|
) * ao_factor(k)
|
|
|
|
! Grad_z
|
|
ao_vgl(k,ielec,iwalk,4) = ( &
|
|
poly_vgl(4,il) * shell_vgl(ishell,ielec,iwalk,1) + &
|
|
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,4) &
|
|
) * ao_factor(k)
|
|
|
|
! Lapl_z
|
|
ao_vgl(k,ielec,iwalk,5) = ( &
|
|
poly_vgl(5,il) * shell_vgl(ishell,ielec,iwalk,1) + &
|
|
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,5) + &
|
|
2.d0 * ( &
|
|
poly_vgl(2,il) * shell_vgl(ishell,ielec,iwalk,2) + &
|
|
poly_vgl(3,il) * shell_vgl(ishell,ielec,iwalk,3) + &
|
|
poly_vgl(4,il) * shell_vgl(ishell,ielec,iwalk,4) ) &
|
|
) * ao_factor(k)
|
|
|
|
k = k+1
|
|
end do
|
|
end do
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
deallocate(poly_vgl, powers)
|
|
end function qmckl_compute_ao_vgl_f
|
|
#+end_src
|
|
|
|
# #+CALL: generate_c_header(table=qmckl_ao_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_private_func) :comments org
|
|
qmckl_exit_code qmckl_compute_ao_vgl (
|
|
const qmckl_context context,
|
|
const int64_t ao_num,
|
|
const int64_t shell_num,
|
|
const int64_t elec_num,
|
|
const int64_t nucl_num,
|
|
const int64_t walk_num,
|
|
const double* elec_coord,
|
|
const double* nucl_coord,
|
|
const int64_t* nucleus_index,
|
|
const int64_t* nucleus_shell_num,
|
|
const double* nucleus_range,
|
|
const int32_t* nucleus_max_ang_mom,
|
|
const int32_t* shell_ang_mom,
|
|
const double* ao_factor,
|
|
const double* shell_vgl,
|
|
double* const ao_vgl );
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_ao_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_ao_vgl &
|
|
(context, &
|
|
ao_num, &
|
|
shell_num, &
|
|
elec_num, &
|
|
nucl_num, &
|
|
walk_num, &
|
|
elec_coord, &
|
|
nucl_coord, &
|
|
nucleus_index, &
|
|
nucleus_shell_num, &
|
|
nucleus_range, &
|
|
nucleus_max_ang_mom, &
|
|
shell_ang_mom, &
|
|
ao_factor, &
|
|
shell_vgl, &
|
|
ao_vgl) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: ao_num
|
|
integer (c_int64_t) , intent(in) , value :: shell_num
|
|
integer (c_int64_t) , intent(in) , value :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num)
|
|
real (c_double ) , intent(in) :: nucl_coord(nucl_num,3)
|
|
integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num)
|
|
integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num)
|
|
real (c_double ) , intent(in) :: nucleus_range(nucl_num)
|
|
integer (c_int32_t) , intent(in) :: nucleus_max_ang_mom(nucl_num)
|
|
integer (c_int32_t) , intent(in) :: shell_ang_mom(shell_num)
|
|
real (c_double ) , intent(in) :: ao_factor(ao_num)
|
|
real (c_double ) , intent(in) :: shell_vgl(shell_num,elec_num,walk_num,5)
|
|
real (c_double ) , intent(out) :: ao_vgl(ao_num,elec_num,walk_num,5)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_ao_vgl_f
|
|
info = qmckl_compute_ao_vgl_f &
|
|
(context, &
|
|
ao_num, &
|
|
shell_num, &
|
|
elec_num, &
|
|
nucl_num, &
|
|
walk_num, &
|
|
elec_coord, &
|
|
nucl_coord, &
|
|
nucleus_index, &
|
|
nucleus_shell_num, &
|
|
nucleus_range, &
|
|
nucleus_max_ang_mom, &
|
|
shell_ang_mom, &
|
|
ao_factor, &
|
|
shell_vgl, &
|
|
ao_vgl)
|
|
|
|
end function qmckl_compute_ao_vgl
|
|
#+end_src
|
|
|
|
#+begin_src python :results output :exports none
|
|
import numpy as np
|
|
|
|
def f(a,x,y):
|
|
return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )
|
|
|
|
def df(a,x,y,n):
|
|
h0 = 1.e-6
|
|
if n == 1: h = np.array([h0,0.,0.])
|
|
elif n == 2: h = np.array([0.,h0,0.])
|
|
elif n == 3: h = np.array([0.,0.,h0])
|
|
return ( f(a,x+h,y) - f(a,x-h,y) ) / (2.*h0)
|
|
|
|
def d2f(a,x,y,n):
|
|
h0 = 1.e-6
|
|
if n == 1: h = np.array([h0,0.,0.])
|
|
elif n == 2: h = np.array([0.,h0,0.])
|
|
elif n == 3: h = np.array([0.,0.,h0])
|
|
return ( f(a,x+h,y) - 2.*f(a,x,y) + f(a,x-h,y) ) / h0**2
|
|
|
|
def lf(a,x,y):
|
|
return d2f(a,x,y,1) + d2f(a,x,y,2) + d2f(a,x,y,3)
|
|
|
|
elec_26_w1 = np.array( [ 1.49050402641, 2.90106987953, -1.05920815468 ] )
|
|
elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] )
|
|
nucl_1 = np.array( [ -2.302574592081335e+00, -3.542027060505035e-01, -5.334129934317614e-02] )
|
|
|
|
#double prim_vgl[prim_num][5][walk_num][elec_num];
|
|
x = elec_26_w1 ; y = nucl_1
|
|
a = [( 403.830000, 0.001473 * 5.9876577632594533e+04),
|
|
( 121.170000, 0.012672 * 7.2836806319891484e+03),
|
|
( 46.345000, 0.058045 * 1.3549226646722386e+03),
|
|
( 19.721000, 0.170510 * 3.0376315094739988e+02),
|
|
( 8.862400, 0.318596 * 7.4924579607137730e+01),
|
|
( 3.996200, 0.384502 * 1.8590543353806009e+01),
|
|
( 1.763600, 0.273774 * 4.4423176930919421e+00),
|
|
( 0.706190, 0.074397 * 8.9541051939952665e-01)]
|
|
|
|
print ( "[0][0][26][219] : %25.15e"%(f(a,x,y) * (x[0] - y[0])**2) )
|
|
print ( "[1][0][26][219] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + 2.*f(a,x,y) * (x[0] - y[0])) )
|
|
|
|
print ( "[0][0][26][220] : %25.15e"%(f(a,x,y) * (x[0] - y[0]) * (x[1] - y[1])) )
|
|
print ( "[1][0][26][220] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + f(a,x,y) * (x[1] - y[1])) )
|
|
|
|
print ( "[0][0][26][221] : %25.15e"%(f(a,x,y) * (x[0] - y[0]) * (x[2] - y[2])) )
|
|
print ( "[1][0][26][221] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[2] - y[2]) + f(a,x,y) * (x[2] - y[2])) )
|
|
|
|
print ( "[0][0][26][222] : %25.15e"%(f(a,x,y) * (x[1] - y[1]) * (x[1] - y[1])) )
|
|
print ( "[1][0][26][222] : %25.15e"%(df(a,x,y,1)* (x[1] - y[1]) * (x[1] - y[1])) )
|
|
|
|
print ( "[0][0][26][223] : %25.15e"%(f(a,x,y) * (x[1] - y[1]) * (x[2] - y[2])) )
|
|
print ( "[1][0][26][223] : %25.15e"%(df(a,x,y,1)* (x[1] - y[1]) * (x[2] - y[2])) )
|
|
|
|
print ( "[0][0][26][224] : %25.15e"%(f(a,x,y) * (x[2] - y[2]) * (x[2] - y[2])) )
|
|
print ( "[1][0][26][224] : %25.15e"%(df(a,x,y,1)* (x[2] - y[2]) * (x[2] - y[2])) )
|
|
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
#+begin_example
|
|
[0][0][26][219] : 1.020302912653649e-08
|
|
[1][0][26][219] : -4.153046808203204e-08
|
|
[0][0][26][220] : 8.756380857379661e-09
|
|
[1][0][26][220] : -4.460176677299534e-08
|
|
[0][0][26][221] : -2.705688401075445e-09
|
|
[1][0][26][221] : 1.378177639720419e-08
|
|
[0][0][26][222] : 7.514847283937212e-09
|
|
[1][0][26][222] : -4.025905373647693e-08
|
|
[0][0][26][223] : -2.322059246071533e-09
|
|
[1][0][26][223] : 1.243989457599443e-08
|
|
[0][0][26][224] : 7.175074806631758e-10
|
|
[1][0][26][224] : -3.843880138733679e-09
|
|
#+end_example
|
|
|
|
*** Test
|
|
|
|
#+begin_src c :tangle (eval c_test) :exports none
|
|
{
|
|
#define walk_num chbrclf_walk_num
|
|
#define elec_num chbrclf_elec_num
|
|
#define shell_num chbrclf_shell_num
|
|
#define ao_num chbrclf_ao_num
|
|
|
|
int64_t elec_up_num = chbrclf_elec_up_num;
|
|
int64_t elec_dn_num = chbrclf_elec_dn_num;
|
|
double* elec_coord = &(chbrclf_elec_coord[0][0][0]);
|
|
|
|
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_electron_walk_num (context, walk_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
assert(qmckl_electron_provided(context));
|
|
|
|
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
|
|
double ao_vgl[5][walk_num][elec_num][ao_num];
|
|
|
|
rc = qmckl_get_ao_vgl(context, &(ao_vgl[0][0][0][0]));
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
printf("\n");
|
|
printf(" ao_vgl ao_vgl[0][0][26][219] %25.15e\n", ao_vgl[0][0][26][219]);
|
|
printf(" ao_vgl ao_vgl[1][0][26][219] %25.15e\n", ao_vgl[1][0][26][219]);
|
|
printf(" ao_vgl ao_vgl[0][0][26][220] %25.15e\n", ao_vgl[0][0][26][220]);
|
|
printf(" ao_vgl ao_vgl[1][0][26][220] %25.15e\n", ao_vgl[1][0][26][220]);
|
|
printf(" ao_vgl ao_vgl[0][0][26][221] %25.15e\n", ao_vgl[0][0][26][221]);
|
|
printf(" ao_vgl ao_vgl[1][0][26][221] %25.15e\n", ao_vgl[1][0][26][221]);
|
|
printf(" ao_vgl ao_vgl[0][0][26][222] %25.15e\n", ao_vgl[0][0][26][222]);
|
|
printf(" ao_vgl ao_vgl[1][0][26][222] %25.15e\n", ao_vgl[1][0][26][222]);
|
|
printf(" ao_vgl ao_vgl[0][0][26][223] %25.15e\n", ao_vgl[0][0][26][223]);
|
|
printf(" ao_vgl ao_vgl[1][0][26][223] %25.15e\n", ao_vgl[1][0][26][223]);
|
|
printf(" ao_vgl ao_vgl[0][0][26][224] %25.15e\n", ao_vgl[0][0][26][224]);
|
|
printf(" ao_vgl ao_vgl[1][0][26][224] %25.15e\n", ao_vgl[1][0][26][224]);
|
|
printf("\n");
|
|
|
|
assert( fabs(ao_vgl[0][0][26][219] - ( 1.020298798341620e-08)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[1][0][26][219] - ( -4.928035238010602e-08)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[0][0][26][220] - ( 8.756345547784206e-09)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[1][0][26][220] - ( -4.460158690983819e-08)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[0][0][26][221] - ( -2.705677490544664e-09)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[1][0][26][221] - ( 1.378172082017231e-08)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[0][0][26][222] - ( 7.514816980753531e-09)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[1][0][26][222] - ( -4.025889138635182e-08)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[0][0][26][223] - ( -2.322049882502961e-09)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[1][0][26][223] - ( 1.243984441042288e-08)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[0][0][26][224] - ( 7.175045873560788e-10)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[1][0][26][224] - ( -3.843864637762753e-09)) < 1.e-14 );
|
|
}
|
|
|
|
#+end_src
|
|
|
|
* End of files :noexport:
|
|
|
|
#+begin_src c :tangle (eval h_private_type)
|
|
#endif
|
|
#+end_src
|
|
|
|
*** Test
|
|
#+begin_src c :tangle (eval c_test)
|
|
rc = qmckl_context_destroy(context);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
return 0;
|
|
}
|
|
#+end_src
|
|
|
|
**✸ Compute file names
|
|
#+begin_src emacs-lisp
|
|
; The following is required to compute the file names
|
|
|
|
(setq pwd (file-name-directory buffer-file-name))
|
|
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
|
|
(setq f (concat pwd name "_f.f90"))
|
|
(setq fh (concat pwd name "_fh.f90"))
|
|
(setq c (concat pwd name ".c"))
|
|
(setq h (concat name ".h"))
|
|
(setq h_private (concat name "_private.h"))
|
|
(setq c_test (concat pwd "test_" name ".c"))
|
|
(setq f_test (concat pwd "test_" name "_f.f90"))
|
|
|
|
; Minted
|
|
(require 'ox-latex)
|
|
(setq org-latex-listings 'minted)
|
|
(add-to-list 'org-latex-packages-alist '("" "listings"))
|
|
(add-to-list 'org-latex-packages-alist '("" "color"))
|
|
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
| | color |
|
|
| | listings |
|
|
|
|
|
|
# -*- mode: org -*-
|
|
# vim: syntax=c
|