mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-12-23 12:57:20 +01:00
6106 lines
207 KiB
Org Mode
6106 lines
207 KiB
Org Mode
B
|
|
#+TITLE: Atomic Orbitals
|
|
#+SETUPFILE: ../tools/theme.setup
|
|
#+INCLUDE: ../tools/lib.org
|
|
|
|
* Introduction
|
|
|
|
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
|
|
/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 (integrate) 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 and $P$ are the
|
|
generating functions of the given angular momentum $\eta(i)$. 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_func)
|
|
#ifndef QMCKL_AO_HPF
|
|
#define QMCKL_AO_HPF
|
|
|
|
#+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 f90 :tangle (eval f) :noweb yes
|
|
#ifdef HAVE_CONFIG_H
|
|
#include "config.h"
|
|
#endif
|
|
#+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 <string.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 <stdio.h>
|
|
#include <math.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
|
|
|
|
** Constant data
|
|
|
|
The following arrays are stored in the context, and need to be set
|
|
when initializing the library:
|
|
|
|
#+NAME: constant_data
|
|
|---------------------+----------------------+----------------------------------------------------------------------|
|
|
| Variable | Type | Description |
|
|
|---------------------+----------------------+----------------------------------------------------------------------|
|
|
| ~type~ | ~char~ | Gaussian (~'G'~) or Slater (~'S'~) |
|
|
| ~shell_num~ | ~int64_t~ | Number of shells |
|
|
| ~prim_num~ | ~int64_t~ | Total number of primitives |
|
|
| ~nucleus_index~ | ~int64_t[nucl_num]~ | Index of the first shell of each nucleus |
|
|
| ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | Number of shells per nucleus |
|
|
| ~shell_ang_mom~ | ~int32_t[shell_num]~ | Angular momentum of each shell |
|
|
| ~shell_prim_num~ | ~int64_t[shell_num]~ | Number of primitives in each shell |
|
|
| ~shell_prim_index~ | ~int64_t[shell_num]~ | Address of the first primitive of each shell in the ~EXPONENT~ array |
|
|
| ~shell_factor~ | ~double[shell_num]~ | Normalization factor for each shell |
|
|
| ~exponent~ | ~double[prim_num]~ | Array of exponents |
|
|
| ~coefficient~ | ~double[prim_num]~ | Array of coefficients |
|
|
| ~prim_factor~ | ~double[prim_num]~ | Normalization factors of the primtives |
|
|
| ~ao_num~ | ~int64_t~ | Number of AOs |
|
|
| ~ao_cartesian~ | ~bool~ | If true, use polynomials. Otherwise, use spherical harmonics |
|
|
| ~ao_factor~ | ~double[ao_num]~ | Normalization factor of the AO |
|
|
|
|
For H_2 with the following basis set,
|
|
|
|
#+NAME: basis
|
|
#+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:
|
|
|
|
#+NAME: params
|
|
#+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
|
|
|
|
A scalar variable ~$V$~ present in this table can be set or get by
|
|
calling the functions:
|
|
|
|
#+NAME: template_scalar_c
|
|
#+begin_src C :tangle no
|
|
qmckl_exit_code qmckl_set_ao_basis_$V$ ( qmckl_context context,
|
|
const $type_of_V$ $V$);
|
|
|
|
qmckl_exit_code qmckl_get_ao_basis_$V$ ( const qmckl_context context,
|
|
$type_of_V$ const $V$);
|
|
#+end_src
|
|
|
|
#+begin_src f90 :tangle no
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_$V$ (context, $V$) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
$f_type_of_V$ , intent(in) , value :: $V$
|
|
end function qmckl_set_ao_basis_$V$
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_ao_basis_$V$ (context, $V$) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
$f_type_of_V$ , intent(out) :: $V$
|
|
end function qmckl_get_ao_basis_$V$
|
|
end interface
|
|
#+end_src
|
|
For array variables, use the rule:
|
|
|
|
#+NAME: template_array_c
|
|
#+begin_src C :tangle no
|
|
qmckl_exit_code qmckl_set_ao_basis_$V$ ( qmckl_context context,
|
|
const $type_of_V$ $V$,
|
|
const int64_t size_max);
|
|
|
|
qmckl_exit_code qmckl_get_ao_basis_$V$ ( const qmckl_context context,
|
|
$type_of_V$ const $V$,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src f90 :tangle no
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_$V$ (context, &
|
|
$V$, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
$f_type_of_V$ , intent(in) , value :: $V$
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_set_ao_basis_$V$
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_ao_basis_$V$ (context, &
|
|
$V$, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
$f_type_of_V$ , intent(out) :: $V$
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_get_ao_basis_$V$
|
|
end interface
|
|
#+end_src
|
|
|
|
*** Data structure :noexport:
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_type) :noweb yes
|
|
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;
|
|
uint64_t primitive_vgl_date;
|
|
double * shell_vgl;
|
|
uint64_t shell_vgl_date;
|
|
double * ao_vgl;
|
|
uint64_t ao_vgl_date;
|
|
|
|
int32_t uninitialized;
|
|
bool provided;
|
|
bool ao_cartesian;
|
|
char type;
|
|
#ifdef HAVE_HPC
|
|
<<HPC_struct>>
|
|
#endif
|
|
} 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 qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_init_ao_basis",
|
|
NULL);
|
|
}
|
|
|
|
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
|
|
|
|
*** Initialization functions ~size_max~ is the dimension of the input array, which should be
|
|
equal of larger than the value given in the table of section [[Context]].
|
|
|
|
**** C interface
|
|
|
|
#+NAME:pre2
|
|
#+begin_src c :exports none
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_*",
|
|
NULL);
|
|
}
|
|
|
|
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
|
|
|
|
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 basis_type);
|
|
#+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 basis_type)
|
|
{
|
|
<<pre2>>
|
|
|
|
if (basis_type != 'G' && basis_type != 'S') {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_ao_basis_type",
|
|
NULL);
|
|
}
|
|
|
|
int32_t mask = 1;
|
|
ctx->ao_basis.type = basis_type;
|
|
|
|
<<post2>>
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_shell_num (qmckl_context context,
|
|
const int64_t shell_num);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
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");
|
|
}
|
|
|
|
const int64_t prim_num = ctx->ao_basis.prim_num;
|
|
|
|
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>>
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_prim_num (qmckl_context context,
|
|
const int64_t prim_num);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
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");
|
|
}
|
|
|
|
const int64_t shell_num = ctx->ao_basis.shell_num;
|
|
|
|
if (shell_num <= 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_shell_num",
|
|
"shell_num is not set");
|
|
}
|
|
|
|
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>>
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_ao_num (qmckl_context context,
|
|
const int64_t ao_num);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
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");
|
|
}
|
|
|
|
const int64_t shell_num = ctx->ao_basis.shell_num;
|
|
if (shell_num <= 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_ao_basis_shell_num",
|
|
"shell_num is not set");
|
|
}
|
|
|
|
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>>
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_nucleus_index (qmckl_context context,
|
|
const int64_t* nucleus_index,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_nucleus_index (qmckl_context context,
|
|
const int64_t* nucleus_index,
|
|
const int64_t size_max)
|
|
{
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 4;
|
|
|
|
const int64_t nucl_num = ctx->nucleus.num;
|
|
|
|
if (nucl_num <= 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_nucleus_index",
|
|
"nucl_num is not set");
|
|
}
|
|
|
|
if (size_max < nucl_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_nucleus_index",
|
|
"input array too small");
|
|
}
|
|
|
|
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 = nucl_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>>
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_nucleus_shell_num (qmckl_context context,
|
|
const int64_t* nucleus_shell_num,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_nucleus_shell_num (qmckl_context context,
|
|
const int64_t* nucleus_shell_num,
|
|
const int64_t size_max)
|
|
{
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 3;
|
|
|
|
const int64_t nucl_num = ctx->nucleus.num;
|
|
|
|
if (nucl_num <= 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_nucleus_shell_num",
|
|
"shell_num is not set");
|
|
}
|
|
|
|
if (size_max < nucl_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_set_ao_basis_nucleus_shell_num",
|
|
"input array too small");
|
|
}
|
|
|
|
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 = nucl_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>>
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_shell_ang_mom (qmckl_context context,
|
|
const int32_t* shell_ang_mom,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_shell_ang_mom (qmckl_context context,
|
|
const int32_t* shell_ang_mom,
|
|
const int64_t size_max)
|
|
{
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 5;
|
|
|
|
const int64_t shell_num = ctx->ao_basis.shell_num;
|
|
|
|
if (shell_num == 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_shell_ang_mom",
|
|
"shell_num is not set");
|
|
}
|
|
|
|
if (size_max < shell_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_set_ao_basis_shell_ang_mom",
|
|
"input array too small");
|
|
}
|
|
|
|
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>>
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_shell_prim_num (qmckl_context context,
|
|
const int64_t* shell_prim_num,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_shell_prim_num (qmckl_context context,
|
|
const int64_t* shell_prim_num,
|
|
const int64_t size_max)
|
|
{
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 6;
|
|
|
|
const int64_t shell_num = ctx->ao_basis.shell_num;
|
|
|
|
if (shell_num <= 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_shell_prim_num",
|
|
"shell_num is not set");
|
|
}
|
|
|
|
if (size_max < shell_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_set_ao_basis_shell_prim_num",
|
|
"input array too small");
|
|
}
|
|
|
|
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>>
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_shell_prim_index (qmckl_context context,
|
|
const int64_t* shell_prim_index,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_shell_prim_index (qmckl_context context,
|
|
const int64_t* shell_prim_index,
|
|
const int64_t size_max)
|
|
{
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 7;
|
|
|
|
const int64_t shell_num = ctx->ao_basis.shell_num;
|
|
|
|
if (shell_num <= 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_shell_prim_index",
|
|
"shell_num is not set");
|
|
}
|
|
|
|
if (size_max < shell_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_set_ao_basis_shell_prim_index",
|
|
"input array too small");
|
|
}
|
|
|
|
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>>
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_shell_factor (qmckl_context context,
|
|
const double* shell_factor,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_shell_factor (qmckl_context context,
|
|
const double* shell_factor,
|
|
const int64_t size_max)
|
|
{
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 8;
|
|
|
|
const int64_t shell_num = ctx->ao_basis.shell_num;
|
|
|
|
if (shell_num <= 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_shell_factor",
|
|
"shell_num is not set");
|
|
}
|
|
|
|
if (size_max < shell_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_set_ao_basis_shell_factor",
|
|
"input array too small");
|
|
}
|
|
|
|
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>>
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_exponent (qmckl_context context,
|
|
const double* exponent,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_exponent (qmckl_context context,
|
|
const double* exponent,
|
|
const int64_t size_max)
|
|
{
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 9;
|
|
|
|
const int64_t prim_num = ctx->ao_basis.prim_num;
|
|
|
|
if (prim_num <= 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_exponent",
|
|
"prim_num is not set");
|
|
}
|
|
|
|
if (size_max < prim_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_set_ao_basis_exponent",
|
|
"input array too small");
|
|
}
|
|
|
|
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>>
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_coefficient (qmckl_context context,
|
|
const double* coefficient,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_coefficient (qmckl_context context,
|
|
const double* coefficient,
|
|
const int64_t size_max)
|
|
{
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 10;
|
|
|
|
const int64_t prim_num = ctx->ao_basis.prim_num;
|
|
|
|
if (prim_num <= 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_coefficient",
|
|
"prim_num is not set");
|
|
}
|
|
|
|
if (size_max < prim_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_set_ao_basis_coefficient",
|
|
"input array too small");
|
|
}
|
|
|
|
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
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_prim_factor (qmckl_context context,
|
|
const double* prim_factor,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_prim_factor (qmckl_context context,
|
|
const double* prim_factor,
|
|
const int64_t size_max)
|
|
{
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 11;
|
|
|
|
const int64_t prim_num = ctx->ao_basis.prim_num;
|
|
|
|
if (prim_num <= 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_prim_factor",
|
|
"prim_num is not set");
|
|
}
|
|
|
|
if (size_max < prim_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_set_ao_basis_prim_factor",
|
|
"input array too small");
|
|
}
|
|
|
|
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>>
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_ao_factor (qmckl_context context,
|
|
const double* ao_factor,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_ao_factor (qmckl_context context,
|
|
const double* ao_factor,
|
|
const int64_t size_max)
|
|
{
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 13;
|
|
|
|
const int64_t ao_num = ctx->ao_basis.ao_num;
|
|
|
|
if (ao_num <= 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_ao_basis_ao_factor",
|
|
"ao_num is not set");
|
|
}
|
|
|
|
if (size_max < ao_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_set_ao_basis_ao_factor",
|
|
"input array too small");
|
|
}
|
|
|
|
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>>
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_cartesian (qmckl_context context,
|
|
const bool cartesian);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_set_ao_basis_cartesian (qmckl_context context,
|
|
const bool cartesian)
|
|
{
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1;
|
|
ctx->ao_basis.ao_cartesian = cartesian;
|
|
|
|
<<post2>>
|
|
}
|
|
#+end_src
|
|
|
|
**** Fortran interface
|
|
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_type (context, &
|
|
basis_type) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
character(c_char) , intent(in) , value :: basis_type
|
|
end function qmckl_set_ao_basis_type
|
|
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 qmckl_set_ao_basis_shell_num
|
|
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 qmckl_set_ao_basis_prim_num
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_nucleus_index(context, &
|
|
idx, size_max) 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(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_set_ao_basis_nucleus_index
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_nucleus_shell_num(context, &
|
|
shell_num, size_max) 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(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_set_ao_basis_nucleus_shell_num
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_shell_ang_mom(context, &
|
|
shell_ang_mom, size_max) 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(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_set_ao_basis_shell_ang_mom
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_shell_prim_num(context, &
|
|
shell_prim_num, size_max) 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(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_set_ao_basis_shell_prim_num
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_shell_prim_index(context, &
|
|
shell_prim_index, size_max) 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(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_set_ao_basis_shell_prim_index
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_shell_factor(context, &
|
|
shell_factor, size_max) 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(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_set_ao_basis_shell_factor
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_exponent(context, &
|
|
exponent, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double) , intent(in) :: exponent(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_set_ao_basis_exponent
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_coefficient(context, &
|
|
coefficient, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double) , intent(in) :: coefficient(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_set_ao_basis_coefficient
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_prim_factor(context, &
|
|
prim_factor, size_max) 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(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_set_ao_basis_prim_factor
|
|
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 qmckl_set_ao_basis_ao_num
|
|
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 qmckl_set_ao_basis_cartesian
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_ao_basis_ao_factor(context, &
|
|
ao_factor, size_max) 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(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_set_ao_basis_ao_factor
|
|
end interface
|
|
#+end_src
|
|
|
|
*** Access functions ~size_max~ is the dimension of the input array, which should be
|
|
equal of larger than the value given in the table of section [[Context]].
|
|
|
|
**** C interface
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_type (const qmckl_context context,
|
|
char* const basis_type);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_type (const qmckl_context context,
|
|
char* const basis_type)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_type",
|
|
NULL);
|
|
}
|
|
|
|
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 qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_get_ao_basis_type",
|
|
NULL);
|
|
|
|
}
|
|
|
|
assert (ctx->ao_basis.type != (char) 0);
|
|
|
|
basis_type[0] = ctx->ao_basis.type;
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_shell_num (const qmckl_context context,
|
|
int64_t* const shell_num);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_shell_num (const qmckl_context context,
|
|
int64_t* const shell_num)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_shell_factor",
|
|
NULL);
|
|
}
|
|
|
|
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 qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_get_ao_basis_shell_num",
|
|
NULL);
|
|
}
|
|
|
|
assert (ctx->ao_basis.shell_num > (int64_t) 0);
|
|
,*shell_num = ctx->ao_basis.shell_num;
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_prim_num (const qmckl_context context,
|
|
int64_t* const prim_num);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_prim_num (const qmckl_context context,
|
|
int64_t* const prim_num)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_prim_num",
|
|
NULL);
|
|
}
|
|
|
|
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 qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_get_ao_basis_prim_num",
|
|
NULL);
|
|
}
|
|
|
|
assert (ctx->ao_basis.prim_num > (int64_t) 0);
|
|
|
|
,*prim_num = ctx->ao_basis.prim_num;
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_nucleus_shell_num (const qmckl_context context,
|
|
int64_t* const nucleus_shell_num,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_nucleus_shell_num (const qmckl_context context,
|
|
int64_t* const nucleus_shell_num,
|
|
const int64_t size_max)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_nucleus_shell_num",
|
|
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 qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_get_ao_basis_nucleus_shell_num",
|
|
NULL);
|
|
}
|
|
|
|
if (nucleus_shell_num == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_ao_basis_nucleus_shell_num",
|
|
"NULL pointer");
|
|
}
|
|
|
|
if (size_max < ctx->nucleus.num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_get_ao_basis_nucleus_shell_num",
|
|
"Array too small. Expected nucl_num");
|
|
}
|
|
|
|
assert (ctx->ao_basis.nucleus_shell_num != NULL);
|
|
memcpy(nucleus_shell_num, ctx->ao_basis.nucleus_shell_num,
|
|
(size_t) ctx->nucleus.num * sizeof(int64_t));
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_nucleus_index (const qmckl_context context,
|
|
int64_t* const nucleus_index,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_nucleus_index (const qmckl_context context,
|
|
int64_t* const nucleus_index,
|
|
const int64_t size_max)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_nucleus_index",
|
|
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 qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_get_ao_basis_nucleus_index",
|
|
NULL);
|
|
}
|
|
|
|
if (nucleus_index == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_ao_basis_nucleus_index",
|
|
"NULL pointer");
|
|
}
|
|
|
|
if (size_max < ctx->nucleus.num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_get_ao_basis_nucleus_index",
|
|
"Array too small. Expected shell_num");
|
|
}
|
|
|
|
assert (ctx->ao_basis.nucleus_index != NULL);
|
|
memcpy(nucleus_index, ctx->ao_basis.nucleus_index,
|
|
(size_t) ctx->nucleus.num * sizeof(int64_t));
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context,
|
|
int32_t* const shell_ang_mom,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context,
|
|
int32_t* const shell_ang_mom,
|
|
const int64_t size_max)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_shell_ang_mom",
|
|
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 qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_get_ao_basis_shell_ang_mom",
|
|
NULL);
|
|
}
|
|
|
|
if (shell_ang_mom == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_ao_basis_shell_ang_mom",
|
|
"NULL pointer");
|
|
}
|
|
|
|
if (size_max < ctx->ao_basis.shell_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_get_ao_basis_shell_ang_mom",
|
|
"Array too small. Expected shell_num");
|
|
}
|
|
|
|
assert (ctx->ao_basis.shell_ang_mom != NULL);
|
|
memcpy(shell_ang_mom, ctx->ao_basis.shell_ang_mom,
|
|
(size_t) ctx->ao_basis.shell_num * sizeof(int32_t));
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_shell_prim_num (const qmckl_context context,
|
|
int64_t* const shell_prim_num,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_shell_prim_num (const qmckl_context context,
|
|
int64_t* const shell_prim_num,
|
|
const int64_t size_max)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_shell_prim_num",
|
|
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 qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_get_ao_basis_shell_prim_num",
|
|
NULL);
|
|
}
|
|
|
|
if (shell_prim_num == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_ao_basis_shell_prim_num",
|
|
"NULL pointer");
|
|
}
|
|
|
|
if (size_max < ctx->ao_basis.shell_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_get_ao_basis_shell_prim_num",
|
|
"Array too small. Expected shell_num");
|
|
}
|
|
|
|
assert (ctx->ao_basis.shell_prim_num != NULL);
|
|
memcpy(shell_prim_num, ctx->ao_basis.shell_prim_num,
|
|
(size_t) ctx->ao_basis.shell_num * sizeof(int64_t));
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_shell_prim_index (const qmckl_context context,
|
|
int64_t* const shell_prim_index,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_shell_prim_index (const qmckl_context context,
|
|
int64_t* const shell_prim_index,
|
|
const int64_t size_max)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_shell_prim_index",
|
|
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 qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_get_ao_basis_shell_prim_index",
|
|
NULL);
|
|
}
|
|
|
|
if (shell_prim_index == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_ao_basis_shell_prim_index",
|
|
"NULL pointer");
|
|
}
|
|
|
|
if (size_max < ctx->ao_basis.shell_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_get_ao_basis_shell_prim_index",
|
|
"Array too small. Expected shell_num");
|
|
}
|
|
|
|
assert (ctx->ao_basis.shell_prim_index != NULL);
|
|
memcpy(shell_prim_index, ctx->ao_basis.shell_prim_index,
|
|
(size_t) ctx->ao_basis.shell_num * sizeof(int64_t));
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_shell_factor (const qmckl_context context,
|
|
double* const shell_factor,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_shell_factor (const qmckl_context context,
|
|
double* const shell_factor,
|
|
const int64_t size_max)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_shell_factor",
|
|
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 qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_get_ao_basis_shell_factor",
|
|
NULL);
|
|
}
|
|
|
|
if (shell_factor == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_ao_basis_shell_factor",
|
|
"NULL pointer");
|
|
}
|
|
|
|
if (size_max < ctx->ao_basis.shell_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_get_ao_basis_shell_factor",
|
|
"Array too small. Expected shell_num");
|
|
}
|
|
|
|
assert (ctx->ao_basis.shell_factor != NULL);
|
|
memcpy(shell_factor, ctx->ao_basis.shell_factor,
|
|
(size_t) ctx->ao_basis.shell_num * sizeof(double));
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_exponent (const qmckl_context context,
|
|
double* const exponent,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_exponent (const qmckl_context context,
|
|
double* const exponent,
|
|
const int64_t size_max)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_exponent",
|
|
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 qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_get_ao_basis_exponent",
|
|
NULL);
|
|
}
|
|
|
|
if (exponent == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_ao_basis_exponent",
|
|
"NULL pointer");
|
|
}
|
|
|
|
if (size_max < ctx->ao_basis.prim_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_get_ao_basis_exponent",
|
|
"Array too small. Expected prim_num");
|
|
}
|
|
|
|
assert (ctx->ao_basis.exponent != NULL);
|
|
memcpy(exponent, ctx->ao_basis.exponent,
|
|
(size_t) ctx->ao_basis.prim_num * sizeof(double));
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_coefficient (const qmckl_context context,
|
|
double* const coefficient,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_coefficient (const qmckl_context context,
|
|
double* const coefficient,
|
|
const int64_t size_max)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_coefficient",
|
|
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 qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_get_ao_basis_coefficient",
|
|
NULL);
|
|
}
|
|
|
|
if (coefficient == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_ao_basis_coefficient",
|
|
"NULL pointer");
|
|
}
|
|
|
|
if (size_max < ctx->ao_basis.prim_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_get_ao_basis_coefficient",
|
|
"Array too small. Expected prim_num");
|
|
}
|
|
assert (ctx->ao_basis.coefficient != NULL);
|
|
memcpy(coefficient, ctx->ao_basis.coefficient,
|
|
(size_t) ctx->ao_basis.prim_num * sizeof(double));
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_prim_factor (const qmckl_context context,
|
|
double* const prim_factor,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_prim_factor (const qmckl_context context,
|
|
double* const prim_factor,
|
|
const int64_t size_max)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_prim_factor",
|
|
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 qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_get_ao_basis_prim_factor",
|
|
NULL);
|
|
}
|
|
|
|
if (prim_factor == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_ao_basis_prim_factor",
|
|
"NULL pointer");
|
|
}
|
|
|
|
if (size_max < ctx->ao_basis.prim_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_get_ao_basis_prim_factor",
|
|
"Array too small. Expected prim_num");
|
|
}
|
|
|
|
assert (ctx->ao_basis.prim_factor != NULL);
|
|
memcpy(prim_factor, ctx->ao_basis.prim_factor,
|
|
(size_t) ctx->ao_basis.prim_num * sizeof(double));
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_ao_num (const qmckl_context context,
|
|
int64_t* const ao_num);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_ao_num (const qmckl_context context,
|
|
int64_t* const ao_num)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_ao_num",
|
|
NULL);
|
|
}
|
|
|
|
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 qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_get_ao_basis_ao_num",
|
|
NULL);
|
|
}
|
|
|
|
assert (ctx->ao_basis.ao_num > (int64_t) 0);
|
|
|
|
*ao_num = ctx->ao_basis.ao_num;
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_ao_factor (const qmckl_context context,
|
|
double* const ao_factor,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_ao_factor (const qmckl_context context,
|
|
double* const ao_factor,
|
|
const int64_t size_max)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_ao_factor",
|
|
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 qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_get_ao_basis_ao_factor",
|
|
NULL);
|
|
}
|
|
|
|
if (ao_factor == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_ao_basis_ao_factor",
|
|
"NULL pointer");
|
|
}
|
|
|
|
if (size_max < ctx->ao_basis.ao_num) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_get_ao_basis_ao_factor",
|
|
"Array too small. Expected ao_num");
|
|
}
|
|
|
|
assert (ctx->ao_basis.ao_factor != NULL);
|
|
memcpy(ao_factor, ctx->ao_basis.ao_factor, ctx->ao_basis.ao_num * sizeof(double));
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+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
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
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
|
|
|
|
**** Fortran interface
|
|
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_ao_basis_type (context, &
|
|
basis_type) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
character(c_char) , intent(out) :: basis_type
|
|
end function qmckl_get_ao_basis_type
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_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(out) :: num
|
|
end function qmckl_get_ao_basis_shell_num
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_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(out) :: num
|
|
end function qmckl_get_ao_basis_prim_num
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_ao_basis_nucleus_shell_num(context, &
|
|
shell_num, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(out) :: shell_num(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_get_ao_basis_nucleus_shell_num
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_ao_basis_nucleus_index(context, &
|
|
idx, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(out) :: idx(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_get_ao_basis_nucleus_index
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_ao_basis_shell_ang_mom(context, &
|
|
shell_ang_mom, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int32_t) , intent(out) :: shell_ang_mom(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_get_ao_basis_shell_ang_mom
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_ao_basis_shell_prim_num(context, &
|
|
shell_prim_num, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(out) :: shell_prim_num(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_get_ao_basis_shell_prim_num
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_ao_basis_shell_prim_index(context, &
|
|
shell_prim_index, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(out) :: shell_prim_index(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_get_ao_basis_shell_prim_index
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_ao_basis_shell_factor(context, &
|
|
shell_factor, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double) , intent(out) :: shell_factor(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_get_ao_basis_shell_factor
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_ao_basis_exponent(context, &
|
|
exponent, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double) , intent(out) :: exponent(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_get_ao_basis_exponent
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_ao_basis_coefficient(context, &
|
|
coefficient, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double) , intent(out) :: coefficient(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_get_ao_basis_coefficient
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_ao_basis_prim_factor(context, &
|
|
prim_factor, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double) , intent(out) :: prim_factor(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_get_ao_basis_prim_factor
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_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(out) :: num
|
|
end function qmckl_get_ao_basis_ao_num
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_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(out) :: cartesian
|
|
end function qmckl_get_ao_basis_cartesian
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_ao_basis_ao_factor(context, &
|
|
ao_factor, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double) , intent(out) :: ao_factor(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_get_ao_basis_ao_factor
|
|
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]), 3*nucl_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
|
|
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, nucl_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, nucl_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, shell_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_shell_factor (context, shell_factor, shell_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, shell_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, shell_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_exponent (context, exponent, prim_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_coefficient (context, coefficient, prim_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_ao_basis_provided(context));
|
|
|
|
rc = qmckl_set_ao_basis_prim_factor (context, prim_factor, prim_num);
|
|
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, ao_num);
|
|
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 ;
|
|
|
|
|
|
rc = qmckl_get_ao_basis_type (context, &typ_test);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
assert(typ == typ_test);
|
|
|
|
rc = qmckl_get_ao_basis_shell_num (context, &shell_num_test);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
assert(shell_num == shell_num_test);
|
|
|
|
rc = qmckl_get_ao_basis_prim_num (context, &prim_num_test);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
assert(prim_num == prim_num_test);
|
|
|
|
nucleus_index_test = (int64_t*) malloc (nucl_num * sizeof(int64_t));
|
|
rc = qmckl_get_ao_basis_nucleus_index (context, nucleus_index_test, nucl_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
for (int64_t i=0 ; i < nucl_num ; ++i) {
|
|
assert(nucleus_index_test[i] == nucleus_index[i]);
|
|
}
|
|
free(nucleus_index_test);
|
|
|
|
nucleus_shell_num_test = (int64_t*) malloc ( nucl_num * sizeof(int64_t));
|
|
rc = qmckl_get_ao_basis_nucleus_shell_num (context, nucleus_shell_num_test, nucl_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
for (int64_t i=0 ; i < nucl_num ; ++i) {
|
|
assert(nucleus_shell_num_test[i] == nucleus_shell_num[i]);
|
|
}
|
|
|
|
shell_ang_mom_test = (int32_t*) malloc ( shell_num * sizeof(int32_t));
|
|
rc = qmckl_get_ao_basis_shell_ang_mom (context, shell_ang_mom_test, shell_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
for (int64_t i=0 ; i < shell_num ; ++i) {
|
|
assert(shell_ang_mom_test[i] == shell_ang_mom[i]);
|
|
}
|
|
free(shell_ang_mom_test);
|
|
|
|
shell_factor_test = (double*) malloc ( shell_num * sizeof(double));
|
|
rc = qmckl_get_ao_basis_shell_factor (context, shell_factor_test, shell_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
for (int64_t i=0 ; i < shell_num ; ++i) {
|
|
assert(shell_factor_test[i] == shell_factor[i]);
|
|
}
|
|
free(shell_factor_test);
|
|
|
|
shell_prim_num_test = (int64_t*) malloc ( shell_num * sizeof(int64_t));
|
|
rc = qmckl_get_ao_basis_shell_prim_num (context, shell_prim_num_test, shell_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
for (int64_t i=0 ; i < shell_num ; ++i) {
|
|
assert(shell_prim_num_test[i] == shell_prim_num[i]);
|
|
}
|
|
free(shell_prim_num_test);
|
|
|
|
shell_prim_index_test = (int64_t*) malloc ( shell_num * sizeof(int64_t));
|
|
rc = qmckl_get_ao_basis_shell_prim_index (context, shell_prim_index_test, shell_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
for (int64_t i=0 ; i < shell_num ; ++i) {
|
|
assert(shell_prim_index_test[i] == shell_prim_index[i]);
|
|
}
|
|
free(shell_prim_index_test);
|
|
|
|
exponent_test = (double*) malloc ( prim_num * sizeof(double));
|
|
rc = qmckl_get_ao_basis_exponent(context, exponent_test, prim_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
for (int64_t i=0 ; i < prim_num ; ++i) {
|
|
assert(exponent_test[i] == exponent[i]);
|
|
}
|
|
free(exponent_test);
|
|
|
|
coefficient_test = (double*) malloc ( prim_num * sizeof(double));
|
|
rc = qmckl_get_ao_basis_coefficient(context, coefficient_test, prim_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
for (int64_t i=0 ; i < prim_num ; ++i) {
|
|
assert(coefficient_test[i] == coefficient[i]);
|
|
}
|
|
free(coefficient_test);
|
|
|
|
prim_factor_test = (double*) malloc ( prim_num * sizeof(double));
|
|
rc = qmckl_get_ao_basis_prim_factor (context, prim_factor_test, prim_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
for (int64_t i=0 ; i < prim_num ; ++i) {
|
|
assert(prim_factor_test[i] == prim_factor[i]);
|
|
}
|
|
free(prim_factor_test);
|
|
|
|
rc = qmckl_get_ao_basis_ao_num(context, &ao_num_test);
|
|
assert(ao_num == ao_num_test);
|
|
|
|
ao_factor_test = (double*) malloc ( ao_num * sizeof(double));
|
|
rc = qmckl_get_ao_basis_ao_factor (context, ao_factor_test, ao_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
for (int64_t i=0 ; i < ao_num ; ++i) {
|
|
assert(ao_factor_test[i] == ao_factor[i]);
|
|
}
|
|
free(ao_factor_test);
|
|
|
|
#+end_src
|
|
|
|
** Computed data
|
|
|
|
The following data is computed as described in the next sections:
|
|
|
|
|----------------------+-----------------------------------+----------------------------------------------------------------------------------------------|
|
|
| Variable | Type | Description |
|
|
|----------------------+-----------------------------------+----------------------------------------------------------------------------------------------|
|
|
| ~primitive_vgl~ | ~double[point_num][5][prim_num]~ | Value, gradients, Laplacian of the primitives at current positions |
|
|
| ~primitive_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the primitives at current positions |
|
|
| ~shell_vgl~ | ~double[point_num][5][shell_num]~ | Value, gradients, Laplacian of the primitives at current positions |
|
|
| ~shell_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the AOs at current positions |
|
|
| ~ao_vgl~ | ~double[point_num][5][ao_num]~ | Value, gradients, Laplacian of the primitives at current positions |
|
|
| ~ao_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the AOs at current positions |
|
|
|
|
*** After initialization
|
|
|
|
When the basis set is completely entered, extra data structures may be
|
|
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_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_finalize_basis",
|
|
NULL);
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int64_t nucl_num = 0;
|
|
|
|
qmckl_exit_code 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;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
rc = QMCKL_SUCCESS;
|
|
#ifdef HAVE_HPC
|
|
rc = qmckl_finalize_basis_hpc(context);
|
|
#endif
|
|
|
|
return rc;
|
|
}
|
|
#+end_src
|
|
|
|
*** HPC-specific data structures
|
|
|
|
For faster access, we provide extra arrays for the shell information as:
|
|
- $C_{psa}$ = =coef_per_nucleus (iprim, ishell, inucl)=
|
|
- $\gamma_{pa}$ =expo_per_nucleus (iprim, inucl)=
|
|
|
|
such that the computation of the radial parts can be vectorized.
|
|
|
|
Exponents are sorted in increasing order to exit loops quickly,
|
|
and only unique exponents are kept. This also allows to pack the
|
|
exponents to enable vectorization of exponentials.
|
|
|
|
The computation of the radial part is made as
|
|
\[
|
|
R_{sa} = \sum_p C_{psa} \gamma_{pa}
|
|
\]
|
|
which is a matrix-vector product.
|
|
|
|
#+NAME:HPC_struct
|
|
#+begin_src c :comments org :exports none
|
|
/* HPC specific data structures */
|
|
qmckl_tensor coef_per_nucleus;
|
|
qmckl_matrix expo_per_nucleus;
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :exports none
|
|
#ifdef HAVE_HPC
|
|
qmckl_exit_code qmckl_finalize_basis_hpc (qmckl_context context);
|
|
#endif
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :exports none
|
|
#ifdef HAVE_HPC
|
|
qmckl_exit_code qmckl_finalize_basis_hpc (qmckl_context context)
|
|
{
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
|
|
int64_t shell_max = 0;
|
|
int64_t prim_max = 0;
|
|
int64_t nucl_num = ctx->nucleus.num;
|
|
for (int inucl=0 ; inucl < nucl_num ; ++inucl) {
|
|
shell_max = ctx->ao_basis.nucleus_shell_num[inucl] > shell_max ?
|
|
ctx->ao_basis.nucleus_shell_num[inucl] : shell_max;
|
|
|
|
int64_t prim_num = 0;
|
|
const int64_t ishell_start = ctx->ao_basis.nucleus_index[inucl];
|
|
const int64_t ishell_end = ctx->ao_basis.nucleus_index[inucl] + ctx->ao_basis.nucleus_shell_num[inucl];
|
|
for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {
|
|
prim_num += ctx->ao_basis.shell_prim_num[ishell];
|
|
}
|
|
prim_max = prim_num > prim_max ?
|
|
prim_num : prim_max;
|
|
}
|
|
|
|
int64_t size[3] = { prim_max, shell_max, nucl_num };
|
|
ctx->ao_basis.coef_per_nucleus = qmckl_tensor_alloc( context, 3, size );
|
|
ctx->ao_basis.coef_per_nucleus = qmckl_tensor_set(ctx->ao_basis.coef_per_nucleus, 0.);
|
|
|
|
ctx->ao_basis.expo_per_nucleus = qmckl_matrix_alloc( context, prim_max, nucl_num );
|
|
ctx->ao_basis.expo_per_nucleus = qmckl_matrix_set(ctx->ao_basis.expo_per_nucleus, 0.);
|
|
|
|
double expo[prim_max];
|
|
double coef[shell_max][prim_max];
|
|
for (int inucl=0 ; inucl < nucl_num ; ++inucl) {
|
|
memset(expo, 0, sizeof(expo));
|
|
memset(coef, 0, sizeof(expo));
|
|
|
|
int64_t idx = 0;
|
|
const int64_t ishell_start = ctx->ao_basis.nucleus_index[inucl];
|
|
const int64_t ishell_end = ctx->ao_basis.nucleus_index[inucl] + ctx->ao_basis.nucleus_shell_num[inucl];
|
|
for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {
|
|
|
|
const int64_t iprim_start = ctx->ao_basis.shell_prim_index[ishell];
|
|
const int64_t iprim_end = ctx->ao_basis.shell_prim_index[ishell] + ctx->ao_basis.shell_prim_num[ishell];
|
|
for (int64_t iprim = iprim_start ; iprim < iprim_end ; ++iprim) {
|
|
expo[idx] = ctx->ao_basis.coefficient_normalized[iprim];
|
|
coef[ishell - ishell_start][idx] = ctx->ao_basis.coefficient_normalized[iprim];
|
|
idx += 1;
|
|
}
|
|
}
|
|
|
|
for (int64_t i=0 ; i<prim_max ; ++i) {
|
|
qmckl_mat( ctx->ao_basis.expo_per_nucleus, i, inucl ) = expo[i];
|
|
}
|
|
|
|
for (int64_t j=0 ; j<shell_max ; ++j) {
|
|
for (int64_t i=0 ; i<prim_max ; ++i) {
|
|
qmckl_ten3( ctx->ao_basis.coef_per_nucleus, i, j, inucl ) = coef[j][i];
|
|
}
|
|
}
|
|
|
|
}
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#endif
|
|
|
|
#+end_src
|
|
|
|
*** Access functions
|
|
|
|
#+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,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
Returns the array of values, gradients an Laplacian of primitive
|
|
basis functions evaluated at the current coordinates.
|
|
See section [[Computation of primitives]].
|
|
|
|
#+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,
|
|
const int64_t size_max)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_primitive_vgl",
|
|
NULL);
|
|
}
|
|
|
|
if (size_max <= 0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_get_ao_basis_primitive_vgl",
|
|
"size_max <= 0");
|
|
}
|
|
|
|
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);
|
|
|
|
int64_t sze = ctx->ao_basis.prim_num * 5 * ctx->point.num;
|
|
if (size_max < sze) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_get_ao_basis_primitive_vgl",
|
|
"input array too small");
|
|
}
|
|
memcpy(primitive_vgl, ctx->ao_basis.primitive_vgl, (size_t) 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_primitive_vgl &
|
|
(context, primitive_vgl, size_max) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
double precision, intent(out) :: primitive_vgl(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function
|
|
end interface
|
|
#+end_src
|
|
|
|
|
|
#+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,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
Returns the array of values, gradients an Laplacian of contracted shells
|
|
evaluated at the current coordinates. See section [[Computation of shells]].
|
|
|
|
#+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,
|
|
const int64_t size_max)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_shell_vgl",
|
|
NULL);
|
|
}
|
|
|
|
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);
|
|
|
|
int64_t sze = ctx->ao_basis.shell_num * 5 * ctx->point.num;
|
|
if (size_max < sze) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_get_ao_basis_shell_vgl",
|
|
"input array too small");
|
|
}
|
|
memcpy(shell_vgl, ctx->ao_basis.shell_vgl, (size_t)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, size_max) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
double precision, intent(out) :: shell_vgl(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function
|
|
end interface
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_ao_vgl (qmckl_context context,
|
|
double* const ao_vgl,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
Returns the array of values, gradients an Laplacian of the atomic orbitals
|
|
evaluated at the current coordinates.
|
|
See section [[Combining radial and polynomial parts]].
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_ao_vgl (qmckl_context context,
|
|
double* const ao_vgl,
|
|
const int64_t size_max)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_ao_vgl",
|
|
NULL);
|
|
}
|
|
|
|
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);
|
|
|
|
int64_t sze = ctx->ao_basis.ao_num * 5 * ctx->point.num;
|
|
if (size_max < sze) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_get_ao_basis_ao_vgl",
|
|
"input array too small");
|
|
}
|
|
memcpy(ao_vgl, ctx->ao_basis.ao_vgl, (size_t) 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_ao_vgl (context, &
|
|
ao_vgl, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
double precision, intent(out) :: ao_vgl(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_get_ao_basis_ao_vgl
|
|
end interface
|
|
#+end_src
|
|
|
|
Uses the give array to compute the VGL.
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_ao_vgl_inplace (qmckl_context context,
|
|
double* const ao_vgl,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_ao_basis_ao_vgl_inplace (qmckl_context context,
|
|
double* const ao_vgl,
|
|
const int64_t size_max)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_ao_basis_ao_vgl",
|
|
NULL);
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int64_t sze = ctx->ao_basis.ao_num * 5 * ctx->point.num;
|
|
if (size_max < sze) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_get_ao_basis_ao_vgl",
|
|
"input array too small");
|
|
}
|
|
|
|
rc = qmckl_context_touch(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
double* old_array = ctx->ao_basis.ao_vgl;
|
|
|
|
ctx->ao_basis.ao_vgl = ao_vgl;
|
|
|
|
rc = qmckl_provide_ao_vgl(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
ctx->ao_basis.ao_vgl = old_array;
|
|
|
|
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_ao_vgl_inplace (context, &
|
|
ao_vgl, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
double precision, intent(out) :: ao_vgl(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function qmckl_get_ao_basis_ao_vgl_inplace
|
|
end interface
|
|
#+end_src
|
|
|
|
* Radial part
|
|
|
|
** General functions for Gaussian basis functions ~qmckl_ao_gaussian_vgl~ computes the values, gradients and
|
|
Laplacians at a given point of ~n~ Gaussian functions centered at
|
|
the same point:
|
|
|
|
\[ v_i = \exp(-a_i |X-R|^2) \]
|
|
\[ \nabla_x v_i = -2 a_i (X_x - R_x) v_i \]
|
|
\[ \nabla_y v_i = -2 a_i (X_y - R_y) v_i \]
|
|
\[ \nabla_z v_i = -2 a_i (X_z - R_z) v_i \]
|
|
\[ \Delta v_i = a_i (4 |X-R|^2 a_i - 6) v_i \]
|
|
|
|
|--------------+------------------+------------------------------------------------------|
|
|
| Variable | Type | Description |
|
|
|--------------+------------------+------------------------------------------------------|
|
|
| ~context~ | ~qmckl_context~ | Global state |
|
|
| ~X(3)~ | ~double[3]~ | Array containing the coordinates of the points |
|
|
| ~R(3)~ | ~double[3]~ | Array containing the x,y,z coordinates of the center |
|
|
| ~n~ | ~int64_t~ | Number of computed Gaussians |
|
|
| ~A(n)~ | ~double[n]~ | Exponents of the Gaussians |
|
|
| ~VGL(ldv,5)~ | ~double[5][ldv]~ | Value, gradients and Laplacian of the Gaussians |
|
|
| ~ldv~ | ~int64_t~ | Leading dimension of array ~VGL~ |
|
|
|--------------+------------------+------------------------------------------------------|
|
|
|
|
Requirements:
|
|
|
|
- ~context~ \ne 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 :noexport:
|
|
|
|
#+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, z
|
|
double precision, allocatable :: VGL(:,:), A(:)
|
|
double precision :: epsilon
|
|
|
|
epsilon = 3.d0 * 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
|
|
z = dabs(1.d0 - VGL(i,1) / (dexp(-A(i) * r2)) )
|
|
if ( z > epsilon ) then
|
|
print *, z, epsilon
|
|
return
|
|
end if
|
|
|
|
test_qmckl_ao_gaussian_vgl = -12
|
|
z = dabs(1.d0 - VGL(i,2) / (&
|
|
-2.d0 * A(i) * Y(1) * dexp(-A(i) * r2) ))
|
|
if ( z > epsilon ) then
|
|
print *, z, epsilon
|
|
return
|
|
end if
|
|
|
|
test_qmckl_ao_gaussian_vgl = -13
|
|
z = dabs(1.d0 - VGL(i,3) / (&
|
|
-2.d0 * A(i) * Y(2) * dexp(-A(i) * r2) ))
|
|
if ( z > epsilon ) then
|
|
print *, z, epsilon
|
|
return
|
|
end if
|
|
|
|
test_qmckl_ao_gaussian_vgl = -14
|
|
z = dabs(1.d0 - VGL(i,4) / (&
|
|
-2.d0 * A(i) * Y(3) * dexp(-A(i) * r2) ))
|
|
if ( z > epsilon ) then
|
|
print *, z, epsilon
|
|
return
|
|
end if
|
|
|
|
test_qmckl_ao_gaussian_vgl = -15
|
|
z = dabs(1.d0 - VGL(i,5) / (&
|
|
A(i) * (4.d0*r2*A(i) - 6.d0) * dexp(-A(i) * r2) ))
|
|
if ( z > epsilon ) then
|
|
print *, z, epsilon
|
|
return
|
|
end if
|
|
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 General functions for Slater basis functions :noexport:
|
|
** TODO General functions for Radial functions on a grid :noexport:
|
|
** TODO Helper functions to accelerate calculations :noexport:
|
|
|
|
|--------------------------+---------------------+--------------------------------------------------------|
|
|
| Variable | Type | Description |
|
|
|--------------------------+---------------------+--------------------------------------------------------|
|
|
| ~coefficient_normalized~ | ~double[prim_num]~ | Normalized primitive coefficients |
|
|
| ~nucleus_prim_index~ | ~int64_t[nucl_num]~ | Index of the first primitive for each nucleus |
|
|
| ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | Maximum angular momentum for each nucleus |
|
|
| ~nucleus_range~ | ~double[nucl_num]~ | Distance beyond which all the AOs are zero |
|
|
|--------------------------+---------------------+--------------------------------------------------------|
|
|
| ~nucl_shell_index~ | ~int64_t[nucl_num]~ | Index of the first shell for each nucleus |
|
|
| ~exponent_sorted~ | ~double[prim_num]~ | Array of exponents for sorted primitives |
|
|
| ~coeff_norm_sorted~ | ~double[prim_num]~ | Array of normalized coefficients for sorted primitives |
|
|
| ~prim_factor_sorted~ | ~double[prim_num]~ | Normalization factors of the sorted primtives |
|
|
|
|
** Computation of primitives
|
|
: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
|
|
| Variable | Type | In/Out | Description |
|
|
|----------------------+----------------------------------+--------+--------------------------------------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~prim_num~ | ~int64_t~ | in | Number of primitives |
|
|
| ~point_num~ | ~int64_t~ | in | Number of points considered |
|
|
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
|
|
| ~nucleus_prim_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st primitive of each nucleus |
|
|
| ~coord~ | ~double[3][point_num]~ | in | Coordinates |
|
|
| ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates |
|
|
| ~expo~ | ~double[prim_num]~ | in | Exponents of the primitives |
|
|
| ~primitive_vgl~ | ~double[point_num][5][prim_num]~ | out | Value, gradients and Laplacian of the primitives |
|
|
|
|
#+CALL: generate_c_header(table=qmckl_ao_basis_primitive_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_basis_primitive_gaussian_vgl"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl (
|
|
const qmckl_context context,
|
|
const int64_t prim_num,
|
|
const int64_t point_num,
|
|
const int64_t nucl_num,
|
|
const int64_t* nucleus_prim_index,
|
|
const double* coord,
|
|
const double* nucl_coord,
|
|
const double* expo,
|
|
double* const primitive_vgl );
|
|
#+end_src
|
|
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f( &
|
|
context, prim_num, point_num, nucl_num, &
|
|
nucleus_prim_index, 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) :: point_num
|
|
integer*8 , intent(in) :: nucleus_prim_index(nucl_num+1)
|
|
double precision , intent(in) :: coord(point_num,3)
|
|
double precision , intent(in) :: nucl_coord(nucl_num,3)
|
|
double precision , intent(in) :: expo(prim_num)
|
|
double precision , intent(out) :: primitive_vgl(prim_num,5,point_num)
|
|
|
|
integer*8 :: inucl, iprim, ipoint
|
|
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 ipoint = 1, point_num
|
|
x = coord(ipoint,1) - nucl_coord(inucl,1)
|
|
y = coord(ipoint,2) - nucl_coord(inucl,2)
|
|
z = coord(ipoint,3) - 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, 1, ipoint) = v
|
|
primitive_vgl(iprim, 2, ipoint) = two_a * x
|
|
primitive_vgl(iprim, 3, ipoint) = two_a * y
|
|
primitive_vgl(iprim, 4, ipoint) = two_a * z
|
|
primitive_vgl(iprim, 5, ipoint) = two_a * (3.d0 - 2.d0*ar2)
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end function qmckl_compute_ao_basis_primitive_gaussian_vgl_f
|
|
#+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, &
|
|
point_num, &
|
|
nucl_num, &
|
|
nucleus_prim_index, &
|
|
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 :: point_num
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
integer (c_int64_t) , intent(in) :: nucleus_prim_index(nucl_num)
|
|
real (c_double ) , intent(in) :: coord(point_num,3)
|
|
real (c_double ) , intent(in) :: nucl_coord(nucl_num,3)
|
|
real (c_double ) , intent(in) :: expo(prim_num)
|
|
real (c_double ) , intent(out) :: primitive_vgl(prim_num,5,point_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, &
|
|
point_num, &
|
|
nucl_num, &
|
|
nucleus_prim_index, &
|
|
coord, &
|
|
nucl_coord, &
|
|
expo, &
|
|
primitive_vgl)
|
|
|
|
end function qmckl_compute_ao_basis_primitive_gaussian_vgl
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
|
|
#+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_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_provide_get_ao_basis_primitive_vgl",
|
|
NULL);
|
|
}
|
|
|
|
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->point.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->point.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->point.num,
|
|
ctx->nucleus.num,
|
|
ctx->ao_basis.nucleus_prim_index,
|
|
ctx->point.coord.data,
|
|
ctx->nucleus.coord.data,
|
|
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
|
|
|
|
*** Test :noexport:
|
|
|
|
#+begin_src python :results output :exports none :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 ] )
|
|
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][elec_num];
|
|
a = 0.9059; x = elec_26_w1 ; y = nucl_1
|
|
print ( "[7][0][26] : %e"% f(a,x,y))
|
|
print ( "[7][1][26] : %e"% df(a,x,y,1))
|
|
print ( "[7][2][26] : %e"% df(a,x,y,2))
|
|
print ( "[7][3][26] : %e"% df(a,x,y,3))
|
|
print ( "[7][4][26] : %e"% lf(a,x,y))
|
|
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: [7][0][26] : 1.050157e-03
|
|
: [7][1][26] : -7.501497e-04
|
|
: [7][2][26] : -3.825069e-03
|
|
: [7][3][26] : 3.495056e-03
|
|
: [7][4][26] : 2.040013e-02
|
|
|
|
#+begin_src c :tangle (eval c_test) :exports none
|
|
{
|
|
#define walk_num 1 // 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, walk_num*elec_num*3);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
|
|
|
|
double prim_vgl[elec_num*walk_num][5][prim_num];
|
|
|
|
rc = qmckl_get_ao_basis_primitive_vgl(context, &(prim_vgl[0][0][0]),
|
|
(int64_t) 5*elec_num*walk_num*prim_num );
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
assert( fabs(prim_vgl[26][0][7] - ( 1.0501570432064878E-003)) < 1.e-14 );
|
|
assert( fabs(prim_vgl[26][1][7] - (-7.5014974095310560E-004)) < 1.e-14 );
|
|
assert( fabs(prim_vgl[26][2][7] - (-3.8250692897610380E-003)) < 1.e-14 );
|
|
assert( fabs(prim_vgl[26][3][7] - ( 3.4950559194080275E-003)) < 1.e-14 );
|
|
assert( fabs(prim_vgl[26][4][7] - ( 2.0392163767356572E-002)) < 1.e-14 );
|
|
|
|
}
|
|
|
|
|
|
#+end_src
|
|
|
|
*** Ideas for improvement :noexport:
|
|
|
|
#+begin_src c
|
|
// j : electrons
|
|
// l : primitives
|
|
|
|
k=0;
|
|
for (j=0 ; j<point_num ; ++j) {
|
|
for (i=0 ; i<nucl_num ; ++i) {
|
|
|
|
r2 = nucl_point_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
|
|
: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
|
|
| Variable | Type | In/Out | Description |
|
|
|---------------------+-----------------------------------+--------+----------------------------------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~prim_num~ | ~int64_t~ | in | Number of primitives |
|
|
| ~shell_num~ | ~int64_t~ | in | Number of shells |
|
|
| ~point_num~ | ~int64_t~ | in | Number of points |
|
|
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
|
|
| ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells for each nucleus |
|
|
| ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus |
|
|
| ~nucleus_range~ | ~double[nucl_num]~ | in | Range of the nucleus |
|
|
| ~shell_prim_index~ | ~int64_t[shell_num]~ | in | Index of the 1st primitive of each shell |
|
|
| ~shell_prim_num~ | ~int64_t[shell_num]~ | in | Number of primitives per shell |
|
|
| ~coord~ | ~double[3][point_num]~ | in | Coordinates |
|
|
| ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates |
|
|
| ~expo~ | ~double[prim_num]~ | in | Exponents of the primitives |
|
|
| ~coef_normalized~ | ~double[prim_num]~ | in | Coefficients of the primitives |
|
|
| ~shell_vgl~ | ~double[point_num][5][shell_num]~ | out | Value, gradients and Laplacian of the shells |
|
|
|
|
#+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_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 point_num,
|
|
const int64_t nucl_num,
|
|
const int64_t* nucleus_shell_num,
|
|
const int64_t* nucleus_index,
|
|
const double* nucleus_range,
|
|
const int64_t* shell_prim_index,
|
|
const int64_t* shell_prim_num,
|
|
const double* coord,
|
|
const double* nucl_coord,
|
|
const double* expo,
|
|
const double* coef_normalized,
|
|
double* const shell_vgl );
|
|
#+end_src
|
|
|
|
#+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, point_num, nucl_num, &
|
|
nucleus_shell_num, nucleus_index, nucleus_range, &
|
|
shell_prim_index, shell_prim_num, 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) :: point_num
|
|
integer*8 , intent(in) :: nucleus_shell_num(nucl_num)
|
|
integer*8 , intent(in) :: nucleus_index(nucl_num)
|
|
double precision , intent(in) :: nucleus_range(nucl_num)
|
|
integer*8 , intent(in) :: shell_prim_index(shell_num)
|
|
integer*8 , intent(in) :: shell_prim_num(shell_num)
|
|
double precision , intent(in) :: coord(point_num,3)
|
|
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,5,point_num)
|
|
|
|
integer*8 :: inucl, iprim, ipoint, ishell
|
|
integer*8 :: ishell_start, ishell_end
|
|
integer*8 :: iprim_start , iprim_end
|
|
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-12)
|
|
|
|
do ipoint = 1, point_num
|
|
|
|
do inucl=1,nucl_num
|
|
|
|
x = coord(ipoint,1) - nucl_coord(inucl,1)
|
|
y = coord(ipoint,2) - nucl_coord(inucl,2)
|
|
z = coord(ipoint,3) - nucl_coord(inucl,3)
|
|
|
|
r2 = x*x + y*y + z*z
|
|
|
|
if (r2 > cutoff*nucleus_range(inucl)) then
|
|
cycle
|
|
end if
|
|
|
|
! C is zero-based, so shift bounds by one
|
|
ishell_start = nucleus_index(inucl) + 1
|
|
ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl)
|
|
|
|
do ishell=ishell_start, ishell_end
|
|
|
|
shell_vgl(ishell, 1, ipoint) = 0.d0
|
|
shell_vgl(ishell, 2, ipoint) = 0.d0
|
|
shell_vgl(ishell, 3, ipoint) = 0.d0
|
|
shell_vgl(ishell, 4, ipoint) = 0.d0
|
|
shell_vgl(ishell, 5, ipoint) = 0.d0
|
|
|
|
iprim_start = shell_prim_index(ishell) + 1
|
|
iprim_end = shell_prim_index(ishell) + shell_prim_num(ishell)
|
|
|
|
do iprim = iprim_start, iprim_end
|
|
|
|
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, 1, ipoint) = &
|
|
shell_vgl(ishell, 1, ipoint) + v
|
|
|
|
shell_vgl(ishell, 2, ipoint) = &
|
|
shell_vgl(ishell, 2, ipoint) + two_a * x
|
|
|
|
shell_vgl(ishell, 3, ipoint) = &
|
|
shell_vgl(ishell, 3, ipoint) + two_a * y
|
|
|
|
shell_vgl(ishell, 4, ipoint) = &
|
|
shell_vgl(ishell, 4, ipoint) + two_a * z
|
|
|
|
shell_vgl(ishell, 5, ipoint) = &
|
|
shell_vgl(ishell, 5, ipoint) + two_a * (3.d0 - 2.d0*ar2)
|
|
|
|
end do
|
|
|
|
end do
|
|
end do
|
|
|
|
end do
|
|
|
|
end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
|
|
#+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, &
|
|
point_num, &
|
|
nucl_num, &
|
|
nucleus_shell_num, &
|
|
nucleus_index, &
|
|
nucleus_range, &
|
|
shell_prim_index, &
|
|
shell_prim_num, &
|
|
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 :: point_num
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num)
|
|
integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num)
|
|
real (c_double ) , intent(in) :: nucleus_range(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) :: coord(point_num,3)
|
|
real (c_double ) , intent(in) :: nucl_coord(nucl_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,5,point_num)
|
|
|
|
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, &
|
|
point_num, &
|
|
nucl_num, &
|
|
nucleus_shell_num, &
|
|
nucleus_index, &
|
|
nucleus_range, &
|
|
shell_prim_index, &
|
|
shell_prim_num, &
|
|
coord, &
|
|
nucl_coord, &
|
|
expo, &
|
|
coef_normalized, &
|
|
shell_vgl)
|
|
|
|
end function qmckl_compute_ao_basis_shell_gaussian_vgl
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
|
|
#+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_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_provide_ao_basis_shell_vgl",
|
|
NULL);
|
|
}
|
|
|
|
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_provide_ao_basis_shell_vgl",
|
|
NULL);
|
|
}
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->point.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->point.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->point.num,
|
|
ctx->nucleus.num,
|
|
ctx->ao_basis.nucleus_shell_num,
|
|
ctx->ao_basis.nucleus_index,
|
|
ctx->ao_basis.nucleus_range,
|
|
ctx->ao_basis.shell_prim_index,
|
|
ctx->ao_basis.shell_prim_num,
|
|
ctx->point.coord.data,
|
|
ctx->nucleus.coord.data,
|
|
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
|
|
|
|
*** Test :noexport:
|
|
|
|
#+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][point_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][26] : %25.15e"% f(a,x,y))
|
|
print ( "[1][1][26] : %25.15e"% df(a,x,y,1))
|
|
print ( "[1][2][26] : %25.15e"% df(a,x,y,2))
|
|
print ( "[1][3][26] : %25.15e"% df(a,x,y,3))
|
|
print ( "[1][4][26] : %25.15e"% lf(a,x,y))
|
|
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: [1][0][26] : 3.564393437193867e-02
|
|
: [1][1][26] : -6.030177988891605e-03
|
|
: [1][2][26] : -3.074832579871845e-02
|
|
: [1][3][26] : 2.809546963133958e-02
|
|
: [1][4][26] : 1.903338597841753e-02
|
|
|
|
#+begin_src c :tangle (eval c_test) :exports none
|
|
{
|
|
#define walk_num 1 // 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, walk_num*elec_num*3);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
|
|
double shell_vgl[elec_num][5][shell_num];
|
|
|
|
rc = qmckl_get_ao_basis_shell_vgl(context, &(shell_vgl[0][0][0]),
|
|
(int64_t) 5*elec_num*shell_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
printf(" shell_vgl[26][0][1] %25.15e\n", shell_vgl[26][0][1]);
|
|
printf(" shell_vgl[26][1][1] %25.15e\n", shell_vgl[26][1][1]);
|
|
printf(" shell_vgl[26][2][1] %25.15e\n", shell_vgl[26][2][1]);
|
|
printf(" shell_vgl[26][3][1] %25.15e\n", shell_vgl[26][3][1]);
|
|
printf(" shell_vgl[26][4][1] %25.15e\n", shell_vgl[26][4][1]);
|
|
|
|
assert( fabs(shell_vgl[26][0][1] - ( 3.564393437193868e-02)) < 1.e-14 );
|
|
assert( fabs(shell_vgl[26][1][1] - (-6.030177987072189e-03)) < 1.e-14 );
|
|
assert( fabs(shell_vgl[26][2][1] - (-3.074832579537582e-02)) < 1.e-14 );
|
|
assert( fabs(shell_vgl[26][3][1] - ( 2.809546963519935e-02)) < 1.e-14 );
|
|
assert( fabs(shell_vgl[26][4][1] - ( 1.896046117183968e-02)) < 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
|
|
| Variable | Type | In/Out | Description |
|
|
|-----------+-----------------+--------+---------------------------------------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~n~ | int64_t | in | Number of values |
|
|
| ~X~ | double[n] | in | Array containing the input values |
|
|
| ~LMAX~ | int32_t[n] | in | Array containing the maximum power for each value |
|
|
| ~P~ | double[n][ldp] | out | Array containing all the powers of ~X~ |
|
|
| ~ldp~ | int64_t | 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]~
|
|
|
|
#+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
|
|
|
|
#+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
|
|
|
|
#+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
|
|
|
|
#+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
|
|
|
|
*** Test :noexport:
|
|
|
|
#+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
|
|
|
|
** 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
|
|
| Variable | Type | In/Out | Description |
|
|
|-----------+-------------------+--------+------------------------------------------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~X~ | ~double[3]~ | in | Array containing the coordinates of the points |
|
|
| ~R~ | ~double[3]~ | in | Array containing the x,y,z coordinates of the center |
|
|
| ~lmax~ | ~int32_t~ | in | Maximum angular momentum |
|
|
| ~n~ | ~int64_t~ | inout | Number of computed polynomials |
|
|
| ~L~ | ~int32_t[n][ldl]~ | out | Contains a,b,c for all ~n~ results |
|
|
| ~ldl~ | ~int64_t~ | in | Leading dimension of ~L~ |
|
|
| ~VGL~ | ~double[n][ldv]~ | out | Value, gradients and Laplacian of the polynomials |
|
|
| ~ldv~ | ~int64_t~ | in | Leading dimension of array ~VGL~ |
|
|
|
|
Requirements:
|
|
|
|
- ~context~ \ne ~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"
|
|
|
|
#+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
|
|
|
|
#+CALL: generate_c_header(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_vgl_doc")
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_ao_polynomial_vgl_doc (
|
|
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
|
|
|
|
#+begin_src c :tangle (eval c) :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 )
|
|
{
|
|
#ifdef HAVE_HPC
|
|
//return qmckl_ao_polynomial_vgl_hpc (context, X, R, lmax, n, L, ldl, VGL, ldv);
|
|
return qmckl_ao_polynomial_vgl_doc (context, X, R, lmax, n, L, ldl, VGL, ldv);
|
|
#else
|
|
return qmckl_ao_polynomial_vgl_doc (context, X, R, lmax, n, L, ldl, VGL, ldv);
|
|
#endif
|
|
}
|
|
#+end_src
|
|
|
|
#+begin_src f90 :tangle (eval f)
|
|
integer function qmckl_ao_polynomial_vgl_doc_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)
|
|
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
|
|
|
|
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_doc_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_vgl_doc" )
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_ao_polynomial_vgl_doc &
|
|
(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_doc_f
|
|
info = qmckl_ao_polynomial_vgl_doc_f &
|
|
(context, X, R, lmax, n, L, ldl, VGL, ldv)
|
|
|
|
end function qmckl_ao_polynomial_vgl_doc
|
|
#+end_src
|
|
|
|
#+CALL: generate_f_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("FRetType"),fname="qmckl_ao_polynomial_vgl_doc" )
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(c_int32_t) function qmckl_ao_polynomial_vgl_doc &
|
|
(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_doc
|
|
end interface
|
|
#+end_src
|
|
|
|
#+CALL: generate_f_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("FRetType"),fname="qmckl_ao_polynomial_vgl" )
|
|
|
|
#+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
|
|
|
|
#+CALL: generate_c_header(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_transp_vgl")
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_ao_polynomial_transp_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
|
|
|
|
#+CALL: generate_c_header(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_transp_vgl_doc")
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_ao_polynomial_transp_vgl_doc (
|
|
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
|
|
|
|
#+CALL: generate_c_header(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_transp_vgl_hpc")
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_ao_polynomial_transp_vgl_hpc (
|
|
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
|
|
|
|
#+begin_src c :tangle (eval c) :comments org
|
|
qmckl_exit_code
|
|
qmckl_ao_polynomial_transp_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 )
|
|
{
|
|
#ifdef HAVE_HPC
|
|
return qmckl_ao_polynomial_transp_vgl_hpc (context, X, R, lmax, n, L, ldl, VGL, ldv);
|
|
#else
|
|
return qmckl_ao_polynomial_transp_vgl_doc (context, X, R, lmax, n, L, ldl, VGL, ldv);
|
|
#endif
|
|
}
|
|
#+end_src
|
|
|
|
#+begin_src f90 :tangle (eval f)
|
|
integer function qmckl_ao_polynomial_transp_vgl_doc_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,5)
|
|
integer*8 , intent(in) :: ldv
|
|
|
|
integer*8 :: i,j
|
|
integer :: a,b,c,d
|
|
real*8 :: Y(3)
|
|
real*8 :: pows(-2:21,3) ! lmax < 22
|
|
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 < (lmax+1)*(lmax+2)*(lmax+3)/6) then
|
|
info = QMCKL_INVALID_ARG_9
|
|
return
|
|
endif
|
|
|
|
|
|
if (lmax > 0) then
|
|
|
|
do i=1,3
|
|
Y(i) = X(i) - R(i)
|
|
end do
|
|
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
|
|
|
|
l (1:3,1:4) = 0
|
|
VGL(1:4,1:5) = 0.d0
|
|
|
|
VGL(1 ,1 ) = 1.d0
|
|
|
|
l (1,2) = 1
|
|
VGL(2,1) = Y(1)
|
|
VGL(2,2) = 1.d0
|
|
|
|
l (2,3) = 1
|
|
VGL(3,1) = Y(2)
|
|
VGL(3,3) = 1.d0
|
|
|
|
l (3,4) = 1
|
|
VGL(4,1) = Y(3)
|
|
VGL(4,4) = 1.d0
|
|
|
|
n=4
|
|
else
|
|
VGL(1,1) = 1.d0
|
|
VGL(1,2:5) = 0.d0
|
|
l(1:3,1) = 0
|
|
n=1
|
|
return
|
|
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
|
|
|
|
xy = pows(a,1) * pows(b,2)
|
|
yz = pows(b,2) * pows(c,3)
|
|
xz = pows(a,1) * pows(c,3)
|
|
|
|
l(1,n) = a
|
|
l(2,n) = b
|
|
l(3,n) = c
|
|
|
|
VGL(n,1) = xy * pows(c,3)
|
|
|
|
xy = dc * xy
|
|
xz = db * xz
|
|
yz = da * yz
|
|
|
|
VGL(n,2) = pows(a-1,1) * yz
|
|
VGL(n,3) = pows(b-1,2) * xz
|
|
VGL(n,4) = pows(c-1,3) * xy
|
|
|
|
VGL(n,5) = &
|
|
(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_transp_vgl_doc_f
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval c) :comments org
|
|
qmckl_exit_code
|
|
qmckl_ao_polynomial_transp_vgl_hpc (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 )
|
|
{
|
|
const qmckl_context_struct* ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL && X != NULL && R != NULL && n != NULL && L != NULL && VGL != NULL);
|
|
if (lmax < 0) return QMCKL_INVALID_ARG_4;
|
|
if (ldl < 3) return QMCKL_INVALID_ARG_7;
|
|
|
|
int32_t m;
|
|
|
|
switch (lmax) {
|
|
case 0:
|
|
{
|
|
if (ldv < 1) return QMCKL_INVALID_ARG_9;
|
|
L[0] = 0; L[1] = 0; L[2] = 0;
|
|
|
|
VGL[0 ] = 1.0;
|
|
VGL[ldv ] = 0.0;
|
|
VGL[ldv<<1 ] = 0.0;
|
|
VGL[(ldv<<1)+ldv] = 0.0;
|
|
VGL[ldv<<2 ] = 0.0;
|
|
m=1;
|
|
break;
|
|
}
|
|
case 1:
|
|
{
|
|
if (ldv < 4) return QMCKL_INVALID_ARG_9;
|
|
|
|
if (ldl == 3) {
|
|
|
|
const int32_t ltmp[12] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
|
|
for (int i=0 ; i<12 ; ++i)
|
|
L[i] = ltmp[i];
|
|
|
|
} else {
|
|
|
|
int32_t* const l[4] = {L, L+ldl, L+(ldl<<1), L+ldl+(ldl<<1)};
|
|
l[0][0] = 0; l[0][1] = 0; l[0][2] = 0;
|
|
l[1][0] = 1; l[1][1] = 0; l[1][2] = 0;
|
|
l[2][0] = 0; l[2][1] = 1; l[2][2] = 0;
|
|
l[3][0] = 0; l[3][1] = 0; l[3][2] = 1;
|
|
}
|
|
|
|
|
|
if (ldv == 4) {
|
|
|
|
const double tmp[20] = {1.0, X[0]-R[0], X[1]-R[1], X[2]-R[2],
|
|
0.0, 1.0, 0.0, 0.0,
|
|
0.0, 0.0, 1.0, 0.0,
|
|
0.0, 0.0, 0.0, 1.0,
|
|
0.0, 0.0, 0.0, 0.0};
|
|
|
|
for (int i=0 ; i<20 ; ++i)
|
|
VGL[i] = tmp[i];
|
|
|
|
} else {
|
|
|
|
double* const vgl1 = VGL;
|
|
double* const vgl2 = VGL + ldv;
|
|
double* const vgl3 = VGL + (ldv << 1);
|
|
double* const vgl4 = VGL + ldv + (ldv << 1);
|
|
double* const vgl5 = VGL + (ldv << 2);
|
|
|
|
for (int32_t k=0 ; k<4 ; ++k) {
|
|
vgl2[k] = 0.0;
|
|
vgl3[k] = 0.0;
|
|
vgl4[k] = 0.0;
|
|
vgl5[k] = 0.0;
|
|
}
|
|
vgl1[0] = 1.0;
|
|
vgl1[1] = X[0]-R[0];
|
|
vgl1[2] = X[1]-R[1];
|
|
vgl1[3] = X[2]-R[2];
|
|
vgl2[1] = 1.0;
|
|
vgl3[2] = 1.0;
|
|
vgl4[3] = 1.0;
|
|
}
|
|
m=4;
|
|
break;
|
|
}
|
|
case 2:
|
|
{
|
|
if (ldv < 10) return QMCKL_INVALID_ARG_9;
|
|
if (ldl == 3) {
|
|
const int32_t ltmp[30] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1,
|
|
2, 0, 0, 1, 1, 0, 1, 0, 1, 0, 2, 0,
|
|
0, 1, 1, 0, 0, 2};
|
|
for (int i=0 ; i<30 ; ++i)
|
|
L[i] = ltmp[i];
|
|
|
|
} else {
|
|
int32_t* l[10];
|
|
for (int32_t i=0 ; i<10 ; ++i) {
|
|
l[i] = L + i*ldl;
|
|
}
|
|
|
|
l[0][0] = 0; l[0][1] = 0; l[0][2] = 0;
|
|
l[1][0] = 1; l[1][1] = 0; l[1][2] = 0;
|
|
l[2][0] = 0; l[2][1] = 1; l[2][2] = 0;
|
|
l[3][0] = 0; l[3][1] = 0; l[3][2] = 1;
|
|
l[4][0] = 2; l[4][1] = 0; l[4][2] = 0;
|
|
l[5][0] = 1; l[5][1] = 1; l[5][2] = 0;
|
|
l[6][0] = 1; l[6][1] = 0; l[6][2] = 1;
|
|
l[7][0] = 0; l[7][1] = 2; l[7][2] = 0;
|
|
l[8][0] = 0; l[8][1] = 1; l[8][2] = 1;
|
|
l[9][0] = 0; l[9][1] = 0; l[9][2] = 2;
|
|
}
|
|
|
|
const double Y[3] = { X[0]-R[0],
|
|
X[1]-R[1],
|
|
X[2]-R[2] };
|
|
|
|
if (ldv == 50) {
|
|
const double tmp[50] = {
|
|
1.0, Y[0], Y[1], Y[2], Y[0] * Y[0],
|
|
Y[0] * Y[1], Y[0] * Y[2], Y[1] * Y[1], Y[1] * Y[2], Y[2] * Y[2],
|
|
0.0, 1.0, 0.0, 0.0, Y[0] + Y[0],
|
|
Y[1], Y[2], 0.0, 0.0, 0.0,
|
|
0.0, 0.0, 1.0, 0.0, 0.0,
|
|
Y[0], 0.0, Y[1] + Y[1], Y[2], 0.0,
|
|
0.0, 0.0, 0.0, 1.0, 0.0,
|
|
0.0, Y[0], 0.0, Y[1], Y[2] + Y[2],
|
|
0.0, 0.0, 0.0, 0.0, 2.0,
|
|
0.0, 0.0, 2.0, 0., 2.0 };
|
|
for (int i=0 ; i<50 ; ++i)
|
|
VGL[i] = tmp[i];
|
|
} else {
|
|
double* const vgl1 = VGL;
|
|
double* const vgl2 = VGL + ldv;
|
|
double* const vgl3 = VGL + (ldv << 1);
|
|
double* const vgl4 = VGL + 3*ldv;
|
|
double* const vgl5 = VGL + (ldv << 2);
|
|
|
|
vgl1[0] = 1.0 ; vgl1[1] = Y[0] ; vgl1[2] = Y[1];
|
|
vgl1[3] = Y[2] ; vgl1[4] = Y[0]*Y[0]; vgl1[5] = Y[0]*Y[1];
|
|
vgl1[6] = Y[0]*Y[2]; vgl1[7] = Y[1]*Y[1]; vgl1[8] = Y[1]*Y[2];
|
|
vgl1[9] = Y[2]*Y[2];
|
|
|
|
vgl2[0] = 0.0 ; vgl2[1] = 1.0 ; vgl2[2] = 0.0 ;
|
|
vgl2[3] = 0.0 ; vgl2[4] = Y[0]+Y[0]; vgl2[5] = Y[1];
|
|
vgl2[6] = Y[2]; vgl2[7] = 0.0 ; vgl2[8] = 0.0 ;
|
|
vgl2[9] = 0.0 ;
|
|
|
|
vgl3[0] = 0.0; vgl3[1] = 0.0 ; vgl3[2] = 1.0 ;
|
|
vgl3[3] = 0.0; vgl3[4] = 0.0 ; vgl3[5] = Y[0];
|
|
vgl3[6] = 0.0; vgl3[7] = Y[1]+Y[1]; vgl3[8] = Y[2];
|
|
vgl3[9] = 0.0;
|
|
|
|
vgl4[0] = 0.0 ; vgl4[1] = 0.0; vgl4[2] = 0.0 ;
|
|
vgl4[3] = 1.0 ; vgl4[4] = 0.0; vgl4[5] = 0.0 ;
|
|
vgl4[6] = Y[0] ; vgl4[7] = 0.0; vgl4[8] = Y[1];
|
|
vgl4[9] = Y[2]+Y[2];
|
|
|
|
vgl5[0] = 0.0; vgl5[1] = 0.0; vgl5[2] = 0.0;
|
|
vgl5[3] = 0.0; vgl5[4] = 2.0; vgl5[5] = 0.0;
|
|
vgl5[6] = 0.0; vgl5[7] = 2.0; vgl5[8] = 0.0;
|
|
vgl5[9] = 2.0;
|
|
}
|
|
m=10;
|
|
break;
|
|
}
|
|
default:
|
|
{
|
|
const int32_t size_max = (lmax+1)*(lmax+2)*(lmax+3)/6;
|
|
if (ldv < size_max) return QMCKL_INVALID_ARG_9;
|
|
|
|
double* const vgl1 = VGL;
|
|
double* const vgl2 = VGL + ldv;
|
|
double* const vgl3 = VGL + (ldv<<1);
|
|
double* const vgl4 = VGL + ldv + (ldv<<1);
|
|
double* const vgl5 = VGL + (ldv<<2);
|
|
|
|
const double Y[3] = { X[0]-R[0],
|
|
X[1]-R[1],
|
|
X[2]-R[2] };
|
|
|
|
assert(size_max > lmax+3);
|
|
double pows[3][size_max];
|
|
|
|
for (int32_t i=0 ; i<3 ; ++i) {
|
|
pows[0][i] = 1.0;
|
|
pows[1][i] = 1.0;
|
|
pows[2][i] = 1.0;
|
|
}
|
|
|
|
for (int32_t i=3 ; i<=lmax+2 ; ++i) {
|
|
pows[0][i] = pows[0][i-1] * Y[0];
|
|
pows[1][i] = pows[1][i-1] * Y[1];
|
|
pows[2][i] = pows[2][i-1] * Y[2];
|
|
}
|
|
|
|
int32_t* l[size_max];
|
|
for (int32_t i=0 ; i<size_max ; ++i) {
|
|
l[i] = &(L[i*ldl]);
|
|
}
|
|
|
|
for (int32_t i=0 ; i<4 ; ++i) {
|
|
l[i][0] = 0;
|
|
l[i][1] = 0;
|
|
l[i][2] = 0;
|
|
}
|
|
l[1][0] = 1;
|
|
l[2][1] = 1;
|
|
l[3][2] = 1;
|
|
|
|
for (int32_t k=0 ; k<4 ; ++k) {
|
|
vgl2[k] = 0.0;
|
|
vgl3[k] = 0.0;
|
|
vgl4[k] = 0.0;
|
|
vgl5[k] = 0.0;
|
|
}
|
|
vgl1[0] = 1.0;
|
|
vgl1[1] = Y[0];
|
|
vgl1[2] = Y[1];
|
|
vgl1[3] = Y[2];
|
|
|
|
vgl2[1] = 1.0;
|
|
vgl3[2] = 1.0;
|
|
vgl4[3] = 1.0;
|
|
m=4;
|
|
|
|
double dd = 2.0;
|
|
|
|
for (int32_t d=2 ; d<= lmax ; ++d) {
|
|
double da = dd;
|
|
|
|
for (int32_t a=d ; a>=0 ; --a) {
|
|
double db = dd-da;
|
|
|
|
for (int32_t b=d-a ; b>=0 ; --b) {
|
|
const int32_t c = d - a - b;
|
|
const double dc = dd - da - db;
|
|
|
|
double xy = pows[0][a+2] * pows[1][b+2];
|
|
double yz = pows[1][b+2] * pows[2][c+2];
|
|
double xz = pows[0][a+2] * pows[2][c+2];
|
|
|
|
l[m][0] = a;
|
|
l[m][1] = b;
|
|
l[m][2] = c;
|
|
|
|
vgl1[m] = xy * pows[2][c+2];
|
|
|
|
xy *= dc;
|
|
xz *= db;
|
|
yz *= da;
|
|
|
|
vgl2[m] = pows[0][a+1] * yz;
|
|
vgl3[m] = pows[1][b+1] * xz;
|
|
vgl4[m] = pows[2][c+1] * xy;
|
|
|
|
vgl5[m] = (da-1.) * pows[0][a] * yz +
|
|
(db-1.) * pows[1][b] * xz +
|
|
(dc-1.) * pows[2][c] * xy;
|
|
|
|
db -= 1.0;
|
|
++m;
|
|
}
|
|
da -= 1.0;
|
|
}
|
|
dd += 1.0;
|
|
}
|
|
}
|
|
}
|
|
,*n = m;
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_transp_vgl_doc")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_ao_polynomial_transp_vgl_doc &
|
|
(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_transp_vgl_doc_f
|
|
info = qmckl_ao_polynomial_transp_vgl_doc_f &
|
|
(context, X, R, lmax, n, L, ldl, VGL, ldv)
|
|
|
|
end function qmckl_ao_polynomial_transp_vgl_doc
|
|
#+end_src
|
|
|
|
#+CALL: generate_f_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("FRetType"),fname="qmckl_ao_polynomial_transp_vgl_doc")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(c_int32_t) function qmckl_ao_polynomial_transp_vgl_doc &
|
|
(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_transp_vgl_doc
|
|
end interface
|
|
#+end_src
|
|
|
|
#+CALL: generate_f_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("FRetType"),fname="qmckl_ao_polynomial_transp_vgl")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(c_int32_t) function qmckl_ao_polynomial_transp_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_transp_vgl
|
|
end interface
|
|
#+end_src
|
|
|
|
*** Test :noexport:
|
|
|
|
#+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));
|
|
|
|
double X[3] = { 1.1, 2.2, 3.3 };
|
|
double R[3] = { 0.2, 1.1, 3.0 };
|
|
int32_t ldv[8] = {1, 4, 10, 20, 35, 56, 84, 120};
|
|
for (int32_t ldl=3 ; ldl<=5 ; ++ldl) {
|
|
int64_t n;
|
|
int32_t L0[200][ldl];
|
|
int32_t L1[200][ldl];
|
|
printf("ldl=%d\n", ldl);
|
|
for (int32_t lmax=0 ; lmax<=7 ; lmax++) {
|
|
double VGL0[5][ldv[lmax]];
|
|
double VGL1[5][ldv[lmax]];
|
|
memset(&L0[0][0], 0, sizeof(L0));
|
|
memset(&L1[0][0], 0, sizeof(L1));
|
|
memset(&VGL0[0][0], 0, sizeof(VGL0));
|
|
memset(&VGL1[0][0], 0, sizeof(VGL1));
|
|
rc = qmckl_ao_polynomial_transp_vgl_doc (context, X, R, lmax, &n, &(L0[0][0]), ldl, &(VGL0[0][0]), ldv[lmax]);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
rc = qmckl_ao_polynomial_transp_vgl_hpc (context, X, R, lmax, &n, &(L1[0][0]), ldl, &(VGL1[0][0]), ldv[lmax]);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
printf("lmax=%d\n", lmax);
|
|
for (int32_t l=0 ; l<n ; ++l) {
|
|
for (int32_t k=0 ; k<3 ; ++k) {
|
|
printf("L[%d][%d] = %d %d\n", l, k, L0[l][k], L1[l][k]);
|
|
assert( L0[l][k] == L1[l][k] );
|
|
}
|
|
}
|
|
|
|
for (int32_t k=0 ; k<5 ; ++k) {
|
|
for (int32_t l=0 ; l<n ; ++l) {
|
|
printf("VGL[%d][%d] = %e %e %e\n", k, l, VGL0[k][l], VGL1[k][l], VGL0[k][l]-VGL1[k][l]);
|
|
assert( fabs(1.-(VGL0[k][l]+1.e-100)/(VGL1[k][l]+1.e-100)) < 1.e-15 );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
#+end_src
|
|
|
|
* Combining radial and polynomial parts
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_ao_vgl
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
** Unoptimized version
|
|
#+NAME: qmckl_ao_vgl_args_doc
|
|
| Variable | Type | In/Out | Description |
|
|
|-----------------------+-----------------------------------+--------+----------------------------------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~ao_num~ | ~int64_t~ | in | Number of AOs |
|
|
| ~shell_num~ | ~int64_t~ | in | Number of shells |
|
|
| ~point_num~ | ~int64_t~ | in | Number of points |
|
|
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
|
|
| ~coord~ | ~double[3][point_num]~ | in | Coordinates |
|
|
| ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates |
|
|
| ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus |
|
|
| ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus |
|
|
| ~nucleus_range~ | ~double[nucl_num]~ | in | Range beyond which all is zero |
|
|
| ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus |
|
|
| ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell |
|
|
| ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs |
|
|
| ~shell_vgl~ | ~double[point_num][5][shell_num]~ | in | Value, gradients and Laplacian of the shells |
|
|
| ~ao_vgl~ | ~double[point_num][5][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_doc_f(context, &
|
|
ao_num, shell_num, point_num, nucl_num, &
|
|
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) :: point_num
|
|
integer*8 , intent(in) :: nucl_num
|
|
double precision , intent(in) :: coord(point_num,3)
|
|
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,5,point_num)
|
|
double precision , intent(out) :: ao_vgl(ao_num,5,point_num)
|
|
|
|
double precision :: e_coord(3), n_coord(3)
|
|
integer*8 :: n_poly
|
|
integer :: l, il, k
|
|
integer*8 :: ipoint, inucl, ishell
|
|
integer*8 :: ishell_start, ishell_end
|
|
integer :: lstart(0:20)
|
|
double precision :: x, y, z, r2
|
|
double precision :: cutoff
|
|
integer, external :: qmckl_ao_polynomial_vgl_doc_f
|
|
|
|
double precision, allocatable :: poly_vgl(:,:)
|
|
integer , allocatable :: powers(:,:), ao_index(:)
|
|
|
|
allocate(poly_vgl(5,ao_num), powers(3,ao_num), ao_index(ao_num))
|
|
|
|
! Pre-computed data
|
|
do l=0,20
|
|
lstart(l) = l*(l+1)*(l+2)/6 +1
|
|
end do
|
|
|
|
k=1
|
|
do inucl=1,nucl_num
|
|
ishell_start = nucleus_index(inucl) + 1
|
|
ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl)
|
|
do ishell = ishell_start, ishell_end
|
|
l = shell_ang_mom(ishell)
|
|
ao_index(ishell) = k
|
|
k = k + lstart(l+1) - lstart(l)
|
|
end do
|
|
end do
|
|
info = QMCKL_SUCCESS
|
|
|
|
! Don't compute polynomials when the radial part is zero.
|
|
cutoff = -dlog(1.d-12)
|
|
|
|
do ipoint = 1, point_num
|
|
e_coord(1) = coord(ipoint,1)
|
|
e_coord(2) = coord(ipoint,2)
|
|
e_coord(3) = coord(ipoint,3)
|
|
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 point 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 + y*y + z*z
|
|
|
|
if (r2 > cutoff*nucleus_range(inucl)) then
|
|
cycle
|
|
end if
|
|
|
|
! Compute polynomials
|
|
info = qmckl_ao_polynomial_vgl_doc_f(context, e_coord, n_coord, &
|
|
nucleus_max_ang_mom(inucl), n_poly, powers, 3_8, &
|
|
poly_vgl, 5_8)
|
|
|
|
! Loop over shells
|
|
ishell_start = nucleus_index(inucl) + 1
|
|
ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl)
|
|
do ishell = ishell_start, ishell_end
|
|
k = ao_index(ishell)
|
|
l = shell_ang_mom(ishell)
|
|
do il = lstart(l), lstart(l+1)-1
|
|
! Value
|
|
ao_vgl(k,1,ipoint) = &
|
|
poly_vgl(1,il) * shell_vgl(ishell,1,ipoint) * ao_factor(k)
|
|
|
|
! Grad_x
|
|
ao_vgl(k,2,ipoint) = ( &
|
|
poly_vgl(2,il) * shell_vgl(ishell,1,ipoint) + &
|
|
poly_vgl(1,il) * shell_vgl(ishell,2,ipoint) &
|
|
) * ao_factor(k)
|
|
|
|
! Grad_y
|
|
ao_vgl(k,3,ipoint) = ( &
|
|
poly_vgl(3,il) * shell_vgl(ishell,1,ipoint) + &
|
|
poly_vgl(1,il) * shell_vgl(ishell,3,ipoint) &
|
|
) * ao_factor(k)
|
|
|
|
! Grad_z
|
|
ao_vgl(k,4,ipoint) = ( &
|
|
poly_vgl(4,il) * shell_vgl(ishell,1,ipoint) + &
|
|
poly_vgl(1,il) * shell_vgl(ishell,4,ipoint) &
|
|
) * ao_factor(k)
|
|
|
|
! Lapl_z
|
|
ao_vgl(k,5,ipoint) = ( &
|
|
poly_vgl(5,il) * shell_vgl(ishell,1,ipoint) + &
|
|
poly_vgl(1,il) * shell_vgl(ishell,5,ipoint) + &
|
|
2.d0 * ( &
|
|
poly_vgl(2,il) * shell_vgl(ishell,2,ipoint) + &
|
|
poly_vgl(3,il) * shell_vgl(ishell,3,ipoint) + &
|
|
poly_vgl(4,il) * shell_vgl(ishell,4,ipoint) ) &
|
|
) * ao_factor(k)
|
|
|
|
k = k+1
|
|
end do
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
deallocate(poly_vgl, powers)
|
|
end function qmckl_compute_ao_vgl_doc_f
|
|
#+end_src
|
|
|
|
** HPC version
|
|
#+NAME: qmckl_ao_vgl_args_hpc_gaussian
|
|
| Variable | Type | In/Out | Description |
|
|
|-----------------------+--------------------------------+--------+----------------------------------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~ao_num~ | ~int64_t~ | in | Number of AOs |
|
|
| ~shell_num~ | ~int64_t~ | in | Number of shells |
|
|
| ~prim_num~ | ~int64_t~ | in | Number of primitives |
|
|
| ~point_num~ | ~int64_t~ | in | Number of points |
|
|
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
|
|
| ~coord~ | ~double[3][point_num]~ | in | Coordinates |
|
|
| ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates |
|
|
| ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus |
|
|
| ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus |
|
|
| ~nucleus_range~ | ~double[nucl_num]~ | in | Range beyond which all is zero |
|
|
| ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus |
|
|
| ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell |
|
|
| ~shell_prim_index~ | ~int64_t[shell_num]~ | in | Index of the 1st primitive of each shell |
|
|
| ~shell_prim_num~ | ~int64_t[shell_num]~ | in | Number of primitives per shell |
|
|
| ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs |
|
|
| ~ao_expo~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells |
|
|
| ~coef_normalized~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells |
|
|
| ~ao_vgl~ | ~double[point_num][5][ao_num]~ | out | Value, gradients and Laplacian of the AOs |
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
#ifdef HAVE_HPC
|
|
qmckl_exit_code
|
|
qmckl_compute_ao_vgl_hpc_gaussian (
|
|
const qmckl_context context,
|
|
const int64_t ao_num,
|
|
const int64_t shell_num,
|
|
const int64_t prim_num,
|
|
const int64_t point_num,
|
|
const int64_t nucl_num,
|
|
const double* 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 int64_t* shell_prim_index,
|
|
const int64_t* shell_prim_num,
|
|
const double* ao_factor,
|
|
const double* expo,
|
|
const double* coef_normalized,
|
|
double* const ao_vgl )
|
|
{
|
|
int32_t lstart[32];
|
|
for (int32_t l=0 ; l<32 ; ++l) {
|
|
lstart[l] = l*(l+1)*(l+2)/6;
|
|
}
|
|
|
|
int64_t ao_index[shell_num+1];
|
|
int64_t size_max = 0;
|
|
{
|
|
int64_t k=0;
|
|
for (int inucl=0 ; inucl < nucl_num ; ++inucl) {
|
|
const int64_t ishell_start = nucleus_index[inucl];
|
|
const int64_t ishell_end = nucleus_index[inucl] + nucleus_shell_num[inucl];
|
|
for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {
|
|
const int l = shell_ang_mom[ishell];
|
|
ao_index[ishell] = k;
|
|
k += lstart[l+1] - lstart[l];
|
|
size_max = size_max < lstart[l+1] ? lstart[l+1] : size_max;
|
|
}
|
|
}
|
|
ao_index[shell_num] = ao_num+1;
|
|
}
|
|
|
|
/* Don't compute polynomials when the radial part is zero. */
|
|
double cutoff = -log(1.e-12);
|
|
|
|
#ifdef HAVE_OPENMP
|
|
#pragma omp parallel
|
|
#endif
|
|
{
|
|
qmckl_exit_code rc;
|
|
double c_[prim_num];
|
|
double expo_[prim_num];
|
|
double ar2[prim_num];
|
|
int32_t powers[prim_num];
|
|
double poly_vgl_l1[4][4] = {{1.0, 0.0, 0.0, 0.0},
|
|
{0.0, 1.0, 0.0, 0.0},
|
|
{0.0, 0.0, 1.0, 0.0},
|
|
{0.0, 0.0, 0.0, 1.0}};
|
|
double poly_vgl_l2[5][10] = {{1., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
|
|
{0., 1., 0., 0., 0., 0., 0., 0., 0., 0.},
|
|
{0., 0., 1., 0., 0., 0., 0., 0., 0., 0.},
|
|
{0., 0., 0., 1., 0., 0., 0., 0., 0., 0.},
|
|
{0., 0., 0., 0., 2., 0., 0., 2., 0., 2.}};
|
|
double poly_vgl[5][size_max];
|
|
|
|
#ifdef HAVE_OPENMP
|
|
#pragma omp for
|
|
#endif
|
|
for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
|
|
const double e_coord[3] = { coord[ipoint],
|
|
coord[ipoint + point_num],
|
|
coord[ipoint + 2*point_num] };
|
|
|
|
for (int64_t inucl=0 ; inucl < nucl_num ; ++inucl) {
|
|
const double n_coord[3] = { nucl_coord[inucl],
|
|
nucl_coord[inucl + nucl_num],
|
|
nucl_coord[inucl + 2*nucl_num] };
|
|
|
|
/* Test if the point is in the range of the nucleus */
|
|
const double x = e_coord[0] - n_coord[0];
|
|
const double y = e_coord[1] - n_coord[1];
|
|
const double z = e_coord[2] - n_coord[2];
|
|
|
|
const double r2 = x*x + y*y + z*z;
|
|
|
|
if (r2 > cutoff * nucleus_range[inucl]) {
|
|
continue;
|
|
}
|
|
|
|
int64_t n_poly;
|
|
switch (nucleus_max_ang_mom[inucl]) {
|
|
case 0:
|
|
break;
|
|
|
|
case 1:
|
|
poly_vgl_l1[0][1] = x;
|
|
poly_vgl_l1[0][2] = y;
|
|
poly_vgl_l1[0][3] = z;
|
|
break;
|
|
|
|
case 2:
|
|
poly_vgl_l2[0][1] = x;
|
|
poly_vgl_l2[0][2] = y;
|
|
poly_vgl_l2[0][3] = z;
|
|
poly_vgl_l2[0][4] = x*x;
|
|
poly_vgl_l2[0][5] = x*y;
|
|
poly_vgl_l2[0][6] = x*z;
|
|
poly_vgl_l2[0][7] = y*y;
|
|
poly_vgl_l2[0][8] = y*z;
|
|
poly_vgl_l2[0][9] = z*z;
|
|
poly_vgl_l2[1][4] = x+x;
|
|
poly_vgl_l2[1][5] = y;
|
|
poly_vgl_l2[1][6] = z;
|
|
poly_vgl_l2[2][5] = x;
|
|
poly_vgl_l2[2][7] = y+y;
|
|
poly_vgl_l2[2][8] = z;
|
|
poly_vgl_l2[3][6] = x;
|
|
poly_vgl_l2[3][8] = y;
|
|
poly_vgl_l2[3][9] = z+z;
|
|
break;
|
|
|
|
default:
|
|
rc = qmckl_ao_polynomial_transp_vgl_hpc(context, e_coord, n_coord,
|
|
nucleus_max_ang_mom[inucl],
|
|
&n_poly, powers, (int64_t) 3,
|
|
&(poly_vgl[0][0]), size_max);
|
|
|
|
assert (rc == QMCKL_SUCCESS);
|
|
break;
|
|
}
|
|
|
|
const int64_t ishell_start = nucleus_index[inucl];
|
|
const int64_t ishell_end = nucleus_index[inucl] + nucleus_shell_num[inucl];
|
|
for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {
|
|
|
|
int64_t nidx = 0;
|
|
const int64_t iprim_start = shell_prim_index[ishell];
|
|
const int64_t iprim_end = shell_prim_index[ishell] + shell_prim_num[ishell];
|
|
for (int64_t iprim = iprim_start ; iprim < iprim_end ; ++iprim) {
|
|
const double v = expo[iprim]*r2;
|
|
if (v <= cutoff) {
|
|
ar2[nidx] = v;
|
|
c_[nidx] = coef_normalized[iprim];
|
|
expo_[nidx] = expo[iprim];
|
|
++nidx;
|
|
}
|
|
}
|
|
|
|
double s1_ = 0.;
|
|
double s5_ = 0.;
|
|
double s6_ = 0.;
|
|
for (int idx=0 ; idx<nidx ; ++idx) {
|
|
const double v = c_[idx] * exp(-ar2[idx]);
|
|
const double t = v*expo_[idx];
|
|
s1_ += v;
|
|
s6_ -= t;
|
|
s5_ += t*(4.0*ar2[idx]-6.0);
|
|
}
|
|
s6_ += s6_;
|
|
const double s1 = s1_;
|
|
const double s2 = s6_*x;
|
|
const double s3 = s6_*y;
|
|
const double s4 = s6_*z;
|
|
const double s5 = s5_;
|
|
|
|
const int64_t k = ao_index[ishell];
|
|
double* __restrict__ const ao_vgl_1 = ao_vgl + ipoint*5*ao_num + k;
|
|
|
|
const int32_t l = shell_ang_mom[ishell];
|
|
const int32_t n = lstart[l+1]-lstart[l];
|
|
|
|
double* __restrict__ const ao_vgl_2 = ao_vgl_1 + ao_num;
|
|
double* __restrict__ const ao_vgl_3 = ao_vgl_1 + (ao_num<<1);
|
|
double* __restrict__ const ao_vgl_4 = ao_vgl_1 + (ao_num<<1) + ao_num;
|
|
double* __restrict__ const ao_vgl_5 = ao_vgl_1 + (ao_num<<2);
|
|
|
|
double* __restrict__ poly_vgl_1;
|
|
double* __restrict__ poly_vgl_2;
|
|
double* __restrict__ poly_vgl_3;
|
|
double* __restrict__ poly_vgl_4;
|
|
double* __restrict__ poly_vgl_5;
|
|
if (nidx > 0) {
|
|
const double* __restrict__ f = ao_factor + k;
|
|
const int64_t idx = lstart[l];
|
|
|
|
switch (nucleus_max_ang_mom[inucl]) {
|
|
case 0:
|
|
break;
|
|
case 1:
|
|
poly_vgl_1 = &(poly_vgl_l1[0][idx]);
|
|
poly_vgl_2 = &(poly_vgl_l1[1][idx]);
|
|
poly_vgl_3 = &(poly_vgl_l1[2][idx]);
|
|
poly_vgl_4 = &(poly_vgl_l1[3][idx]);
|
|
break;
|
|
case 2:
|
|
poly_vgl_1 = &(poly_vgl_l2[0][idx]);
|
|
poly_vgl_2 = &(poly_vgl_l2[1][idx]);
|
|
poly_vgl_3 = &(poly_vgl_l2[2][idx]);
|
|
poly_vgl_4 = &(poly_vgl_l2[3][idx]);
|
|
break;
|
|
default:
|
|
poly_vgl_1 = &(poly_vgl[0][idx]);
|
|
poly_vgl_2 = &(poly_vgl[1][idx]);
|
|
poly_vgl_3 = &(poly_vgl[2][idx]);
|
|
poly_vgl_4 = &(poly_vgl[3][idx]);
|
|
poly_vgl_5 = &(poly_vgl[4][idx]);
|
|
}
|
|
switch (n) {
|
|
case(1):
|
|
ao_vgl_1[0] = s1 * f[0];
|
|
ao_vgl_2[0] = s2 * f[0];
|
|
ao_vgl_3[0] = s3 * f[0];
|
|
ao_vgl_4[0] = s4 * f[0];
|
|
ao_vgl_5[0] = s5;
|
|
break;
|
|
case (3):
|
|
#ifdef HAVE_OPENMP
|
|
#pragma omp simd
|
|
#endif
|
|
for (int il=0 ; il<3 ; ++il) {
|
|
ao_vgl_1[il] = poly_vgl_1[il] * s1 * f[il];
|
|
ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
|
|
ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il];
|
|
ao_vgl_4[il] = (poly_vgl_4[il] * s1 + poly_vgl_1[il] * s4) * f[il];
|
|
ao_vgl_5[il] = (poly_vgl_1[il] * s5 +
|
|
2.0*(poly_vgl_2[il] * s2 +
|
|
poly_vgl_3[il] * s3 +
|
|
poly_vgl_4[il] * s4 )) * f[il];
|
|
}
|
|
break;
|
|
case(5):
|
|
#ifdef HAVE_OPENMP
|
|
#pragma omp simd
|
|
#endif
|
|
for (int il=0 ; il<5 ; ++il) {
|
|
ao_vgl_1[il] = poly_vgl_1[il] * s1 * f[il];
|
|
ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
|
|
ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il];
|
|
ao_vgl_4[il] = (poly_vgl_4[il] * s1 + poly_vgl_1[il] * s4) * f[il];
|
|
ao_vgl_5[il] = (poly_vgl_1[il] * s5 +
|
|
2.0*(poly_vgl_2[il] * s2 +
|
|
poly_vgl_3[il] * s3 +
|
|
poly_vgl_4[il] * s4 )) * f[il];
|
|
}
|
|
break;
|
|
default:
|
|
#ifdef HAVE_OPENMP
|
|
#pragma omp simd simdlen(8)
|
|
#endif
|
|
for (int il=0 ; il<n ; ++il) {
|
|
ao_vgl_1[il] = poly_vgl_1[il] * s1 * f[il];
|
|
ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
|
|
ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il];
|
|
ao_vgl_4[il] = (poly_vgl_4[il] * s1 + poly_vgl_1[il] * s4) * f[il];
|
|
ao_vgl_5[il] = (poly_vgl_5[il] * s1 + poly_vgl_1[il] * s5 +
|
|
2.0*(poly_vgl_2[il] * s2 +
|
|
poly_vgl_3[il] * s3 +
|
|
poly_vgl_4[il] * s4 )) * f[il];
|
|
}
|
|
break;
|
|
}
|
|
} else {
|
|
for (int64_t il=0 ; il<n ; ++il) {
|
|
ao_vgl_1[il] = 0.0;
|
|
ao_vgl_2[il] = 0.0;
|
|
ao_vgl_3[il] = 0.0;
|
|
ao_vgl_4[il] = 0.0;
|
|
ao_vgl_5[il] = 0.0;
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
}
|
|
}
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#endif
|
|
#+end_src
|
|
|
|
** Interfaces
|
|
# #+CALL: generate_c_header(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl"))
|
|
# (Commented because the header needs to go into h_private_func)
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_private_func) :comments org
|
|
qmckl_exit_code qmckl_compute_ao_vgl_doc (
|
|
const qmckl_context context,
|
|
const int64_t ao_num,
|
|
const int64_t shell_num,
|
|
const int64_t point_num,
|
|
const int64_t nucl_num,
|
|
const double* 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
|
|
#+begin_src c :tangle (eval h_private_func) :comments org
|
|
#ifdef HAVE_HPC
|
|
qmckl_exit_code qmckl_compute_ao_vgl_hpc_gaussian (
|
|
const qmckl_context context,
|
|
const int64_t ao_num,
|
|
const int64_t shell_num,
|
|
const int64_t prim_num,
|
|
const int64_t point_num,
|
|
const int64_t nucl_num,
|
|
const double* 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 int64_t* shell_prim_index,
|
|
const int64_t* shell_prim_num,
|
|
const double* ao_factor,
|
|
const double* expo,
|
|
const double* coef_normalized,
|
|
double* const ao_vgl );
|
|
#endif
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl_doc"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_ao_vgl_doc &
|
|
(context, &
|
|
ao_num, &
|
|
shell_num, &
|
|
point_num, &
|
|
nucl_num, &
|
|
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 :: point_num
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
real (c_double ) , intent(in) :: coord(point_num,3)
|
|
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,5,point_num)
|
|
real (c_double ) , intent(out) :: ao_vgl(ao_num,5,point_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_ao_vgl_doc_f
|
|
info = qmckl_compute_ao_vgl_doc_f &
|
|
(context, &
|
|
ao_num, &
|
|
shell_num, &
|
|
point_num, &
|
|
nucl_num, &
|
|
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_doc
|
|
#+end_src
|
|
|
|
|
|
*** Provide :noexport:
|
|
|
|
#+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_failwith( context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_provide_ao_vgl",
|
|
NULL);
|
|
}
|
|
|
|
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);
|
|
}
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->point.date > ctx->ao_basis.ao_vgl_date) {
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
/* Provide required data */
|
|
#ifndef HAVE_HPC
|
|
rc = qmckl_provide_ao_basis_shell_vgl(context);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_shell_vgl", NULL);
|
|
}
|
|
#endif
|
|
|
|
/* 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->point.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;
|
|
}
|
|
#ifdef HAVE_HPC
|
|
if (ctx->ao_basis.type == 'G') {
|
|
rc = qmckl_compute_ao_vgl_hpc_gaussian(context,
|
|
ctx->ao_basis.ao_num,
|
|
ctx->ao_basis.shell_num,
|
|
ctx->ao_basis.prim_num,
|
|
ctx->point.num,
|
|
ctx->nucleus.num,
|
|
ctx->point.coord.data,
|
|
ctx->nucleus.coord.data,
|
|
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.shell_prim_index,
|
|
ctx->ao_basis.shell_prim_num,
|
|
ctx->ao_basis.ao_factor,
|
|
ctx->ao_basis.exponent,
|
|
ctx->ao_basis.coefficient_normalized,
|
|
ctx->ao_basis.ao_vgl);
|
|
/*
|
|
} else if (ctx->ao_basis.type == 'S') {
|
|
rc = qmck_compute_ao_vgl_hpc_slater(context,
|
|
ctx->ao_basis.ao_num,
|
|
ctx->ao_basis.shell_num,
|
|
ctx->ao_basis.prim_num,
|
|
ctx->point.num,
|
|
ctx->nucleus.num,
|
|
ctx->point.coord.data,
|
|
ctx->nucleus.coord.data,
|
|
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.shell_prim_index,
|
|
ctx->ao_basis.shell_prim_num,
|
|
ctx->ao_basis.ao_factor,
|
|
ctx->ao_basis.exponent,
|
|
ctx->ao_basis.coefficient_normalized,
|
|
ctx->ao_basis.ao_vgl);
|
|
,*/
|
|
} else {
|
|
rc = qmckl_compute_ao_vgl_doc(context,
|
|
ctx->ao_basis.ao_num,
|
|
ctx->ao_basis.shell_num,
|
|
ctx->point.num,
|
|
ctx->nucleus.num,
|
|
ctx->point.coord.data,
|
|
ctx->nucleus.coord.data,
|
|
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);
|
|
}
|
|
#else
|
|
rc = qmckl_compute_ao_vgl_doc(context,
|
|
ctx->ao_basis.ao_num,
|
|
ctx->ao_basis.shell_num,
|
|
ctx->point.num,
|
|
ctx->nucleus.num,
|
|
ctx->point.coord.data,
|
|
ctx->nucleus.coord.data,
|
|
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);
|
|
#endif
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->ao_basis.ao_vgl_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Test :noexport:
|
|
|
|
#+begin_src python :results output :exports none
|
|
import numpy as np
|
|
from math import sqrt
|
|
|
|
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 ao_vgl[prim_num][5][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)]
|
|
|
|
norm = sqrt(3.)
|
|
print ( "[0][26][219] : %25.15e"%(f(a,x,y) * (x[0] - y[0])**2) )
|
|
print ( "[1][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][26][220] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[1] - y[1]) ))
|
|
print ( "[1][26][220] : %25.15e"%(norm*df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + norm*f(a,x,y) * (x[1] - y[1])) )
|
|
|
|
print ( "[0][26][221] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[2] - y[2])) )
|
|
print ( "[1][26][221] : %25.15e"%(norm*df(a,x,y,1)* (x[0] - y[0]) * (x[2] - y[2]) + norm*f(a,x,y) * (x[2] - y[2])) )
|
|
|
|
print ( "[0][26][222] : %25.15e"%(f(a,x,y) * (x[1] - y[1]) * (x[1] - y[1])) )
|
|
print ( "[1][26][222] : %25.15e"%(df(a,x,y,1)* (x[1] - y[1]) * (x[1] - y[1])) )
|
|
|
|
print ( "[0][26][223] : %25.15e"%(norm*f(a,x,y) * (x[1] - y[1]) * (x[2] - y[2])) )
|
|
print ( "[1][26][223] : %25.15e"%(norm*df(a,x,y,1)* (x[1] - y[1]) * (x[2] - y[2])) )
|
|
|
|
print ( "[0][26][224] : %25.15e"%(f(a,x,y) * (x[2] - y[2]) * (x[2] - y[2])) )
|
|
print ( "[1][26][224] : %25.15e"%(df(a,x,y,1)* (x[2] - y[2]) * (x[2] - y[2])) )
|
|
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
#+begin_example
|
|
[0][26][219] : 1.020302912653649e-08
|
|
[1][26][219] : -4.153046808203204e-08
|
|
[0][26][220] : 1.516649653540510e-08
|
|
[1][26][220] : -7.725252615816528e-08
|
|
[0][26][221] : -4.686389780112468e-09
|
|
[1][26][221] : 2.387073693851122e-08
|
|
[0][26][222] : 7.514847283937212e-09
|
|
[1][26][222] : -4.025905373647693e-08
|
|
[0][26][223] : -4.021924592380977e-09
|
|
[1][26][223] : 2.154652944642284e-08
|
|
[0][26][224] : 7.175074806631758e-10
|
|
[1][26][224] : -3.843880138733679e-09
|
|
#+end_example
|
|
|
|
#+begin_src c :tangle (eval c_test) :exports none
|
|
{
|
|
#define walk_num 1 // 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, walk_num*elec_num*3);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
|
|
double ao_vgl[elec_num][5][ao_num];
|
|
|
|
rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]),
|
|
(int64_t) 5*elec_num*ao_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
printf("\n");
|
|
printf(" ao_vgl ao_vgl[0][26][219] %25.15e\n", ao_vgl[26][0][219]);
|
|
printf(" ao_vgl ao_vgl[1][26][219] %25.15e\n", ao_vgl[26][1][219]);
|
|
printf(" ao_vgl ao_vgl[0][26][220] %25.15e\n", ao_vgl[26][0][220]);
|
|
printf(" ao_vgl ao_vgl[1][26][220] %25.15e\n", ao_vgl[26][1][220]);
|
|
printf(" ao_vgl ao_vgl[0][26][221] %25.15e\n", ao_vgl[26][0][221]);
|
|
printf(" ao_vgl ao_vgl[1][26][221] %25.15e\n", ao_vgl[26][1][221]);
|
|
printf(" ao_vgl ao_vgl[0][26][222] %25.15e\n", ao_vgl[26][0][222]);
|
|
printf(" ao_vgl ao_vgl[1][26][222] %25.15e\n", ao_vgl[26][1][222]);
|
|
printf(" ao_vgl ao_vgl[0][26][223] %25.15e\n", ao_vgl[26][0][223]);
|
|
printf(" ao_vgl ao_vgl[1][26][223] %25.15e\n", ao_vgl[26][1][223]);
|
|
printf(" ao_vgl ao_vgl[0][26][224] %25.15e\n", ao_vgl[26][0][224]);
|
|
printf(" ao_vgl ao_vgl[1][26][224] %25.15e\n", ao_vgl[26][1][224]);
|
|
printf("\n");
|
|
|
|
assert( fabs(ao_vgl[26][0][219] - ( 1.020298798341620e-08)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[26][1][219] - (-4.928035238010602e-08)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[26][0][220] - ( 1.516643537739178e-08)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[26][1][220] - (-7.725221462603871e-08)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[26][0][221] - (-4.686370882518819e-09)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[26][1][221] - ( 2.387064067626827e-08)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[26][0][222] - ( 7.514816980753531e-09)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[26][1][222] - (-4.025889138635182e-08)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[26][0][223] - (-4.021908374204471e-09)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[26][1][223] - ( 2.154644255710413e-08)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[26][0][224] - ( 7.175045873560788e-10)) < 1.e-14 );
|
|
assert( fabs(ao_vgl[26][1][224] - (-3.843864637762753e-09)) < 1.e-14 );
|
|
}
|
|
|
|
#+end_src
|
|
|
|
* End of files :noexport:
|
|
|
|
#+begin_src c :tangle (eval h_private_type)
|
|
#endif
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval h_private_func)
|
|
#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
|
|
|
|
|
|
* TODO [0/1] Missing features :noexport:
|
|
- [ ] Error messages to tell what is missing when initializing
|