mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-11-04 21:24:11 +01:00
3018 lines
102 KiB
Org Mode
3018 lines
102 KiB
Org Mode
#+TITLE: Electrons
|
|
#+SETUPFILE: ../tools/theme.setup
|
|
#+INCLUDE: ../tools/lib.org
|
|
|
|
In conventional QMC simulations, up-spin and down-spin electrons are
|
|
different. The ~electron~ data structure contains the number of
|
|
up-spin and down-spin electrons, and the electron coordinates.
|
|
|
|
* Headers :noexport:
|
|
#+begin_src elisp :noexport :results none
|
|
(org-babel-lob-ingest "../tools/lib.org")
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :tangle (eval h_private_type)
|
|
#ifndef QMCKL_ELECTRON_HPT
|
|
#define QMCKL_ELECTRON_HPT
|
|
#include <stdbool.h>
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval h_private_func)
|
|
#ifndef QMCKL_ELECTRON_HPF
|
|
#define QMCKL_ELECTRON_HPF
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval c_test) :noweb yes
|
|
#include "qmckl.h"
|
|
#include <assert.h>
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
#ifdef HAVE_CONFIG_H
|
|
#include "config.h"
|
|
#endif
|
|
|
|
#include "chbrclf.h"
|
|
|
|
int main() {
|
|
qmckl_context context;
|
|
context = qmckl_context_create();
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval c)
|
|
#ifdef HAVE_CONFIG_H
|
|
#include "config.h"
|
|
#endif
|
|
|
|
#ifdef HAVE_STDINT_H
|
|
#include <stdint.h>
|
|
#elif HAVE_INTTYPES_H
|
|
#include <inttypes.h>
|
|
#endif
|
|
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
#include <stdbool.h>
|
|
#include <assert.h>
|
|
#include <math.h>
|
|
|
|
#include <stdio.h>
|
|
|
|
#include "qmckl.h"
|
|
#include "qmckl_context_private_type.h"
|
|
#include "qmckl_memory_private_type.h"
|
|
#include "qmckl_memory_private_func.h"
|
|
#include "qmckl_electron_private_func.h"
|
|
#+end_src
|
|
|
|
* Context
|
|
|
|
The following data stored in the context:
|
|
|
|
| Variable | Type | Description |
|
|
|---------------------------+----------------+---------------------------------------------------------------|
|
|
| ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data |
|
|
| ~num~ | ~int64_t~ | Total number of electrons |
|
|
| ~up_num~ | ~int64_t~ | Number of up-spin electrons |
|
|
| ~down_num~ | ~int64_t~ | Number of down-spin electrons |
|
|
| ~walk_num~ | ~int64_t~ | Number of walkers |
|
|
| ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor |
|
|
| ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor |
|
|
| ~provided~ | ~bool~ | If true, ~electron~ is valid |
|
|
| ~coord_new~ | ~qmckl_matrix~ | Current set of electron coordinates. Pointer to ~ctx->points~ |
|
|
| ~coord_old~ | ~qmckl_matrix~ | Old set of electron coordinates |
|
|
| ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates |
|
|
|
|
Computed data:
|
|
|
|
| Variable | Type | Description |
|
|
|-------------------------------------+--------------------------------------+----------------------------------------------------------------------|
|
|
| ~ee_distance~ | ~double[walk_num][num][num]~ | Electron-electron distances |
|
|
| ~ee_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
|
|
| ~en_distance~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances |
|
|
| ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
|
|
| ~ee_distance_rescaled~ | ~double[walk_num][num][num]~ | Electron-electron rescaled distances |
|
|
| ~ee_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
|
|
| ~ee_distance_rescaled_deriv_e~ | ~double[walk_num][4][num][num]~ | Electron-electron rescaled distances derivatives |
|
|
| ~ee_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives |
|
|
| ~ee_pot~ | ~double[walk_num]~ | Electron-electron rescaled distances derivatives |
|
|
| ~ee_pot_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives |
|
|
| ~en_pot~ | double[walk_num] | Electron-nucleus potential energy |
|
|
| ~en_pot_date~ | int64_t | Date when the electron-nucleus potential energy was computed |
|
|
| ~en_distance_rescaled~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances |
|
|
| ~en_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
|
|
| ~en_distance_rescaled_deriv_e~ | ~double[walk_num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives |
|
|
| ~en_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives |
|
|
|
|
** Data structure
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_type)
|
|
typedef struct qmckl_electron_struct {
|
|
int64_t num;
|
|
int64_t up_num;
|
|
int64_t down_num;
|
|
int64_t walk_num;
|
|
double rescale_factor_kappa_ee;
|
|
double rescale_factor_kappa_en;
|
|
int64_t coord_new_date;
|
|
int64_t ee_distance_date;
|
|
int64_t en_distance_date;
|
|
int64_t ee_pot_date;
|
|
int64_t en_pot_date;
|
|
int64_t ee_distance_rescaled_date;
|
|
int64_t ee_distance_rescaled_deriv_e_date;
|
|
int64_t en_distance_rescaled_date;
|
|
int64_t en_distance_rescaled_deriv_e_date;
|
|
qmckl_matrix coord_new;
|
|
qmckl_matrix coord_old;
|
|
double* ee_distance;
|
|
double* en_distance;
|
|
double* ee_pot;
|
|
double* en_pot;
|
|
double* ee_distance_rescaled;
|
|
double* ee_distance_rescaled_deriv_e;
|
|
double* en_distance_rescaled;
|
|
double* en_distance_rescaled_deriv_e;
|
|
int32_t uninitialized;
|
|
bool provided;
|
|
} qmckl_electron_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_electron(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c)
|
|
qmckl_exit_code qmckl_init_electron(qmckl_context context) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return false;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
ctx->electron.uninitialized = (1 << 2) - 1;
|
|
|
|
/* Default values */
|
|
ctx->electron.rescale_factor_kappa_ee = 1.0;
|
|
ctx->electron.rescale_factor_kappa_en = 1.0;
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
bool qmckl_electron_provided (const qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
bool qmckl_electron_provided(const qmckl_context context) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return false;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
return ctx->electron.provided;
|
|
}
|
|
#+end_src
|
|
|
|
** Access functions
|
|
|
|
Access functions return ~QMCKL_SUCCESS~ when the data has been
|
|
successfully retrieved. It returnes ~QMCKL_INVALID_CONTEXT~ when
|
|
the context is not a valid context, and ~QMCKL_NOT_PROVIDED~ when
|
|
the data has not been provided. If the function returns
|
|
successfully, the variable pointed by the pointer given in argument
|
|
contains the requested data. Otherwise, this variable is untouched.
|
|
|
|
#+NAME:post
|
|
#+begin_src c :exports none
|
|
if ( (ctx->electron.uninitialized & mask) != 0) {
|
|
return NULL;
|
|
}
|
|
#+end_src
|
|
|
|
*** Number of electrons
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :exports none
|
|
qmckl_exit_code qmckl_get_electron_num (const qmckl_context context, int64_t* const num);
|
|
qmckl_exit_code qmckl_get_electron_up_num (const qmckl_context context, int64_t* const up_num);
|
|
qmckl_exit_code qmckl_get_electron_down_num (const qmckl_context context, int64_t* const down_num);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_electron_num (const qmckl_context context, int64_t* const num) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_INVALID_CONTEXT;
|
|
}
|
|
|
|
if (num == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_electron_num",
|
|
"num is a null pointer");
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 0;
|
|
|
|
if ( (ctx->electron.uninitialized & mask) != 0) {
|
|
return QMCKL_NOT_PROVIDED;
|
|
}
|
|
|
|
assert (ctx->electron.num > (int64_t) 0);
|
|
,*num = ctx->electron.num;
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
|
|
qmckl_exit_code
|
|
qmckl_get_electron_up_num (const qmckl_context context, int64_t* const up_num) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_INVALID_CONTEXT;
|
|
}
|
|
|
|
if (up_num == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_electron_up_num",
|
|
"up_num is a null pointer");
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 0;
|
|
|
|
if ( (ctx->electron.uninitialized & mask) != 0) {
|
|
return QMCKL_NOT_PROVIDED;
|
|
}
|
|
|
|
assert (ctx->electron.up_num > (int64_t) 0);
|
|
,*up_num = ctx->electron.up_num;
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
|
|
qmckl_exit_code
|
|
qmckl_get_electron_down_num (const qmckl_context context, int64_t* const down_num) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_INVALID_CONTEXT;
|
|
}
|
|
|
|
if (down_num == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_electron_down_num",
|
|
"down_num is a null pointer");
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 0;
|
|
|
|
if ( (ctx->electron.uninitialized & mask) != 0) {
|
|
return QMCKL_NOT_PROVIDED;
|
|
}
|
|
|
|
assert (ctx->electron.down_num >= (int64_t) 0);
|
|
,*down_num = ctx->electron.down_num;
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
#+end_src
|
|
|
|
*** Number of walkers
|
|
|
|
A walker is a set of electron coordinates that are arguments of
|
|
the wave function. ~walk_num~ is the number of walkers.
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :exports none
|
|
qmckl_exit_code qmckl_get_electron_walk_num (const qmckl_context context, int64_t* const walk_num);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_electron_walk_num (const qmckl_context context, int64_t* const walk_num) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_INVALID_CONTEXT;
|
|
}
|
|
|
|
if (walk_num == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_electron_walk_num",
|
|
"walk_num is a null pointer");
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 1;
|
|
|
|
if ( (ctx->electron.uninitialized & mask) != 0) {
|
|
return QMCKL_NOT_PROVIDED;
|
|
}
|
|
|
|
assert (ctx->electron.walk_num > (int64_t) 0);
|
|
,*walk_num = ctx->electron.walk_num;
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Scaling factors Kappa
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :exports none
|
|
qmckl_exit_code qmckl_get_electron_rescale_factor_ee (const qmckl_context context, double* const rescale_factor_kappa_ee);
|
|
qmckl_exit_code qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const rescale_factor_kappa_en);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_electron_rescale_factor_ee (const qmckl_context context, double* const rescale_factor_kappa_ee) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_INVALID_CONTEXT;
|
|
}
|
|
|
|
if (rescale_factor_kappa_ee == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_electron_rescale_factor_ee",
|
|
"rescale_factor_kappa_ee is a null pointer");
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
assert (ctx->electron.rescale_factor_kappa_ee > 0.0);
|
|
|
|
*rescale_factor_kappa_ee = ctx->electron.rescale_factor_kappa_ee;
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
|
|
qmckl_exit_code
|
|
qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const rescale_factor_kappa_en) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_INVALID_CONTEXT;
|
|
}
|
|
|
|
if (rescale_factor_kappa_en == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_electron_rescale_factor_en",
|
|
"rescale_factor_kappa_en is a null pointer");
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
assert (ctx->electron.rescale_factor_kappa_en > 0.0);
|
|
*rescale_factor_kappa_en = ctx->electron.rescale_factor_kappa_en;
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Electron coordinates
|
|
|
|
Returns the current electron coordinates. The pointer is assumed
|
|
to point on a memory block of size ~size_max~ \ge ~3 * elec_num * walk_num~.
|
|
The order of the indices is:
|
|
|
|
| | Normal | Transposed |
|
|
|---------+--------------------------+--------------------------|
|
|
| C | ~[walk_num*elec_num][3]~ | ~[3][walk_num*elec_num]~ |
|
|
| Fortran | ~(3,walk_num*elec_num)~ | ~(walk_num*elec_num, 3)~ |
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_electron_coord (const qmckl_context context,
|
|
const char transp,
|
|
double* const coord,
|
|
const int64_t size_max);
|
|
#+end_src
|
|
|
|
As the ~coord_new~ attribute is a pointer equal to ~points~,
|
|
returning the current electron coordinates is equivalent to
|
|
returning the current points.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_get_electron_coord (const qmckl_context context,
|
|
const char transp,
|
|
double* const coord,
|
|
const int64_t size_max)
|
|
{
|
|
if (transp != 'N' && transp != 'T') {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_electron_coord",
|
|
"transp should be 'N' or 'T'");
|
|
}
|
|
|
|
if (coord == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_get_electron_coord",
|
|
"coord is a null pointer");
|
|
}
|
|
|
|
if (size_max <= 0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_4,
|
|
"qmckl_get_electron_coord",
|
|
"size_max should be > 0");
|
|
}
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_INVALID_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
if (!ctx->electron.provided) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_get_electron_coord",
|
|
NULL);
|
|
}
|
|
|
|
assert (ctx->point.coord.data == ctx->electron.coord_new.data);
|
|
assert (ctx->point.coord.size[0] == ctx->electron.coord_new.size[0]);
|
|
assert (ctx->point.coord.size[1] == ctx->electron.coord_new.size[1]);
|
|
|
|
return qmckl_get_point(context, transp, coord, size_max);
|
|
}
|
|
|
|
#+end_src
|
|
|
|
** Initialization functions
|
|
|
|
To set the data relative to the electrons in the context, the
|
|
following functions need to be called. When the data structure is
|
|
initialized, the internal ~coord_new~ and ~coord_old~ arrays are
|
|
both allocated.
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code qmckl_set_electron_num (qmckl_context context, const int64_t up_num, const int64_t down_num);
|
|
qmckl_exit_code qmckl_set_electron_walk_num (qmckl_context context, const int64_t walk_num);
|
|
qmckl_exit_code qmckl_set_electron_coord (qmckl_context context, const char transp, const double* coord, const int64_t size_max);
|
|
|
|
qmckl_exit_code qmckl_set_electron_rescale_factor_ee (qmckl_context context, const double kappa_ee);
|
|
qmckl_exit_code qmckl_set_electron_rescale_factor_en (qmckl_context context, const double kappa_en);
|
|
#+end_src
|
|
|
|
#+NAME:pre2
|
|
#+begin_src c :exports none
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
#+end_src
|
|
|
|
#+NAME:post2
|
|
#+begin_src c :exports none
|
|
ctx->electron.uninitialized &= ~mask;
|
|
ctx->electron.provided = (ctx->electron.uninitialized == 0);
|
|
|
|
if (ctx->electron.provided) {
|
|
if (ctx->electron.coord_new.data != NULL) {
|
|
const qmckl_exit_code rc = qmckl_matrix_free(context, &(ctx->electron.coord_new));
|
|
assert (rc == QMCKL_SUCCESS);
|
|
}
|
|
if (ctx->electron.coord_old.data != NULL) {
|
|
const qmckl_exit_code rc = qmckl_matrix_free(context, &(ctx->electron.coord_old));
|
|
assert (rc == QMCKL_SUCCESS);
|
|
}
|
|
|
|
const int64_t point_num = ctx->electron.walk_num * ctx->electron.num;
|
|
|
|
qmckl_matrix coord_new = qmckl_matrix_alloc(context, point_num, 3);
|
|
|
|
if (coord_new.data == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_set_electron_num",
|
|
NULL);
|
|
}
|
|
ctx->electron.coord_new = coord_new;
|
|
|
|
qmckl_matrix coord_old = qmckl_matrix_alloc(context, point_num, 3);
|
|
|
|
if (coord_old.data == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_set_electron_num",
|
|
NULL);
|
|
}
|
|
ctx->electron.coord_old = coord_old;
|
|
|
|
ctx->point.num = point_num;
|
|
ctx->point.coord = coord_new;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
#+end_src
|
|
|
|
To set the number of electrons, we give the number of up-spin and
|
|
down-spin electrons to the context and we set the number of walkers.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_set_electron_num(qmckl_context context,
|
|
const int64_t up_num,
|
|
const int64_t down_num) {
|
|
<<pre2>>
|
|
|
|
if (up_num <= 0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_electron_num",
|
|
"up_num <= 0");
|
|
}
|
|
|
|
if (down_num < 0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_set_electron_num",
|
|
"down_num < 0");
|
|
}
|
|
|
|
int32_t mask = 1 << 0;
|
|
|
|
ctx->electron.up_num = up_num;
|
|
ctx->electron.down_num = down_num;
|
|
ctx->electron.num = up_num + down_num;
|
|
|
|
<<post2>>
|
|
}
|
|
#+end_src
|
|
|
|
The following function sets the number of walkers.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_set_electron_walk_num(qmckl_context context, const int64_t walk_num) {
|
|
|
|
<<pre2>>
|
|
|
|
if (walk_num <= 0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_electron_walk_num",
|
|
"walk_num <= 0");
|
|
}
|
|
|
|
int32_t mask = 1 << 1;
|
|
ctx->electron.walk_num = walk_num;
|
|
|
|
<<post2>>
|
|
}
|
|
#+end_src
|
|
|
|
Next we set the rescale parameter for the rescaled distance metric.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_set_electron_rescale_factor_ee(qmckl_context context,
|
|
const double rescale_factor_kappa_ee) {
|
|
<<pre2>>
|
|
|
|
if (rescale_factor_kappa_ee <= 0.0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_electron_rescale_factor_ee",
|
|
"rescale_factor_kappa_ee <= 0.0");
|
|
}
|
|
|
|
ctx->electron.rescale_factor_kappa_ee = rescale_factor_kappa_ee;
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
qmckl_exit_code
|
|
qmckl_set_electron_rescale_factor_en(qmckl_context context,
|
|
const double rescale_factor_kappa_en) {
|
|
<<pre2>>
|
|
|
|
if (rescale_factor_kappa_en <= 0.0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_electron_rescale_factor_en",
|
|
"rescale_factor_kappa_en <= 0.0");
|
|
}
|
|
|
|
ctx->electron.rescale_factor_kappa_en = rescale_factor_kappa_en;
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_electron_num(context, alpha, beta) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: alpha
|
|
integer (c_int64_t) , intent(in) , value :: beta
|
|
end function
|
|
end interface
|
|
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_electron_walk_num(context, walk_num) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
end function
|
|
end interface
|
|
#+end_src
|
|
|
|
The following function sets the electron coordinates of all the
|
|
walkers. When this is done, the pointers to the old and new sets
|
|
of coordinates are swapped, and the new coordinates are
|
|
overwritten. This can be done only when the data relative to
|
|
electrons have been set. ~size_max~ should be equal to ~elec_num * walk_num * 3~, to be symmetric
|
|
with ~qmckl_get_electron_coord~.
|
|
|
|
Important: changing the electron coordinates increments the date
|
|
in the context.
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code
|
|
qmckl_set_electron_coord(qmckl_context context,
|
|
const char transp,
|
|
const double* coord,
|
|
const int64_t size_max)
|
|
{
|
|
|
|
<<pre2>>
|
|
|
|
if (transp != 'N' && transp != 'T') {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_electron_coord",
|
|
"transp should be 'N' or 'T'");
|
|
}
|
|
|
|
if (coord == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_3,
|
|
"qmckl_set_electron_coord",
|
|
"coord is a null pointer");
|
|
}
|
|
|
|
const int64_t elec_num = ctx->electron.num;
|
|
|
|
if (elec_num == 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_electron_coord",
|
|
"elec_num is not set");
|
|
}
|
|
|
|
const int64_t walk_num = ctx->electron.walk_num;
|
|
|
|
const int64_t point_num = ctx->point.num ;
|
|
|
|
if (walk_num == 0L) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_electron_coord",
|
|
"walk_num is not set");
|
|
}
|
|
|
|
assert(point_num == elec_num * walk_num);
|
|
|
|
/* Swap pointers */
|
|
ctx->point.coord = ctx->electron.coord_old;
|
|
ctx->electron.coord_old = ctx->electron.coord_new ;
|
|
|
|
qmckl_exit_code rc;
|
|
rc = qmckl_set_point(context, transp, coord, size_max/3);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
ctx->electron.coord_new = ctx->point.coord ;
|
|
|
|
ctx->electron.coord_new_date = ctx->date;
|
|
|
|
return QMCKL_SUCCESS;
|
|
|
|
}
|
|
#+end_src
|
|
|
|
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
|
|
interface
|
|
integer(c_int32_t) function qmckl_set_electron_coord(context, transp, coord, size_max) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
character , intent(in) , value :: transp
|
|
double precision , intent(in) :: coord(*)
|
|
integer (c_int64_t) , intent(in) , value :: size_max
|
|
end function
|
|
end interface
|
|
#+end_src
|
|
|
|
** Test
|
|
|
|
#+begin_src python :results output :exports none
|
|
import numpy as np
|
|
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
/* Reference input data */
|
|
int64_t walk_num = chbrclf_walk_num;
|
|
int64_t elec_num = chbrclf_elec_num;
|
|
int64_t elec_up_num = chbrclf_elec_up_num;
|
|
int64_t elec_dn_num = chbrclf_elec_dn_num;
|
|
double rescale_factor_kappa_ee = 1.0;
|
|
double rescale_factor_kappa_en = 1.0;
|
|
double nucl_rescale_factor_kappa = 1.0;
|
|
double* elec_coord = &(chbrclf_elec_coord[0][0][0]);
|
|
|
|
int64_t nucl_num = chbrclf_nucl_num;
|
|
double* charge = chbrclf_charge;
|
|
double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
|
|
|
|
/* --- */
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
assert(!qmckl_electron_provided(context));
|
|
|
|
int64_t n;
|
|
rc = qmckl_get_electron_num (context, &n);
|
|
assert(rc == QMCKL_NOT_PROVIDED);
|
|
|
|
rc = qmckl_get_electron_up_num (context, &n);
|
|
assert(rc == QMCKL_NOT_PROVIDED);
|
|
|
|
rc = qmckl_get_electron_down_num (context, &n);
|
|
assert(rc == QMCKL_NOT_PROVIDED);
|
|
|
|
|
|
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_electron_provided(context));
|
|
|
|
rc = qmckl_get_electron_up_num (context, &n);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(n == elec_up_num);
|
|
|
|
rc = qmckl_get_electron_down_num (context, &n);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(n == elec_dn_num);
|
|
|
|
rc = qmckl_get_electron_num (context, &n);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(n == elec_num);
|
|
|
|
double k_ee = 0.;
|
|
double k_en = 0.;
|
|
rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(k_ee == 1.0);
|
|
|
|
rc = qmckl_get_electron_rescale_factor_en (context, &k_en);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(k_en == 1.0);
|
|
|
|
rc = qmckl_set_electron_rescale_factor_en(context, rescale_factor_kappa_en);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_electron_rescale_factor_ee(context, rescale_factor_kappa_ee);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(k_ee == rescale_factor_kappa_ee);
|
|
|
|
rc = qmckl_get_electron_rescale_factor_en (context, &k_en);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(k_en == rescale_factor_kappa_en);
|
|
|
|
|
|
int64_t w;
|
|
rc = qmckl_get_electron_walk_num (context, &w);
|
|
assert(rc == QMCKL_NOT_PROVIDED);
|
|
|
|
|
|
rc = qmckl_set_electron_walk_num (context, walk_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_get_electron_walk_num (context, &w);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(w == walk_num);
|
|
|
|
assert(qmckl_electron_provided(context));
|
|
|
|
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
double elec_coord2[walk_num*3*elec_num];
|
|
|
|
rc = qmckl_get_electron_coord (context, 'N', elec_coord2, walk_num*3*elec_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
for (int64_t i=0 ; i<3*elec_num*walk_num ; ++i) {
|
|
printf("%f %f\n", elec_coord[i], elec_coord2[i]);
|
|
assert( elec_coord[i] == elec_coord2[i] );
|
|
}
|
|
|
|
#+end_src
|
|
|
|
* Computation
|
|
|
|
The computed data is stored in the context so that it can be reused
|
|
by different kernels. To ensure that the data is valid, for each
|
|
computed data the date of the context is stored when it is computed.
|
|
To know if some data needs to be recomputed, we check if the date of
|
|
the dependencies are more recent than the date of the data to
|
|
compute. If it is the case, then the data is recomputed and the
|
|
current date is stored.
|
|
|
|
** Electron-electron distances
|
|
|
|
*** Get
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_electron_ee_distance(qmckl_context context, double* const distance);
|
|
#+end_src
|
|
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_electron_ee_distance(context, distance) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double ) , intent(out) :: distance(*)
|
|
end function
|
|
end interface
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_electron_ee_distance(qmckl_context context, double* const distance)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_ee_distance(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num;
|
|
memcpy(distance, ctx->electron.ee_distance, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->electron.coord_new_date > ctx->electron.ee_distance_date) {
|
|
|
|
/* Allocate array */
|
|
if (ctx->electron.ee_distance == NULL) {
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = ctx->electron.num * ctx->electron.num *
|
|
ctx->electron.walk_num * sizeof(double);
|
|
double* ee_distance = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (ee_distance == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_ee_distance",
|
|
NULL);
|
|
}
|
|
ctx->electron.ee_distance = ee_distance;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_ee_distance(context,
|
|
ctx->electron.num,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.coord_new.data,
|
|
ctx->electron.ee_distance);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->electron.ee_distance_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_ee_distance
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_ee_distance_args
|
|
| Variable | Type | In/Out | Description |
|
|
|---------------+----------------------------------------+--------+-----------------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
|
|
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
|
|
| ~coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates |
|
|
| ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | out | Electron-electron distances |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_ee_distance_f(context, elec_num, walk_num, coord, ee_distance) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: elec_num
|
|
integer*8 , intent(in) :: walk_num
|
|
double precision , intent(in) :: coord(elec_num,walk_num,3)
|
|
double precision , intent(out) :: ee_distance(elec_num,elec_num,walk_num)
|
|
|
|
integer*8 :: k, i, j
|
|
double precision :: x, y, z
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
do k=1,walk_num
|
|
info = qmckl_distance(context, 'T', 'T', elec_num, elec_num, &
|
|
coord(1,k,1), elec_num * walk_num, &
|
|
coord(1,k,1), elec_num * walk_num, &
|
|
ee_distance(1,1,k), elec_num)
|
|
if (info /= QMCKL_SUCCESS) then
|
|
exit
|
|
endif
|
|
end do
|
|
|
|
end function qmckl_compute_ee_distance_f
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval h_private_func) :comments org :exports none
|
|
qmckl_exit_code qmckl_compute_ee_distance (
|
|
const qmckl_context context,
|
|
const int64_t elec_num,
|
|
const int64_t walk_num,
|
|
const double* coord,
|
|
double* const ee_distance );
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_ee_distance_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_ee_distance &
|
|
(context, elec_num, walk_num, coord, ee_distance) &
|
|
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 :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
real (c_double ) , intent(in) :: coord(elec_num,3,walk_num)
|
|
real (c_double ) , intent(out) :: ee_distance(elec_num,elec_num,walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_ee_distance_f
|
|
info = qmckl_compute_ee_distance_f &
|
|
(context, elec_num, walk_num, coord, ee_distance)
|
|
|
|
end function qmckl_compute_ee_distance
|
|
#+end_src
|
|
|
|
*** Test
|
|
|
|
#+begin_src python :results output :exports none
|
|
import numpy as np
|
|
|
|
elec_1_w1 = np.array( [ -2.26995253563, -5.15737533569, -2.22940072417 ])
|
|
elec_2_w1 = np.array( [ 3.51983380318, -1.08717381954, -1.19617708027 ])
|
|
elec_1_w2 = np.array( [ -2.34410619736, -3.20016115904, -1.53496759012 ])
|
|
elec_2_w2 = np.array( [ 3.17996025085, -1.40260577202, 1.49473607540 ])
|
|
|
|
print ( "[0][0][0] : ", np.linalg.norm(elec_1_w1-elec_1_w1) )
|
|
print ( "[0][1][0] : ", np.linalg.norm(elec_1_w1-elec_2_w1) )
|
|
print ( "[1][0][0] : ", np.linalg.norm(elec_2_w1-elec_1_w1) )
|
|
print ( "[0][0][1] : ", np.linalg.norm(elec_1_w2-elec_1_w2) )
|
|
print ( "[0][1][1] : ", np.linalg.norm(elec_1_w2-elec_2_w2) )
|
|
print ( "[1][0][1] : ", np.linalg.norm(elec_2_w2-elec_1_w2) )
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: [0][0][0] : 0.0
|
|
: [0][1][0] : 7.152322512964209
|
|
: [1][0][0] : 7.152322512964209
|
|
: [0][0][1] : 0.0
|
|
: [0][1][1] : 6.5517646321055665
|
|
: [1][0][1] : 6.5517646321055665
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
assert(qmckl_electron_provided(context));
|
|
|
|
|
|
double ee_distance[walk_num * elec_num * elec_num];
|
|
rc = qmckl_get_electron_ee_distance(context, ee_distance);
|
|
|
|
// (e1,e2,w)
|
|
// (0,0,0) == 0.
|
|
assert(ee_distance[0] == 0.);
|
|
|
|
// (1,0,0) == (0,1,0)
|
|
assert(ee_distance[1] == ee_distance[elec_num]);
|
|
|
|
// value of (1,0,0)
|
|
assert(fabs(ee_distance[1]-7.152322512964209) < 1.e-12);
|
|
|
|
// (0,0,1) == 0.
|
|
assert(ee_distance[elec_num*elec_num] == 0.);
|
|
|
|
// (1,0,1) == (0,1,1)
|
|
assert(ee_distance[elec_num*elec_num+1] == ee_distance[elec_num*elec_num+elec_num]);
|
|
|
|
// value of (1,0,1)
|
|
assert(fabs(ee_distance[elec_num*elec_num+1]-6.5517646321055665) < 1.e-12);
|
|
|
|
#+end_src
|
|
|
|
** Electron-electron rescaled distances ~ee_distance_rescaled~ stores the matrix of the rescaled distances between all
|
|
pairs of electrons:
|
|
|
|
\[
|
|
C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa
|
|
\]
|
|
|
|
where \(C_{ij}\) is the matrix of electron-electron distances.
|
|
|
|
*** Get
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_electron_ee_distance_rescaled(qmckl_context context, double* const distance_rescaled);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_electron_ee_distance_rescaled(qmckl_context context, double* const distance_rescaled)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_ee_distance_rescaled(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num;
|
|
memcpy(distance_rescaled, ctx->electron.ee_distance_rescaled, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->electron.coord_new_date > ctx->electron.ee_distance_rescaled_date) {
|
|
|
|
/* Allocate array */
|
|
if (ctx->electron.ee_distance_rescaled == NULL) {
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = ctx->electron.num * ctx->electron.num *
|
|
ctx->electron.walk_num * sizeof(double);
|
|
double* ee_distance_rescaled = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (ee_distance_rescaled == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_ee_distance_rescaled",
|
|
NULL);
|
|
}
|
|
ctx->electron.ee_distance_rescaled = ee_distance_rescaled;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_ee_distance_rescaled(context,
|
|
ctx->electron.num,
|
|
ctx->electron.rescale_factor_kappa_en,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.coord_new.data,
|
|
ctx->electron.ee_distance_rescaled);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->electron.ee_distance_rescaled_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_ee_distance_rescaled
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_ee_distance_rescaled_args
|
|
| Variable | Type | In/Out | Description |
|
|
|---------------------------+----------------------------------------+--------+--------------------------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
|
|
| ~rescale_factor_kappa_ee~ | ~double~ | in | Factor to rescale ee distances |
|
|
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
|
|
| ~coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates |
|
|
| ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | out | Electron-electron rescaled distances |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_ee_distance_rescaled_f(context, elec_num, rescale_factor_kappa_ee, walk_num, &
|
|
coord, ee_distance_rescaled) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: elec_num
|
|
double precision , intent(in) :: rescale_factor_kappa_ee
|
|
integer*8 , intent(in) :: walk_num
|
|
double precision , intent(in) :: coord(elec_num,3,walk_num)
|
|
double precision , intent(out) :: ee_distance_rescaled(elec_num,elec_num,walk_num)
|
|
|
|
integer*8 :: k
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
do k=1,walk_num
|
|
info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, elec_num, &
|
|
coord(1,k,1), elec_num * walk_num, &
|
|
coord(1,k,1), elec_num * walk_num, &
|
|
ee_distance_rescaled(1,1,k), elec_num, rescale_factor_kappa_ee)
|
|
if (info /= QMCKL_SUCCESS) then
|
|
exit
|
|
endif
|
|
end do
|
|
|
|
end function qmckl_compute_ee_distance_rescaled_f
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval h_private_func) :comments org :exports none
|
|
qmckl_exit_code qmckl_compute_ee_distance_rescaled (
|
|
const qmckl_context context,
|
|
const int64_t elec_num,
|
|
const double rescale_factor_kappa_ee,
|
|
const int64_t walk_num,
|
|
const double* coord,
|
|
double* const ee_distance_rescaled );
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_ee_distance_rescaled_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_ee_distance_rescaled &
|
|
(context, elec_num, rescale_factor_kappa_ee, walk_num, coord, ee_distance_rescaled) &
|
|
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 :: elec_num
|
|
real (c_double ) , intent(in) , value :: rescale_factor_kappa_ee
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
real (c_double ) , intent(in) :: coord(elec_num,3,walk_num)
|
|
real (c_double ) , intent(out) :: ee_distance_rescaled(elec_num,elec_num,walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_ee_distance_rescaled_f
|
|
info = qmckl_compute_ee_distance_rescaled_f &
|
|
(context, elec_num, rescale_factor_kappa_ee, walk_num, coord, ee_distance_rescaled)
|
|
|
|
end function qmckl_compute_ee_distance_rescaled
|
|
#+end_src
|
|
|
|
*** Test
|
|
|
|
#+begin_src python :results output :exports none
|
|
import numpy as np
|
|
|
|
kappa = 1.0
|
|
|
|
elec_1_w1 = np.array( [ -2.26995253563, -5.15737533569, -2.22940072417 ])
|
|
elec_2_w1 = np.array( [ 3.51983380318, -1.08717381954, -1.19617708027 ])
|
|
elec_1_w2 = np.array( [ -2.34410619736, -3.20016115904, -1.53496759012 ])
|
|
elec_2_w2 = np.array( [ 3.17996025085, -1.40260577202, 1.49473607540 ])
|
|
|
|
print ( "[0][0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-elec_1_w1)) )/kappa )
|
|
print ( "[0][1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-elec_2_w1)) )/kappa )
|
|
print ( "[1][0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w1-elec_1_w1)) )/kappa )
|
|
print ( "[0][0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w2-elec_1_w2)) )/kappa )
|
|
print ( "[0][1][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w2-elec_2_w2)) )/kappa )
|
|
print ( "[1][0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w2-elec_1_w2)) )/kappa )
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: [0][0][0] : 0.0
|
|
: [0][1][0] : 0.9992169566605263
|
|
: [1][0][0] : 0.9992169566605263
|
|
: [0][0][1] : 0.0
|
|
: [0][1][1] : 0.9985724058042633
|
|
: [1][0][1] : 0.9985724058042633
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
assert(qmckl_electron_provided(context));
|
|
|
|
|
|
double ee_distance_rescaled[walk_num * elec_num * elec_num];
|
|
rc = qmckl_get_electron_ee_distance_rescaled(context, ee_distance_rescaled);
|
|
|
|
// (e1,e2,w)
|
|
// (0,0,0) == 0.
|
|
assert(ee_distance_rescaled[0] == 0.);
|
|
|
|
// (1,0,0) == (0,1,0)
|
|
assert(ee_distance_rescaled[1] == ee_distance_rescaled[elec_num]);
|
|
|
|
// value of (1,0,0)
|
|
assert(fabs(ee_distance_rescaled[1]-0.9992169566605263) < 1.e-12);
|
|
|
|
// (0,0,1) == 0.
|
|
assert(ee_distance_rescaled[elec_num*elec_num] == 0.);
|
|
|
|
// (1,0,1) == (0,1,1)
|
|
assert(ee_distance_rescaled[elec_num*elec_num+1] == ee_distance_rescaled[elec_num*elec_num+elec_num]);
|
|
|
|
// value of (1,0,1)
|
|
assert(fabs(ee_distance_rescaled[elec_num*elec_num+1]-0.9985724058042633) < 1.e-12);
|
|
|
|
#+end_src
|
|
|
|
** Electron-electron rescaled distance gradients and laplacian with respect to electron coords
|
|
|
|
The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$
|
|
needs to be perturbed with respect to the electorn coordinates.
|
|
This data is stored in the ~ee_distance_rescaled_deriv_e~ tensor. The
|
|
The first three elements of this three index tensor ~[4][num][num]~ gives the
|
|
derivatives in the x, y, and z directions $dx, dy, dz$ and the last index
|
|
gives the Laplacian $\partial x^2 + \partial y^2 + \partial z^2$.
|
|
|
|
*** Get
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_electron_ee_distance_rescaled_deriv_e(qmckl_context context, double* const distance_rescaled_deriv_e);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_electron_ee_distance_rescaled_deriv_e(qmckl_context context, double* const distance_rescaled_deriv_e)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_ee_distance_rescaled_deriv_e(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = 4 * ctx->electron.num * ctx->electron.num * ctx->electron.walk_num;
|
|
memcpy(distance_rescaled_deriv_e, ctx->electron.ee_distance_rescaled_deriv_e, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->electron.coord_new_date > ctx->electron.ee_distance_rescaled_deriv_e_date) {
|
|
|
|
/* Allocate array */
|
|
if (ctx->electron.ee_distance_rescaled_deriv_e == NULL) {
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = 4 * ctx->electron.num * ctx->electron.num *
|
|
ctx->electron.walk_num * sizeof(double);
|
|
double* ee_distance_rescaled_deriv_e = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (ee_distance_rescaled_deriv_e == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_ee_distance_rescaled_deriv_e",
|
|
NULL);
|
|
}
|
|
ctx->electron.ee_distance_rescaled_deriv_e = ee_distance_rescaled_deriv_e;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_ee_distance_rescaled_deriv_e(context,
|
|
ctx->electron.num,
|
|
ctx->electron.rescale_factor_kappa_en,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.coord_new.data,
|
|
ctx->electron.ee_distance_rescaled_deriv_e);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->electron.ee_distance_rescaled_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_ee_distance_rescaled_deriv_e
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_ee_distance_rescaled_deriv_e_args
|
|
| Variable | Type | In/Out | Description |
|
|
|---------------------------+-------------------------------------------+--------+-------------------------------------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
|
|
| ~rescale_factor_kappa_ee~ | ~double~ | in | Factor to rescale ee distances |
|
|
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
|
|
| ~coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates |
|
|
| ~ee_distance_deriv_e~ | ~double[walk_num][4][elec_num][elec_num]~ | out | Electron-electron rescaled distance derivatives |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_ee_distance_rescaled_deriv_e_f(context, elec_num, rescale_factor_kappa_ee, walk_num, &
|
|
coord, ee_distance_rescaled_deriv_e) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: elec_num
|
|
double precision , intent(in) :: rescale_factor_kappa_ee
|
|
integer*8 , intent(in) :: walk_num
|
|
double precision , intent(in) :: coord(elec_num,3,walk_num)
|
|
double precision , intent(out) :: ee_distance_rescaled_deriv_e(4,elec_num,elec_num,walk_num)
|
|
|
|
integer*8 :: k
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
do k=1,walk_num
|
|
info = qmckl_distance_rescaled_deriv_e(context, 'T', 'T', elec_num, elec_num, &
|
|
coord(1,1,k), elec_num, &
|
|
coord(1,1,k), elec_num, &
|
|
ee_distance_rescaled_deriv_e(1,1,1,k), elec_num, rescale_factor_kappa_ee)
|
|
if (info /= QMCKL_SUCCESS) then
|
|
exit
|
|
endif
|
|
end do
|
|
|
|
end function qmckl_compute_ee_distance_rescaled_deriv_e_f
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval h_private_func) :comments org :exports none
|
|
qmckl_exit_code qmckl_compute_ee_distance_rescaled_deriv_e (
|
|
const qmckl_context context,
|
|
const int64_t elec_num,
|
|
const double rescale_factor_kappa_ee,
|
|
const int64_t walk_num,
|
|
const double* coord,
|
|
double* const ee_distance_rescaled_deriv_e );
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_ee_distance_rescaled_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_ee_distance_rescaled_deriv_e &
|
|
(context, elec_num, rescale_factor_kappa_ee, walk_num, coord, ee_distance_rescaled_deriv_e) &
|
|
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 :: elec_num
|
|
real (c_double ) , intent(in) , value :: rescale_factor_kappa_ee
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
real (c_double ) , intent(in) :: coord(elec_num,3,walk_num)
|
|
real (c_double ) , intent(out) :: ee_distance_rescaled_deriv_e(4,elec_num,elec_num,walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_ee_distance_rescaled_deriv_e_f
|
|
info = qmckl_compute_ee_distance_rescaled_deriv_e_f &
|
|
(context, elec_num, rescale_factor_kappa_ee, walk_num, coord, ee_distance_rescaled_deriv_e)
|
|
|
|
end function qmckl_compute_ee_distance_rescaled_deriv_e
|
|
#+end_src
|
|
|
|
*** Test
|
|
|
|
#+begin_src python :results output :exports none
|
|
import numpy as np
|
|
|
|
# TODO
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
assert(qmckl_electron_provided(context));
|
|
|
|
|
|
double ee_distance_rescaled_deriv_e[4 * walk_num * elec_num * elec_num];
|
|
rc = qmckl_get_electron_ee_distance_rescaled_deriv_e(context, ee_distance_rescaled_deriv_e);
|
|
|
|
// TODO: Get exact values
|
|
//// (e1,e2,w)
|
|
//// (0,0,0) == 0.
|
|
//assert(ee_distance[0] == 0.);
|
|
//
|
|
//// (1,0,0) == (0,1,0)
|
|
//assert(ee_distance[1] == ee_distance[elec_num]);
|
|
//
|
|
//// value of (1,0,0)
|
|
//assert(fabs(ee_distance[1]-7.152322512964209) < 1.e-12);
|
|
//
|
|
//// (0,0,1) == 0.
|
|
//assert(ee_distance[elec_num*elec_num] == 0.);
|
|
//
|
|
//// (1,0,1) == (0,1,1)
|
|
//assert(ee_distance[elec_num*elec_num+1] == ee_distance[elec_num*elec_num+elec_num]);
|
|
//
|
|
//// value of (1,0,1)
|
|
//assert(fabs(ee_distance[elec_num*elec_num+1]-6.5517646321055665) < 1.e-12);
|
|
|
|
#+end_src
|
|
|
|
** Electron-electron potential ~ee_pot~ calculates the ~ee~ potential energy.
|
|
|
|
\[
|
|
\mathcal{V}_{ee} = \sum_{i=1}^{N_e}\sum_{j>i}^{N_e}\frac{1}{r_{ij}}
|
|
\]
|
|
|
|
where \(\mathcal{V}_{ee}\) is the ~ee~ potential and \[r_{ij}\] the ~ee~
|
|
distance.
|
|
|
|
*** Get
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* const ee_pot);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* const ee_pot)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_ee_potential(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = ctx->electron.walk_num * sizeof(double);
|
|
memcpy(ee_pot, ctx->electron.ee_pot, sze);
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
if (!ctx->electron.provided) return QMCKL_NOT_PROVIDED;
|
|
|
|
qmckl_exit_code rc = qmckl_provide_ee_distance(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->electron.coord_new_date > ctx->electron.ee_pot_date) {
|
|
|
|
/* Allocate array */
|
|
if (ctx->electron.ee_pot == NULL) {
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = ctx->electron.walk_num * sizeof(double);
|
|
double* ee_pot = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (ee_pot == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_ee_potential",
|
|
NULL);
|
|
}
|
|
ctx->electron.ee_pot = ee_pot;
|
|
}
|
|
|
|
rc = qmckl_compute_ee_potential(context,
|
|
ctx->electron.num,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.ee_distance,
|
|
ctx->electron.ee_pot);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->electron.ee_pot_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_ee_potential
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_ee_potential_args
|
|
| Variable | Type | In/Out | Description |
|
|
|---------------+----------------------------------------+--------+--------------------------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
|
|
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
|
|
| ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | in | Electron-electron rescaled distances |
|
|
| ~ee_pot~ | ~double[walk_num]~ | out | Electron-electron potential |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_ee_potential_f(context, elec_num, walk_num, &
|
|
ee_distance, ee_pot) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: elec_num
|
|
integer*8 , intent(in) :: walk_num
|
|
double precision , intent(in) :: ee_distance(elec_num,elec_num,walk_num)
|
|
double precision , intent(out) :: ee_pot(walk_num)
|
|
|
|
integer*8 :: nw, i, j
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
ee_pot = 0.0d0
|
|
do nw=1,walk_num
|
|
do j=2,elec_num
|
|
do i=1,j-1
|
|
if (dabs(ee_distance(i,j,nw)) > 1e-5) then
|
|
ee_pot(nw) = ee_pot(nw) + 1.0d0/(ee_distance(i,j,nw))
|
|
endif
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end function qmckl_compute_ee_potential_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_header(table=qmckl_ee_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_compute_ee_potential (
|
|
const qmckl_context context,
|
|
const int64_t elec_num,
|
|
const int64_t walk_num,
|
|
const double* ee_distance,
|
|
double* const ee_pot );
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_ee_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_ee_potential &
|
|
(context, elec_num, walk_num, ee_distance, ee_pot) &
|
|
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 :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
real (c_double ) , intent(in) :: ee_distance(elec_num,elec_num,walk_num)
|
|
real (c_double ) , intent(out) :: ee_pot(walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_ee_potential_f
|
|
info = qmckl_compute_ee_potential_f &
|
|
(context, elec_num, walk_num, ee_distance, ee_pot)
|
|
|
|
end function qmckl_compute_ee_potential
|
|
#+end_src
|
|
|
|
*** Test
|
|
#+begin_src c :tangle (eval c_test)
|
|
double ee_pot[walk_num];
|
|
|
|
rc = qmckl_get_electron_ee_potential(context, &(ee_pot[0]));
|
|
assert (rc == QMCKL_SUCCESS);
|
|
#+end_src
|
|
** Electron-nucleus distances
|
|
|
|
*** Get
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_electron_en_distance(qmckl_context context, double* distance);
|
|
#+end_src
|
|
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(c_int32_t) function qmckl_get_electron_en_distance(context, distance) &
|
|
bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
implicit none
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
real (c_double ) , intent(out) :: distance(*)
|
|
end function
|
|
end interface
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_electron_en_distance(qmckl_context context, double* distance)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_en_distance(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num;
|
|
memcpy(distance, ctx->electron.en_distance, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_en_distance(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
if (!(ctx->nucleus.provided)) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_provide_en_distance",
|
|
NULL);
|
|
}
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->electron.coord_new_date > ctx->electron.en_distance_date) {
|
|
|
|
/* Allocate array */
|
|
if (ctx->electron.en_distance == NULL) {
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = ctx->electron.num * ctx->nucleus.num *
|
|
ctx->electron.walk_num * sizeof(double);
|
|
double* en_distance = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (en_distance == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_en_distance",
|
|
NULL);
|
|
}
|
|
ctx->electron.en_distance = en_distance;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_en_distance(context,
|
|
ctx->electron.num,
|
|
ctx->nucleus.num,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.coord_new.data,
|
|
ctx->nucleus.coord.data,
|
|
ctx->electron.en_distance);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->electron.en_distance_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_en_distance
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_en_distance_args
|
|
| Variable | Type | In/Out | Description |
|
|
|---------------+----------------------------------------+--------+----------------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
|
|
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
|
|
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
|
|
| ~elec_coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates |
|
|
| ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates |
|
|
| ~en_distance~ | ~double[walk_num][nucl_num][elec_num]~ | out | Electron-nucleus distances |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_num, elec_coord, nucl_coord, en_distance) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: elec_num
|
|
integer*8 , intent(in) :: nucl_num
|
|
integer*8 , intent(in) :: walk_num
|
|
double precision , intent(in) :: elec_coord(elec_num,walk_num,3)
|
|
double precision , intent(in) :: nucl_coord(nucl_num,3)
|
|
double precision , intent(out) :: en_distance(elec_num,nucl_num,walk_num)
|
|
|
|
integer*8 :: k
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (nucl_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_4
|
|
return
|
|
endif
|
|
|
|
do k=1,walk_num
|
|
info = qmckl_distance(context, 'T', 'T', elec_num, nucl_num, &
|
|
elec_coord(1,k,1), elec_num * walk_num, &
|
|
nucl_coord, nucl_num, &
|
|
en_distance(1,1,k), elec_num)
|
|
if (info /= QMCKL_SUCCESS) then
|
|
exit
|
|
endif
|
|
end do
|
|
|
|
end function qmckl_compute_en_distance_f
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval h_private_func) :comments org :exports none
|
|
qmckl_exit_code qmckl_compute_en_distance (
|
|
const qmckl_context context,
|
|
const int64_t elec_num,
|
|
const int64_t nucl_num,
|
|
const int64_t walk_num,
|
|
const double* elec_coord,
|
|
const double* nucl_coord,
|
|
double* const en_distance );
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_en_distance_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_en_distance &
|
|
(context, elec_num, nucl_num, walk_num, elec_coord, nucl_coord, en_distance) &
|
|
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 :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3)
|
|
real (c_double ) , intent(in) :: nucl_coord(elec_num,3)
|
|
real (c_double ) , intent(out) :: en_distance(elec_num,nucl_num,walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_en_distance_f
|
|
info = qmckl_compute_en_distance_f &
|
|
(context, elec_num, nucl_num, walk_num, elec_coord, nucl_coord, en_distance)
|
|
|
|
end function qmckl_compute_en_distance
|
|
#+end_src
|
|
|
|
*** Test
|
|
|
|
#+begin_src python :results output :exports none
|
|
import numpy as np
|
|
|
|
elec_1_w1 = np.array( [ -2.26995253563, -5.15737533569, -2.22940072417 ])
|
|
elec_2_w1 = np.array( [ 3.51983380318, -1.08717381954, -1.19617708027 ])
|
|
elec_1_w2 = np.array( [ -2.34410619736, -3.20016115904, -1.53496759012 ])
|
|
elec_2_w2 = np.array( [ 3.17996025085, -1.40260577202, 1.49473607540 ])
|
|
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 ] )
|
|
|
|
print ( "[0][0][0] : ", np.linalg.norm(elec_1_w1-nucl_1) )
|
|
print ( "[0][1][0] : ", np.linalg.norm(elec_1_w1-nucl_2) )
|
|
print ( "[0][0][1] : ", np.linalg.norm(elec_2_w1-nucl_1) )
|
|
print ( "[1][0][0] : ", np.linalg.norm(elec_1_w2-nucl_1) )
|
|
print ( "[1][1][0] : ", np.linalg.norm(elec_1_w2-nucl_2) )
|
|
print ( "[1][0][1] : ", np.linalg.norm(elec_2_w2-nucl_1) )
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: [0][0][0] : 7.546738741619978
|
|
: [0][1][0] : 8.77102435246984
|
|
: [0][0][1] : 3.698922010513608
|
|
: [1][0][0] : 5.824059436060509
|
|
: [1][1][0] : 7.080482110317645
|
|
: [1][0][1] : 3.1804527583077356
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
|
|
assert(!qmckl_nucleus_provided(context));
|
|
assert(qmckl_electron_provided(context));
|
|
|
|
rc = qmckl_set_nucleus_num (context, nucl_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_nucleus_charge (context, charge, nucl_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord, 3*nucl_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
assert(qmckl_nucleus_provided(context));
|
|
|
|
double en_distance[walk_num][nucl_num][elec_num];
|
|
|
|
rc = qmckl_get_electron_en_distance(context, &(en_distance[0][0][0]));
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
// (e,n,w) in Fortran notation
|
|
// (1,1,1)
|
|
assert(fabs(en_distance[0][0][0] - 7.546738741619978) < 1.e-12);
|
|
|
|
// (1,2,1)
|
|
assert(fabs(en_distance[0][1][0] - 8.77102435246984) < 1.e-12);
|
|
|
|
// (2,1,1)
|
|
assert(fabs(en_distance[0][0][1] - 3.698922010513608) < 1.e-12);
|
|
|
|
// (1,1,2)
|
|
assert(fabs(en_distance[1][0][0] - 5.824059436060509) < 1.e-12);
|
|
|
|
// (1,2,2)
|
|
assert(fabs(en_distance[1][1][0] - 7.080482110317645) < 1.e-12);
|
|
|
|
// (2,1,2)
|
|
assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12);
|
|
|
|
#+end_src
|
|
|
|
** Electron-nucleus rescaled distances ~en_distance_rescaled~ stores the matrix of the rescaled distances between
|
|
electrons and nuclei.
|
|
|
|
\[
|
|
C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa
|
|
\]
|
|
|
|
where \(C_{ij}\) is the matrix of electron-nucleus distances.
|
|
|
|
*** Get
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled);
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_en_distance_rescaled(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num;
|
|
memcpy(distance_rescaled, ctx->electron.en_distance_rescaled, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
if (!(ctx->nucleus.provided)) {
|
|
return QMCKL_NOT_PROVIDED;
|
|
}
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->electron.coord_new_date > ctx->electron.en_distance_rescaled_date) {
|
|
|
|
/* Allocate array */
|
|
if (ctx->electron.en_distance_rescaled == NULL) {
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = ctx->electron.num * ctx->nucleus.num *
|
|
ctx->electron.walk_num * sizeof(double);
|
|
double* en_distance_rescaled = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (en_distance_rescaled == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_en_distance_rescaled",
|
|
NULL);
|
|
}
|
|
ctx->electron.en_distance_rescaled = en_distance_rescaled;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_en_distance_rescaled(context,
|
|
ctx->electron.num,
|
|
ctx->nucleus.num,
|
|
ctx->electron.rescale_factor_kappa_en,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.coord_new.data,
|
|
ctx->nucleus.coord.data,
|
|
ctx->electron.en_distance_rescaled);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->electron.en_distance_rescaled_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_en_distance_rescaled
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_en_distance_rescaled_args
|
|
| Variable | Type | In/Out | Description |
|
|
|---------------------------+----------------------------------------+--------+-----------------------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
|
|
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
|
|
| ~rescale_factor_kappa_en~ | ~double~ | in | The factor for rescaled distances |
|
|
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
|
|
| ~elec_coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates |
|
|
| ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates |
|
|
| ~en_distance_rescaled~ | ~double[walk_num][nucl_num][elec_num]~ | out | Electron-nucleus distances |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_en_distance_rescaled_f(context, elec_num, nucl_num, rescale_factor_kappa_en, walk_num, elec_coord, &
|
|
nucl_coord, en_distance_rescaled) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: elec_num
|
|
integer*8 , intent(in) :: nucl_num
|
|
double precision , intent(in) :: rescale_factor_kappa_en
|
|
integer*8 , intent(in) :: walk_num
|
|
double precision , intent(in) :: elec_coord(elec_num,walk_num,3)
|
|
double precision , intent(in) :: nucl_coord(nucl_num,3)
|
|
double precision , intent(out) :: en_distance_rescaled(elec_num,nucl_num,walk_num)
|
|
|
|
integer*8 :: k
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (nucl_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
! TODO: comparison with 0
|
|
!if (rescale_factor_kappa_en <= 0) then
|
|
! info = QMCKL_INVALID_ARG_4
|
|
! return
|
|
!endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_5
|
|
return
|
|
endif
|
|
|
|
do k=1,walk_num
|
|
info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, nucl_num, &
|
|
elec_coord(1,k,1), elec_num*walk_num, &
|
|
nucl_coord, nucl_num, &
|
|
en_distance_rescaled(1,1,k), elec_num, rescale_factor_kappa_en)
|
|
if (info /= QMCKL_SUCCESS) then
|
|
exit
|
|
endif
|
|
end do
|
|
|
|
end function qmckl_compute_en_distance_rescaled_f
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval h_private_func) :comments org :exports none
|
|
qmckl_exit_code qmckl_compute_en_distance_rescaled (
|
|
const qmckl_context context,
|
|
const int64_t elec_num,
|
|
const int64_t nucl_num,
|
|
const double rescale_factor_kappa_en,
|
|
const int64_t walk_num,
|
|
const double* elec_coord,
|
|
const double* nucl_coord,
|
|
double* const en_distance_rescaled );
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_en_distance_rescaled_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_en_distance_rescaled &
|
|
(context, elec_num, nucl_num, rescale_factor_kappa_en, walk_num, elec_coord, nucl_coord, en_distance_rescaled) &
|
|
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 :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
real (c_double ) , intent(in) , value :: rescale_factor_kappa_en
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3)
|
|
real (c_double ) , intent(in) :: nucl_coord(elec_num,3)
|
|
real (c_double ) , intent(out) :: en_distance_rescaled(elec_num,nucl_num,walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_en_distance_rescaled_f
|
|
info = qmckl_compute_en_distance_rescaled_f &
|
|
(context, elec_num, nucl_num, rescale_factor_kappa_en, walk_num, elec_coord, nucl_coord, en_distance_rescaled)
|
|
|
|
end function qmckl_compute_en_distance_rescaled
|
|
#+end_src
|
|
|
|
*** Test
|
|
|
|
#+begin_src python :results output :exports none
|
|
import numpy as np
|
|
|
|
kappa = 1.0
|
|
|
|
elec_1_w1 = np.array( [ -2.26995253563, -5.15737533569, -2.22940072417 ])
|
|
elec_2_w1 = np.array( [ 3.51983380318, -1.08717381954, -1.19617708027 ])
|
|
elec_1_w2 = np.array( [ -2.34410619736, -3.20016115904, -1.53496759012 ])
|
|
elec_2_w2 = np.array( [ 3.17996025085, -1.40260577202, 1.49473607540 ])
|
|
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 ] )
|
|
|
|
print ( "[0][0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-nucl_1)) )/kappa )
|
|
print ( "[0][1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-nucl_2)) )/kappa )
|
|
print ( "[0][0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w1-nucl_1)) )/kappa )
|
|
print ( "[1][0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w2-nucl_1)) )/kappa )
|
|
print ( "[1][1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w2-nucl_2)) )/kappa )
|
|
print ( "[1][0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w2-nucl_1)) )/kappa )
|
|
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: [0][0][0] : 0.9994721712909764
|
|
: [0][1][0] : 0.9998448354439821
|
|
: [0][0][1] : 0.9752498074577688
|
|
: [1][0][0] : 0.9970444172399963
|
|
: [1][1][0] : 0.9991586325813303
|
|
: [1][0][1] : 0.9584331688679852
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
|
|
assert(qmckl_electron_provided(context));
|
|
|
|
rc = qmckl_set_nucleus_num (context, nucl_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_nucleus_charge (context, charge, nucl_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord, 3*nucl_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
assert(qmckl_nucleus_provided(context));
|
|
|
|
double en_distance_rescaled[walk_num][nucl_num][elec_num];
|
|
|
|
rc = qmckl_get_electron_en_distance_rescaled(context, &(en_distance_rescaled[0][0][0]));
|
|
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
// (e,n,w) in Fortran notation
|
|
// (1,1,1)
|
|
assert(fabs(en_distance_rescaled[0][0][0] - 0.9994721712909764) < 1.e-12);
|
|
|
|
// (1,2,1)
|
|
assert(fabs(en_distance_rescaled[0][1][0] - 0.9998448354439821) < 1.e-12);
|
|
|
|
// (2,1,1)
|
|
assert(fabs(en_distance_rescaled[0][0][1] - 0.9752498074577688) < 1.e-12);
|
|
|
|
// (1,1,2)
|
|
assert(fabs(en_distance_rescaled[1][0][0] - 0.9970444172399963) < 1.e-12);
|
|
|
|
// (1,2,2)
|
|
assert(fabs(en_distance_rescaled[1][1][0] - 0.9991586325813303) < 1.e-12);
|
|
|
|
// (2,1,2)
|
|
assert(fabs(en_distance_rescaled[1][0][1] - 0.9584331688679852) < 1.e-12);
|
|
|
|
#+end_src
|
|
|
|
** Electron-nucleus rescaled distance gradients and laplacian with respect to electron coords
|
|
|
|
The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$
|
|
needs to be perturbed with respect to the nuclear coordinates.
|
|
This data is stored in the ~en_distance_rescaled_deriv_e~ tensor. The
|
|
The first three elements of this three index tensor ~[4][nucl_num][elec_num]~ gives the
|
|
derivatives in the x, y, and z directions $dx, dy, dz$ and the last index
|
|
gives the Laplacian $\partial x^2 + \partial y^2 + \partial z^2$.
|
|
|
|
*** Get
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_electron_en_distance_rescaled_deriv_e(qmckl_context context, double* distance_rescaled_deriv_e);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_electron_en_distance_rescaled_deriv_e(qmckl_context context, double* distance_rescaled_deriv_e)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_en_distance_rescaled_deriv_e(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num;
|
|
memcpy(distance_rescaled_deriv_e, ctx->electron.en_distance_rescaled_deriv_e, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
if (!(ctx->nucleus.provided)) {
|
|
return QMCKL_NOT_PROVIDED;
|
|
}
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->electron.coord_new_date > ctx->electron.en_distance_rescaled_deriv_e_date) {
|
|
|
|
/* Allocate array */
|
|
if (ctx->electron.en_distance_rescaled_deriv_e == NULL) {
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = 4 * ctx->electron.num * ctx->nucleus.num *
|
|
ctx->electron.walk_num * sizeof(double);
|
|
double* en_distance_rescaled_deriv_e = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (en_distance_rescaled_deriv_e == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_en_distance_rescaled_deriv_e",
|
|
NULL);
|
|
}
|
|
ctx->electron.en_distance_rescaled_deriv_e = en_distance_rescaled_deriv_e;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_en_distance_rescaled_deriv_e(context,
|
|
ctx->electron.num,
|
|
ctx->nucleus.num,
|
|
ctx->electron.rescale_factor_kappa_en,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.coord_new.data,
|
|
ctx->nucleus.coord.data,
|
|
ctx->electron.en_distance_rescaled_deriv_e);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->electron.en_distance_rescaled_deriv_e_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_en_distance_rescaled_deriv_e
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_en_distance_rescaled_deriv_e_args
|
|
| Variable | Type | In/Out | Description |
|
|
|--------------------------------+-------------------------------------------+--------+---------------------------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
|
|
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
|
|
| ~rescale_factor_kappa_en~ | ~double~ | in | The factor for rescaled distances |
|
|
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
|
|
| ~elec_coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates |
|
|
| ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates |
|
|
| ~en_distance_rescaled_deriv_e~ | ~double[walk_num][4][nucl_num][elec_num]~ | out | Electron-nucleus distance derivatives |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_en_distance_rescaled_deriv_e_f(context, elec_num, nucl_num, &
|
|
rescale_factor_kappa_en, walk_num, elec_coord, &
|
|
nucl_coord, en_distance_rescaled_deriv_e) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: elec_num
|
|
integer*8 , intent(in) :: nucl_num
|
|
double precision , intent(in) :: rescale_factor_kappa_en
|
|
integer*8 , intent(in) :: walk_num
|
|
double precision , intent(in) :: elec_coord(elec_num,walk_num,3)
|
|
double precision , intent(in) :: nucl_coord(nucl_num,3)
|
|
double precision , intent(out) :: en_distance_rescaled_deriv_e(4,elec_num,nucl_num,walk_num)
|
|
|
|
integer*8 :: k
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (nucl_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
! TODO: comparison with 0
|
|
!if (rescale_factor_kappa_en <= 0) then
|
|
! info = QMCKL_INVALID_ARG_4
|
|
! return
|
|
!endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_5
|
|
return
|
|
endif
|
|
|
|
do k=1,walk_num
|
|
info = qmckl_distance_rescaled_deriv_e(context, 'T', 'T', elec_num, nucl_num, &
|
|
elec_coord(1,k,1), elec_num*walk_num, &
|
|
nucl_coord, nucl_num, &
|
|
en_distance_rescaled_deriv_e(1,1,1,k), elec_num, rescale_factor_kappa_en)
|
|
if (info /= QMCKL_SUCCESS) then
|
|
exit
|
|
endif
|
|
end do
|
|
|
|
end function qmckl_compute_en_distance_rescaled_deriv_e_f
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval h_private_func) :comments org :exports none
|
|
qmckl_exit_code qmckl_compute_en_distance_rescaled_deriv_e (
|
|
const qmckl_context context,
|
|
const int64_t elec_num,
|
|
const int64_t nucl_num,
|
|
const double rescale_factor_kappa_en,
|
|
const int64_t walk_num,
|
|
const double* elec_coord,
|
|
const double* nucl_coord,
|
|
double* const en_distance_rescaled_deriv_e );
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_en_distance_rescaled_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_en_distance_rescaled_deriv_e &
|
|
(context, elec_num, nucl_num, rescale_factor_kappa_en, walk_num, elec_coord, nucl_coord, en_distance_rescaled_deriv_e) &
|
|
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 :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
real (c_double ) , intent(in) , value :: rescale_factor_kappa_en
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3)
|
|
real (c_double ) , intent(in) :: nucl_coord(elec_num,3)
|
|
real (c_double ) , intent(out) :: en_distance_rescaled_deriv_e(elec_num,nucl_num,walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_en_distance_rescaled_deriv_e_f
|
|
info = qmckl_compute_en_distance_rescaled_deriv_e_f &
|
|
(context, elec_num, nucl_num, rescale_factor_kappa_en, walk_num, elec_coord, nucl_coord, en_distance_rescaled_deriv_e)
|
|
|
|
end function qmckl_compute_en_distance_rescaled_deriv_e
|
|
#+end_src
|
|
|
|
*** Test
|
|
|
|
#+begin_src python :results output :exports none
|
|
import numpy as np
|
|
|
|
# TODO
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
|
|
assert(qmckl_electron_provided(context));
|
|
|
|
rc = qmckl_set_nucleus_num (context, nucl_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_nucleus_rescale_factor (context, nucl_rescale_factor_kappa);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_nucleus_charge (context, charge, nucl_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord, 3*nucl_num);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
assert(qmckl_nucleus_provided(context));
|
|
|
|
double en_distance_rescaled_deriv_e[walk_num][4][nucl_num][elec_num];
|
|
|
|
rc = qmckl_get_electron_en_distance_rescaled_deriv_e(context, &(en_distance_rescaled_deriv_e[0][0][0][0]));
|
|
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
// TODO: check exact values
|
|
//// (e,n,w) in Fortran notation
|
|
//// (1,1,1)
|
|
//assert(fabs(en_distance_rescaled[0][0][0] - 7.546738741619978) < 1.e-12);
|
|
//
|
|
//// (1,2,1)
|
|
//assert(fabs(en_distance_rescaled[0][1][0] - 8.77102435246984) < 1.e-12);
|
|
//
|
|
//// (2,1,1)
|
|
//assert(fabs(en_distance_rescaled[0][0][1] - 3.698922010513608) < 1.e-12);
|
|
//
|
|
//// (1,1,2)
|
|
//assert(fabs(en_distance_rescaled[1][0][0] - 5.824059436060509) < 1.e-12);
|
|
//
|
|
//// (1,2,2)
|
|
//assert(fabs(en_distance_rescaled[1][1][0] - 7.080482110317645) < 1.e-12);
|
|
//
|
|
//// (2,1,2)
|
|
//assert(fabs(en_distance_rescaled[1][0][1] - 3.1804527583077356) < 1.e-12);
|
|
|
|
#+end_src
|
|
|
|
** Electron-nucleus potential ~en_potential~ stores the ~en~ potential energy
|
|
|
|
\[
|
|
\mathcal{V}_{en} = -\sum_{i=1}^{N_e}\sum_{A=1}^{N_n}\frac{Z_A}{r_{iA}}
|
|
\]
|
|
|
|
where \(\mathcal{V}_{en}\) is the ~en~ potential, \[r_{iA}\] the ~en~
|
|
distance and \[Z_A\] is the nuclear charge.
|
|
|
|
*** Get
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* const en_pot);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* const en_pot)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_en_potential(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = ctx->electron.walk_num * sizeof(double);
|
|
memcpy(en_pot, ctx->electron.en_pot, sze);
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_en_potential(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_en_potential(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
if (!ctx->electron.provided) return QMCKL_NOT_PROVIDED;
|
|
if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED;
|
|
|
|
qmckl_exit_code rc = qmckl_provide_en_distance(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->electron.coord_new_date > ctx->electron.en_pot_date) {
|
|
|
|
/* Allocate array */
|
|
if (ctx->electron.en_pot == NULL) {
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = ctx->electron.walk_num * sizeof(double);
|
|
double* en_pot = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (en_pot == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_en_potential",
|
|
NULL);
|
|
}
|
|
ctx->electron.en_pot = en_pot;
|
|
}
|
|
|
|
rc = qmckl_compute_en_potential(context,
|
|
ctx->electron.num,
|
|
ctx->nucleus.num,
|
|
ctx->electron.walk_num,
|
|
ctx->nucleus.charge.data,
|
|
ctx->electron.en_distance,
|
|
ctx->electron.en_pot);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->electron.en_pot_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_en_potential
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_en_potential_args
|
|
| Variable | Type | In/Out | Description |
|
|
|---------------+----------------------------------------+--------+--------------------------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
|
|
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
|
|
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
|
|
| ~charge~ | ~double[nucl_num]~ | in | charge of nucleus |
|
|
| ~en_distance~ | ~double[walk_num][nucl_num][elec_num]~ | in | Electron-electron rescaled distances |
|
|
| ~en_pot~ | ~double[walk_num]~ | out | Electron-electron potential |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_en_potential_f(context, elec_num, nucl_num, walk_num, &
|
|
charge, en_distance, en_pot) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: elec_num
|
|
integer*8 , intent(in) :: nucl_num
|
|
integer*8 , intent(in) :: walk_num
|
|
double precision , intent(in) :: charge(nucl_num)
|
|
double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num)
|
|
double precision , intent(out) :: en_pot(walk_num)
|
|
|
|
integer*8 :: nw, i, j
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
en_pot = 0.0d0
|
|
do nw=1,walk_num
|
|
do j=1,nucl_num
|
|
do i=1,elec_num
|
|
if (dabs(en_distance(i,j,nw)) > 1e-5) then
|
|
en_pot(nw) = en_pot(nw) - charge(j)/(en_distance(i,j,nw))
|
|
endif
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end function qmckl_compute_en_potential_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_header(table=qmckl_en_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_compute_en_potential (
|
|
const qmckl_context context,
|
|
const int64_t elec_num,
|
|
const int64_t nucl_num,
|
|
const int64_t walk_num,
|
|
const double* charge,
|
|
const double* en_distance,
|
|
double* const en_pot );
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_en_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_en_potential &
|
|
(context, elec_num, nucl_num, walk_num, charge, en_distance, en_pot) &
|
|
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 :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
real (c_double ) , intent(in) :: charge(nucl_num)
|
|
real (c_double ) , intent(in) :: en_distance(elec_num,nucl_num,walk_num)
|
|
real (c_double ) , intent(out) :: en_pot(walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_en_potential_f
|
|
info = qmckl_compute_en_potential_f &
|
|
(context, elec_num, nucl_num, walk_num, charge, en_distance, en_pot)
|
|
|
|
end function qmckl_compute_en_potential
|
|
#+end_src
|
|
|
|
*** Test
|
|
#+begin_src c :tangle (eval c_test)
|
|
double en_pot[walk_num];
|
|
|
|
rc = qmckl_get_electron_en_potential(context, &(en_pot[0]));
|
|
assert (rc == QMCKL_SUCCESS);
|
|
#+end_src
|
|
|
|
** Generate initial coordinates
|
|
|
|
*** Compute :noexport:
|
|
|
|
# begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
subroutine draw_init_points
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Place randomly electrons around nuclei
|
|
END_DOC
|
|
integer :: iwalk
|
|
logical, allocatable :: do_elec(:)
|
|
integer :: acc_num
|
|
|
|
real, allocatable :: xmin(:,:)
|
|
|
|
integer :: i, j, k, l, kk
|
|
|
|
real :: norm
|
|
allocate (do_elec(elec_num), xmin(3,elec_num))
|
|
xmin = -huge(1.)
|
|
norm = 0.
|
|
do i=1,elec_alpha_num
|
|
do j=1,ao_num
|
|
norm += mo_coef_transp(i,j)*mo_coef_transp(i,j)
|
|
enddo
|
|
enddo
|
|
norm = sqrt(norm/float(elec_alpha_num))
|
|
call rinfo( irp_here, 'Norm : ', norm )
|
|
call rinfo( irp_here, 'mo_scale: ' , mo_scale )
|
|
mo_coef_transp = mo_coef_transp/norm
|
|
|
|
double precision :: qmc_ranf
|
|
real :: mo_max
|
|
do i=1,elec_alpha_num
|
|
l=1
|
|
xmin(1,i) = mo_coef_transp(i,1)*mo_coef_transp(i,1) - 0.001*qmc_ranf()
|
|
do j=2,ao_num
|
|
xmin(2,i) = mo_coef_transp(i,j)*mo_coef_transp(i,j) - 0.001*qmc_ranf()
|
|
if (xmin(2,i) > xmin(1,i) ) then
|
|
xmin(1,i) = xmin(2,i)
|
|
l = ao_nucl(j)
|
|
endif
|
|
enddo
|
|
xmin(1,i) = nucl_coord(l,1)
|
|
xmin(2,i) = nucl_coord(l,2)
|
|
xmin(3,i) = nucl_coord(l,3)
|
|
enddo
|
|
|
|
call iinfo(irp_here, 'Det num = ', det_num )
|
|
do k=1,elec_beta_num
|
|
i = k+elec_alpha_num
|
|
l=1
|
|
xmin(1,i) = mo_coef_transp(k,1)*mo_coef_transp(k,1) - 0.001*qmc_ranf()
|
|
do j=2,ao_num
|
|
xmin(2,i) = mo_coef_transp(k,j)*mo_coef_transp(k,j) - 0.001*qmc_ranf()
|
|
if (xmin(2,i) > xmin(1,i) ) then
|
|
xmin(1,i) = xmin(2,i)
|
|
l = ao_nucl(j)
|
|
endif
|
|
enddo
|
|
xmin(1,i) = nucl_coord(l,1)
|
|
xmin(2,i) = nucl_coord(l,2)
|
|
xmin(3,i) = nucl_coord(l,3)
|
|
enddo
|
|
|
|
call rinfo( irp_here, 'time step =', time_step )
|
|
do iwalk=1,walk_num
|
|
print *, 'Generating initial positions for walker', iwalk
|
|
acc_num = 0
|
|
do_elec = .True.
|
|
integer :: iter
|
|
do iter = 1,10000
|
|
if (acc_num >= elec_num) then
|
|
exit
|
|
endif
|
|
double precision :: gauss
|
|
real :: re_compute
|
|
re_compute = 0.
|
|
do while (re_compute < 1.e-6)
|
|
do i=1,elec_num
|
|
if (do_elec(i)) then
|
|
do l=1,3
|
|
elec_coord(i,l) = xmin(l,i) + 1.5*(0.5-qmc_ranf())
|
|
enddo
|
|
endif
|
|
enddo
|
|
TOUCH elec_coord
|
|
re_compute = minval(nucl_elec_dist(1:nucl_num,1:elec_num))
|
|
enddo
|
|
|
|
do i=1,elec_alpha_num
|
|
if (do_elec(i)) then
|
|
if ( mo_value_transp(i,i)**2 >= qmc_ranf()) then
|
|
acc_num += 1
|
|
do_elec(i) = .False.
|
|
endif
|
|
endif
|
|
enddo
|
|
|
|
do i=1,elec_beta_num
|
|
if (do_elec(i+elec_alpha_num)) then
|
|
if ( mo_value_transp(i,i+elec_alpha_num)**2 >= qmc_ranf()) then
|
|
acc_num += 1
|
|
do_elec(i+elec_alpha_num) = .False.
|
|
endif
|
|
endif
|
|
enddo
|
|
|
|
enddo
|
|
|
|
do l=1,3
|
|
do i=1,elec_num+1
|
|
elec_coord_full(i,l,iwalk) = elec_coord(i,l)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
if (.not.is_worker) then
|
|
call ezfio_set_electrons_elec_coord_pool_size(walk_num)
|
|
call ezfio_set_electrons_elec_coord_pool(elec_coord_full)
|
|
endif
|
|
SOFT_TOUCH elec_coord elec_coord_full
|
|
deallocate (do_elec, xmin)
|
|
|
|
end
|
|
# 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)
|
|
if (qmckl_context_destroy(context) != QMCKL_SUCCESS)
|
|
return QMCKL_FAILURE;
|
|
return 0;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute file names
|
|
#+begin_src emacs-lisp
|
|
; The following is required to compute the file names
|
|
|
|
(setq pwd (file-name-directory buffer-file-name))
|
|
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
|
|
(setq f (concat pwd name "_f.f90"))
|
|
(setq fh (concat pwd name "_fh.f90"))
|
|
(setq c (concat pwd name ".c"))
|
|
(setq h (concat name ".h"))
|
|
(setq h_private (concat name "_private.h"))
|
|
(setq c_test (concat pwd "test_" name ".c"))
|
|
(setq f_test (concat pwd "test_" name "_f.f90"))
|
|
|
|
; Minted
|
|
(require 'ox-latex)
|
|
(setq org-latex-listings 'minted)
|
|
(add-to-list 'org-latex-packages-alist '("" "listings"))
|
|
(add-to-list 'org-latex-packages-alist '("" "color"))
|
|
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
| | color |
|
|
| | listings |
|
|
|
|
|
|
# -*- mode: org -*-
|
|
# vim: syntax=c
|
|
|
|
|