mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-08-19 03:41:44 +02:00
1713 lines
50 KiB
Org Mode
1713 lines
50 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}
|
|
\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. As this normalization requires
|
|
the ability to compute overlap integrals, it should be written in the
|
|
file to ensure that the file is self-contained and does not require
|
|
the client program to have the ability to compute such integrals.
|
|
|
|
Atomic orbitals (AOs) are defined as
|
|
|
|
\[
|
|
\chi_i (\mathbf{r}) = 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.
|
|
|
|
In this section we describe 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
|
|
|
|
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"
|
|
#+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 |
|
|
| ~shell_center~ | ~[shell_num]~ | Id of the nucleus on which each shell is centered |
|
|
| ~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 |
|
|
|
|
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
|
|
shell_center = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2]
|
|
shell_ang_mom = ['S', 'S', 'S', 'P', 'P', 'D', 'S', 'S', 'S', 'P', 'P', 'D']
|
|
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 = [1, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 20]
|
|
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]
|
|
#+END_EXAMPLE
|
|
|
|
** Data structure
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_type)
|
|
typedef struct qmckl_ao_basis_struct {
|
|
int32_t uninitialized;
|
|
int64_t shell_num;
|
|
int64_t prim_num;
|
|
int64_t * shell_center;
|
|
char * shell_ang_mom;
|
|
int64_t * shell_prim_num;
|
|
int64_t * shell_prim_index;
|
|
double * shell_factor;
|
|
double * exponent ;
|
|
double * coefficient ;
|
|
bool provided;
|
|
char type;
|
|
} qmckl_ao_basis_struct;
|
|
#+end_src
|
|
|
|
The ~uninitialized~ integer contains one bit set to one for each
|
|
initialization function which has not bee called. It becomes equal
|
|
to zero after all initialization functions have been called. The
|
|
struct is then initialized and ~provided == true~.
|
|
|
|
** 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_shell_center (const qmckl_context context);
|
|
char* 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);
|
|
#+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_shell_center (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.shell_center != NULL);
|
|
return ctx->ao_basis.shell_center;
|
|
}
|
|
|
|
|
|
char* 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 << 4;
|
|
|
|
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 << 5;
|
|
|
|
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 << 6;
|
|
|
|
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 << 7;
|
|
|
|
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 << 8;
|
|
|
|
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 << 9;
|
|
|
|
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
|
|
return NULL;
|
|
}
|
|
|
|
assert (ctx->ao_basis.coefficient != NULL);
|
|
return ctx->ao_basis.coefficient;
|
|
}
|
|
|
|
|
|
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. When
|
|
|
|
#+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_shell_prim_index (qmckl_context context, const int64_t * shell_prim_index);
|
|
qmckl_exit_code qmckl_set_ao_basis_shell_center (qmckl_context context, const int64_t * shell_center);
|
|
qmckl_exit_code qmckl_set_ao_basis_shell_ang_mom (qmckl_context context, const char * 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_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);
|
|
#+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);
|
|
|
|
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_shell_center(qmckl_context context, const int64_t* shell_center) {
|
|
<<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_shell_center",
|
|
"shell_num is not set");
|
|
}
|
|
|
|
if (ctx->ao_basis.shell_center != NULL) {
|
|
qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_center);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc,
|
|
"qmckl_set_ao_basis_shell_center",
|
|
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_center",
|
|
NULL);
|
|
}
|
|
|
|
memcpy(new_array, shell_center, mem_info.size);
|
|
|
|
ctx->ao_basis.shell_center = new_array;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
|
|
qmckl_exit_code qmckl_set_ao_basis_shell_ang_mom(qmckl_context context, const char* shell_ang_mom) {
|
|
<<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_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(char);
|
|
char* new_array = (char*) 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 << 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_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 << 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_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 << 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_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 << 8;
|
|
|
|
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 << 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_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>>
|
|
}
|
|
|
|
#+end_src
|
|
|
|
** TODO Fortran interfaces
|
|
|
|
** Test :noexport:
|
|
|
|
#+begin_src c :tangle (eval c_test) :exports none :exports none
|
|
/* Reference input data */
|
|
|
|
char typ = 'G';
|
|
#define shell_num ((int64_t) 12)
|
|
#define prim_num ((int64_t) 20)
|
|
|
|
int64_t shell_center [shell_num] =
|
|
{ 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2 };
|
|
|
|
char shell_ang_mom [shell_num] =
|
|
{ 'S', 'S', 'S', 'P', 'P', 'D', 'S', 'S', 'S', 'P', 'P', 'D' };
|
|
|
|
double shell_factor [shell_num] =
|
|
{ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1. };
|
|
|
|
int64_t shell_prim_num [shell_num] =
|
|
{5, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1};
|
|
|
|
int64_t shell_prim_index [shell_num] =
|
|
{1, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 20};
|
|
|
|
double exponent [prim_num] =
|
|
{ 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 };
|
|
|
|
double coefficient [prim_num] =
|
|
{ 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 };
|
|
/* --- */
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
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_shell_center (context, shell_center);
|
|
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_center (context, shell_prim_num);
|
|
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));
|
|
|
|
#+end_src
|
|
|
|
* Polynomial part
|
|
** 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
|
|
implicit none
|
|
|
|
integer(qmckl_context), intent(in), value :: context
|
|
|
|
integer*8 :: n, LDP
|
|
integer, allocatable :: LMAX(:)
|
|
double precision, allocatable :: X(:), P(:,:)
|
|
integer*8 :: i,j
|
|
double precision :: epsilon
|
|
|
|
epsilon = qmckl_get_numprec_epsilon(context)
|
|
|
|
n = 100;
|
|
LDP = 10;
|
|
|
|
allocate(X(n), P(LDP,n), LMAX(n))
|
|
|
|
do j=1,n
|
|
X(j) = -5.d0 + 0.1d0 * dble(j)
|
|
LMAX(j) = 1 + int(mod(j, 5),4)
|
|
end do
|
|
|
|
test_qmckl_ao_power = qmckl_ao_power(context, n, X, LMAX, P, LDP)
|
|
if (test_qmckl_ao_power /= QMCKL_SUCCESS) return
|
|
|
|
test_qmckl_ao_power = QMCKL_FAILURE
|
|
|
|
do j=1,n
|
|
do i=1,LMAX(j)
|
|
if ( X(j)**i == 0.d0 ) then
|
|
if ( P(i,j) /= 0.d0) return
|
|
else
|
|
if ( dabs(1.d0 - P(i,j) / (X(j)**i)) > epsilon ) return
|
|
end if
|
|
end do
|
|
end do
|
|
|
|
test_qmckl_ao_power = QMCKL_SUCCESS
|
|
deallocate(X,P,LMAX)
|
|
end function test_qmckl_ao_power
|
|
#+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
|
|
|
|
** 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)
|
|
integer, external :: qmckl_ao_power_f
|
|
double precision :: xy, yz, xz
|
|
double precision :: da, db, dc, dd
|
|
|
|
info = 0
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (lmax < 0) then
|
|
info = QMCKL_INVALID_ARG_4
|
|
return
|
|
endif
|
|
|
|
if (ldl < 3) then
|
|
info = QMCKL_INVALID_ARG_7
|
|
return
|
|
endif
|
|
|
|
if (ldv < 5) then
|
|
info = QMCKL_INVALID_ARG_9
|
|
return
|
|
endif
|
|
|
|
|
|
do i=1,3
|
|
Y(i) = X(i) - R(i)
|
|
end do
|
|
|
|
lmax_array(1:3) = lmax
|
|
if (lmax == 0) then
|
|
VGL(1,1) = 1.d0
|
|
vgL(2:5,1) = 0.d0
|
|
l(1:3,1) = 0
|
|
n=1
|
|
else if (lmax > 0) then
|
|
pows(-2:0,1:3) = 1.d0
|
|
do i=1,lmax
|
|
pows(i,1) = pows(i-1,1) * Y(1)
|
|
pows(i,2) = pows(i-1,2) * Y(2)
|
|
pows(i,3) = pows(i-1,3) * Y(3)
|
|
end do
|
|
|
|
VGL(1:5,1:4) = 0.d0
|
|
l (1:3,1:4) = 0
|
|
|
|
VGL(1 ,1 ) = 1.d0
|
|
vgl(1:5,2:4) = 0.d0
|
|
|
|
l (1,2) = 1
|
|
vgl(1,2) = pows(1,1)
|
|
vgL(2,2) = 1.d0
|
|
|
|
l (2,3) = 1
|
|
vgl(1,3) = pows(1,2)
|
|
vgL(3,3) = 1.d0
|
|
|
|
l (3,4) = 1
|
|
vgl(1,4) = pows(1,3)
|
|
vgL(4,4) = 1.d0
|
|
|
|
n=4
|
|
endif
|
|
|
|
! l>=2
|
|
dd = 2.d0
|
|
do d=2,lmax
|
|
da = dd
|
|
do a=d,0,-1
|
|
db = dd-da
|
|
do b=d-a,0,-1
|
|
c = d - a - b
|
|
dc = dd - da - db
|
|
n = n+1
|
|
|
|
l(1,n) = a
|
|
l(2,n) = b
|
|
l(3,n) = c
|
|
|
|
xy = pows(a,1) * pows(b,2)
|
|
yz = pows(b,2) * pows(c,3)
|
|
xz = pows(a,1) * pows(c,3)
|
|
|
|
vgl(1,n) = xy * pows(c,3)
|
|
|
|
xy = dc * xy
|
|
xz = db * xz
|
|
yz = da * yz
|
|
|
|
vgl(2,n) = pows(a-1,1) * yz
|
|
vgl(3,n) = pows(b-1,2) * xz
|
|
vgl(4,n) = pows(c-1,3) * xy
|
|
|
|
vgl(5,n) = &
|
|
(da-1.d0) * pows(a-2,1) * yz + &
|
|
(db-1.d0) * pows(b-2,2) * xz + &
|
|
(dc-1.d0) * pows(c-2,3) * xy
|
|
|
|
db = db - 1.d0
|
|
end do
|
|
da = da - 1.d0
|
|
end do
|
|
dd = dd + 1.d0
|
|
end do
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
end function qmckl_ao_polynomial_vgl_f
|
|
#+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
|
|
implicit none
|
|
|
|
integer(c_int64_t), intent(in), value :: context
|
|
|
|
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)
|
|
|
|
if (test_qmckl_ao_polynomial_vgl /= QMCKL_SUCCESS) return
|
|
if (n /= d) return
|
|
|
|
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
|
|
|
|
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
|
|
|
|
* Radial part
|
|
** 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
|
|
implicit none
|
|
|
|
integer(c_int64_t), intent(in), value :: context
|
|
|
|
integer*8 :: n, ldv, j, i
|
|
double precision :: X(3), R(3), Y(3), r2
|
|
double precision, allocatable :: VGL(:,:), A(:)
|
|
double precision :: epsilon
|
|
|
|
epsilon = qmckl_get_numprec_epsilon(context)
|
|
|
|
X = (/ 1.1 , 2.2 , 3.3 /)
|
|
R = (/ 0.1 , 1.2 , -2.3 /)
|
|
Y(:) = X(:) - R(:)
|
|
r2 = Y(1)**2 + Y(2)**2 + Y(3)**2
|
|
|
|
n = 10;
|
|
ldv = 100;
|
|
|
|
allocate (A(n), VGL(ldv,5))
|
|
do i=1,n
|
|
A(i) = 0.0013 * dble(ishft(1,i))
|
|
end do
|
|
|
|
|
|
test_qmckl_ao_gaussian_vgl = &
|
|
qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv)
|
|
if (test_qmckl_ao_gaussian_vgl /= 0) return
|
|
|
|
test_qmckl_ao_gaussian_vgl = -1
|
|
|
|
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
|
|
|
|
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 Slater basis functions
|
|
|
|
** TODO Radial functions on a grid
|
|
* Combining radial and polynomial parts
|
|
* 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
|
|
|
|
|