mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-11-03 20:54:09 +01:00
3548 lines
125 KiB
Org Mode
3548 lines
125 KiB
Org Mode
#+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];
|
|
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.");
|
|
}
|
|
|
|
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);
|
|
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
|
|
IVDEP
|
|
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
|
|
IVDEP
|
|
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
|
|
IVDEP
|
|
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
|
|
IVDEP
|
|
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
|
|
IVDEP
|
|
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
|
|
IVDEP
|
|
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
|
|
|
|
! info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, 1.d0, &
|
|
! coefficient_t, int(size(coefficient_t,1),8), &
|
|
! ao_vgl, int(size(ao_vgl,1),8), 0.d0, &
|
|
! mo_vgl, int(size(mo_vgl,1),8))
|
|
|
|
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
|
|
IVDEP
|
|
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
|
|
IVDEP
|
|
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
|
|
IVDEP
|
|
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
|
|
IVDEP
|
|
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
|
|
IVDEP
|
|
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
|
|
IVDEP
|
|
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
|
|
IVDEP
|
|
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
|
|
IVDEP
|
|
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) cycle
|
|
inucl = ao_nucl(k)+1
|
|
if ( (en_distance(inucl,j) < r_cusp(inucl)) .and. (ao_ang_mom(k) == 0) ) cycle
|
|
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 do
|
|
|
|
! Cusp adjustment
|
|
do inucl=1,nucl_num
|
|
r = en_distance(inucl,j)
|
|
if (r > r_cusp(inucl)) cycle
|
|
|
|
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
|
|
|
|
enddo
|
|
enddo ! 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
|
|
IVDEP
|
|
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
|
|
IVDEP
|
|
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) {
|
|
if (ria[inucl] < r_cusp[inucl]) {
|
|
const double r = ria[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;
|
|
|
|
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) )));
|
|
|
|
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
|
|
|
|
|