TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-12-22 20:36:01 +01:00
Code Issues Projects Releases Wiki Activity
26f913518a
qmckl/org/qmckl_mo.org
Anthony Scemama 4638b40724 Fix cusp
2024-07-11 13:53:10 +02:00

3555 lines
126 KiB
Org Mode
Raw Blame History

#+TITLE: Molecular Orbitals
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
The molecular orbitals (MOs) are defined in the basis of AOs along with a AO to MO
coefficient matrix $C$. Using these coefficients (e.g. from Hartree Fock method)
the MOs are defined as follows:
\[
\phi_i(\mathbf{r}) = C_i * \chi_i (\mathbf{r})
\]
By default, all the MOs are computed. A subset of MOs can be selected
for computation, for example to remove computation of the useless
virtual MOs present in a MO coefficient matrix.
For each nucleus, a radius $r_{\text{cusp}}$ can be given such that if
the electron-nucleus distance is smaller than $r_{\text{cusp}}$, all
the MOs are replaced by functions with the correct cusp condition and
such that the values and the gradients of the MOs are continuous at
$r_{\text{cusp}}$.
Molecular orbitals (MOs) are defined in the basis of atomic orbitals
(AOs) using a coefficient matrix $C$, which determines how the AOs are
combined to form the MOs.
The equation for the MOs $\phi_i(\mathbf{r})$ is given by:
$$\phi_i(\mathbf{r}) = \sum_{\mu} C_{\mu i} \chi_{\mu}(\mathbf{r})$$
where $i$ labels the MO, $\mu$ labels the AO, $C_{\mu i}$ is the
coefficient of AO $\mu$ in MO $i$, and $\chi_{\mu}(\mathbf{r})$ is the
AO itself.
In some cases, it may be desirable to only compute a subset of the
MOs, for example to exclude virtual MOs that do not contribute
to the wave function. This can be achieved by selecting a subset of
columns from the coefficient matrix $C$.
The exact wave function must have a cusp at the positions of the
nuclei to ensure that the kinetic energy diverges and cancels out the
divergence of the potential, resulting in a finite total energy.
To ensure that the cusp condition is satisfied, a modification can be
made to the molecular orbitals (MOs) when the electron-nucleus
distance is smaller than a certain radius $r_{\text{cusp}}$. At
distances smaller than $r_{\text{cusp}}$, the MOs are replaced by
functions that have the correct electron-nucleus cusp condition and
that ensure that the values and the gradients of the MOs are
continuous at $r_{\text{cusp}}$.
A radius $r_{\text{cusp}\, A}$ is given for each nucleus $A$, default is
zero. If an electron is closer to the nucleus $A$ than
$r_{\text{cusp}\, A}$, the MOs are locally replaced as
\[
\phi_{\text{cusp}\, i}(\mathbf{r}) = \phi_i(\mathbf{r}) - \phi_{s_A i}(\mathbf{r}) + \sum_{k=0}^{3} f_k\, |\mathbf{r}-\mathbf{R}_A|^k
\]
where $\phi_{s_A i}$ are the contributions of the $s$ AOs centered at
$A$ to MO $i$.
The coefficients $f_k$ are such that the value and gradient of the
MO are continuous at $r_{\text{cusp}}$, and the electron-nucleus
cusp is exact.
* 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_MO_HPF
#define QMCKL_MO_HPF
#include "qmckl_blas_private_type.h"
#+end_src
#+begin_src c :tangle (eval h_private_type)
#ifndef QMCKL_MO_HPT
#define QMCKL_MO_HPT
#include <stdbool.h>
#include "qmckl_blas_private_type.h"
#+end_src
#+begin_src c :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include "assert.h"
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <stdio.h>
#include <math.h>
#include "chbrclf.h"
#include "qmckl_electron_private_func.h"
#include "qmckl_ao_private_func.h"
#include "qmckl_mo_private_func.h"
int main() {
qmckl_context context;
context = qmckl_context_create();
qmckl_exit_code rc;
#+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 <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"
#include "qmckl_mo_private_type.h"
#include "qmckl_mo_private_func.h"
#+end_src
* Context
The following arrays are stored in the context:
|-----------------+--------------------+----------------------------------------------|
| ~mo_num~ | | Number of MOs |
| ~coefficient~ | ~[mo_num][ao_num]~ | MO coefficients |
| ~coefficient_t~ | ~[ao_num][mo_num]~ | Transposed of the Orbital coefficients |
| ~r_cusp~ | ~[nucl_num]~ | Radius of the functions for Cusp adjustments |
|-----------------+--------------------+----------------------------------------------|
Computed data:
|--------------+--------------------------+-----------------------------------------------------------|
| ~cusp_param~ | ~[nucl_num][4][mo_num]~ | Parameters of the functions for Cusp adjustments |
| ~mo_value~ | ~[point_num][mo_num]~ | Value of the MOs at point positions |
| ~mo_vgl~ | ~[point_num][5][mo_num]~ | Value, gradients, Laplacian of the MOs at point positions |
|--------------+--------------------------+-----------------------------------------------------------|
** Data structure
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_mo_basis_struct {
int64_t mo_num;
double * restrict coefficient;
double * restrict coefficient_t;
double * restrict r_cusp;
double * restrict mo_vgl;
double * restrict mo_value;
qmckl_tensor cusp_param;
uint64_t mo_vgl_date;
uint64_t mo_value_date;
int32_t uninitialized;
bool provided;
} qmckl_mo_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_mo_basis(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code qmckl_init_mo_basis(qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
ctx->mo_basis.r_cusp = NULL;
ctx->mo_basis.uninitialized = (1 << 2) - 1;
return QMCKL_SUCCESS;
}
#+end_src
** Initialization functions
To set the basis set, all the following functions need to be
called.
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_set_mo_basis_mo_num (qmckl_context context, const int64_t mo_num);
qmckl_exit_code qmckl_set_mo_basis_coefficient (qmckl_context context, const double * coefficient, const int64_t size_max);
qmckl_exit_code qmckl_set_mo_basis_r_cusp (qmckl_context context, const double * r_cusp, const int64_t size_max);
#+end_src
#+NAME:pre
#+begin_src c :exports none
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
if (mask != 0 && !(ctx->mo_basis.uninitialized & mask)) {
return qmckl_failwith( context,
QMCKL_ALREADY_SET,
"qmckl_set_mo_*",
NULL);
}
#+end_src
#+NAME:post
#+begin_src c :exports none
ctx->mo_basis.uninitialized &= ~mask;
ctx->mo_basis.provided = (ctx->mo_basis.uninitialized == 0);
if (ctx->mo_basis.provided) {
qmckl_exit_code rc_ = qmckl_finalize_mo_basis(context);
if (rc_ != QMCKL_SUCCESS) return rc_;
}
return QMCKL_SUCCESS;
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_set_mo_basis_mo_num(qmckl_context context, const int64_t mo_num) {
int32_t mask = 1 ;
<<pre>>
if (mo_num <= 0) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_set_mo_basis_mo_num",
"mo_num <= 0");
}
ctx->mo_basis.mo_num = mo_num;
<<post>>
}
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_set_mo_basis_coefficient(qmckl_context context,
const double* coefficient,
const int64_t size_max)
{
int32_t mask = 1 << 1;
<<pre>>
if (ctx->mo_basis.coefficient != NULL) {
qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_set_mo_basis_coefficient",
NULL);
}
}
if (size_max < ctx->ao_basis.ao_num * ctx->mo_basis.mo_num) {
return qmckl_failwith( context, QMCKL_INVALID_ARG_3,
"qmckl_set_mo_basis_coefficient",
"Array too small");
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.ao_num * ctx->mo_basis.mo_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_mo_basis_coefficient",
NULL);
}
memcpy(new_array, coefficient, mem_info.size);
ctx->mo_basis.coefficient = new_array;
<<post>>
}
#+end_src
When the basis set is completely entered, other data structures are
computed to accelerate the calculations.
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_finalize_mo_basis",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (ctx->mo_basis.coefficient_t != NULL) {
qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient_t);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_finalize_mo_basis",
NULL);
}
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(double);
double* new_array = (double*) qmckl_malloc(context, mem_info);
if (new_array == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_finalize_mo_basis",
NULL);
}
assert (ctx->mo_basis.coefficient != NULL);
for (int64_t i=0 ; i<ctx->ao_basis.ao_num ; ++i) {
for (int64_t j=0 ; j<ctx->mo_basis.mo_num ; ++j) {
new_array[i*ctx->mo_basis.mo_num + j] = ctx->mo_basis.coefficient[j*ctx->ao_basis.ao_num + i];
}
}
ctx->mo_basis.coefficient_t = new_array;
qmckl_exit_code rc;
if (ctx->mo_basis.mo_vgl != NULL) {
rc = qmckl_free(context, ctx->mo_basis.mo_vgl);
if (rc != QMCKL_SUCCESS) return rc;
ctx->mo_basis.mo_vgl = NULL;
ctx->mo_basis.mo_vgl_date = 0;
}
if (ctx->mo_basis.mo_value != NULL) {
rc = qmckl_free(context, ctx->mo_basis.mo_value);
if (rc != QMCKL_SUCCESS) return rc;
ctx->mo_basis.mo_value = NULL;
ctx->mo_basis.mo_value_date = 0;
}
return qmckl_context_touch(context);
}
#+end_src
** Cusp adjsutment functions
To activate the cusp adjustment, the user must enter the radius of
the fitting function for each atom.
This function requires the computation of the value and gradients
of the $s$ AOs at the distance equal to the radius, and the values
of the non-$s$ AOs at the center.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_mo_basis_r_cusp(qmckl_context context,
const double* r_cusp,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
if (r_cusp == NULL) {
return qmckl_failwith( context, QMCKL_INVALID_ARG_2,
"qmckl_set_mo_basis_r_cusp",
"r_cusp: Null pointer");
}
if (size_max < ctx->nucleus.num) {
return qmckl_failwith( context, QMCKL_INVALID_ARG_3,
"qmckl_set_mo_basis_r_cusp",
"Array too small");
}
// Nullify r_cusp
if (ctx->mo_basis.r_cusp != NULL) {
rc = qmckl_free(context, ctx->mo_basis.r_cusp);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->mo_basis.r_cusp = NULL;
}
// Save old points
int64_t old_point_num = ctx->point.num;
double* old_coord = NULL;
if (old_point_num > 0) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = old_point_num * 3 * sizeof(double);
old_coord = (double*) qmckl_malloc(context, mem_info);
rc = qmckl_get_point(context, 'T', old_coord, (old_point_num * 3));
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_FAILURE,
"qmckl_set_mo_basis_r_cusp",
"Unable to get old coordinates");
}
}
double* coord = (double*) malloc(ctx->nucleus.num * 3 * sizeof(double));
// Set r_cusp
{
assert (ctx->mo_basis.r_cusp == NULL);
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->nucleus.num * sizeof(double);
ctx->mo_basis.r_cusp = (double*) qmckl_malloc(context, mem_info);
if (ctx->mo_basis.r_cusp == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_mo_basis_r_cusp",
NULL);
}
memcpy(ctx->mo_basis.r_cusp, r_cusp, mem_info.size);
}
// Allocate cusp parameters and set them to zero
{
if (ctx->mo_basis.cusp_param.size[0] != 0) {
rc = qmckl_tensor_free(context, &(ctx->mo_basis.cusp_param));
if (rc != QMCKL_SUCCESS) return rc;
}
int64_t sze[3] = { ctx->mo_basis.mo_num, 4, ctx->nucleus.num };
ctx->mo_basis.cusp_param = qmckl_tensor_alloc(context, 3, &(sze[0]));
ctx->mo_basis.cusp_param = qmckl_tensor_set(ctx->mo_basis.cusp_param, 0.);
}
// Evaluate MO value at nucleus without s components
qmckl_matrix mo_value_at_nucl_no_s;
{
mo_value_at_nucl_no_s = qmckl_matrix_alloc(context, ctx->mo_basis.mo_num, ctx->nucleus.num);
rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, ctx->nucleus.num * 3);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_set_point(context, 'T', ctx->nucleus.num, coord, (ctx->nucleus.num * 3));
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_FAILURE,
"qmckl_set_mo_basis_r_cusp",
"Unable to set coordinates at the nuclei");
}
rc = qmckl_get_mo_basis_mo_value(context,
&(qmckl_mat(mo_value_at_nucl_no_s,0,0)),
ctx->mo_basis.mo_num * ctx->nucleus.num);
if (rc != QMCKL_SUCCESS) return rc;
}
// Evaluate MO vgl at r_cusp without s components
qmckl_tensor mo_vgl_at_r_cusp_s;
{
int64_t sze[3] = { ctx->mo_basis.mo_num, 3, ctx->nucleus.num };
mo_vgl_at_r_cusp_s = qmckl_tensor_alloc(context, 3, &(sze[0]));
}
{
qmckl_tensor ao_vgl_at_r_cusp_s;
int64_t sze[3] = { ctx->ao_basis.ao_num, 5, ctx->nucleus.num };
ao_vgl_at_r_cusp_s = qmckl_tensor_alloc(context, 3, &(sze[0]));
rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, ctx->nucleus.num * 3);
if (rc != QMCKL_SUCCESS) return rc;
for (int64_t i=0 ; i<ctx->nucleus.num ; ++i) {
coord[2*ctx->nucleus.num + i] += r_cusp[i];
}
rc = qmckl_set_point(context, 'T', ctx->nucleus.num, coord, (ctx->nucleus.num * 3));
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_FAILURE,
"qmckl_set_mo_basis_r_cusp",
"Unable to set coordinates at r_cusp");
}
rc = qmckl_get_ao_basis_ao_vgl(context,
&(qmckl_ten3(ao_vgl_at_r_cusp_s,0,0,0)),
ctx->ao_basis.ao_num * 5 * ctx->point.num);
for (int64_t inucl=0 ; inucl<ctx->nucleus.num ; ++inucl) {
for (int64_t i=0 ; i<ctx->mo_basis.mo_num ; ++i) {
qmckl_ten3(mo_vgl_at_r_cusp_s,i,0,inucl) = 0.;
qmckl_ten3(mo_vgl_at_r_cusp_s,i,1,inucl) = 0.;
qmckl_ten3(mo_vgl_at_r_cusp_s,i,2,inucl) = 0.;
for (int64_t k=0 ; k<ctx->ao_basis.ao_num ; ++k) {
if ( ctx->ao_basis.ao_nucl[k] == inucl && ctx->ao_basis.ao_ang_mom[k] == 0) {
const double ck = ctx->mo_basis.coefficient[k + i*ctx->ao_basis.ao_num];
qmckl_ten3(mo_vgl_at_r_cusp_s,i,0,inucl) = qmckl_ten3(mo_vgl_at_r_cusp_s,i,0,inucl) + ck * qmckl_ten3(ao_vgl_at_r_cusp_s,k,0,inucl);
qmckl_ten3(mo_vgl_at_r_cusp_s,i,1,inucl) = qmckl_ten3(mo_vgl_at_r_cusp_s,i,1,inucl) + ck * qmckl_ten3(ao_vgl_at_r_cusp_s,k,3,inucl);
qmckl_ten3(mo_vgl_at_r_cusp_s,i,2,inucl) = qmckl_ten3(mo_vgl_at_r_cusp_s,i,2,inucl) + ck * qmckl_ten3(ao_vgl_at_r_cusp_s,k,4,inucl);
}
}
}
}
rc = qmckl_tensor_free(context,&ao_vgl_at_r_cusp_s);
if (rc != QMCKL_SUCCESS) return rc;
}
// Compute parameters
{
for (int64_t inucl=0 ; inucl < ctx->nucleus.num ; ++inucl) {
const double Z = qmckl_vec(ctx->nucleus.charge,inucl);
if (Z < 0.1) continue; // Avoid dummy atoms
const double R = r_cusp[inucl];
assert (R != 0.);
assert (fabs(2.*R*Z-6.) > 1.e-6);
assert (fabs(R*R*Z-3.*R) > 1.e-6);
assert (fabs(R*R*R*Z-3.*R*R) > 1.e-6);
for (int64_t i=0 ; i<ctx->mo_basis.mo_num ; ++i) {
const double phi = qmckl_ten3(mo_vgl_at_r_cusp_s,i,0,inucl);
const double grad_phi = qmckl_ten3(mo_vgl_at_r_cusp_s,i,1,inucl);
const double lap_phi = qmckl_ten3(mo_vgl_at_r_cusp_s,i,2,inucl);
const double eta = qmckl_mat(mo_value_at_nucl_no_s,i,inucl);
qmckl_ten3(ctx->mo_basis.cusp_param,i,0,inucl) =
-(R*(2.*eta*Z-6.*grad_phi)+lap_phi*R*R+6.*phi)/(2.*R*Z-6.);
qmckl_ten3(ctx->mo_basis.cusp_param,i,1,inucl) =
(lap_phi*R*R*Z-6.*grad_phi*R*Z+6.*phi*Z+6.*eta*Z)/(2.*R*Z-6.);
qmckl_ten3(ctx->mo_basis.cusp_param,i,2,inucl) =
-(R*(-5.*grad_phi*Z-1.5*lap_phi)+lap_phi*R*R*Z+3.*phi*Z+3.*eta*Z+6.*grad_phi)/(R*R*Z-3.*R);
qmckl_ten3(ctx->mo_basis.cusp_param,i,3,inucl) =
(R*(-2.*grad_phi*Z-lap_phi)+0.5*lap_phi*R*R*Z+phi*Z+eta*Z+3.*grad_phi)/(R*R*R*Z-3.*R*R);
}
}
}
free(coord);
qmckl_matrix_free(context, &mo_value_at_nucl_no_s);
qmckl_tensor_free(context, &mo_vgl_at_r_cusp_s);
// Restore old points
if (old_point_num > 0) {
rc = qmckl_set_point(context, 'T', old_point_num, old_coord, (old_point_num * 3));
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_FAILURE,
"qmckl_set_mo_basis_r_cusp",
"Unable to set old coordinates");
}
rc = qmckl_free(context, old_coord);
if (rc != QMCKL_SUCCESS) return rc;
old_coord = NULL;
}
return QMCKL_SUCCESS;
}
#+end_src
** Access functions
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_mo_basis_mo_num (const qmckl_context context,
int64_t* mo_num);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_mo_basis_mo_num (const qmckl_context context,
int64_t* mo_num)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_mo_basis_mo_num",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1;
if ( (ctx->mo_basis.uninitialized & mask) != 0) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_get_mo_basis_mo_num",
NULL);
}
assert (ctx->mo_basis.mo_num > (int64_t) 0);
,*mo_num = ctx->mo_basis.mo_num;
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_mo_basis_coefficient (const qmckl_context context,
double* const coefficient,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code
qmckl_get_mo_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_mo_basis_coefficient",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 1;
if ( (ctx->mo_basis.uninitialized & mask) != 0) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_get_mo_basis_coefficient",
NULL);
}
if (coefficient == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_get_mo_basis_coefficient",
"NULL pointer");
}
if (size_max < ctx->ao_basis.ao_num * ctx->mo_basis.mo_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_mo_basis_coefficient",
"Array too small: expected mo_num * ao_num.");
}
assert (ctx->mo_basis.coefficient != NULL);
memcpy(coefficient, ctx->mo_basis.coefficient,
ctx->ao_basis.ao_num * ctx->mo_basis.mo_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_mo_basis_provided (const qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :exports none
bool qmckl_mo_basis_provided(const qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
return ctx->mo_basis.provided;
}
#+end_src
*** Fortran interfaces
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(qmckl_exit_code) function qmckl_get_mo_basis_mo_num (context, &
mo_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) :: mo_num
end function qmckl_get_mo_basis_mo_num
end interface
interface
integer(qmckl_exit_code) function qmckl_get_mo_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_mo_basis_coefficient
end interface
interface
integer(qmckl_exit_code) function qmckl_set_mo_basis_r_cusp(context, &
r_cusp, 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) :: r_cusp(*)
integer (c_int64_t) , intent(in), value :: size_max
end function qmckl_set_mo_basis_r_cusp
end interface
#+end_src
** Update
It may be desirable to remove certain molecular orbitals (MOs) that
do not significantly contribute to the wave function. In
particular, in a single determinant calculation, the virtual MOs
can be removed as they do not participate in the ground state
configuration.
To select a subset of MOs that will be kept, an array of integers of
size ~mo_num~ can be created. If the integer corresponding to an MO is
zero, that MO is dropped and will not be included in the
calculation. If the integer is non-zero, the MO will be kept.
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_mo_basis_select_mo (const qmckl_context context,
const int32_t* keep,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code
qmckl_mo_basis_select_mo (const qmckl_context context,
const int32_t* keep,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_mo_basis_select_mo",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if ( !(qmckl_mo_basis_provided(context)) ) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_get_mo_basis_select_mo",
NULL);
}
if (keep == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_get_mo_basis_select_mo",
"NULL pointer");
}
const int64_t mo_num = ctx->mo_basis.mo_num;
const int64_t ao_num = ctx->ao_basis.ao_num;
if (size_max < mo_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_mo_basis_select_mo",
"Array too small: expected mo_num.");
}
/*
if (ctx->mo_basis.r_cusp != NULL) {
return qmckl_failwith( context,
QMCKL_ALREADY_SET,
"qmckl_get_mo_basis_select_mo",
"r_cusp is already set. Please set it only after selecting MOs.");
}
*/
int64_t mo_num_new = 0;
for (int64_t i=0 ; i<mo_num ; ++i) {
if (keep[i] != 0) ++mo_num_new;
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ao_num * mo_num_new * sizeof(double);
double* restrict coefficient = (double*) qmckl_malloc(context, mem_info);
int64_t k = 0;
for (int64_t i=0 ; i<mo_num ; ++i) {
if (keep[i] != 0) {
memcpy(&(coefficient[k*ao_num]), &(ctx->mo_basis.coefficient[i*ao_num]), ao_num*sizeof(double));
++k;
}
}
qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient);
if (rc != QMCKL_SUCCESS) return rc;
ctx->mo_basis.coefficient = coefficient;
ctx->mo_basis.mo_num = mo_num_new;
rc = qmckl_finalize_mo_basis(context);
if (ctx->mo_basis.r_cusp != NULL) {
double * tmp = (double*) calloc (ctx->nucleus.num, sizeof(double));
assert(tmp != NULL);
memcpy(tmp, ctx->mo_basis.r_cusp, ctx->nucleus.num * sizeof(double));
rc = qmckl_set_mo_basis_r_cusp(context, tmp,
mo_num_new * 4 * ctx->nucleus.num);
free(tmp);
}
return rc;
}
#+end_src
*** Fortran interface
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(qmckl_exit_code) function qmckl_mo_basis_select_mo (context, &
keep, 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) :: keep(*)
integer (c_int64_t) , intent(in), value :: size_max
end function qmckl_mo_basis_select_mo
end interface
#+end_src
* Computation
** Parameters of the cusp-correction functions
** Computation of MOs: values only
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_mo_basis_mo_value(qmckl_context context,
double* const mo_value,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_mo_basis_mo_value(qmckl_context context,
double* const mo_value,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_mo_basis_mo_value(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->point.num * ctx->mo_basis.mo_num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_mo_basis_mo_value",
"input array too small");
}
memcpy(mo_value, ctx->mo_basis.mo_value, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(qmckl_exit_code) function qmckl_get_mo_basis_mo_value (context, &
mo_value, 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) :: mo_value(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_mo_basis_mo_value
end interface
#+end_src
Uses the given array to compute the values.
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
double* const mo_value,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
double* const mo_value,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_mo_basis_mo_value",
NULL);
}
qmckl_exit_code rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->mo_basis.mo_num * ctx->point.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_mo_basis_mo_value",
"input array too small");
}
rc = qmckl_context_touch(context);
if (rc != QMCKL_SUCCESS) return rc;
double* old_array = ctx->mo_basis.mo_value;
ctx->mo_basis.mo_value = mo_value;
rc = qmckl_provide_mo_basis_mo_value(context);
if (rc != QMCKL_SUCCESS) return rc;
ctx->mo_basis.mo_value = old_array;
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(qmckl_exit_code) function qmckl_get_mo_basis_mo_value_inplace (context, &
mo_value, 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) :: mo_value(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_mo_basis_mo_value_inplace
end interface
#+end_src
*** Provide
#+CALL: write_provider_header( group="mo_basis", data="mo_value" )
#+RESULTS:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :export none
qmckl_exit_code qmckl_provide_mo_basis_mo_value(qmckl_context context);
#+end_src
#+CALL: write_provider_pre( group="mo_basis", data="mo_value", dimension="ctx->mo_basis.mo_num * ctx->point.num")
#+RESULTS:
#+begin_src c :comments org :tangle (eval c) :noweb yes :export none
qmckl_exit_code qmckl_provide_mo_basis_mo_value(qmckl_context context)
{
qmckl_exit_code rc = QMCKL_SUCCESS;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_provide_mo_basis_mo_value",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->mo_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_provide_mo_basis_mo_value",
NULL);
}
/* Compute if necessary */
if (ctx->point.date > ctx->mo_basis.mo_value_date) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->mo_basis.mo_num * ctx->point.num * sizeof(double);
if (ctx->mo_basis.mo_value != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->mo_basis.mo_value, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->mo_basis.mo_value);
assert (rc == QMCKL_SUCCESS);
ctx->mo_basis.mo_value = NULL;
}
}
/* Allocate array */
if (ctx->mo_basis.mo_value == NULL) {
double* mo_value = (double*) qmckl_malloc(context, mem_info);
if (mo_value == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_mo_basis_mo_value",
NULL);
}
ctx->mo_basis.mo_value = mo_value;
}
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
if (ctx->mo_basis.mo_vgl_date == ctx->point.date) {
// mo_vgl has been computed at this step: Just copy the data.
double * v = &(ctx->mo_basis.mo_value[0]);
double * vgl = &(ctx->mo_basis.mo_vgl[0]);
for (int i=0 ; i<ctx->point.num ; ++i) {
for (int k=0 ; k<ctx->mo_basis.mo_num ; ++k) {
v[k] = vgl[k];
}
v += ctx->mo_basis.mo_num;
vgl += ctx->mo_basis.mo_num * 5;
}
} else {
rc = qmckl_provide_ao_basis_ao_value(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_value",
NULL);
}
if (ctx->mo_basis.r_cusp == NULL) {
/* No cusp correction */
rc = qmckl_compute_mo_basis_mo_value(context,
ctx->ao_basis.ao_num,
ctx->mo_basis.mo_num,
ctx->point.num,
ctx->mo_basis.coefficient_t,
ctx->ao_basis.ao_value,
ctx->mo_basis.mo_value);
} else {
rc = qmckl_provide_en_distance(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_provide_mo_basis_mo_value",
"en_distance");
}
rc = qmckl_compute_mo_basis_mo_value_cusp(context,
ctx->nucleus.num,
ctx->ao_basis.ao_num,
ctx->mo_basis.mo_num,
ctx->point.num,
ctx->ao_basis.ao_nucl,
ctx->ao_basis.ao_ang_mom,
ctx->electron.en_distance,
ctx->mo_basis.r_cusp,
ctx->mo_basis.cusp_param,
ctx->mo_basis.coefficient_t,
ctx->ao_basis.ao_value,
ctx->mo_basis.mo_value);
}
}
#+end_src
#+CALL: write_provider_post( group="mo_basis", data="mo_value" )
#+RESULTS:
#+begin_src c :comments org :tangle (eval c) :noweb yes :export none
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->mo_basis.mo_value_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
*** Compute
:PROPERTIES:
:Name: qmckl_compute_mo_basis_mo_value
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_mo_basis_mo_value_args
| Variable | Type | In/Out | Description |
|-----------------+-----------------------------+--------+-------------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~ao_num~ | ~int64_t~ | in | Number of AOs |
| ~mo_num~ | ~int64_t~ | in | Number of MOs |
| ~point_num~ | ~int64_t~ | in | Number of points |
| ~coefficient_t~ | ~double[mo_num][ao_num]~ | in | Transpose of the AO to MO transformation matrix |
| ~ao_value~ | ~double[point_num][ao_num]~ | in | Value of the AOs |
| ~mo_value~ | ~double[point_num][mo_num]~ | out | Value of the MOs |
The matrix of AO values is very sparse, so we use a sparse-dense
matrix multiplication instead of a dgemm, as exposed in
https://dx.doi.org/10.1007/978-3-642-38718-0_14.
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_mo_basis_mo_value_doc_f(context, &
ao_num, mo_num, point_num, &
coefficient_t, ao_value, mo_value) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: ao_num, mo_num
integer*8 , intent(in) :: point_num
double precision , intent(in) :: ao_value(ao_num,point_num)
double precision , intent(in) :: coefficient_t(mo_num,ao_num)
double precision , intent(out) :: mo_value(mo_num,point_num)
integer*8 :: j,k
info = QMCKL_SUCCESS
do j=1,point_num
mo_value(:,j) = 0.d0
do k=1,ao_num
if (ao_value(k,j) == 0.d0) cycle
mo_value(:,j) = mo_value(:,j) + coefficient_t(:,k) * ao_value(k,j)
end do
end do
end function qmckl_compute_mo_basis_mo_value_doc_f
#+end_src
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value"))
#+RESULTS:
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_value (
const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* coefficient_t,
const double* ao_value,
double* const mo_value );
#+end_src
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc"))
#+RESULTS:
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_value_doc (
const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* coefficient_t,
const double* ao_value,
double* const mo_value );
#+end_src
#+CALL: generate_c_interface(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_mo_basis_mo_value_doc &
(context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value) &
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 :: mo_num
integer (c_int64_t) , intent(in) , value :: point_num
real (c_double ) , intent(in) :: coefficient_t(ao_num,mo_num)
real (c_double ) , intent(in) :: ao_value(ao_num,point_num)
real (c_double ) , intent(out) :: mo_value(mo_num,point_num)
integer(c_int32_t), external :: qmckl_compute_mo_basis_mo_value_doc_f
info = qmckl_compute_mo_basis_mo_value_doc_f &
(context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value)
end function qmckl_compute_mo_basis_mo_value_doc
#+end_src
#+begin_src c :tangle (eval c) :comments org
qmckl_exit_code
qmckl_compute_mo_basis_mo_value (const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* coefficient_t,
const double* ao_value,
double* const mo_value )
{
#ifdef HAVE_HPC
return qmckl_compute_mo_basis_mo_value_hpc (context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value);
#else
return qmckl_compute_mo_basis_mo_value_doc (context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value);
#endif
}
#+end_src
*** HPC version
#+begin_src c :tangle (eval h_private_func) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* coefficient_t,
const double* ao_value,
double* const mo_value );
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_hpc_sp (const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* coefficient_t,
const double* ao_value,
double* const mo_value );
#endif
#+end_src
**** Single-precision
#+begin_src c :tangle (eval c) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_hpc_sp (const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* restrict coefficient_t,
const double* restrict ao_value,
double* restrict const mo_value )
{
assert (context != QMCKL_NULL_CONTEXT);
float* __attribute__((aligned(64))) coefficient_t_sp = calloc(ao_num*mo_num, sizeof(float));
if (coefficient_t_sp == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_compute_mo_basis_mo_value_hpc_sp",
"coefficient_t_sp");
};
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num*ao_num ; ++i) {
coefficient_t_sp[i] = (float) coefficient_t[i];
}
#ifdef HAVE_OPENMP
#pragma omp parallel
#endif
{
int64_t* __attribute__((aligned(64))) idx = calloc((size_t) ao_num, sizeof(int64_t));
float* __attribute__((aligned(64))) av1 = calloc((size_t) ao_num, sizeof(float));
assert (idx != NULL);
assert (av1 != NULL);
#ifdef HAVE_OPENMP
#pragma omp for
#endif
for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
double* restrict const vgl1 = &(mo_value[ipoint*mo_num]);
const double* restrict avgl1 = &(ao_value[ipoint*ao_num]);
float* __attribute__((aligned(64))) vgl_sp = calloc((size_t) mo_num, sizeof(float));
assert (vgl_sp != NULL);
int64_t nidx=0;
for (int64_t k=0 ; k<ao_num ; ++k) {
if (avgl1[k] != 0.) {
idx[nidx] = k;
av1[nidx] = (float) avgl1[k];
++nidx;
}
}
int64_t n=0;
for (n=0 ; n < nidx-8 ; n+=8) {
const float* restrict ck1 = coefficient_t_sp + idx[n ]*mo_num;
const float* restrict ck2 = coefficient_t_sp + idx[n+1]*mo_num;
const float* restrict ck3 = coefficient_t_sp + idx[n+2]*mo_num;
const float* restrict ck4 = coefficient_t_sp + idx[n+3]*mo_num;
const float* restrict ck5 = coefficient_t_sp + idx[n+4]*mo_num;
const float* restrict ck6 = coefficient_t_sp + idx[n+5]*mo_num;
const float* restrict ck7 = coefficient_t_sp + idx[n+6]*mo_num;
const float* restrict ck8 = coefficient_t_sp + idx[n+7]*mo_num;
const float a11 = av1[n ];
const float a21 = av1[n+1];
const float a31 = av1[n+2];
const float a41 = av1[n+3];
const float a51 = av1[n+4];
const float a61 = av1[n+5];
const float a71 = av1[n+6];
const float a81 = av1[n+7];
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl_sp[i] = vgl_sp[i] +
ck1[i] * a11 + ck2[i] * a21 + ck3[i] * a31 + ck4[i] * a41 +
ck5[i] * a51 + ck6[i] * a61 + ck7[i] * a71 + ck8[i] * a81;
}
}
for (int64_t m=n ; m < nidx ; m+=1) {
const float* restrict ck = coefficient_t_sp + idx[m]*mo_num;
const float a1 = av1[m];
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl_sp[i] = vgl_sp[i] + ck[i] * a1;
}
}
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] = (double) vgl_sp[i];
}
free(vgl_sp);
}
free(av1);
free(idx);
}
free(coefficient_t_sp);
return QMCKL_SUCCESS;
}
#endif
#+end_src
**** Double-precision
#+begin_src c :tangle (eval c) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* restrict coefficient_t,
const double* restrict ao_value,
double* restrict const mo_value )
{
assert (context != QMCKL_NULL_CONTEXT);
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
/* Don't compute polynomials when the radial part is zero. */
const int precision = ctx->numprec.precision;
const bool single_precision = precision <= 24;
if (single_precision) {
return qmckl_compute_mo_basis_mo_value_hpc_sp (context,
ao_num,
mo_num,
point_num,
coefficient_t,
ao_value,
mo_value );
}
#ifdef HAVE_OPENMP
#pragma omp parallel
#endif
{
int64_t* __attribute__((aligned(64))) idx = calloc(ao_num, sizeof(int64_t));
double* __attribute__((aligned(64))) av1 = calloc(ao_num, sizeof(double));
assert (idx != NULL);
assert (av1 != NULL);
#ifdef HAVE_OPENMP
#pragma omp for
#endif
for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
double* restrict const vgl1 = &(mo_value[ipoint*mo_num]);
const double* restrict avgl1 = &(ao_value[ipoint*ao_num]);
memset(vgl1, 0, mo_num*sizeof(double));
int64_t nidx=0;
int64_t idx[ao_num];
double av1[ao_num];
for (int64_t k=0 ; k<ao_num ; ++k) {
if (avgl1[k] != 0.) {
idx[nidx] = k;
av1[nidx] = avgl1[k];
++nidx;
}
}
int64_t n=0;
for (n=0 ; n < nidx-8 ; n+=8) {
const double* restrict ck1 = coefficient_t + idx[n ]*mo_num;
const double* restrict ck2 = coefficient_t + idx[n+1]*mo_num;
const double* restrict ck3 = coefficient_t + idx[n+2]*mo_num;
const double* restrict ck4 = coefficient_t + idx[n+3]*mo_num;
const double* restrict ck5 = coefficient_t + idx[n+4]*mo_num;
const double* restrict ck6 = coefficient_t + idx[n+5]*mo_num;
const double* restrict ck7 = coefficient_t + idx[n+6]*mo_num;
const double* restrict ck8 = coefficient_t + idx[n+7]*mo_num;
const double a11 = av1[n ];
const double a21 = av1[n+1];
const double a31 = av1[n+2];
const double a41 = av1[n+3];
const double a51 = av1[n+4];
const double a61 = av1[n+5];
const double a71 = av1[n+6];
const double a81 = av1[n+7];
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] = vgl1[i] +
ck1[i] * a11 + ck2[i] * a21 + ck3[i] * a31 + ck4[i] * a41 +
ck5[i] * a51 + ck6[i] * a61 + ck7[i] * a71 + ck8[i] * a81;
}
}
for (int64_t m=n ; m < nidx ; m+=1) {
const double* restrict ck = coefficient_t + idx[m]*mo_num;
const double a1 = av1[m];
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] = vgl1[i] + ck[i] * a1;
}
}
}
free(av1);
free(idx);
}
return QMCKL_SUCCESS;
}
#endif
#+end_src
** Computation of MOs: values, gradient, Laplacian
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_mo_basis_mo_vgl(qmckl_context context,
double* const mo_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_mo_basis_mo_vgl(qmckl_context context,
double* const mo_vgl,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->point.num * 5 * ctx->mo_basis.mo_num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_mo_basis_mo_vgl",
"input array too small");
}
memcpy(mo_vgl, ctx->mo_basis.mo_vgl, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(qmckl_exit_code) function qmckl_get_mo_basis_mo_vgl (context, &
mo_vgl, 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) :: mo_vgl(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_mo_basis_mo_vgl
end interface
#+end_src
Uses the given array to compute the VGL.
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_mo_basis_mo_vgl_inplace (qmckl_context context,
double* const mo_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_mo_basis_mo_vgl_inplace (qmckl_context context,
double* const mo_vgl,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_mo_basis_mo_vgl",
NULL);
}
qmckl_exit_code rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->mo_basis.mo_num * 5 * ctx->point.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_mo_basis_mo_vgl",
"input array too small");
}
rc = qmckl_context_touch(context);
if (rc != QMCKL_SUCCESS) return rc;
double* old_array = ctx->mo_basis.mo_vgl;
ctx->mo_basis.mo_vgl = mo_vgl;
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
ctx->mo_basis.mo_vgl = old_array;
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(qmckl_exit_code) function qmckl_get_mo_basis_mo_vgl_inplace (context, &
mo_vgl, 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) :: mo_vgl(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_mo_basis_mo_vgl_inplace
end interface
#+end_src
*** Provide
#+CALL: write_provider_header( group="mo_basis", data="mo_vgl" )
#+RESULTS:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :export none
qmckl_exit_code qmckl_provide_mo_basis_mo_vgl(qmckl_context context);
#+end_src
#+CALL: write_provider_pre( group="mo_basis", data="mo_vgl", dimension="5 * ctx->mo_basis.mo_num * ctx->point.num")
#+RESULTS:
#+begin_src c :comments org :tangle (eval c) :noweb yes :export none
qmckl_exit_code qmckl_provide_mo_basis_mo_vgl(qmckl_context context)
{
qmckl_exit_code rc = QMCKL_SUCCESS;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_provide_mo_basis_mo_vgl",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->mo_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_provide_mo_basis_mo_vgl",
NULL);
}
/* Compute if necessary */
if (ctx->point.date > ctx->mo_basis.mo_vgl_date) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 5 * ctx->mo_basis.mo_num * ctx->point.num * sizeof(double);
if (ctx->mo_basis.mo_vgl != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->mo_basis.mo_vgl, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->mo_basis.mo_vgl);
assert (rc == QMCKL_SUCCESS);
ctx->mo_basis.mo_vgl = NULL;
}
}
/* Allocate array */
if (ctx->mo_basis.mo_vgl == NULL) {
double* mo_vgl = (double*) qmckl_malloc(context, mem_info);
if (mo_vgl == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_mo_basis_mo_vgl",
NULL);
}
ctx->mo_basis.mo_vgl = mo_vgl;
}
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_basis",
NULL);
}
if (ctx->mo_basis.r_cusp == NULL) {
/* No cusp correction */
rc = qmckl_compute_mo_basis_mo_vgl(context,
ctx->ao_basis.ao_num,
ctx->mo_basis.mo_num,
ctx->point.num,
ctx->mo_basis.coefficient_t,
ctx->ao_basis.ao_vgl,
ctx->mo_basis.mo_vgl);
} else {
rc = qmckl_provide_en_distance(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron_en_distance",
NULL);
}
rc = qmckl_compute_mo_basis_mo_vgl_cusp(context,
ctx->nucleus.num,
ctx->ao_basis.ao_num,
ctx->mo_basis.mo_num,
ctx->point.num,
ctx->ao_basis.ao_nucl,
ctx->ao_basis.ao_ang_mom,
ctx->electron.en_distance,
ctx->nucleus.coord,
ctx->point.coord,
ctx->mo_basis.r_cusp,
ctx->mo_basis.cusp_param,
ctx->mo_basis.coefficient_t,
ctx->ao_basis.ao_vgl,
ctx->mo_basis.mo_vgl);
}
#+end_src
#+CALL: write_provider_post( group="mo_basis", data="mo_vgl" )
#+RESULTS:
#+begin_src c :comments org :tangle (eval c) :noweb yes :export none
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->mo_basis.mo_vgl_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
*** Compute
:PROPERTIES:
:Name: qmckl_compute_mo_basis_mo_vgl
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_mo_basis_mo_vgl_args
| Variable | Type | In/Out | Description |
|---------------------+--------------------------------+--------+-------------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~ao_num~ | ~int64_t~ | in | Number of AOs |
| ~mo_num~ | ~int64_t~ | in | Number of MOs |
| ~point_num~ | ~int64_t~ | in | Number of points |
| ~coefficient_t~ | ~double[mo_num][ao_num]~ | in | Transpose of the AO to MO transformation matrix |
| ~ao_vgl~ | ~double[point_num][5][ao_num]~ | in | Value, gradients and Laplacian of the AOs |
| ~mo_vgl~ | ~double[point_num][5][mo_num]~ | out | Value, gradients and Laplacian of the MOs |
The matrix of AO values is very sparse, so we use a sparse-dense
matrix multiplication instead of a dgemm, as exposed in
https://dx.doi.org/10.1007/978-3-642-38718-0_14.
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_mo_basis_mo_vgl_doc_f(context, &
ao_num, mo_num, point_num, &
coefficient_t, ao_vgl, mo_vgl) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: ao_num, mo_num
integer*8 , intent(in) :: point_num
double precision , intent(in) :: ao_vgl(ao_num,5,point_num)
double precision , intent(in) :: coefficient_t(mo_num,ao_num)
double precision , intent(out) :: mo_vgl(mo_num,5,point_num)
integer*8 :: i,j,k
double precision :: c1, c2, c3, c4, c5
info = QMCKL_SUCCESS
do j=1,point_num
mo_vgl(:,:,j) = 0.d0
do k=1,ao_num
if (ao_vgl(k,1,j) /= 0.d0) then
c1 = ao_vgl(k,1,j)
c2 = ao_vgl(k,2,j)
c3 = ao_vgl(k,3,j)
c4 = ao_vgl(k,4,j)
c5 = ao_vgl(k,5,j)
do i=1,mo_num
mo_vgl(i,1,j) = mo_vgl(i,1,j) + coefficient_t(i,k) * c1
mo_vgl(i,2,j) = mo_vgl(i,2,j) + coefficient_t(i,k) * c2
mo_vgl(i,3,j) = mo_vgl(i,3,j) + coefficient_t(i,k) * c3
mo_vgl(i,4,j) = mo_vgl(i,4,j) + coefficient_t(i,k) * c4
mo_vgl(i,5,j) = mo_vgl(i,5,j) + coefficient_t(i,k) * c5
end do
end if
end do
end do
end function qmckl_compute_mo_basis_mo_vgl_doc_f
#+end_src
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl"))
#+RESULTS:
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_vgl (
const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* coefficient_t,
const double* ao_vgl,
double* const mo_vgl );
#+end_src
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_doc"))
#+RESULTS:
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_doc (
const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* coefficient_t,
const double* ao_vgl,
double* const mo_vgl );
#+end_src
#+CALL: generate_c_interface(table=qmckl_mo_basis_mo_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_doc"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_mo_basis_mo_vgl_doc &
(context, ao_num, mo_num, point_num, coefficient_t, ao_vgl, mo_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 :: mo_num
integer (c_int64_t) , intent(in) , value :: point_num
real (c_double ) , intent(in) :: coefficient_t(ao_num,mo_num)
real (c_double ) , intent(in) :: ao_vgl(ao_num,5,point_num)
real (c_double ) , intent(out) :: mo_vgl(mo_num,5,point_num)
integer(c_int32_t), external :: qmckl_compute_mo_basis_mo_vgl_doc_f
info = qmckl_compute_mo_basis_mo_vgl_doc_f &
(context, ao_num, mo_num, point_num, coefficient_t, ao_vgl, mo_vgl)
end function qmckl_compute_mo_basis_mo_vgl_doc
#+end_src
#+begin_src c :tangle (eval c) :comments org
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl (const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* coefficient_t,
const double* ao_vgl,
double* const mo_vgl )
{
#ifdef HAVE_HPC
return qmckl_compute_mo_basis_mo_vgl_hpc (context, ao_num, mo_num, point_num, coefficient_t, ao_vgl, mo_vgl);
#else
return qmckl_compute_mo_basis_mo_vgl_doc (context, ao_num, mo_num, point_num, coefficient_t, ao_vgl, mo_vgl);
#endif
}
#+end_src
*** HPC version :noexport:
**** Double precision
#+begin_src c :tangle (eval h_private_func) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* coefficient_t,
const double* ao_vgl,
double* const mo_vgl );
#endif
#+end_src
#+begin_src c :tangle (eval c) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* restrict coefficient_t,
const double* restrict ao_vgl,
double* restrict const mo_vgl )
{
assert (context != QMCKL_NULL_CONTEXT);
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
/* Don't compute polynomials when the radial part is zero. */
const int precision = ctx->numprec.precision;
const bool single_precision = precision <= 26;
if (single_precision) {
return qmckl_compute_mo_basis_mo_vgl_hpc_sp (context,
ao_num,
mo_num,
point_num,
coefficient_t,
ao_vgl,
mo_vgl );
}
#ifdef HAVE_OPENMP
#pragma omp parallel
#endif
{
int64_t* __attribute__((aligned(64))) idx = calloc(ao_num, sizeof(int64_t));
double* __attribute__((aligned(64))) av1 = calloc(ao_num, sizeof(double));
double* __attribute__((aligned(64))) av2 = calloc(ao_num, sizeof(double));
double* __attribute__((aligned(64))) av3 = calloc(ao_num, sizeof(double));
double* __attribute__((aligned(64))) av4 = calloc(ao_num, sizeof(double));
double* __attribute__((aligned(64))) av5 = calloc(ao_num, sizeof(double));
assert (idx != NULL);
assert (av1 != NULL);
assert (av2 != NULL);
assert (av3 != NULL);
assert (av4 != NULL);
assert (av5 != NULL);
#ifdef HAVE_OPENMP
#pragma omp for
#endif
for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
double* restrict const vgl1 = &(mo_vgl[ipoint*5*mo_num]);
double* restrict const vgl2 = vgl1 + mo_num;
double* restrict const vgl3 = vgl1 + (mo_num << 1);
double* restrict const vgl4 = vgl1 + (mo_num << 1) + mo_num;
double* restrict const vgl5 = vgl1 + (mo_num << 2);
const double* restrict avgl1 = &(ao_vgl[ipoint*5*ao_num]);
const double* restrict avgl2 = avgl1 + ao_num;
const double* restrict avgl3 = avgl1 + (ao_num << 1);
const double* restrict avgl4 = avgl1 + (ao_num << 1) + ao_num;
const double* restrict avgl5 = avgl1 + (ao_num << 2);
memset(vgl1,0,5*mo_num*sizeof(double));
int64_t nidx=0;
for (int64_t k=0 ; k<ao_num ; ++k) {
if (avgl1[k] != 0.) {
idx[nidx] = k;
av1[nidx] = avgl1[k];
av2[nidx] = avgl2[k];
av3[nidx] = avgl3[k];
av4[nidx] = avgl4[k];
av5[nidx] = avgl5[k];
++nidx;
}
}
int64_t n=0;
for (n=0 ; n < nidx-4 ; n+=4) {
const double* restrict ck1 = coefficient_t + idx[n ]*mo_num;
const double* restrict ck2 = coefficient_t + idx[n+1]*mo_num;
const double* restrict ck3 = coefficient_t + idx[n+2]*mo_num;
const double* restrict ck4 = coefficient_t + idx[n+3]*mo_num;
const double a11 = av1[n ];
const double a21 = av1[n+1];
const double a31 = av1[n+2];
const double a41 = av1[n+3];
const double a12 = av2[n ];
const double a22 = av2[n+1];
const double a32 = av2[n+2];
const double a42 = av2[n+3];
const double a13 = av3[n ];
const double a23 = av3[n+1];
const double a33 = av3[n+2];
const double a43 = av3[n+3];
const double a14 = av4[n ];
const double a24 = av4[n+1];
const double a34 = av4[n+2];
const double a44 = av4[n+3];
const double a15 = av5[n ];
const double a25 = av5[n+1];
const double a35 = av5[n+2];
const double a45 = av5[n+3];
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] = vgl1[i] + ck1[i] * a11 + ck2[i] * a21 + ck3[i] * a31 + ck4[i] * a41;
vgl2[i] = vgl2[i] + ck1[i] * a12 + ck2[i] * a22 + ck3[i] * a32 + ck4[i] * a42;
vgl3[i] = vgl3[i] + ck1[i] * a13 + ck2[i] * a23 + ck3[i] * a33 + ck4[i] * a43;
vgl4[i] = vgl4[i] + ck1[i] * a14 + ck2[i] * a24 + ck3[i] * a34 + ck4[i] * a44;
vgl5[i] = vgl5[i] + ck1[i] * a15 + ck2[i] * a25 + ck3[i] * a35 + ck4[i] * a45;
}
}
for (int64_t m=n ; m < nidx ; m+=1) {
const double* restrict ck = coefficient_t + idx[m]*mo_num;
const double a1 = av1[m];
const double a2 = av2[m];
const double a3 = av3[m];
const double a4 = av4[m];
const double a5 = av5[m];
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] = vgl1[i] + ck[i] * a1;
vgl2[i] = vgl2[i] + ck[i] * a2;
vgl3[i] = vgl3[i] + ck[i] * a3;
vgl4[i] = vgl4[i] + ck[i] * a4;
vgl5[i] = vgl5[i] + ck[i] * a5;
}
}
}
free(idx);
free(av1);
free(av2);
free(av3);
free(av4);
free(av5);
}
return QMCKL_SUCCESS;
}
#endif
#+end_src
**** Single precision
#+begin_src c :tangle (eval h_private_func) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl_hpc_sp (const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* coefficient_t,
const double* ao_vgl,
double* const mo_vgl );
#endif
#+end_src
#+begin_src c :tangle (eval c) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl_hpc_sp (const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* restrict coefficient_t,
const double* restrict ao_vgl,
double* restrict const mo_vgl )
{
assert (context != QMCKL_NULL_CONTEXT);
float* __attribute__((aligned(64))) coefficient_t_sp = calloc(ao_num*mo_num, sizeof(float));
if (coefficient_t_sp == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_compute_mo_basis_mo_vgl_hpc_sp",
"coefficient_t_sp");
};
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num*ao_num ; ++i) {
coefficient_t_sp[i] = (float) coefficient_t[i];
}
#ifdef HAVE_OPENMP
#pragma omp parallel
#endif
{
int64_t* __attribute__((aligned(64))) idx = calloc(ao_num, sizeof(int64_t));
float* __attribute__((aligned(64))) av1 = calloc(ao_num, sizeof(double));
float* __attribute__((aligned(64))) av2 = calloc(ao_num, sizeof(double));
float* __attribute__((aligned(64))) av3 = calloc(ao_num, sizeof(double));
float* __attribute__((aligned(64))) av4 = calloc(ao_num, sizeof(double));
float* __attribute__((aligned(64))) av5 = calloc(ao_num, sizeof(double));
assert (idx != NULL);
assert (av1 != NULL);
assert (av2 != NULL);
assert (av3 != NULL);
assert (av4 != NULL);
assert (av5 != NULL);
#ifdef HAVE_OPENMP
#pragma omp for
#endif
for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
float* __attribute__((aligned(64))) vgl_sp1 = calloc(mo_num, sizeof(float));
float* __attribute__((aligned(64))) vgl_sp2 = calloc(mo_num, sizeof(float));
float* __attribute__((aligned(64))) vgl_sp3 = calloc(mo_num, sizeof(float));
float* __attribute__((aligned(64))) vgl_sp4 = calloc(mo_num, sizeof(float));
float* __attribute__((aligned(64))) vgl_sp5 = calloc(mo_num, sizeof(float));
assert (vgl_sp1 != NULL);
assert (vgl_sp2 != NULL);
assert (vgl_sp3 != NULL);
assert (vgl_sp4 != NULL);
assert (vgl_sp5 != NULL);
double* restrict const vgl1 = &(mo_vgl[ipoint*5*mo_num]);
double* restrict const vgl2 = vgl1 + mo_num;
double* restrict const vgl3 = vgl1 + (mo_num << 1);
double* restrict const vgl4 = vgl1 + (mo_num << 1) + mo_num;
double* restrict const vgl5 = vgl1 + (mo_num << 2);
const double* restrict avgl1 = &(ao_vgl[ipoint*5*ao_num]);
const double* restrict avgl2 = avgl1 + ao_num;
const double* restrict avgl3 = avgl1 + (ao_num << 1);
const double* restrict avgl4 = avgl1 + (ao_num << 1) + ao_num;
const double* restrict avgl5 = avgl1 + (ao_num << 2);
int64_t nidx=0;
for (int64_t k=0 ; k<ao_num ; ++k) {
if (avgl1[k] != 0.) {
idx[nidx] = k;
av1[nidx] = (float) avgl1[k];
av2[nidx] = (float) avgl2[k];
av3[nidx] = (float) avgl3[k];
av4[nidx] = (float) avgl4[k];
av5[nidx] = (float) avgl5[k];
++nidx;
}
}
int64_t n=0;
for (n=0 ; n < nidx-4 ; n+=4) {
const float* restrict ck1 = coefficient_t_sp + idx[n ]*mo_num;
const float* restrict ck2 = coefficient_t_sp + idx[n+1]*mo_num;
const float* restrict ck3 = coefficient_t_sp + idx[n+2]*mo_num;
const float* restrict ck4 = coefficient_t_sp + idx[n+3]*mo_num;
const float a11 = av1[n ];
const float a21 = av1[n+1];
const float a31 = av1[n+2];
const float a41 = av1[n+3];
const float a12 = av2[n ];
const float a22 = av2[n+1];
const float a32 = av2[n+2];
const float a42 = av2[n+3];
const float a13 = av3[n ];
const float a23 = av3[n+1];
const float a33 = av3[n+2];
const float a43 = av3[n+3];
const float a14 = av4[n ];
const float a24 = av4[n+1];
const float a34 = av4[n+2];
const float a44 = av4[n+3];
const float a15 = av5[n ];
const float a25 = av5[n+1];
const float a35 = av5[n+2];
const float a45 = av5[n+3];
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl_sp1[i] = vgl_sp1[i] + ck1[i] * a11 + ck2[i] * a21 + ck3[i] * a31 + ck4[i] * a41;
vgl_sp2[i] = vgl_sp2[i] + ck1[i] * a12 + ck2[i] * a22 + ck3[i] * a32 + ck4[i] * a42;
vgl_sp3[i] = vgl_sp3[i] + ck1[i] * a13 + ck2[i] * a23 + ck3[i] * a33 + ck4[i] * a43;
vgl_sp4[i] = vgl_sp4[i] + ck1[i] * a14 + ck2[i] * a24 + ck3[i] * a34 + ck4[i] * a44;
vgl_sp5[i] = vgl_sp5[i] + ck1[i] * a15 + ck2[i] * a25 + ck3[i] * a35 + ck4[i] * a45;
}
}
for (int64_t m=n ; m < nidx ; m+=1) {
const float* restrict ck = coefficient_t_sp + idx[m]*mo_num;
const float a1 = av1[m];
const float a2 = av2[m];
const float a3 = av3[m];
const float a4 = av4[m];
const float a5 = av5[m];
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl_sp1[i] = vgl_sp1[i] + ck[i] * a1;
vgl_sp2[i] = vgl_sp2[i] + ck[i] * a2;
vgl_sp3[i] = vgl_sp3[i] + ck[i] * a3;
vgl_sp4[i] = vgl_sp4[i] + ck[i] * a4;
vgl_sp5[i] = vgl_sp5[i] + ck[i] * a5;
}
}
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] = (double) vgl_sp1[i];
vgl2[i] = (double) vgl_sp2[i];
vgl3[i] = (double) vgl_sp3[i];
vgl4[i] = (double) vgl_sp4[i];
vgl5[i] = (double) vgl_sp5[i];
}
free(vgl_sp1);
free(vgl_sp2);
free(vgl_sp3);
free(vgl_sp4);
free(vgl_sp5);
}
free(idx);
free(av1);
free(av2);
free(av3);
free(av4);
free(av5);
}
free(coefficient_t_sp);
return QMCKL_SUCCESS;
}
#endif
#+end_src
** Computation of cusp-corrected MOs: values only
*** Compute
:PROPERTIES:
:Name: qmckl_compute_mo_basis_mo_value_cusp
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_mo_basis_mo_value_cusp_args
| Variable | Type | In/Out | Description |
|-----------------+-------------------------------+--------+-------------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~ao_num~ | ~int64_t~ | in | Number of AOs |
| ~mo_num~ | ~int64_t~ | in | Number of MOs |
| ~point_num~ | ~int64_t~ | in | Number of points |
| ~ao_nucl~ | ~int64_t[ao_num]~ | in | Nucleus on which the AO is centered |
| ~ao_ang_mom~ | ~int32_t[ao_num]~ | in | Angular momentum of the shell |
| ~en_distance~ | ~double[point_num][nucl_num]~ | in | Electron-nucleus distances |
| ~r_cusp~ | ~double[nucl_num]~ | in | Cusp-adjustment radius |
| ~cusp_param~ | ~double[nucl_num][4][mo_num]~ | in | Cusp-adjustment parameters |
| ~coefficient_t~ | ~double[mo_num][ao_num]~ | in | Transpose of the AO to MO transformation matrix |
| ~ao_value~ | ~double[point_num][ao_num]~ | in | Value of the AOs |
| ~mo_value~ | ~double[point_num][mo_num]~ | out | Cusp correction for the values of the MOs |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_mo_basis_mo_value_cusp_doc_f(context, &
nucl_num, ao_num, mo_num, point_num, ao_nucl, ao_ang_mom, en_distance, &
r_cusp, cusp_param, coefficient_t, ao_value, mo_value) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: nucl_num, ao_num, mo_num, point_num
integer*8 , intent(in) :: ao_nucl(ao_num)
integer*4 , intent(in) :: ao_ang_mom(ao_num)
double precision , intent(in) :: en_distance(nucl_num, point_num)
double precision , intent(in) :: r_cusp(nucl_num)
double precision , intent(in) :: cusp_param(mo_num, 4, nucl_num)
double precision , intent(in) :: coefficient_t(mo_num, ao_num)
double precision , intent(in) :: ao_value(ao_num, point_num)
double precision , intent(out) :: mo_value(mo_num, point_num)
integer*8 :: i, j, k, inucl
double precision :: r
info = QMCKL_SUCCESS
do i=1,point_num
mo_value(:,i) = 0.d0
do k=1,ao_num
if (ao_value(k,i) == 0.d0) cycle
inucl = ao_nucl(k)+1
if ( (en_distance(inucl,i) < r_cusp(inucl)) .and. (ao_ang_mom(k) == 0) ) cycle
mo_value(:,i) = mo_value(:,i) + coefficient_t(:,k) * ao_value(k,i)
end do ! k
do inucl=1,nucl_num
r = en_distance(inucl,i)
if (r > r_cusp(inucl)) cycle
do j=1,mo_num
mo_value(j,i) = mo_value(j,i) + &
cusp_param(j,1,inucl) + r*(cusp_param(j,2,inucl) + r*( &
cusp_param(j,3,inucl) + r* cusp_param(j,4,inucl) ))
enddo
enddo ! inucl
enddo ! i
end function qmckl_compute_mo_basis_mo_value_cusp_doc_f
#+end_src
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_cusp"))
#+RESULTS:
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_value_cusp (
const qmckl_context context,
const int64_t nucl_num,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const int64_t* ao_nucl,
const int32_t* ao_ang_mom,
const double* en_distance,
const double* r_cusp,
const qmckl_tensor cusp_param,
const double* coefficient_t,
const double* ao_value,
double* const mo_value );
#+end_src
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_cusp_doc"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_value_cusp_doc (
const qmckl_context context,
const int64_t nucl_num,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const int64_t* ao_nucl,
const int32_t* ao_ang_mom,
const double* en_distance,
const double* r_cusp,
const double* cusp_param,
const double* coefficient_t,
const double* ao_value,
double* const mo_value );
#+end_src
#+CALL: generate_c_interface(table=qmckl_mo_basis_mo_value_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_cusp_doc"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_mo_basis_mo_value_cusp_doc &
(context, &
nucl_num, &
ao_num, &
mo_num, &
point_num, &
ao_nucl, &
ao_ang_mom, &
en_distance, &
r_cusp, &
cusp_param, &
coefficient_t, &
ao_value, &
mo_value) &
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 :: nucl_num
integer (c_int64_t) , intent(in) , value :: ao_num
integer (c_int64_t) , intent(in) , value :: mo_num
integer (c_int64_t) , intent(in) , value :: point_num
integer (c_int64_t) , intent(in) :: ao_nucl(ao_num)
integer (c_int32_t) , intent(in) :: ao_ang_mom(ao_num)
real (c_double ) , intent(in) :: en_distance(nucl_num,point_num)
real (c_double ) , intent(in) :: r_cusp(nucl_num)
real (c_double ) , intent(in) :: cusp_param(mo_num,4,nucl_num)
real (c_double ) , intent(in) :: coefficient_t(ao_num,mo_num)
real (c_double ) , intent(in) :: ao_value(ao_num,point_num)
real (c_double ) , intent(out) :: mo_value(mo_num,point_num)
integer(c_int32_t), external :: qmckl_compute_mo_basis_mo_value_cusp_doc_f
info = qmckl_compute_mo_basis_mo_value_cusp_doc_f &
(context, &
nucl_num, &
ao_num, &
mo_num, &
point_num, &
ao_nucl, &
ao_ang_mom, &
en_distance, &
r_cusp, &
cusp_param, &
coefficient_t, &
ao_value, &
mo_value)
end function qmckl_compute_mo_basis_mo_value_cusp_doc
#+end_src
#+begin_src c :tangle (eval c) :comments org
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_cusp (const qmckl_context context,
const int64_t nucl_num,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const int64_t* ao_nucl,
const int32_t* ao_ang_mom,
const double* en_distance,
const double* r_cusp,
const qmckl_tensor cusp_param_tensor,
const double* coefficient_t,
const double* ao_value,
double* const mo_value )
{
qmckl_exit_code rc;
#ifdef HAVE_HPC
rc = qmckl_compute_mo_basis_mo_value_cusp_hpc (context, nucl_num, ao_num, mo_num, point_num,
ao_nucl, ao_ang_mom, en_distance, r_cusp,
cusp_param_tensor, coefficient_t, ao_value, mo_value );
#else
double * cusp_param = qmckl_alloc_double_of_tensor(context, cusp_param_tensor);
rc = qmckl_compute_mo_basis_mo_value_cusp_doc (context, nucl_num, ao_num, mo_num, point_num,
ao_nucl, ao_ang_mom, en_distance, r_cusp,
cusp_param, coefficient_t, ao_value, mo_value );
qmckl_free(context, cusp_param);
#endif
return rc;
}
#+end_src
*** HPC version
#+begin_src c :tangle (eval h_private_func) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_cusp_hpc (const qmckl_context context,
const int64_t nucl_num,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const int64_t* ao_nucl,
const int32_t* ao_ang_mom,
const double* en_distance,
const double* r_cusp,
const qmckl_tensor cusp_param,
const double* coefficient_t,
const double* ao_value,
double* const mo_value );
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_cusp_hpc_sp (const qmckl_context context,
const int64_t nucl_num,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const int64_t* ao_nucl,
const int32_t* ao_ang_mom,
const double* en_distance,
const double* r_cusp,
const qmckl_tensor cusp_param,
const double* coefficient_t,
const double* ao_value,
double* const mo_value );
#endif
#+end_src
#+begin_src c :tangle (eval c) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_cusp_hpc (const qmckl_context context,
const int64_t nucl_num,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const int64_t* ao_nucl,
const int32_t* ao_ang_mom,
const double* en_distance,
const double* r_cusp,
const qmckl_tensor cusp_param,
const double* coefficient_t,
const double* ao_value,
double* const mo_value)
{
assert (context != QMCKL_NULL_CONTEXT);
#ifdef HAVE_OPENMP
#pragma omp parallel for
#endif
for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
double* restrict const vgl1 = &(mo_value[ipoint*mo_num]);
const double* restrict avgl1 = &(ao_value[ipoint*ao_num]);
const double* restrict ria = &(en_distance[ipoint*nucl_num]);
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] = 0.;
}
int64_t nidx=0;
int64_t idx[ao_num];
double av1[ao_num];
for (int64_t k=0 ; k<ao_num ; ++k) {
if (avgl1[k] != 0.) {
const int64_t inucl = ao_nucl[k];
if (ria[inucl] > r_cusp[inucl] || ao_ang_mom[k] > 0) {
idx[nidx] = k;
av1[nidx] = avgl1[k];
++nidx;
}
}
}
int64_t n=0;
for (n=0 ; n < nidx-4 ; n+=4) {
const double* restrict ck1 = coefficient_t + idx[n ]*mo_num;
const double* restrict ck2 = coefficient_t + idx[n+1]*mo_num;
const double* restrict ck3 = coefficient_t + idx[n+2]*mo_num;
const double* restrict ck4 = coefficient_t + idx[n+3]*mo_num;
const double a11 = av1[n ];
const double a21 = av1[n+1];
const double a31 = av1[n+2];
const double a41 = av1[n+3];
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] = vgl1[i] + ck1[i] * a11 + ck2[i] * a21 + ck3[i] * a31 + ck4[i] * a41;
}
}
for (int64_t m=n ; m < nidx ; m+=1) {
const double* restrict ck = coefficient_t + idx[m]*mo_num;
const double a1 = av1[m];
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] = vgl1[i] + ck[i] * a1;
}
}
for (int64_t inucl=0 ; inucl<nucl_num ; ++inucl) {
if (ria[inucl] < r_cusp[inucl]) {
const double r = ria[inucl];
IVDEP
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] = vgl1[i] + qmckl_ten3(cusp_param,i,0,inucl) + r*(
qmckl_ten3(cusp_param,i,1,inucl) + r*(
qmckl_ten3(cusp_param,i,2,inucl) + r*(
qmckl_ten3(cusp_param,i,3,inucl) )));
}
}
}
}
return QMCKL_SUCCESS;
}
#endif
#+end_src
** Computation of cusp-corrected MOs: values, gradient, Laplacian
*** Compute
:PROPERTIES:
:Name: qmckl_compute_mo_basis_mo_vgl
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_mo_basis_mo_vgl_cusp_args
| Variable | Type | In/Out | Description |
|-----------------+--------------------------------+--------+-------------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~ao_num~ | ~int64_t~ | in | Number of AOs |
| ~mo_num~ | ~int64_t~ | in | Number of MOs |
| ~point_num~ | ~int64_t~ | in | Number of points |
| ~ao_nucl~ | ~int64_t[ao_num]~ | in | Nucleus on which the AO is centered |
| ~ao_ang_mom~ | ~int32_t[ao_num]~ | in | Angular momentum of the shell |
| ~en_distance~ | ~double[point_num][nucl_num]~ | in | Electron-nucleus distances |
| ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates |
| ~point_coord~ | ~double[3][point_num]~ | in | Electron coordinates |
| ~r_cusp~ | ~double[nucl_num]~ | in | Cusp-adjustment radius |
| ~cusp_param~ | ~double[nucl_num][4][mo_num]~ | in | Cusp-adjustment parameters |
| ~coefficient_t~ | ~double[mo_num][ao_num]~ | in | Transpose of the AO to MO transformation matrix |
| ~ao_vgl~ | ~double[point_num][5][ao_num]~ | in | Value, gradients and Laplacian of the AOs |
| ~mo_vgl~ | ~double[point_num][5][mo_num]~ | out | Value, gradients and Laplacian of the MOs |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f(context, &
nucl_num, ao_num, mo_num, point_num, ao_nucl, ao_ang_mom, en_distance, &
nucl_coord, point_coord, r_cusp, cusp_param, coefficient_t, ao_vgl, mo_vgl) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: nucl_num, ao_num, mo_num, point_num
integer*8 , intent(in) :: ao_nucl(ao_num)
integer*4 , intent(in) :: ao_ang_mom(ao_num)
double precision , intent(in) :: en_distance(nucl_num, point_num)
double precision , intent(in) :: nucl_coord(nucl_num,3)
double precision , intent(in) :: point_coord(point_num,3)
double precision , intent(in) :: r_cusp(nucl_num)
double precision , intent(in) :: cusp_param(mo_num,4,nucl_num)
double precision , intent(in) :: coefficient_t(mo_num,ao_num)
double precision , intent(in) :: ao_vgl(ao_num,5,point_num)
double precision , intent(out) :: mo_vgl(mo_num,5,point_num)
integer*8 :: i,j,k, inucl
double precision :: c1, c2, c3, c4, c5
double precision :: r, r_inv, r_vec(3)
do j=1,point_num
! Initial contribution of the MO
mo_vgl(:,:,j) = 0.d0
do k=1,ao_num
if (ao_vgl(k,1,j) /= 0.d0) then
inucl = ao_nucl(k)+1
if ( (en_distance(inucl,j) > r_cusp(inucl)) .or. (ao_ang_mom(k) > 0) ) then
c1 = ao_vgl(k,1,j)
c2 = ao_vgl(k,2,j)
c3 = ao_vgl(k,3,j)
c4 = ao_vgl(k,4,j)
c5 = ao_vgl(k,5,j)
do i=1,mo_num
mo_vgl(i,1,j) = mo_vgl(i,1,j) + coefficient_t(i,k) * c1
mo_vgl(i,2,j) = mo_vgl(i,2,j) + coefficient_t(i,k) * c2
mo_vgl(i,3,j) = mo_vgl(i,3,j) + coefficient_t(i,k) * c3
mo_vgl(i,4,j) = mo_vgl(i,4,j) + coefficient_t(i,k) * c4
mo_vgl(i,5,j) = mo_vgl(i,5,j) + coefficient_t(i,k) * c5
end do
end if
end if
end do
! Cusp adjustment
do inucl=1,nucl_num
r = en_distance(inucl,j)
if (r < r_cusp(inucl)) then
r_vec(1:3) = point_coord(j,1:3) - nucl_coord(inucl,1:3)
r_inv = 1.d0/r
do i=1,mo_num
mo_vgl(i,1,j) = mo_vgl(i,1,j) + &
cusp_param(i,1,inucl) + r*( &
cusp_param(i,2,inucl) + r*( &
cusp_param(i,3,inucl) + r*( &
cusp_param(i,4,inucl) )))
c1 = r_inv * cusp_param(i,2,inucl) + 2.d0*cusp_param(i,3,inucl) + &
r * 3.d0 * cusp_param(i,4,inucl)
mo_vgl(i,2,j) = mo_vgl(i,2,j) + r_vec(1) * c1
mo_vgl(i,3,j) = mo_vgl(i,3,j) + r_vec(2) * c1
mo_vgl(i,4,j) = mo_vgl(i,4,j) + r_vec(3) * c1
mo_vgl(i,5,j) = mo_vgl(i,5,j) + &
2.d0*cusp_param(i,2,inucl)*r_inv + &
6.d0*cusp_param(i,3,inucl) + &
12.d0*cusp_param(i,4,inucl)*r
end do
end if
end do ! inucl
end do
info = QMCKL_SUCCESS
end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f
#+end_src
# #+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_cusp"))
#+RESULTS:
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_cusp (
const qmckl_context context,
const int64_t nucl_num,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const int64_t* ao_nucl,
const int32_t* ao_ang_mom,
const double* en_distance,
const qmckl_matrix nucl_coord,
const qmckl_matrix point_coord,
const double* r_cusp,
const qmckl_tensor cusp_param,
const double* coefficient_t,
const double* ao_vgl,
double* const mo_vgl );
#+end_src
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_cusp_doc"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_cusp_doc (
const qmckl_context context,
const int64_t nucl_num,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const int64_t* ao_nucl,
const int32_t* ao_ang_mom,
const double* en_distance,
const double* nucl_coord,
const double* point_coord,
const double* r_cusp,
const double* cusp_param,
const double* coefficient_t,
const double* ao_vgl,
double* const mo_vgl );
#+end_src
#+CALL: generate_c_interface(table=qmckl_mo_basis_mo_vgl_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_cusp_doc"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_mo_basis_mo_vgl_cusp_doc &
(context, &
nucl_num, &
ao_num, &
mo_num, &
point_num, &
ao_nucl, &
ao_ang_mom, &
en_distance, &
nucl_coord, &
point_coord, &
r_cusp, &
cusp_param, &
coefficient_t, &
ao_vgl, &
mo_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 :: nucl_num
integer (c_int64_t) , intent(in) , value :: ao_num
integer (c_int64_t) , intent(in) , value :: mo_num
integer (c_int64_t) , intent(in) , value :: point_num
integer (c_int64_t) , intent(in) :: ao_nucl(ao_num)
integer (c_int32_t) , intent(in) :: ao_ang_mom(ao_num)
real (c_double ) , intent(in) :: en_distance(nucl_num,point_num)
real (c_double ) , intent(in) :: nucl_coord(nucl_num,3)
real (c_double ) , intent(in) :: point_coord(point_num,3)
real (c_double ) , intent(in) :: r_cusp(nucl_num)
real (c_double ) , intent(in) :: cusp_param(mo_num,4,nucl_num)
real (c_double ) , intent(in) :: coefficient_t(ao_num,mo_num)
real (c_double ) , intent(in) :: ao_vgl(ao_num,5,point_num)
real (c_double ) , intent(out) :: mo_vgl(mo_num,5,point_num)
integer(c_int32_t), external :: qmckl_compute_mo_basis_mo_vgl_cusp_doc_f
info = qmckl_compute_mo_basis_mo_vgl_cusp_doc_f &
(context, &
nucl_num, &
ao_num, &
mo_num, &
point_num, &
ao_nucl, &
ao_ang_mom, &
en_distance, &
nucl_coord, &
point_coord, &
r_cusp, &
cusp_param, &
coefficient_t, &
ao_vgl, &
mo_vgl)
end function qmckl_compute_mo_basis_mo_vgl_cusp_doc
#+end_src
#+begin_src c :tangle (eval c) :comments org
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl_cusp (const qmckl_context context,
const int64_t nucl_num,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const int64_t* ao_nucl,
const int32_t* ao_ang_mom,
const double* en_distance,
const qmckl_matrix nucl_coord_matrix,
const qmckl_matrix point_coord_matrix,
const double* r_cusp,
const qmckl_tensor cusp_param_tensor,
const double* coefficient_t,
const double* ao_vgl,
double* const mo_vgl )
{
qmckl_exit_code rc;
#ifdef HAVE_HPC
rc = qmckl_compute_mo_basis_mo_vgl_cusp_hpc (context, nucl_num, ao_num, mo_num, point_num,
ao_nucl, ao_ang_mom, en_distance, nucl_coord_matrix,
point_coord_matrix, r_cusp, cusp_param_tensor,
coefficient_t, ao_vgl, mo_vgl );
#else
double * nucl_coord = qmckl_alloc_double_of_matrix(context, nucl_coord_matrix);
double * point_coord = qmckl_alloc_double_of_matrix(context, point_coord_matrix);
double * cusp_param = qmckl_alloc_double_of_tensor(context, cusp_param_tensor);
rc = qmckl_compute_mo_basis_mo_vgl_cusp_doc (context, nucl_num, ao_num, mo_num, point_num,
ao_nucl, ao_ang_mom, en_distance, nucl_coord,
point_coord, r_cusp, cusp_param, coefficient_t,
ao_vgl, mo_vgl );
qmckl_free(context, nucl_coord);
qmckl_free(context, point_coord);
qmckl_free(context, cusp_param);
#endif
return rc;
}
#+end_src
*** HPC version
#+begin_src c :tangle (eval h_private_func) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl_cusp_hpc (const qmckl_context context,
const int64_t nucl_num,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const int64_t* ao_nucl,
const int32_t* ao_ang_mom,
const double* en_distance,
const qmckl_matrix nucl_coord,
const qmckl_matrix point_coord,
const double* r_cusp,
const qmckl_tensor cusp_param,
const double* coefficient_t,
const double* ao_vgl,
double* const mo_vgl );
#endif
#+end_src
#+begin_src c :tangle (eval c) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl_cusp_hpc (const qmckl_context context,
const int64_t nucl_num,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const int64_t* ao_nucl,
const int32_t* ao_ang_mom,
const double* en_distance,
const qmckl_matrix nucl_coord,
const qmckl_matrix point_coord,
const double* r_cusp,
const qmckl_tensor cusp_param,
const double* coefficient_t,
const double* ao_vgl,
double* const mo_vgl )
{
assert (context != QMCKL_NULL_CONTEXT);
#ifdef HAVE_OPENMP
#pragma omp parallel for
#endif
for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
double* restrict const vgl1 = &(mo_vgl[ipoint*5*mo_num]);
double* restrict const vgl2 = vgl1 + mo_num;
double* restrict const vgl3 = vgl1 + (mo_num << 1);
double* restrict const vgl4 = vgl1 + (mo_num << 1) + mo_num;
double* restrict const vgl5 = vgl1 + (mo_num << 2);
const double* restrict avgl1 = &(ao_vgl[ipoint*5*ao_num]);
const double* restrict avgl2 = avgl1 + ao_num;
const double* restrict avgl3 = avgl1 + (ao_num << 1);
const double* restrict avgl4 = avgl1 + (ao_num << 1) + ao_num;
const double* restrict avgl5 = avgl1 + (ao_num << 2);
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] = 0.;
vgl2[i] = 0.;
vgl3[i] = 0.;
vgl4[i] = 0.;
vgl5[i] = 0.;
}
const double* restrict ria = &(en_distance[ipoint*nucl_num]);
int64_t nidx=0;
int64_t idx[ao_num];
double av1[ao_num];
double av2[ao_num];
double av3[ao_num];
double av4[ao_num];
double av5[ao_num];
for (int64_t k=0 ; k<ao_num ; ++k) {
if (avgl1[k] != 0.) {
const int64_t inucl = ao_nucl[k];
if (ria[inucl] > r_cusp[inucl] || ao_ang_mom[k] > 0) {
idx[nidx] = k;
av1[nidx] = avgl1[k];
av2[nidx] = avgl2[k];
av3[nidx] = avgl3[k];
av4[nidx] = avgl4[k];
av5[nidx] = avgl5[k];
++nidx;
}
}
}
int64_t n=0;
for (n=0 ; n < nidx-4 ; n+=4) {
const double* restrict ck1 = coefficient_t + idx[n ]*mo_num;
const double* restrict ck2 = coefficient_t + idx[n+1]*mo_num;
const double* restrict ck3 = coefficient_t + idx[n+2]*mo_num;
const double* restrict ck4 = coefficient_t + idx[n+3]*mo_num;
const double a11 = av1[n ];
const double a21 = av1[n+1];
const double a31 = av1[n+2];
const double a41 = av1[n+3];
const double a12 = av2[n ];
const double a22 = av2[n+1];
const double a32 = av2[n+2];
const double a42 = av2[n+3];
const double a13 = av3[n ];
const double a23 = av3[n+1];
const double a33 = av3[n+2];
const double a43 = av3[n+3];
const double a14 = av4[n ];
const double a24 = av4[n+1];
const double a34 = av4[n+2];
const double a44 = av4[n+3];
const double a15 = av5[n ];
const double a25 = av5[n+1];
const double a35 = av5[n+2];
const double a45 = av5[n+3];
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] = vgl1[i] + ck1[i] * a11 + ck2[i] * a21 + ck3[i] * a31 + ck4[i] * a41;
vgl2[i] = vgl2[i] + ck1[i] * a12 + ck2[i] * a22 + ck3[i] * a32 + ck4[i] * a42;
vgl3[i] = vgl3[i] + ck1[i] * a13 + ck2[i] * a23 + ck3[i] * a33 + ck4[i] * a43;
vgl4[i] = vgl4[i] + ck1[i] * a14 + ck2[i] * a24 + ck3[i] * a34 + ck4[i] * a44;
vgl5[i] = vgl5[i] + ck1[i] * a15 + ck2[i] * a25 + ck3[i] * a35 + ck4[i] * a45;
}
}
for (int64_t m=n ; m < nidx ; m+=1) {
const double* restrict ck = coefficient_t + idx[m]*mo_num;
const double a1 = av1[m];
const double a2 = av2[m];
const double a3 = av3[m];
const double a4 = av4[m];
const double a5 = av5[m];
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] = vgl1[i] + ck[i] * a1;
vgl2[i] = vgl2[i] + ck[i] * a2;
vgl3[i] = vgl3[i] + ck[i] * a3;
vgl4[i] = vgl4[i] + ck[i] * a4;
vgl5[i] = vgl5[i] + ck[i] * a5;
}
}
// TODO
for (int64_t inucl=0 ; inucl<nucl_num ; ++inucl) {
const double r = ria[inucl];
if (r < r_cusp[inucl]) {
const double r_vec[3] = {
qmckl_mat(point_coord,ipoint,0) - qmckl_mat(nucl_coord,inucl,0),
qmckl_mat(point_coord,ipoint,1) - qmckl_mat(nucl_coord,inucl,1),
qmckl_mat(point_coord,ipoint,2) - qmckl_mat(nucl_coord,inucl,2) };
const double r_inv = 1./r;
for (int64_t i=0 ; i<mo_num ; ++i) {
vgl1[i] = vgl1[i] +
qmckl_ten3(cusp_param,i,0,inucl) + r*(
qmckl_ten3(cusp_param,i,1,inucl) + r*(
qmckl_ten3(cusp_param,i,2,inucl) + r*(
qmckl_ten3(cusp_param,i,3,inucl) )));
const double c1 = r_inv * qmckl_ten3(cusp_param,i,1,inucl) +
2.0*qmckl_ten3(cusp_param,i,2,inucl) +
r * 3.0 * qmckl_ten3(cusp_param,i,3,inucl);
vgl2[i] = vgl2[i] + r_vec[0] * c1;
vgl3[i] = vgl3[i] + r_vec[1] * c1;
vgl4[i] = vgl4[i] + r_vec[2] * c1;
vgl5[i] = vgl5[i] + 2.0*qmckl_ten3(cusp_param,i,1,inucl)*r_inv +
6.0*qmckl_ten3(cusp_param,i,2,inucl) +
12.0*qmckl_ten3(cusp_param,i,3,inucl)*r;
}
}
}
}
return QMCKL_SUCCESS;
}
#endif
#+end_src
** Rescaling of MO coefficients
When evaluating Slater determinants, the value of the determinants
may get out of the range of double precision. A simple fix is to
rescale the MO coefficients to put back the determinants in the
correct range.
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_mo_basis_rescale(qmckl_context context,
const double scaling_factor);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_mo_basis_rescale(qmckl_context context,
const double scaling_factor)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_NULL_CONTEXT,
"qmckl_mo_basis_rescale",
NULL);
}
if (scaling_factor == 0.) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_mo_basis_rescale",
"scaling factor can't be zero");
}
qmckl_exit_code rc;
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->mo_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_mo_basis_rescale",
NULL);
}
for (int64_t i=0 ; i<ctx->ao_basis.ao_num * ctx->mo_basis.mo_num ; ++i) {
ctx->mo_basis.coefficient[i] *= scaling_factor;
ctx->mo_basis.coefficient_t[i] *= scaling_factor;
}
rc = qmckl_context_touch(context);
return rc;
}
#+end_src
*** Fortran interface
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(qmckl_exit_code) function qmckl_mo_basis_rescale (context, &
scaling_factor) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in), value :: context
real (c_double) , intent(in), value :: scaling_factor
end function qmckl_mo_basis_rescale
end interface
#+end_src
** Test
#+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][0][26] : %25.15e"% f(a,x,y))
print ( "[1][1][0][26] : %25.15e"% df(a,x,y,1))
print ( "[1][2][0][26] : %25.15e"% df(a,x,y,2))
print ( "[1][3][0][26] : %25.15e"% df(a,x,y,3))
print ( "[1][4][0][26] : %25.15e"% lf(a,x,y))
x = elec_15_w2 ; y = nucl_2
a = [(3.387000E+01, 6.068000E-03 *1.0006253235944540e+01),
(5.095000E+00, 4.530800E-02 *2.4169531573445120e+00),
(1.159000E+00, 2.028220E-01 *7.9610924849766440e-01),
(3.258000E-01, 5.039030E-01 *3.0734305383061117e-01),
(1.027000E-01, 3.834210E-01 *1.2929684417481876e-01)]
print ( "[0][1][15][14] : %25.15e"% f(a,x,y))
print ( "[1][1][15][14] : %25.15e"% df(a,x,y,1))
print ( "[2][1][15][14] : %25.15e"% df(a,x,y,2))
print ( "[3][1][15][14] : %25.15e"% df(a,x,y,3))
print ( "[4][1][15][14] : %25.15e"% lf(a,x,y))
#+end_src
#+begin_src c :tangle (eval c_test) :exports none
{
#define walk_num chbrclf_walk_num
#define elec_num chbrclf_elec_num
#define shell_num chbrclf_shell_num
#define ao_num chbrclf_ao_num
int64_t elec_up_num = chbrclf_elec_up_num;
int64_t elec_dn_num = chbrclf_elec_dn_num;
double* elec_coord = &(chbrclf_elec_coord[0][0][0]);
int64_t nucl_num = chbrclf_nucl_num;
const double* nucl_charge = chbrclf_charge;
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
int64_t point_num = walk_num*elec_num;
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
rc = qmckl_set_point(context, 'N', point_num, elec_coord, point_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), nucl_num*3);
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 * 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, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_prim_num (context, chbrclf_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, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_factor (context, shell_factor, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_exponent (context, exponent, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_coefficient (context, coefficient, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_prim_factor (context, prim_factor, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_ao_basis_ao_num(context, chbrclf_ao_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_ao_basis_ao_factor (context, ao_factor, chbrclf_ao_num);
assert(rc == QMCKL_SUCCESS);
assert(qmckl_ao_basis_provided(context));
double ao_vgl[point_num][5][chbrclf_ao_num];
rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]),
(int64_t) 5*point_num*chbrclf_ao_num);
assert (rc == QMCKL_SUCCESS);
/* Set up MO data */
int64_t mo_num = chbrclf_mo_num;
rc = qmckl_set_mo_basis_mo_num(context, mo_num);
assert (rc == QMCKL_SUCCESS);
const double * mo_coefficient = &(chbrclf_mo_coef[0]);
rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient, chbrclf_mo_num*chbrclf_ao_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_mo_basis_provided(context));
rc = qmckl_context_touch(context);
assert (rc == QMCKL_SUCCESS);
double mo_value[point_num][chbrclf_mo_num];
rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), point_num*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
double mo_vgl[point_num][5][chbrclf_mo_num];
rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), point_num*5*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
for (int i=0 ; i< point_num; ++i) {
for (int k=0 ; k< chbrclf_mo_num ; ++k) {
assert(fabs(mo_vgl[i][0][k] - mo_value[i][k]) < 1.e-12) ;
}
}
rc = qmckl_context_touch(context);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), point_num*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
for (int i=0 ; i< point_num; ++i) {
for (int k=0 ; k< chbrclf_mo_num ; ++k) {
assert(fabs(mo_vgl[i][0][k] - mo_value[i][k]) < 1.e-12) ;
}
}
rc = qmckl_mo_basis_rescale(context, 0.);
assert (rc != QMCKL_SUCCESS);
rc = qmckl_mo_basis_rescale(context, 2.);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), point_num*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
for (int i=0 ; i< point_num; ++i) {
for (int k=0 ; k< chbrclf_mo_num ; ++k) {
assert(fabs(2.*mo_vgl[i][0][k] - mo_value[i][k]) < 1.e-12) ;
}
}
rc = qmckl_mo_basis_rescale(context, 0.5);
assert (rc == QMCKL_SUCCESS);
printf("\n");
printf(" mo_vgl mo_vgl[0][26][219] %25.15e\n", mo_vgl[2][0][3]);
printf(" mo_vgl mo_vgl[1][26][219] %25.15e\n", mo_vgl[2][1][3]);
printf(" mo_vgl mo_vgl[0][26][220] %25.15e\n", mo_vgl[2][0][3]);
printf(" mo_vgl mo_vgl[1][26][220] %25.15e\n", mo_vgl[2][1][3]);
printf(" mo_vgl mo_vgl[0][26][221] %25.15e\n", mo_vgl[2][0][3]);
printf(" mo_vgl mo_vgl[1][26][221] %25.15e\n", mo_vgl[2][1][3]);
printf(" mo_vgl mo_vgl[0][26][222] %25.15e\n", mo_vgl[2][0][3]);
printf(" mo_vgl mo_vgl[1][26][222] %25.15e\n", mo_vgl[2][1][3]);
printf(" mo_vgl mo_vgl[0][26][223] %25.15e\n", mo_vgl[2][0][3]);
printf(" mo_vgl mo_vgl[1][26][223] %25.15e\n", mo_vgl[2][1][3]);
printf(" mo_vgl mo_vgl[0][26][224] %25.15e\n", mo_vgl[2][0][3]);
printf(" mo_vgl mo_vgl[1][26][224] %25.15e\n", mo_vgl[2][1][3]);
printf("\n");
/* Check single precision */
rc = qmckl_set_numprec_precision(context,53);
rc = qmckl_context_touch(context);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), point_num*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), point_num*5*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_numprec_precision(context,23);
rc = qmckl_context_touch(context);
assert (rc == QMCKL_SUCCESS);
double mo_value_sp[point_num][chbrclf_mo_num];
rc = qmckl_get_mo_basis_mo_value(context, &(mo_value_sp[0][0]), point_num*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
double mo_vgl_sp[point_num][5][chbrclf_mo_num];
rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl_sp[0][0][0]), point_num*5*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
uint64_t average_prec = 0;
int32_t nbits;
for (int i=0 ; i< point_num; ++i) {
for (int k=0 ; k< chbrclf_mo_num ; ++k) {
nbits = qmckl_test_precision_64(mo_value_sp[i][k], mo_value[i][k]);
// printf("%d %d %25.15e %25.15e %d\n", i, k, mo_value_sp[i][k], mo_value[i][k], nbits);
average_prec += nbits;
}
}
printf("Average precision for %d: %d\n", qmckl_get_numprec_precision(context),
(int) (average_prec/(point_num*chbrclf_mo_num)));
assert(nbits > 12);
fflush(stdout);
average_prec = 0;
for (int i=0 ; i< point_num; ++i) {
for (int k=0 ; k< chbrclf_mo_num ; ++k) {
// printf("%d %d\n", i, k);
nbits = qmckl_test_precision_64(mo_vgl_sp[i][0][k], mo_vgl[i][0][k]);
average_prec += nbits;
// printf("%25.15e %25.15e %d\n", mo_vgl_sp[i][0][k], mo_vgl[i][0][k], nbits);
nbits = qmckl_test_precision_64(mo_vgl_sp[i][1][k], mo_vgl[i][1][k]);
average_prec += nbits;
// printf("%25.15e %25.15e %d\n", mo_vgl_sp[i][1][k], mo_vgl[i][1][k], nbits);
nbits = qmckl_test_precision_64(mo_vgl_sp[i][2][k], mo_vgl[i][2][k]);
average_prec += nbits;
// printf("%25.15e %25.15e %d\n", mo_vgl_sp[i][2][k], mo_vgl[i][2][k], nbits);
nbits = qmckl_test_precision_64(mo_vgl_sp[i][3][k], mo_vgl[i][3][k]);
average_prec += nbits;
// printf("%25.15e %25.15e %d\n", mo_vgl_sp[i][3][k], mo_vgl[i][3][k], nbits);
nbits = qmckl_test_precision_64(mo_vgl_sp[i][4][k], mo_vgl[i][4][k]);
average_prec += nbits;
// printf("%25.15e %25.15e %d\n", mo_vgl_sp[i][4][k], mo_vgl[i][4][k], nbits);
}
}
nbits = (int) (average_prec/(point_num*chbrclf_mo_num*5));
printf("Average precision for %d: %d\n", qmckl_get_numprec_precision(context), nbits);
assert(nbits > 11);
fflush(stdout);
rc = qmckl_set_numprec_precision(context,53);
/* Check selection of MOs */
int32_t keep[mo_num];
for (int i=0 ; i<mo_num ; ++i) {
keep[i] = 0;
}
keep[2] = 1;
keep[5] = 1;
rc = qmckl_mo_basis_select_mo(context, &(keep[0]), mo_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_get_mo_basis_mo_num(context, &mo_num);
printf(" mo_num: %ld\n", (long) mo_num);
assert(mo_num == 2);
double mo_coefficient_new[mo_num][ao_num];
rc = qmckl_get_mo_basis_coefficient (context, &(mo_coefficient_new[0][0]), mo_num*ao_num);
for (int i=0 ; i<ao_num ; ++i) {
assert(mo_coefficient_new[0][i] == mo_coefficient[i + ao_num*2]);
assert(mo_coefficient_new[1][i] == mo_coefficient[i + ao_num*5]);
}
}
#+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
# -*- mode: org -*-
# vim: syntax=c
Reference in New Issue View Git Blame Copy Permalink