TREX
/
qmckl
mirror of https://github.com/TREX-CoE/qmckl.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'master' into add-python-api

This commit is contained in:
q-posev 2022-08-22 12:11:07 +02:00
parent 5a638d4697 5fffb52513
commit 0d6a07705c
23 changed files with 1584 additions and 1042 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
.github/workflows/test-build.yml vendored
View File

@ -40,7 +40,7 @@ jobs:
./autogen.sh
mkdir _build
cd _build
../configure --enable-silent-rules --enable-debug
../configure --enable-silent-rules
make -j 4
sudo make install
@ -76,6 +76,24 @@ jobs:
run: make python-test
working-directory: _build
- name: Configure QMCkl in debug mode
run: |
mkdir _build_debug
cd _build_debug
../configure --enable-debug --enable-silent-rules
make -j2
- name: Run test
run: make -j2 check
working-directory: _build_debug
- name: Archive test log file
if: failure()
uses: actions/upload-artifact@v2
with:
name: test-report-ubuntu-debug
path: _build_debug/test-suite.log
# x86_macos:
#
# runs-on: macos-latest

3
.gitmodules vendored Normal file
View File

@ -0,0 +1,3 @@
[submodule "share/doc/qmckl/html/emacs-htmlize"]
path = share/doc/qmckl/html/emacs-htmlize
url = https://github.com/TREX-CoE/emacs-htmlize

5
Makefile.am
View File

@ -88,6 +88,9 @@ dist_html_DATA = $(HTML_FILES) \
dist_text_DATA = $(TEXT_FILES)
# Include Guix manifest in the source code distribution tarball
tools_qmckl_scm = tools/qmckl.scm
EXTRA_DIST += $(tools_qmckl_scm)
html-local: $(htmlize_el) $(dist_html_DATA)
text: $(htmlize_el) $(dist_text_DATA)
@ -202,7 +205,7 @@ cppcheck.out: $(qmckl_h)
$(qmckl_include_i): $(qmckl_h) $(process_header_py)
$(MKDIR_P) python/src
python $(process_header_py) $(qmckl_h)
python3 $(process_header_py) $(qmckl_h)
mv qmckl_include.i $(qmckl_include_i)

25
README.md
View File

@ -76,6 +76,31 @@ sudo make installcheck
```
## Installation procedure for Guix users
QMCkl can be installed with the [GNU Guix](https://guix.gnu.org) functional package manager.
The [qmckl.scm](https://github.com/TREX-CoE/qmckl/blob/master/tools/qmckl.scm)
Schema file contains the manifest specification for the `qmckl` installations.
It can be installed within the selected `$GUIX_PROFILE` as follows:
```
guix package \
--profile=$GUIX_PROFILE \
--load-path=<path_to_trexio_scm> \
--cores=<n_cores> \
--install-from-file=qmckl.scm
```
where `<path_to_trexio_scm>` should point to a folder, which contains the TREXIO manifest file
[trexio.scm](https://github.com/TREX-CoE/trexio/blob/master/tools/trexio.scm)
(e.g. `~/trexio/tools/` if TREXIO repository was cloned under $HOME).
Installation procedures for both development version (`qmckl-dev`)
and stable releases (`qmckl-hpc`) are provided.
One can switch between them using the return value (last line)
in the `qmckl.scm` file.
## Linking to your program
The `make install` command takes care of installing the QMCkl shared library on the user machine.

4
configure.ac
View File

@ -35,7 +35,7 @@
AC_PREREQ([2.69])
AC_INIT([qmckl],[0.2.1],[https://github.com/TREX-CoE/qmckl/issues],[],[https://trex-coe.github.io/qmckl/index.html])
AC_INIT([qmckl],[0.3.1],[https://github.com/TREX-CoE/qmckl/issues],[],[https://trex-coe.github.io/qmckl/index.html])
AC_CONFIG_AUX_DIR([tools])
AM_INIT_AUTOMAKE([subdir-objects color-tests parallel-tests silent-rules 1.11])
@ -358,8 +358,10 @@ if test "$ok" = "yes"; then
if test "$GCC" = "yes"; then
CFLAGS="$CFLAGS \
-g -Wall -W -Wbad-function-cast -Wcast-qual -Warray-bounds -Wdisabled-optimization \
-fsanitize=address -fno-omit-frame-pointer \
-Wpointer-arith -Wcast-align -Wpedantic -Wextra -Walloc-zero -Werror \
"
LDFLAGS="$LDFLAGS -fsanitize=address"
fi
if test "$GFC" = "yes"; then
FCFLAGS="$FCFLAGS \

368
org/qmckl_ao.org
View File

@ -2423,6 +2423,7 @@ assert (rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i < nucl_num ; ++i) {
assert(nucleus_shell_num_test[i] == nucleus_shell_num[i]);
}
free(nucleus_shell_num_test);
shell_ang_mom_test = (int32_t*) malloc ( shell_num * sizeof(int32_t));
rc = qmckl_get_ao_basis_shell_ang_mom (context, shell_ang_mom_test, shell_num);
@ -3013,7 +3014,7 @@ qmckl_get_ao_basis_ao_vgl (qmckl_context context,
qmckl_exit_code rc;
rc = qmckl_provide_ao_vgl(context);
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
@ -3089,7 +3090,7 @@ qmckl_get_ao_basis_ao_vgl_inplace (qmckl_context context,
ctx->ao_basis.ao_vgl = ao_vgl;
rc = qmckl_provide_ao_vgl(context);
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
ctx->ao_basis.ao_vgl = old_array;
@ -3638,7 +3639,7 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context);
qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
{
qmckl_exit_code rc;
qmckl_exit_code rc = QMCKL_SUCCESS;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
@ -3660,8 +3661,17 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
/* Compute if necessary */
if (ctx->point.date > ctx->ao_basis.primitive_vgl_date) {
if (ctx->point.alloc_date > ctx->ao_basis.primitive_vgl_date) {
if (ctx->ao_basis.primitive_vgl != NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.prim_num * 5 * ctx->point.num * sizeof(double);
if (ctx->ao_basis.primitive_vgl != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->ao_basis.primitive_vgl, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->ao_basis.primitive_vgl);
assert (rc == QMCKL_SUCCESS);
ctx->ao_basis.primitive_vgl = NULL;
@ -3671,8 +3681,6 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
/* Allocate array */
if (ctx->ao_basis.primitive_vgl == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.prim_num * 5 * ctx->point.num * sizeof(double);
double* primitive_vgl = (double*) qmckl_malloc(context, mem_info);
if (primitive_vgl == NULL) {
@ -3779,20 +3787,19 @@ print ( "[7][4][26] : %e"% lf(a,x,y))
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_electron_walk_num (context, walk_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
const int64_t point_num = elec_num;
rc = qmckl_set_point(context, 'N', point_num, elec_coord, point_num*3);
assert(rc == QMCKL_SUCCESS);
double prim_vgl[elec_num*walk_num][5][prim_num];
double prim_vgl[point_num][5][prim_num];
rc = qmckl_get_ao_basis_primitive_vgl(context, &(prim_vgl[0][0][0]),
(int64_t) 5*elec_num*walk_num*prim_num );
(int64_t) 5*point_num*prim_num );
assert (rc == QMCKL_SUCCESS);
printf("prim_vgl[26][0][7] = %e\n",prim_vgl[26][0][7]);
@ -4067,7 +4074,7 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context);
qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
{
qmckl_exit_code rc;
qmckl_exit_code rc = QMCKL_SUCCESS;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
@ -4089,8 +4096,17 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
/* Compute if necessary */
if (ctx->point.date > ctx->ao_basis.shell_vgl_date) {
if (ctx->point.alloc_date > ctx->ao_basis.shell_vgl_date) {
if (ctx->ao_basis.shell_vgl != NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.shell_num * 5 * ctx->point.num * sizeof(double);
if (ctx->ao_basis.shell_vgl != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->ao_basis.shell_vgl, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->ao_basis.shell_vgl);
assert (rc == QMCKL_SUCCESS);
ctx->ao_basis.shell_vgl = NULL;
@ -4100,8 +4116,6 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
/* Allocate array */
if (ctx->ao_basis.shell_vgl == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.shell_num * 5 * ctx->point.num * sizeof(double);
double* shell_vgl = (double*) qmckl_malloc(context, mem_info);
if (shell_vgl == NULL) {
@ -4220,14 +4234,15 @@ print ( "[1][4][26] : %25.15e"% lf(a,x,y))
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
const int64_t point_num = elec_num;
rc = qmckl_set_point(context, 'N', point_num, elec_coord, point_num*3);
assert(rc == QMCKL_SUCCESS);
double shell_vgl[elec_num][5][shell_num];
double shell_vgl[point_num][5][shell_num];
rc = qmckl_get_ao_basis_shell_vgl(context, &(shell_vgl[0][0][0]),
(int64_t) 5*elec_num*shell_num);
(int64_t) 5*point_num*shell_num);
assert (rc == QMCKL_SUCCESS);
printf(" shell_vgl[26][0][1] %25.15e\n", shell_vgl[26][0][1]);
@ -5633,33 +5648,14 @@ end function qmckl_compute_ao_value_doc_f
| ~coef_normalized~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells |
| ~ao_value~ | ~double[point_num][ao_num]~ | out | Values of the AOs |
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
const int64_t ao_num,
const int64_t shell_num,
const int32_t* restrict prim_num_per_nucleus,
const int64_t point_num,
const int64_t nucl_num,
const double* restrict coord,
const double* restrict nucl_coord,
const int64_t* restrict nucleus_index,
const int64_t* restrict nucleus_shell_num,
const double* nucleus_range,
const int32_t* restrict nucleus_max_ang_mom,
const int32_t* restrict shell_ang_mom,
const double* restrict ao_factor,
const qmckl_matrix expo_per_nucleus,
const qmckl_tensor coef_per_nucleus,
double* restrict const ao_value )
{
int32_t lstart[32];
#+NAME:ao_value_constants
#+begin_src c :exports none
int32_t lstart[32] __attribute__((aligned(64)));
for (int32_t l=0 ; l<32 ; ++l) {
lstart[l] = l*(l+1)*(l+2)/6;
}
int64_t ao_index[shell_num+1];
int64_t ao_index[shell_num+1] __attribute__((aligned(64)));
int64_t size_max = 0;
int64_t prim_max = 0;
int64_t shell_max = 0;
@ -5685,17 +5681,42 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
/* Don't compute polynomials when the radial part is zero. */
double cutoff = 27.631021115928547; // -log(1.e-12)
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
const int64_t ao_num,
const int64_t shell_num,
const int32_t* restrict prim_num_per_nucleus,
const int64_t point_num,
const int64_t nucl_num,
const double* restrict coord,
const double* restrict nucl_coord,
const int64_t* restrict nucleus_index,
const int64_t* restrict nucleus_shell_num,
const double* nucleus_range,
const int32_t* restrict nucleus_max_ang_mom,
const int32_t* restrict shell_ang_mom,
const double* restrict ao_factor,
const qmckl_matrix expo_per_nucleus,
const qmckl_tensor coef_per_nucleus,
double* restrict const ao_value )
{
<<ao_value_constants>>
#ifdef HAVE_OPENMP
#pragma omp parallel
#endif
{
qmckl_exit_code rc;
double ar2[prim_max];
int32_t powers[3*size_max];
double poly_vgl[5*size_max];
double ar2[prim_max] __attribute__((aligned(64)));
int32_t powers[3*size_max] __attribute__((aligned(64)));
double poly_vgl[5*size_max] __attribute__((aligned(64)));
double exp_mat[prim_max];
double ce_mat[shell_max];
double exp_mat[prim_max] __attribute__((aligned(64)));
double ce_mat[shell_max] __attribute__((aligned(64)));
double coef_mat[nucl_num][shell_max][prim_max];
for (int i=0 ; i<nucl_num ; ++i) {
@ -5710,15 +5731,17 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
#pragma omp for
#endif
for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
const double e_coord[3] = { coord[ipoint],
coord[ipoint + point_num],
coord[ipoint + 2*point_num] };
const double e_coord[3] __attribute__((aligned(64))) =
{ coord[ipoint],
coord[ipoint + point_num],
coord[ipoint + 2*point_num] };
for (int64_t inucl=0 ; inucl < nucl_num ; ++inucl) {
const double n_coord[3] = { nucl_coord[inucl],
nucl_coord[inucl + nucl_num],
nucl_coord[inucl + 2*nucl_num] };
const double n_coord[3] __attribute__((aligned(64))) =
{ nucl_coord[inucl],
nucl_coord[inucl + nucl_num],
nucl_coord[inucl + 2*nucl_num] };
/* Test if the point is in the range of the nucleus */
const double x = e_coord[0] - n_coord[0];
const double y = e_coord[1] - n_coord[1];
@ -5736,14 +5759,14 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
break;
case 1:
poly_vgl[0] = 0.;
poly_vgl[0] = 1.;
poly_vgl[1] = x;
poly_vgl[2] = y;
poly_vgl[3] = z;
break;
case 2:
poly_vgl[0] = 0.;
poly_vgl[0] = 1.;
poly_vgl[1] = x;
poly_vgl[2] = y;
poly_vgl[3] = z;
@ -5782,10 +5805,8 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
}
for (int i=0 ; i<nucleus_shell_num[inucl] ; ++i) {
ce_mat[i] = 0.;
}
for (int k=0 ; k<nidx; ++k) {
for (int i=0 ; i<nucleus_shell_num[inucl] ; ++i) {
ce_mat[i] = 0.;
for (int k=0 ; k<nidx; ++k) {
ce_mat[i] += coef_mat[inucl][i][k] * exp_mat[k];
}
}
@ -5796,7 +5817,6 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {
const double s1 = ce_mat[ishell-ishell_start];
if (s1 == 0.0) continue;
const int64_t k = ao_index[ishell];
double* restrict const ao_value_1 = ao_value + ipoint*ao_num + k;
@ -5804,6 +5824,13 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
const int32_t l = shell_ang_mom[ishell];
const int32_t n = lstart[l+1]-lstart[l];
if (s1 == 0.0) {
for (int64_t il=0 ; il<n ; ++il) {
ao_value_1[il] = 0.;
}
continue;
}
double* restrict poly_vgl_1 = NULL;
if (nidx > 0) {
const double* restrict f = ao_factor + k;
@ -5812,10 +5839,10 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
poly_vgl_1 = &(poly_vgl[idx]);
switch (n) {
case(1):
case 1:
ao_value_1[0] = s1 * f[0];
break;
case (3):
case 3:
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
@ -5833,7 +5860,7 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
break;
default:
#ifdef HAVE_OPENMP
#pragma omp simd simdlen(8)
#pragma omp simd
#endif
for (int il=0 ; il<n ; ++il) {
ao_value_1[il] = poly_vgl_1[il] * s1 * f[il];
@ -5979,7 +6006,7 @@ qmckl_exit_code qmckl_provide_ao_basis_ao_value(qmckl_context context);
qmckl_exit_code qmckl_provide_ao_basis_ao_value(qmckl_context context)
{
qmckl_exit_code rc;
qmckl_exit_code rc = QMCKL_SUCCESS;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
@ -6001,8 +6028,17 @@ qmckl_exit_code qmckl_provide_ao_basis_ao_value(qmckl_context context)
/* Compute if necessary */
if (ctx->point.date > ctx->ao_basis.ao_value_date) {
if (ctx->point.alloc_date > ctx->ao_basis.ao_value_date) {
if (ctx->ao_basis.ao_value != NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.ao_num * ctx->point.num * sizeof(double);
if (ctx->ao_basis.ao_value != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->ao_basis.ao_value, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->ao_basis.ao_value);
assert (rc == QMCKL_SUCCESS);
ctx->ao_basis.ao_value = NULL;
@ -6012,8 +6048,6 @@ qmckl_exit_code qmckl_provide_ao_basis_ao_value(qmckl_context context)
/* Allocate array */
if (ctx->ao_basis.ao_value == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.ao_num * ctx->point.num * sizeof(double);
double* ao_value = (double*) qmckl_malloc(context, mem_info);
if (ao_value == NULL) {
@ -6210,14 +6244,15 @@ double* elec_coord = &(chbrclf_elec_coord[0][0][0]);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
const int64_t point_num = elec_num;
rc = qmckl_set_point(context, 'N', point_num, elec_coord, point_num*3);
assert(rc == QMCKL_SUCCESS);
double ao_value[elec_num][ao_num];
double ao_value[point_num][ao_num];
rc = qmckl_get_ao_basis_ao_value(context, &(ao_value[0][0]),
(int64_t) elec_num*ao_num);
(int64_t) point_num*ao_num);
assert (rc == QMCKL_SUCCESS);
printf("\n");
@ -6447,36 +6482,7 @@ qmckl_compute_ao_vgl_hpc_gaussian (
const qmckl_tensor coef_per_nucleus,
double* restrict const ao_vgl )
{
int32_t lstart[32];
for (int32_t l=0 ; l<32 ; ++l) {
lstart[l] = l*(l+1)*(l+2)/6;
}
int64_t ao_index[shell_num+1];
int64_t size_max = 0;
int64_t prim_max = 0;
int64_t shell_max = 0;
{
int64_t k=0;
for (int inucl=0 ; inucl < nucl_num ; ++inucl) {
prim_max = prim_num_per_nucleus[inucl] > prim_max ?
prim_num_per_nucleus[inucl] : prim_max;
shell_max = nucleus_shell_num[inucl] > shell_max ?
nucleus_shell_num[inucl] : shell_max;
const int64_t ishell_start = nucleus_index[inucl];
const int64_t ishell_end = nucleus_index[inucl] + nucleus_shell_num[inucl];
for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {
const int l = shell_ang_mom[ishell];
ao_index[ishell] = k;
k += lstart[l+1] - lstart[l];
size_max = size_max < lstart[l+1] ? lstart[l+1] : size_max;
}
}
ao_index[shell_num] = ao_num+1;
}
/* Don't compute when the radial part is zero. */
double cutoff = 27.631021115928547; // -log(1.e-12)
<<ao_value_constants>>
#ifdef HAVE_OPENMP
#pragma omp parallel
@ -6485,18 +6491,6 @@ qmckl_compute_ao_vgl_hpc_gaussian (
qmckl_exit_code rc;
double ar2[prim_max] __attribute__((aligned(64)));
int32_t powers[3*size_max] __attribute__((aligned(64)));
double poly_vgl_l1[4][4] __attribute__((aligned(64))) =
{{1.0, 0.0, 0.0, 0.0},
{0.0, 1.0, 0.0, 0.0},
{0.0, 0.0, 1.0, 0.0},
{0.0, 0.0, 0.0, 1.0}};
double poly_vgl_l2[5][10]__attribute__((aligned(64))) =
{{1., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
{0., 1., 0., 0., 0., 0., 0., 0., 0., 0.},
{0., 0., 1., 0., 0., 0., 0., 0., 0., 0.},
{0., 0., 0., 1., 0., 0., 0., 0., 0., 0.},
{0., 0., 0., 0., 2., 0., 0., 2., 0., 2.}};
double poly_vgl[5][size_max] __attribute__((aligned(64)));
double exp_mat[prim_max][8] __attribute__((aligned(64))) ;
double ce_mat[shell_max][8] __attribute__((aligned(64))) ;
@ -6510,20 +6504,33 @@ qmckl_compute_ao_vgl_hpc_gaussian (
}
}
double poly_vgl_l1[4][4] __attribute__((aligned(64))) =
{{1.0, 0.0, 0.0, 0.0},
{0.0, 1.0, 0.0, 0.0},
{0.0, 0.0, 1.0, 0.0},
{0.0, 0.0, 0.0, 1.0}};
double poly_vgl_l2[5][10]__attribute__((aligned(64))) =
{{1., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
{0., 1., 0., 0., 0., 0., 0., 0., 0., 0.},
{0., 0., 1., 0., 0., 0., 0., 0., 0., 0.},
{0., 0., 0., 1., 0., 0., 0., 0., 0., 0.},
{0., 0., 0., 0., 2., 0., 0., 2., 0., 2.}};
double poly_vgl[5][size_max] __attribute__((aligned(64)));
#ifdef HAVE_OPENMP
#pragma omp for
#endif
for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
const double e_coord[3] __attribute__((aligned(64))) =
{ coord[ipoint],
coord[ipoint + point_num],
coord[ipoint + 2*point_num] };
{ coord[ipoint],
coord[ipoint + point_num],
coord[ipoint + 2*point_num] };
for (int64_t inucl=0 ; inucl < nucl_num ; ++inucl) {
const double n_coord[3] __attribute__((aligned(64))) =
{ nucl_coord[inucl],
nucl_coord[inucl + nucl_num],
nucl_coord[inucl + 2*nucl_num] };
{ nucl_coord[inucl],
nucl_coord[inucl + nucl_num],
nucl_coord[inucl + 2*nucl_num] };
/* Test if the point is in the range of the nucleus */
const double x = e_coord[0] - n_coord[0];
@ -6593,6 +6600,7 @@ qmckl_compute_ao_vgl_hpc_gaussian (
for (int64_t iprim = 0 ; iprim < nidx ; ++iprim) {
exp_mat[iprim][0] = exp(-ar2[iprim]);
}
for (int64_t iprim = 0 ; iprim < nidx ; ++iprim) {
double f = qmckl_mat(expo_per_nucleus, iprim, inucl) * exp_mat[iprim][0];
f = -f-f;
@ -6605,21 +6613,6 @@ qmckl_compute_ao_vgl_hpc_gaussian (
/* --- */
switch (8) {
case(5):
for (int i=0 ; i<nucleus_shell_num[inucl] ; ++i) {
for (int j=0 ; j<5 ; ++j) {
ce_mat[i][j] = 0.;
}
for (int k=0 ; k<nidx; ++k) {
const double cm = coef_mat[inucl][i][k];
for (int j=0 ; j<5 ; ++j) {
ce_mat[i][j] += cm * exp_mat[k][j];
}
}
}
break;
case(8):
for (int i=0 ; i<nucleus_shell_num[inucl] ; ++i) {
@ -6641,6 +6634,21 @@ qmckl_compute_ao_vgl_hpc_gaussian (
}
break;
case(5):
for (int i=0 ; i<nucleus_shell_num[inucl] ; ++i) {
for (int j=0 ; j<5 ; ++j) {
ce_mat[i][j] = 0.;
}
for (int k=0 ; k<nidx; ++k) {
const double cm = coef_mat[inucl][i][k];
for (int j=0 ; j<5 ; ++j) {
ce_mat[i][j] += cm * exp_mat[k][j];
}
}
}
break;
case(512):
for(int i=0; i<nucleus_shell_num[inucl]; ++i){
__m512d cemat_avx512;
@ -6755,14 +6763,14 @@ qmckl_compute_ao_vgl_hpc_gaussian (
poly_vgl_5 = &(poly_vgl[4][idx]);
}
switch (n) {
case(1):
case 1:
ao_vgl_1[0] = s1 * f[0];
ao_vgl_2[0] = s2 * f[0];
ao_vgl_3[0] = s3 * f[0];
ao_vgl_4[0] = s4 * f[0];
ao_vgl_5[0] = s5;
break;
case (3):
case 3:
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
@ -6777,7 +6785,7 @@ qmckl_compute_ao_vgl_hpc_gaussian (
poly_vgl_4[il] * s4 )) * f[il];
}
break;
case(6):
case 6:
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
@ -6937,18 +6945,26 @@ qmckl_compute_ao_vgl_hpc_gaussian (
**** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context);
#+end_src
#+CALL: write_provider_header( group="ao_basis", data="ao_vgl" )
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
#+RESULTS:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :export none
qmckl_exit_code qmckl_provide_ao_basis_ao_vgl(qmckl_context context);
#+end_src
#+CALL: write_provider_pre( group="ao_basis", data="ao_vgl", dimension="ctx->ao_basis.ao_num * 5 * ctx->point.num")
#+RESULTS:
#+begin_src c :comments org :tangle (eval c) :noweb yes :export none
qmckl_exit_code qmckl_provide_ao_basis_ao_vgl(qmckl_context context)
{
qmckl_exit_code rc = QMCKL_SUCCESS;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_provide_ao_vgl",
"qmckl_provide_ao_basis_ao_vgl",
NULL);
}
@ -6958,30 +6974,33 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_vgl",
"qmckl_provide_ao_basis_ao_vgl",
NULL);
}
/* Compute if necessary */
if (ctx->point.date > ctx->ao_basis.ao_vgl_date) {
qmckl_exit_code rc;
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.ao_num * 5 * ctx->point.num * sizeof(double);
/* Provide required data */
#ifdef HAVE_HPC
if (ctx->ao_basis.ao_vgl != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->ao_basis.ao_vgl, &mem_info_test);
#else
rc = qmckl_provide_ao_basis_shell_vgl(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_shell_vgl", NULL);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->ao_basis.ao_vgl);
assert (rc == QMCKL_SUCCESS);
ctx->ao_basis.ao_vgl = NULL;
}
}
#endif
/* Allocate array */
if (ctx->ao_basis.ao_vgl == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.ao_num * 5 * ctx->point.num * sizeof(double);
double* ao_vgl = (double*) qmckl_malloc(context, mem_info);
if (ao_vgl == NULL) {
@ -6992,6 +7011,10 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
}
ctx->ao_basis.ao_vgl = ao_vgl;
}
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
#ifdef HAVE_HPC
if (ctx->ao_basis.type == 'G') {
rc = qmckl_compute_ao_vgl_hpc_gaussian(context,
@ -7068,6 +7091,13 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
ctx->ao_basis.ao_vgl);
,*/
} else {
/* Provide required data */
rc = qmckl_provide_ao_basis_shell_vgl(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_shell_vgl", NULL);
}
rc = qmckl_compute_ao_vgl_doc(context,
ctx->ao_basis.ao_num,
ctx->ao_basis.shell_num,
@ -7085,6 +7115,11 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
ctx->ao_basis.ao_vgl);
}
#else
rc = qmckl_provide_ao_basis_shell_vgl(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_shell_vgl", NULL);
}
rc = qmckl_compute_ao_vgl_doc(context,
ctx->ao_basis.ao_num,
ctx->ao_basis.shell_num,
@ -7101,6 +7136,12 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
ctx->ao_basis.shell_vgl,
ctx->ao_basis.ao_vgl);
#endif
#+end_src
#+CALL: write_provider_post( group="ao_basis", data="ao_vgl" )
#+RESULTS:
#+begin_src c :comments org :tangle (eval c) :noweb yes :export none
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -7110,7 +7151,7 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
return QMCKL_SUCCESS;
}
#+end_src
#+end_src
**** Test :noexport:
@ -7196,14 +7237,15 @@ double* elec_coord = &(chbrclf_elec_coord[0][0][0]);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
const int64_t point_num = elec_num;
rc = qmckl_set_point(context, 'N', point_num, elec_coord, point_num*3);
assert(rc == QMCKL_SUCCESS);
double ao_vgl[elec_num][5][ao_num];
double ao_vgl[point_num][5][ao_num];
rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]),
(int64_t) 5*elec_num*ao_num);
(int64_t) 5*point_num*ao_num);
assert (rc == QMCKL_SUCCESS);
printf("\n");

4
org/qmckl_context.org
View File

@ -212,7 +212,7 @@ qmckl_context_touch(const qmckl_context context)
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
ctx->date += 1UL;
ctx->point.date += 1UL;
ctx->point.date = ctx-> date;
return QMCKL_SUCCESS;
}
#+end_src
@ -520,11 +520,13 @@ assert(qmckl_context_check(context) == context);
assert(qmckl_context_destroy(context) == QMCKL_SUCCESS);
/* Check that context is destroyed */
#ifndef DEBUG
assert(qmckl_context_check(context) != context);
assert(qmckl_context_check(context) == QMCKL_NULL_CONTEXT);
/* Destroy invalid context */
assert(qmckl_context_destroy(QMCKL_NULL_CONTEXT) == QMCKL_INVALID_CONTEXT);
#endif
#+end_src
* End of files :noexport:

234
org/qmckl_determinant.org
View File

@ -99,42 +99,40 @@ int main() {
|------------------+------------------------------------------------+------------------------------------|
| ~type~ | ~char~ | α (~'A'~) or β (~'B'~) determinant |
| ~walk_num~ | ~int64_t~ | Number of walkers |
| ~det_num_alpha~ | ~int64_t~ | Number of determinants per walker |
| ~det_num_beta~ | ~int64_t~ | Number of determinants per walker |
| ~mo_index_alpha~ | ~mo_index[det_num_alpha][walk_num][alpha_num]~ | Index of MOs for each walker |
| ~mo_index_beta~ | ~mo_index[det_num_beta][walk_num][beta_num]~ | Index of MOs for each walker |
| ~mo_index_alpha~ | ~mo_index[det_num_alpha][walker.num][alpha_num]~ | Index of MOs for each walker |
| ~mo_index_beta~ | ~mo_index[det_num_beta][walker.num][beta_num]~ | Index of MOs for each walker |
Computed data:
|-----------------------------+------------------------------------------------------+-------------------------------------------------------------------------------------------|
| ~up_num~ | ~int64_t~ | Number of number of α electrons |
| ~donwn_num~ | ~int64_t~ | Number of number of β electrons |
| ~det_value_alpha~ | ~[det_num_alpha][walk_num]~ | The α slater matrix for each determinant of each walker. |
| ~det_value_alpha_date~ | ~int64_t~ | Date of The α slater matrix for each determinant of each walker. |
| ~det_value_beta~ | ~[det_num_beta][walk_num]~ | The β slater matrix for each determinant of each walker. |
| ~det_value_beta_date~ | ~int64_t~ | Date of The β slater matrix for each determinant of each walker. |
| ~det_adj_matrix_alpha~ | ~[det_num_alpha][walk_num][alpha_num][alpha_num]~ | Adjoint of the α slater matrix for each determinant of each walker. |
| ~det_adj_matrix_alpha_date~ | ~int64_t~ | Date of the Adjoint of the α slater matrix for each determinant of each walker. |
| ~det_adj_matrix_beta~ | ~[det_num_beta][walk_num][beta_num][beta_num]~ | Adjoint of the β slater matrix for each determinant of each walker. |
| ~det_adj_matrix_beta_date~ | ~int64_t~ | Date of the Adjoint of the β slater matrix for each determinant of each walker. |
|-----------------------------+------------------------------------------------------+-------------------------------------------------------------------------------------------|
| ~det_vgl_alpha~ | ~[5][det_num_alpha][walk_num][alpha_num][alpha_num]~ | Value, gradients, Laplacian of Dᵅᵢⱼ(x) at electron positions |
| ~det_vgl_alpha_date~ | ~int64_t~ | Late modification date of Value, gradients, Laplacian of the MOs at electron positions |
| ~det_vgl_beta~ | ~[5][det_num_beta][walk_num][beta_num][beta_num]~ | Value, gradients, Laplacian of Dᵝᵢⱼ(x) at electron positions |
| ~det_vgl_beta_date~ | ~int64_t~ | Late modification date of Value, gradients, Laplacian of the MOs at electron positions |
| ~det_inv_matrix_alpha~ | ~[det_num_alpha][walk_num][alpha_num][alpha_num]~ | Inverse of the α electron slater matrix for each determinant of each walker. |
| ~det_inv_matrix_alpha_date~ | ~int64_t~ | Date for the Inverse of the α electron slater matrix for each determinant of each walker. |
| ~det_inv_matrix_beta~ | ~[det_num_beta][walk_num][beta_num][beta_num]~ | Inverse of the β electron slater matrix for each determinant of each walker. |
| ~det_inv_matrix_beta_date~ | ~int64_t~ | Date for the Inverse of the β electron slater matrix for each determinant of each walker. |
|-----------------------------+------------------------------------------------------+-------------------------------------------------------------------------------------------|
|-----------------------------+--------------------------------------------------------+-------------------------------------------------------------------------------------------|
| ~up_num~ | ~int64_t~ | Number of number of α electrons |
| ~donwn_num~ | ~int64_t~ | Number of number of β electrons |
| ~det_value_alpha~ | ~[det_num_alpha][walker.num]~ | The α slater matrix for each determinant of each walker. |
| ~det_value_alpha_date~ | ~uint64_t~ | Date of The α slater matrix for each determinant of each walker. |
| ~det_value_beta~ | ~[det_num_beta][walker.num]~ | The β slater matrix for each determinant of each walker. |
| ~det_value_beta_date~ | ~uint64_t~ | Date of The β slater matrix for each determinant of each walker. |
| ~det_adj_matrix_alpha~ | ~[det_num_alpha][walker.num][alpha_num][alpha_num]~ | Adjoint of the α slater matrix for each determinant of each walker. |
| ~det_adj_matrix_alpha_date~ | ~uint64_t~ | Date of the Adjoint of the α slater matrix for each determinant of each walker. |
| ~det_adj_matrix_beta~ | ~[det_num_beta][walker.num][beta_num][beta_num]~ | Adjoint of the β slater matrix for each determinant of each walker. |
| ~det_adj_matrix_beta_date~ | ~uint64_t~ | Date of the Adjoint of the β slater matrix for each determinant of each walker. |
|-----------------------------+--------------------------------------------------------+-------------------------------------------------------------------------------------------|
| ~det_vgl_alpha~ | ~[5][det_num_alpha][walker.num][alpha_num][alpha_num]~ | Value, gradients, Laplacian of Dᵅᵢⱼ(x) at electron positions |
| ~det_vgl_alpha_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the MOs at electron positions |
| ~det_vgl_beta~ | ~[5][det_num_beta][walker.num][beta_num][beta_num]~ | Value, gradients, Laplacian of Dᵝᵢⱼ(x) at electron positions |
| ~det_vgl_beta_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the MOs at electron positions |
| ~det_inv_matrix_alpha~ | ~[det_num_alpha][walker.num][alpha_num][alpha_num]~ | Inverse of the α electron slater matrix for each determinant of each walker. |
| ~det_inv_matrix_alpha_date~ | ~uint64_t~ | Date for the Inverse of the α electron slater matrix for each determinant of each walker. |
| ~det_inv_matrix_beta~ | ~[det_num_beta][walker.num][beta_num][beta_num]~ | Inverse of the β electron slater matrix for each determinant of each walker. |
| ~det_inv_matrix_beta_date~ | ~uint64_t~ | Date for the Inverse of the β electron slater matrix for each determinant of each walker. |
|-----------------------------+--------------------------------------------------------+-------------------------------------------------------------------------------------------|
** Data structure
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_determinant_struct {
char type;
int64_t walk_num;
char type;
int64_t det_num_alpha;
int64_t det_num_beta ;
int64_t up_num;
@ -150,14 +148,14 @@ typedef struct qmckl_determinant_struct {
double * det_vgl_beta;
double * det_adj_matrix_beta;
double * det_inv_matrix_beta;
int64_t det_value_alpha_date;
int64_t det_vgl_alpha_date;
int64_t det_adj_matrix_alpha_date;
int64_t det_inv_matrix_alpha_date;
int64_t det_value_beta_date;
int64_t det_vgl_beta_date;
int64_t det_adj_matrix_beta_date;
int64_t det_inv_matrix_beta_date;
uint64_t det_value_alpha_date;
uint64_t det_vgl_alpha_date;
uint64_t det_adj_matrix_alpha_date;
uint64_t det_inv_matrix_alpha_date;
uint64_t det_value_beta_date;
uint64_t det_vgl_beta_date;
uint64_t det_adj_matrix_beta_date;
uint64_t det_inv_matrix_beta_date;
int32_t uninitialized;
bool provided;
@ -185,7 +183,7 @@ qmckl_exit_code qmckl_init_determinant(qmckl_context context) {
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
ctx->det.uninitialized = (1 << 6) - 1;
ctx->det.uninitialized = (1 << 5) - 1;
return QMCKL_SUCCESS;
}
@ -195,7 +193,6 @@ qmckl_exit_code qmckl_init_determinant(qmckl_context context) {
#+begin_src c :comments org :tangle (eval h_private_func) :exports none
char qmckl_get_determinant_type (const qmckl_context context);
int64_t qmckl_get_determinant_walk_num (const qmckl_context context);
int64_t qmckl_get_determinant_det_num_alpha (const qmckl_context context);
int64_t qmckl_get_determinant_det_num_beta (const qmckl_context context);
int64_t* qmckl_get_determinant_mo_index_alpha (const qmckl_context context);
@ -251,24 +248,6 @@ char qmckl_get_determinant_type (const qmckl_context context) {
return ctx->det.type;
}
int64_t qmckl_get_determinant_walk_num (const qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (int64_t) 0;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 1;
if ( (ctx->det.uninitialized & mask) != 0) {
return (int64_t) 0;
}
assert (ctx->det.walk_num > (int64_t) 0);
return ctx->det.walk_num;
}
int64_t qmckl_get_determinant_det_num_alpha (const qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (int64_t) 0;
@ -277,7 +256,7 @@ int64_t qmckl_get_determinant_det_num_alpha (const qmckl_context context) {
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 2;
int32_t mask = 1 << 1;
if ( (ctx->det.uninitialized & mask) != 0) {
return (int64_t) 0;
@ -295,7 +274,7 @@ int64_t qmckl_get_determinant_det_num_beta (const qmckl_context context) {
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 3;
int32_t mask = 1 << 2;
if ( (ctx->det.uninitialized & mask) != 0) {
return (int64_t) 0;
@ -313,7 +292,7 @@ int64_t* qmckl_get_determinant_mo_index_alpha (const qmckl_context context) {
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 4;
int32_t mask = 1 << 3;
if ( (ctx->det.uninitialized & mask) != 0) {
return NULL;
@ -331,7 +310,7 @@ int64_t* qmckl_get_determinant_mo_index_beta (const qmckl_context context) {
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 5;
int32_t mask = 1 << 4;
if ( (ctx->det.uninitialized & mask) != 0) {
return NULL;
@ -350,7 +329,6 @@ int64_t* qmckl_get_determinant_mo_index_beta (const qmckl_context context) {
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_set_determinant_type (const qmckl_context context, const char t);
qmckl_exit_code qmckl_set_determinant_walk_num (const qmckl_context context, const int64_t walk_num);
qmckl_exit_code qmckl_set_determinant_det_num_alpha (const qmckl_context context, const int64_t det_num_alpha);
qmckl_exit_code qmckl_set_determinant_det_num_beta (const qmckl_context context, const int64_t det_num_beta);
qmckl_exit_code qmckl_set_determinant_mo_index_alpha (const qmckl_context context, const int64_t* mo_index_alpha);
@ -395,22 +373,6 @@ qmckl_exit_code qmckl_set_determinant_type(qmckl_context context, const char t)
<<post2>>
}
qmckl_exit_code qmckl_set_determinant_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_determinant_walk_num",
"walk_num <= 0");
}
int32_t mask = 1 << 1;
ctx->det.walk_num = walk_num;
<<post2>>
}
qmckl_exit_code qmckl_set_determinant_det_num_alpha(qmckl_context context, const int64_t det_num_alpha) {
<<pre2>>
@ -421,7 +383,7 @@ qmckl_exit_code qmckl_set_determinant_det_num_alpha(qmckl_context context, const
"det_num_alpha <= 0");
}
int32_t mask = 1 << 2;
int32_t mask = 1 << 1;
ctx->det.det_num_alpha = det_num_alpha;
<<post2>>
@ -437,7 +399,7 @@ qmckl_exit_code qmckl_set_determinant_det_num_beta(qmckl_context context, const
"det_num_beta <= 0");
}
int32_t mask = 1 << 3;
int32_t mask = 1 << 2;
ctx->det.det_num_beta = det_num_beta;
<<post2>>
@ -446,7 +408,7 @@ qmckl_exit_code qmckl_set_determinant_det_num_beta(qmckl_context context, const
qmckl_exit_code qmckl_set_determinant_mo_index_alpha(qmckl_context context, const int64_t* mo_index_alpha) {
<<pre2>>
int32_t mask = 1 << 4;
int32_t mask = 1 << 3;
if (ctx->det.mo_index_alpha != NULL) {
qmckl_exit_code rc = qmckl_free(context, ctx->det.mo_index_alpha);
@ -458,7 +420,7 @@ qmckl_exit_code qmckl_set_determinant_mo_index_alpha(qmckl_context context, con
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->det.walk_num * ctx->det.det_num_alpha *
mem_info.size = ctx->electron.walker.num * ctx->det.det_num_alpha *
ctx->electron.up_num * sizeof(int64_t);
int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);
if (new_array == NULL) {
@ -478,7 +440,7 @@ qmckl_exit_code qmckl_set_determinant_mo_index_alpha(qmckl_context context, con
qmckl_exit_code qmckl_set_determinant_mo_index_beta(qmckl_context context, const int64_t* mo_index_beta) {
<<pre2>>
int32_t mask = 1 << 5;
int32_t mask = 1 << 4;
if (ctx->det.mo_index_beta != NULL) {
qmckl_exit_code rc = qmckl_free(context, ctx->det.mo_index_beta);
@ -490,7 +452,7 @@ qmckl_exit_code qmckl_set_determinant_mo_index_beta(qmckl_context context, cons
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->det.walk_num * ctx->det.det_num_beta *
mem_info.size = ctx->electron.walker.num * ctx->det.det_num_beta *
ctx->electron.down_num * sizeof(int64_t);
int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);
if (new_array == NULL) {
@ -587,10 +549,10 @@ qmckl_exit_code qmckl_get_det_vgl_alpha(qmckl_context context, double * const de
qmckl_exit_code rc;
rc = qmckl_provide_ao_vgl(context);
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_vgl(context);
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_det_vgl_alpha(context);
@ -599,7 +561,7 @@ qmckl_exit_code qmckl_get_det_vgl_alpha(qmckl_context context, double * const de
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = 5 * ctx->det.det_num_alpha * ctx->det.walk_num *
size_t sze = 5 * ctx->det.det_num_alpha * ctx->electron.walker.num *
ctx->electron.up_num * ctx->electron.up_num * sizeof(double);
memcpy(det_vgl_alpha, ctx->det.det_vgl_alpha, sze);
@ -614,10 +576,10 @@ qmckl_exit_code qmckl_get_det_vgl_beta(qmckl_context context, double * const det
qmckl_exit_code rc;
rc = qmckl_provide_ao_vgl(context);
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_vgl(context);
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_det_vgl_beta(context);
@ -626,7 +588,7 @@ qmckl_exit_code qmckl_get_det_vgl_beta(qmckl_context context, double * const det
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = 5 * ctx->det.det_num_beta * ctx->det.walk_num *
size_t sze = 5 * ctx->det.det_num_beta * ctx->electron.walker.num *
ctx->electron.down_num * ctx->electron.down_num * sizeof(double);
memcpy(det_vgl_beta, ctx->det.det_vgl_beta, sze);
@ -679,7 +641,7 @@ qmckl_exit_code qmckl_provide_det_vgl_alpha(qmckl_context context) {
"qmckl_mo_basis",
NULL);
}
rc = qmckl_provide_mo_vgl(context);
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
@ -696,13 +658,18 @@ qmckl_exit_code qmckl_provide_det_vgl_alpha(qmckl_context context) {
}
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->det.det_vgl_alpha_date) {
if (ctx->electron.walker.point.date > ctx->det.det_vgl_alpha_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->det.det_vgl_alpha);
ctx->det.det_vgl_alpha = NULL;
}
/* Allocate array */
if (ctx->det.det_vgl_alpha == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 5 * ctx->det.walk_num * ctx->det.det_num_alpha *
mem_info.size = 5 * ctx->electron.walker.num * ctx->det.det_num_alpha *
ctx->electron.up_num * ctx->electron.up_num * sizeof(double);
double* det_vgl_alpha = (double*) qmckl_malloc(context, mem_info);
@ -718,7 +685,7 @@ qmckl_exit_code qmckl_provide_det_vgl_alpha(qmckl_context context) {
if (ctx->det.type == 'G') {
rc = qmckl_compute_det_vgl_alpha(context,
ctx->det.det_num_alpha,
ctx->det.walk_num,
ctx->electron.walker.num,
ctx->electron.up_num,
ctx->electron.down_num,
ctx->electron.num,
@ -787,13 +754,18 @@ qmckl_exit_code qmckl_provide_det_vgl_beta(qmckl_context context) {
}
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->det.det_vgl_beta_date) {
if (ctx->electron.walker.point.date > ctx->det.det_vgl_beta_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->det.det_vgl_beta);
ctx->det.det_vgl_beta = NULL;
}
/* Allocate array */
if (ctx->det.det_vgl_beta == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 5 * ctx->det.walk_num * ctx->det.det_num_beta *
mem_info.size = 5 * ctx->electron.walker.num * ctx->det.det_num_beta *
ctx->electron.down_num * ctx->electron.down_num * sizeof(double);
double* det_vgl_beta = (double*) qmckl_malloc(context, mem_info);
@ -810,7 +782,7 @@ qmckl_exit_code qmckl_provide_det_vgl_beta(qmckl_context context) {
if (ctx->det.type == 'G') {
rc = qmckl_compute_det_vgl_beta(context,
ctx->det.det_num_beta,
ctx->det.walk_num,
ctx->electron.walker.num,
ctx->electron.up_num,
ctx->electron.down_num,
ctx->electron.num,
@ -1141,12 +1113,9 @@ const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
rc = qmckl_set_electron_num (context, chbrclf_elec_up_num, chbrclf_elec_dn_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_electron_walk_num (context, chbrclf_walk_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord, chbrclf_walk_num*chbrclf_elec_num*3);
rc = qmckl_set_electron_coord (context, 'N', chbrclf_walk_num, elec_coord, chbrclf_walk_num*chbrclf_elec_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num);
@ -1271,9 +1240,6 @@ for(k = 0; k < det_num_beta; ++k)
rc = qmckl_set_determinant_type (context, typ);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_determinant_walk_num (context, chbrclf_walk_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_determinant_det_num_alpha (context, det_num_alpha);
assert (rc == QMCKL_SUCCESS);
@ -1326,10 +1292,10 @@ qmckl_exit_code qmckl_get_det_inv_matrix_alpha(qmckl_context context, double * c
qmckl_exit_code rc;
rc = qmckl_provide_ao_vgl(context);
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_vgl(context);
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_det_vgl_alpha(context);
@ -1341,7 +1307,7 @@ qmckl_exit_code qmckl_get_det_inv_matrix_alpha(qmckl_context context, double * c
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.up_num * ctx->electron.up_num;
size_t sze = ctx->det.det_num_alpha * ctx->electron.walker.num * ctx->electron.up_num * ctx->electron.up_num;
memcpy(det_inv_matrix_alpha, ctx->det.det_inv_matrix_alpha, sze * sizeof(double));
return QMCKL_SUCCESS;
@ -1355,10 +1321,10 @@ qmckl_exit_code qmckl_get_det_inv_matrix_beta(qmckl_context context, double * co
qmckl_exit_code rc;
rc = qmckl_provide_ao_vgl(context);
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_vgl(context);
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_det_vgl_beta(context);
@ -1370,7 +1336,7 @@ qmckl_exit_code qmckl_get_det_inv_matrix_beta(qmckl_context context, double * co
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.down_num * ctx->electron.down_num;
size_t sze = ctx->det.det_num_alpha * ctx->electron.walker.num * ctx->electron.down_num * ctx->electron.down_num;
memcpy(det_inv_matrix_beta, ctx->det.det_inv_matrix_beta, sze * sizeof(double));
return QMCKL_SUCCESS;
@ -1384,10 +1350,10 @@ qmckl_exit_code qmckl_get_det_adj_matrix_alpha(qmckl_context context, double * c
qmckl_exit_code rc;
rc = qmckl_provide_ao_vgl(context);
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_vgl(context);
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_det_vgl_alpha(context);
@ -1399,7 +1365,7 @@ qmckl_exit_code qmckl_get_det_adj_matrix_alpha(qmckl_context context, double * c
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.up_num * ctx->electron.up_num;
size_t sze = ctx->det.det_num_alpha * ctx->electron.walker.num * ctx->electron.up_num * ctx->electron.up_num;
memcpy(det_adj_matrix_alpha, ctx->det.det_adj_matrix_alpha, sze * sizeof(double));
return QMCKL_SUCCESS;
@ -1413,10 +1379,10 @@ qmckl_exit_code qmckl_get_det_adj_matrix_beta(qmckl_context context, double * co
qmckl_exit_code rc;
rc = qmckl_provide_ao_vgl(context);
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_vgl(context);
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_det_vgl_beta(context);
@ -1428,7 +1394,7 @@ qmckl_exit_code qmckl_get_det_adj_matrix_beta(qmckl_context context, double * co
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.down_num * ctx->electron.down_num;
size_t sze = ctx->det.det_num_alpha * ctx->electron.walker.num * ctx->electron.down_num * ctx->electron.down_num;
memcpy(det_adj_matrix_beta, ctx->det.det_adj_matrix_beta, sze * sizeof(double));
return QMCKL_SUCCESS;
@ -1442,10 +1408,10 @@ qmckl_exit_code qmckl_get_det_alpha(qmckl_context context, double * const det_va
qmckl_exit_code rc;
rc = qmckl_provide_ao_vgl(context);
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_vgl(context);
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_det_vgl_alpha(context);
@ -1457,7 +1423,7 @@ qmckl_exit_code qmckl_get_det_alpha(qmckl_context context, double * const det_va
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num;
size_t sze = ctx->det.det_num_alpha * ctx->electron.walker.num;
memcpy(det_value_alpha, ctx->det.det_value_alpha, sze * sizeof(double));
return QMCKL_SUCCESS;
@ -1471,10 +1437,10 @@ qmckl_exit_code qmckl_get_det_beta(qmckl_context context, double * const det_val
qmckl_exit_code rc;
rc = qmckl_provide_ao_vgl(context);
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_vgl(context);
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_det_vgl_beta(context);
@ -1486,7 +1452,7 @@ qmckl_exit_code qmckl_get_det_beta(qmckl_context context, double * const det_val
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num;
size_t sze = ctx->det.det_num_alpha * ctx->electron.walker.num;
memcpy(det_value_beta, ctx->det.det_value_beta, sze * sizeof(double));
return QMCKL_SUCCESS;
@ -1547,13 +1513,18 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_alpha(qmckl_context context) {
}
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->det.det_inv_matrix_alpha_date) {
if (ctx->electron.walker.point.date > ctx->det.det_inv_matrix_alpha_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->det.det_inv_matrix_alpha);
ctx->det.det_inv_matrix_alpha = NULL;
}
/* Allocate array */
if (ctx->det.det_inv_matrix_alpha == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->det.walk_num * ctx->det.det_num_alpha *
mem_info.size = ctx->electron.walker.num * ctx->det.det_num_alpha *
ctx->electron.up_num * ctx->electron.up_num * sizeof(double);
double* det_inv_matrix_alpha = (double*) qmckl_malloc(context, mem_info);
@ -1569,7 +1540,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_alpha(qmckl_context context) {
if (ctx->det.det_adj_matrix_alpha == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->det.walk_num * ctx->det.det_num_alpha *
mem_info.size = ctx->electron.walker.num * ctx->det.det_num_alpha *
ctx->electron.up_num * ctx->electron.up_num * sizeof(double);
double* det_adj_matrix_alpha = (double*) qmckl_malloc(context, mem_info);
@ -1585,7 +1556,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_alpha(qmckl_context context) {
if (ctx->det.det_value_alpha == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->det.walk_num * ctx->det.det_num_alpha * sizeof(double);
mem_info.size = ctx->electron.walker.num * ctx->det.det_num_alpha * sizeof(double);
double* det_value_alpha = (double*) qmckl_malloc(context, mem_info);
if (det_value_alpha == NULL) {
@ -1601,7 +1572,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_alpha(qmckl_context context) {
if (ctx->det.type == 'G') {
rc = qmckl_compute_det_inv_matrix_alpha(context,
ctx->det.det_num_alpha,
ctx->det.walk_num,
ctx->electron.walker.num,
ctx->electron.up_num,
ctx->det.det_vgl_alpha,
ctx->det.det_value_alpha,
@ -1670,13 +1641,18 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_beta(qmckl_context context) {
}
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->det.det_inv_matrix_beta_date) {
if (ctx->electron.walker.point.date > ctx->det.det_inv_matrix_beta_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->det.det_inv_matrix_beta);
ctx->det.det_inv_matrix_beta = NULL;
}
/* Allocate array */
if (ctx->det.det_inv_matrix_beta == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->det.walk_num * ctx->det.det_num_beta *
mem_info.size = ctx->electron.walker.num * ctx->det.det_num_beta *
ctx->electron.down_num * ctx->electron.down_num * sizeof(double);
double* det_inv_matrix_beta = (double*) qmckl_malloc(context, mem_info);
@ -1692,7 +1668,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_beta(qmckl_context context) {
if (ctx->det.det_adj_matrix_beta == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->det.walk_num * ctx->det.det_num_beta *
mem_info.size = ctx->electron.walker.num * ctx->det.det_num_beta *
ctx->electron.down_num * ctx->electron.down_num * sizeof(double);
double* det_adj_matrix_beta = (double*) qmckl_malloc(context, mem_info);
@ -1708,7 +1684,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_beta(qmckl_context context) {
if (ctx->det.det_value_beta == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->det.walk_num * ctx->det.det_num_beta * sizeof(double);
mem_info.size = ctx->electron.walker.num * ctx->det.det_num_beta * sizeof(double);
double* det_value_beta = (double*) qmckl_malloc(context, mem_info);
if (det_value_beta == NULL) {
@ -1724,7 +1700,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_beta(qmckl_context context) {
if (ctx->det.type == 'G') {
rc = qmckl_compute_det_inv_matrix_beta(context,
ctx->det.det_num_beta,
ctx->det.walk_num,
ctx->electron.walker.num,
ctx->electron.down_num,
ctx->det.det_vgl_beta,
ctx->det.det_value_beta,

479
org/qmckl_electron.org
View File

@ -6,6 +6,15 @@ 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.
The electrons are stored as points in the following order: for each walker,
first up-spin electrons and then down-spin electrons.
If the points are set with the ='N'= flag, the order of
the components is =[ (x,y,z), (x,y,z), ... ]=
Using the ='T'= flage, it is =[ (x,x,x,...), (y,y,y,...), (z,z,z,...) ]=
# TODO: replace the qmckl_matrix by qmckl_point data structures.
* Headers :noexport:
#+begin_src elisp :noexport :results none
(org-babel-lob-ingest "../tools/lib.org")
@ -16,6 +25,7 @@ up-spin and down-spin electrons, and the electron coordinates.
#ifndef QMCKL_ELECTRON_HPT
#define QMCKL_ELECTRON_HPT
#include <stdbool.h>
#include "qmckl_point_private_type.h"
#+end_src
#+begin_src c :tangle (eval h_private_func)
@ -69,66 +79,68 @@ int main() {
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 |
| 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 |
| ~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 |
| ~walker~ | ~qmckl_point~ | Current set of walkers |
| ~walker_old~ | ~qmckl_point~ | Previous set of walkers |
Computed data:
| Variable | Type | Description |
|-------------------------------------+--------------------------------------+----------------------------------------------------------------------|
| ~ee_distance~ | ~double[walk_num][num][num]~ | Electron-electron distances |
| ~ee_distance~ | ~double[walker.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~ | ~double[walker.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~ | ~double[walker.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~ | ~double[walker.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 |
| ~ee_potential~ | ~double[walker.num]~ | Electron-electron potential energy |
| ~ee_potential_date~ | ~uint64_t~ | Last modification date of the electron-electron potential |
| ~en_potential~ | ~double[walker.num]~ | Electron-nucleus potential energy |
| ~en_potential_date~ | ~int64_t~ | Date when the electron-nucleus potential energy was computed |
| ~en_distance_rescaled~ | ~double[walker.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~ | ~double[walker.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_walker_struct {
int64_t num;
qmckl_point_struct point;
} qmckl_walker;
typedef struct qmckl_electron_struct {
int64_t num;
int64_t up_num;
int64_t down_num;
int64_t walk_num;
qmckl_walker walker;
qmckl_walker walker_old;
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;
uint64_t ee_distance_date;
uint64_t en_distance_date;
uint64_t ee_potential_date;
uint64_t en_potential_date;
uint64_t ee_distance_rescaled_date;
uint64_t ee_distance_rescaled_deriv_e_date;
uint64_t en_distance_rescaled_date;
uint64_t en_distance_rescaled_deriv_e_date;
double* ee_distance;
double* en_distance;
double* ee_pot;
double* en_pot;
double* ee_potential;
double* en_potential;
double* ee_distance_rescaled;
double* ee_distance_rescaled_deriv_e;
double* en_distance_rescaled;
@ -160,7 +172,7 @@ qmckl_exit_code qmckl_init_electron(qmckl_context context) {
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
ctx->electron.uninitialized = (1 << 2) - 1;
ctx->electron.uninitialized = (1 << 1) - 1;
/* Default values */
ctx->electron.rescale_factor_kappa_ee = 1.0;
@ -326,14 +338,7 @@ qmckl_get_electron_walk_num (const qmckl_context context, int64_t* const walk_nu
qmckl_context_struct* const ctx = (qmckl_context_struct*) 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;
,*walk_num = ctx->electron.walker.num;
return QMCKL_SUCCESS;
}
#+end_src
@ -395,13 +400,13 @@ qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const
*** 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~.
to point on a memory block of size ~size_max~ \ge ~3 * elec_num * walker.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)~ |
| | Normal | Transposed |
|---------+----------------------------+----------------------------|
| C | ~[walker.num*elec_num][3]~ | ~[3][walker.num*elec_num]~ |
| Fortran | ~(3,walker.num*elec_num)~ | ~(walker.num*elec_num, 3)~ |
#+begin_src c :comments org :tangle (eval h_func) :exports none
@ -412,9 +417,9 @@ qmckl_get_electron_coord (const qmckl_context context,
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.
As the ~walker~ attribute is 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
@ -458,9 +463,8 @@ qmckl_get_electron_coord (const qmckl_context context,
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]);
assert (ctx->point.num == ctx->electron.walker.point.num);
assert (ctx->point.coord.data == ctx->electron.walker.point.coord.data);
return qmckl_get_point(context, transp, coord, size_max);
}
@ -472,12 +476,11 @@ qmckl_get_electron_coord (const qmckl_context context,
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.
both not 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_coord (qmckl_context context, const char transp, const int64_t walk_num, 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);
@ -504,42 +507,6 @@ if (mask != 0 && !(ctx->electron.uninitialized & mask)) {
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
@ -577,28 +544,6 @@ qmckl_set_electron_num(qmckl_context context,
}
#+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) {
int32_t mask = 1 << 1;
<<pre2>>
if (walk_num <= 0) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_set_electron_walk_num",
"walk_num <= 0");
}
ctx->electron.walk_num = walk_num;
<<post2>>
}
#+end_src
Next we set the rescale parameter for the rescaled distance metric.
@ -656,27 +601,17 @@ interface
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~.
~size_max~ should be equal equal or geater than ~elec_num *
walker.num * 3~, to be symmetric with ~qmckl_get_electron_coord~.
Important: changing the electron coordinates increments the date
in the context.
@ -685,6 +620,7 @@ end interface
qmckl_exit_code
qmckl_set_electron_coord(qmckl_context context,
const char transp,
const int64_t walk_num,
const double* coord,
const int64_t size_max)
{
@ -700,6 +636,13 @@ qmckl_set_electron_coord(qmckl_context context,
"transp should be 'N' or 'T'");
}
if (walk_num <= 0) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_set_electron_coord",
"walk_num <= 0");
}
if (coord == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
@ -716,30 +659,19 @@ qmckl_set_electron_coord(qmckl_context context,
"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_walker tmp = ctx->electron.walker_old;
ctx->electron.walker_old = ctx->electron.walker;
ctx->electron.walker = tmp;
memcpy(&(ctx->point), &(ctx->electron.walker.point), sizeof(qmckl_point_struct));
qmckl_exit_code rc;
rc = qmckl_set_point(context, transp, size_max/3, coord, size_max);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_point(context, transp, walk_num*elec_num, coord, size_max);
if (rc != QMCKL_SUCCESS) return rc;
ctx->electron.coord_new = ctx->point.coord ;
ctx->electron.coord_new_date = ctx->date;
ctx->electron.walker.num = walk_num;
memcpy(&(ctx->electron.walker.point), &(ctx->point), sizeof(qmckl_point_struct));
return QMCKL_SUCCESS;
@ -748,13 +680,14 @@ qmckl_set_electron_coord(qmckl_context context,
#+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)
integer(c_int32_t) function qmckl_set_electron_coord(context, transp, walk_num, 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
integer (c_int64_t) , intent(in) , value :: walk_num
double precision , intent(in) :: coord(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
@ -802,7 +735,7 @@ 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));
assert(qmckl_electron_provided(context));
rc = qmckl_get_electron_up_num (context, &n);
assert(rc == QMCKL_SUCCESS);
@ -841,23 +774,21 @@ assert(rc == QMCKL_SUCCESS);
assert(k_en == rescale_factor_kappa_en);
int64_t w;
int64_t w = 0;
rc = qmckl_get_electron_walk_num (context, &w);
assert(rc == QMCKL_NOT_PROVIDED);
assert(rc == QMCKL_SUCCESS);
assert(w == 0);
rc = qmckl_set_electron_walk_num (context, walk_num);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', walk_num, elec_coord, walk_num*elec_num*3);
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);
@ -916,7 +847,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance(qmckl_context context, double* co
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num;
size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walker.num;
memcpy(distance, ctx->electron.ee_distance, sze * sizeof(double));
return QMCKL_SUCCESS;
@ -942,14 +873,19 @@ qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context)
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->electron.ee_distance_date) {
if (ctx->electron.walker.point.date > ctx->electron.ee_distance_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->electron.ee_distance);
ctx->electron.ee_distance = NULL;
}
/* 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);
ctx->electron.walker.num * sizeof(double);
double* ee_distance = (double*) qmckl_malloc(context, mem_info);
if (ee_distance == NULL) {
@ -964,8 +900,8 @@ qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context)
qmckl_exit_code rc =
qmckl_compute_ee_distance(context,
ctx->electron.num,
ctx->electron.walk_num,
ctx->electron.coord_new.data,
ctx->electron.walker.num,
ctx->electron.walker.point.coord.data,
ctx->electron.ee_distance);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -1157,7 +1093,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance_rescaled(qmckl_context context, d
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num;
size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walker.num;
memcpy(distance_rescaled, ctx->electron.ee_distance_rescaled, sze * sizeof(double));
return QMCKL_SUCCESS;
@ -1183,14 +1119,19 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context)
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->electron.ee_distance_rescaled_date) {
if (ctx->electron.walker.point.date > ctx->electron.ee_distance_rescaled_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->electron.ee_distance_rescaled);
ctx->electron.ee_distance_rescaled = NULL;
}
/* 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);
ctx->electron.walker.num * sizeof(double);
double* ee_distance_rescaled = (double*) qmckl_malloc(context, mem_info);
if (ee_distance_rescaled == NULL) {
@ -1206,8 +1147,8 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context)
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.walker.num,
ctx->electron.walker.point.coord.data,
ctx->electron.ee_distance_rescaled);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -1403,7 +1344,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance_rescaled_deriv_e(qmckl_context co
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = 4 * ctx->electron.num * ctx->electron.num * ctx->electron.walk_num;
size_t sze = 4 * ctx->electron.num * ctx->electron.num * ctx->electron.walker.num;
memcpy(distance_rescaled_deriv_e, ctx->electron.ee_distance_rescaled_deriv_e, sze * sizeof(double));
return QMCKL_SUCCESS;
@ -1429,14 +1370,19 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->electron.ee_distance_rescaled_deriv_e_date) {
if (ctx->electron.walker.point.date > ctx->electron.ee_distance_rescaled_deriv_e_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->electron.ee_distance_rescaled_deriv_e);
ctx->electron.ee_distance_rescaled_deriv_e = NULL;
}
/* 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);
ctx->electron.walker.num * sizeof(double);
double* ee_distance_rescaled_deriv_e = (double*) qmckl_malloc(context, mem_info);
if (ee_distance_rescaled_deriv_e == NULL) {
@ -1452,8 +1398,8 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context
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.walker.num,
ctx->electron.walker.point.coord.data,
ctx->electron.ee_distance_rescaled_deriv_e);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -1602,7 +1548,7 @@ rc = qmckl_get_electron_ee_distance_rescaled_deriv_e(context, ee_distance_rescal
** Electron-electron potential
~ee_pot~ calculates the ~ee~ potential energy.
~ee_potential~ is given by
\[
\mathcal{V}_{ee} = \sum_{i=1}^{N_e}\sum_{j>i}^{N_e}\frac{1}{r_{ij}}
@ -1614,11 +1560,11 @@ rc = qmckl_get_electron_ee_distance_rescaled_deriv_e(context, ee_distance_rescal
*** 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);
qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* const ee_potential);
#+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)
qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* const ee_potential)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -1632,8 +1578,8 @@ qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* c
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.walk_num * sizeof(double);
memcpy(ee_pot, ctx->electron.ee_pot, sze);
size_t sze = ctx->electron.walker.num * sizeof(double);
memcpy(ee_potential, ctx->electron.ee_potential, sze);
return QMCKL_SUCCESS;
}
@ -1662,34 +1608,39 @@ qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context)
if (rc != QMCKL_SUCCESS) return rc;
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->electron.ee_pot_date) {
if (ctx->electron.walker.point.date > ctx->electron.ee_potential_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->electron.ee_potential);
ctx->electron.ee_distance_rescaled_deriv_e = NULL;
}
/* Allocate array */
if (ctx->electron.ee_pot == NULL) {
if (ctx->electron.ee_potential == 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);
mem_info.size = ctx->electron.walker.num * sizeof(double);
double* ee_potential = (double*) qmckl_malloc(context, mem_info);
if (ee_pot == NULL) {
if (ee_potential == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_ee_potential",
NULL);
}
ctx->electron.ee_pot = ee_pot;
ctx->electron.ee_potential = ee_potential;
}
rc = qmckl_compute_ee_potential(context,
ctx->electron.num,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->electron.ee_distance,
ctx->electron.ee_pot);
ctx->electron.ee_potential);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->electron.ee_pot_date = ctx->date;
ctx->electron.ee_potential_date = ctx->date;
}
return QMCKL_SUCCESS;
@ -1704,17 +1655,17 @@ qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context)
: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 |
| 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_potential~ | ~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) &
ee_distance, ee_potential) &
result(info)
use qmckl
implicit none
@ -1722,7 +1673,7 @@ integer function qmckl_compute_ee_potential_f(context, elec_num, walk_num, &
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)
double precision , intent(out) :: ee_potential(walk_num)
integer*8 :: nw, i, j
@ -1743,12 +1694,12 @@ integer function qmckl_compute_ee_potential_f(context, elec_num, walk_num, &
return
endif
ee_pot = 0.0d0
ee_potential = 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))
ee_potential(nw) = ee_potential(nw) + 1.0d0/(ee_distance(i,j,nw))
endif
end do
end do
@ -1766,7 +1717,7 @@ end function qmckl_compute_ee_potential_f
const int64_t elec_num,
const int64_t walk_num,
const double* ee_distance,
double* const ee_pot );
double* const ee_potential );
#+end_src
#+CALL: generate_c_interface(table=qmckl_ee_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
@ -1774,7 +1725,7 @@ end function qmckl_compute_ee_potential_f
#+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) &
(context, elec_num, walk_num, ee_distance, ee_potential) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
@ -1784,20 +1735,20 @@ end function qmckl_compute_ee_potential_f
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)
real (c_double ) , intent(out) :: ee_potential(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)
(context, elec_num, walk_num, ee_distance, ee_potential)
end function qmckl_compute_ee_potential
#+end_src
*** Test
#+begin_src c :tangle (eval c_test)
double ee_pot[walk_num];
double ee_potential[walk_num];
rc = qmckl_get_electron_ee_potential(context, &(ee_pot[0]));
rc = qmckl_get_electron_ee_potential(context, &(ee_potential[0]));
assert (rc == QMCKL_SUCCESS);
#+end_src
** Electron-nucleus distances
@ -1837,7 +1788,7 @@ qmckl_exit_code qmckl_get_electron_en_distance(qmckl_context context, double* di
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num;
size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num;
memcpy(distance, ctx->electron.en_distance, sze * sizeof(double));
return QMCKL_SUCCESS;
@ -1869,14 +1820,19 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)
}
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->electron.en_distance_date) {
if (ctx->electron.walker.point.date > ctx->electron.en_distance_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->electron.en_distance);
ctx->electron.en_distance = NULL;
}
/* 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);
ctx->electron.walker.num * sizeof(double);
double* en_distance = (double*) qmckl_malloc(context, mem_info);
if (en_distance == NULL) {
@ -1892,8 +1848,8 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)
qmckl_compute_en_distance(context,
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.walk_num,
ctx->electron.coord_new.data,
ctx->electron.walker.num,
ctx->electron.walker.point.coord.data,
ctx->nucleus.coord.data,
ctx->electron.en_distance);
if (rc != QMCKL_SUCCESS) {
@ -2116,7 +2072,7 @@ qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, d
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num;
size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num;
memcpy(distance_rescaled, ctx->electron.en_distance_rescaled, sze * sizeof(double));
return QMCKL_SUCCESS;
@ -2145,14 +2101,19 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context)
}
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->electron.en_distance_rescaled_date) {
if (ctx->electron.walker.point.date > ctx->electron.en_distance_rescaled_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->electron.en_distance_rescaled);
ctx->electron.en_distance_rescaled = NULL;
}
/* 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);
ctx->electron.walker.num * sizeof(double);
double* en_distance_rescaled = (double*) qmckl_malloc(context, mem_info);
if (en_distance_rescaled == NULL) {
@ -2169,8 +2130,8 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context)
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.rescale_factor_kappa_en,
ctx->electron.walk_num,
ctx->electron.coord_new.data,
ctx->electron.walker.num,
ctx->electron.walker.point.coord.data,
ctx->nucleus.coord.data,
ctx->electron.en_distance_rescaled);
if (rc != QMCKL_SUCCESS) {
@ -2394,7 +2355,7 @@ qmckl_exit_code qmckl_get_electron_en_distance_rescaled_deriv_e(qmckl_context co
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num;
size_t sze = 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num;
memcpy(distance_rescaled_deriv_e, ctx->electron.en_distance_rescaled_deriv_e, sze * sizeof(double));
return QMCKL_SUCCESS;
@ -2423,14 +2384,19 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context
}
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->electron.en_distance_rescaled_deriv_e_date) {
if (ctx->electron.walker.point.date > ctx->electron.en_distance_rescaled_deriv_e_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->electron.en_distance_rescaled_deriv_e);
ctx->electron.en_distance_rescaled_deriv_e = NULL;
}
/* 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);
ctx->electron.walker.num * sizeof(double);
double* en_distance_rescaled_deriv_e = (double*) qmckl_malloc(context, mem_info);
if (en_distance_rescaled_deriv_e == NULL) {
@ -2447,8 +2413,8 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.rescale_factor_kappa_en,
ctx->electron.walk_num,
ctx->electron.coord_new.data,
ctx->electron.walker.num,
ctx->electron.walker.point.coord.data,
ctx->nucleus.coord.data,
ctx->electron.en_distance_rescaled_deriv_e);
if (rc != QMCKL_SUCCESS) {
@ -2638,11 +2604,11 @@ assert (rc == QMCKL_SUCCESS);
*** 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);
qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* const en_potential);
#+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)
qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* const en_potential)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
@ -2656,8 +2622,8 @@ qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* c
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.walk_num * sizeof(double);
memcpy(en_pot, ctx->electron.en_pot, sze);
size_t sze = ctx->electron.walker.num * sizeof(double);
memcpy(en_potential, ctx->electron.en_potential, sze);
return QMCKL_SUCCESS;
}
@ -2687,36 +2653,41 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context)
if (rc != QMCKL_SUCCESS) return rc;
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->electron.en_pot_date) {
if (ctx->electron.walker.point.date > ctx->electron.en_potential_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->electron.en_potential);
ctx->electron.en_potential = NULL;
}
/* Allocate array */
if (ctx->electron.en_pot == NULL) {
if (ctx->electron.en_potential == 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);
mem_info.size = ctx->electron.walker.num * sizeof(double);
double* en_potential = (double*) qmckl_malloc(context, mem_info);
if (en_pot == NULL) {
if (en_potential == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_en_potential",
NULL);
}
ctx->electron.en_pot = en_pot;
ctx->electron.en_potential = en_potential;
}
rc = qmckl_compute_en_potential(context,
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->nucleus.charge.data,
ctx->electron.en_distance,
ctx->electron.en_pot);
ctx->electron.en_potential);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->electron.en_pot_date = ctx->date;
ctx->electron.en_potential_date = ctx->date;
}
return QMCKL_SUCCESS;
@ -2731,19 +2702,19 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context)
: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 |
| 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_potential~ | ~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) &
charge, en_distance, en_potential) &
result(info)
use qmckl
implicit none
@ -2753,7 +2724,7 @@ integer function qmckl_compute_en_potential_f(context, elec_num, nucl_num, walk_
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)
double precision , intent(out) :: en_potential(walk_num)
integer*8 :: nw, i, j
@ -2774,12 +2745,12 @@ integer function qmckl_compute_en_potential_f(context, elec_num, nucl_num, walk_
return
endif
en_pot = 0.0d0
en_potential = 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))
en_potential(nw) = en_potential(nw) - charge(j)/(en_distance(i,j,nw))
endif
end do
end do
@ -2799,7 +2770,7 @@ end function qmckl_compute_en_potential_f
const int64_t walk_num,
const double* charge,
const double* en_distance,
double* const en_pot );
double* const en_potential );
#+end_src
#+CALL: generate_c_interface(table=qmckl_en_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
@ -2807,7 +2778,7 @@ end function qmckl_compute_en_potential_f
#+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) &
(context, elec_num, nucl_num, walk_num, charge, en_distance, en_potential) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
@ -2819,20 +2790,20 @@ end function qmckl_compute_en_potential_f
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)
real (c_double ) , intent(out) :: en_potential(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)
(context, elec_num, nucl_num, walk_num, charge, en_distance, en_potential)
end function qmckl_compute_en_potential
#+end_src
*** Test
#+begin_src c :tangle (eval c_test)
double en_pot[walk_num];
double en_potential[walk_num];
rc = qmckl_get_electron_en_potential(context, &(en_pot[0]));
rc = qmckl_get_electron_en_potential(context, &(en_potential[0]));
assert (rc == QMCKL_SUCCESS);
#+end_src

60
org/qmckl_error.org
View File

@ -360,6 +360,7 @@ void qmckl_string_of_error_f(const qmckl_exit_code error, char result[<<MAX_STRI
subroutine qmckl_string_of_error (error, string) bind(C, name='qmckl_string_of_error_f')
use, intrinsic :: iso_c_binding
import
implicit none
integer (qmckl_exit_code), intent(in), value :: error
character, intent(out) :: string(<<MAX_STRING_LENGTH()>>)
end subroutine qmckl_string_of_error
@ -554,6 +555,58 @@ if (x < 0) {
}
#+end_src
* Last error
Returns a string describing the last error, using ~qmckl_get_error~.
# Header
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_last_error(qmckl_context context, char* buffer);
#+end_src
# Source
#+begin_src c :tangle (eval c) :exports none
qmckl_exit_code
qmckl_last_error(qmckl_context context, char* buffer) {
char function_name[QMCKL_MAX_FUN_LEN];
char message[QMCKL_MAX_MSG_LEN];
qmckl_exit_code rc, last_rc;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
strncpy(buffer, "Null context", 13);
return QMCKL_FAILURE;
}
rc = qmckl_get_error(context, &last_rc, function_name, message);
if (rc != QMCKL_SUCCESS) {
return rc;
}
sprintf(buffer, "Error -- %s -- in %s\n%s",
qmckl_string_of_error(last_rc),
function_name, message);
return QMCKL_SUCCESS;
}
#+end_src
** Fortran inteface
#+begin_src f90 :tangle (eval fh_func) :exports none :noweb yes
interface
subroutine qmckl_last_error (context, string) bind(C, name='qmckl_last_error')
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in), value :: context
character, intent(out) :: string(*)
end subroutine qmckl_last_error
end interface
#+end_src
* End of files :noexport:
@ -577,11 +630,14 @@ if (x < 0) {
assert (strcmp(function_name,"qmckl_transpose") == 0);
assert (strcmp(message,"Success") == 0);
exit_code = qmckl_context_destroy(context);
assert(exit_code == QMCKL_SUCCESS);
return 0;
}
#+end_src
# -*- mode: org -*-
# vim: syntax=c
# -*- mode: org -*-
# vim: syntax=c

279
org/qmckl_jastrow.org
View File

@ -145,17 +145,17 @@ int main() {
| ~aord_vector~ | ~double[aord_num + 1][type_nucl_num]~ | in | Order of a polynomial coefficients |
| ~bord_vector~ | ~double[bord_num + 1]~ | in | Order of b polynomial coefficients |
| ~cord_vector~ | ~double[cord_num][type_nucl_num]~ | in | Order of c polynomial coefficients |
| ~factor_ee~ | ~double[walk_num]~ | out | Jastrow factor: electron-electron part |
| ~factor_ee~ | ~double[walker.num]~ | out | Jastrow factor: electron-electron part |
| ~factor_ee_date~ | ~uint64_t~ | out | Jastrow factor: electron-electron part |
| ~factor_en~ | ~double[walk_num]~ | out | Jastrow factor: electron-nucleus part |
| ~factor_en~ | ~double[walker.num]~ | out | Jastrow factor: electron-nucleus part |
| ~factor_en_date~ | ~uint64_t~ | out | Jastrow factor: electron-nucleus part |
| ~factor_een~ | ~double[walk_num]~ | out | Jastrow factor: electron-electron-nucleus part |
| ~factor_een~ | ~double[walker.num]~ | out | Jastrow factor: electron-electron-nucleus part |
| ~factor_een_date~ | ~uint64_t~ | out | Jastrow factor: electron-electron-nucleus part |
| ~factor_ee_deriv_e~ | ~double[4][nelec][walk_num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part |
| ~factor_ee_deriv_e~ | ~double[4][nelec][walker.num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part |
| ~factor_ee_deriv_e_date~ | ~uint64_t~ | out | Keep track of the date for the derivative |
| ~factor_en_deriv_e~ | ~double[4][nelec][walk_num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part |
| ~factor_en_deriv_e~ | ~double[4][nelec][walker.num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part |
| ~factor_en_deriv_e_date~ | ~uint64_t~ | out | Keep track of the date for the en derivative |
| ~factor_een_deriv_e~ | ~double[4][nelec][walk_num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part |
| ~factor_een_deriv_e~ | ~double[4][nelec][walker.num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part |
| ~factor_een_deriv_e_date~ | ~uint64_t~ | out | Keep track of the date for the een derivative |
computed data:
@ -170,14 +170,13 @@ int main() {
| ~cord_vect_full_date~ | ~uint64_t~ | Keep track of changes here |
| ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices |
| ~lkpm_combined_index_date~ | ~uint64_t~ | Transform l,k,p, and m into consecutive indices |
| ~tmp_c~ | ~double[walk_num][cord_num][cord_num+1][nucl_num][elec_num]~ | vector of non-zero coefficients |
| ~dtmp_c~ | ~double[walk_num][elec_num][4][nucl_num][cord_num+1][cord_num]~ | vector of non-zero coefficients |
| ~een_rescaled_n~ | ~double[walk_num][cord_num+1][nucl_num][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord |
| ~tmp_c~ | ~double[walker.num][cord_num][cord_num+1][nucl_num][elec_num]~ | vector of non-zero coefficients |
| ~dtmp_c~ | ~double[walker.num][elec_num][4][nucl_num][cord_num+1][cord_num]~ | vector of non-zero coefficients |
| ~een_rescaled_n~ | ~double[walker.num][cord_num+1][nucl_num][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord |
| ~een_rescaled_n_date~ | ~uint64_t~ | Keep track of the date of creation |
| ~een_rescaled_e_deriv_e~ | ~double[walk_num][cord_num+1][elec_num][4][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons |
| ~een_rescaled_e_deriv_e~ | ~double[walker.num][cord_num+1][elec_num][4][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons |
| ~een_rescaled_e_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation |
| ~een_rescaled_n_deriv_e~ | ~double[walk_num][cord_num+1][nucl_num][4][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons |
| ~een_rescaled_n_deriv_e~ | ~double[walker.num][cord_num+1][nucl_num][4][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons |
| ~een_rescaled_n_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation |
#+NAME: jastrow_data
@ -1166,32 +1165,10 @@ 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);
assert(qmckl_electron_provided(context));
double k_ee = 0.;
double k_en = 0.;
@ -1218,21 +1195,7 @@ 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*3*elec_num);
rc = qmckl_set_electron_coord (context, 'N', walk_num, elec_coord, walk_num*3*elec_num);
assert(rc == QMCKL_SUCCESS);
double elec_coord2[walk_num*3*elec_num];
@ -1248,18 +1211,10 @@ for (int64_t i=0 ; i<3*elec_num ; ++i) {
assert(!qmckl_nucleus_provided(context));
rc = qmckl_get_nucleus_num (context, &n);
assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_nucleus_provided(context));
rc = qmckl_get_nucleus_num (context, &n);
assert(rc == QMCKL_SUCCESS);
assert(n == nucl_num);
double k;
rc = qmckl_get_nucleus_rescale_factor (context, &k);
assert(rc == QMCKL_SUCCESS);
@ -1640,12 +1595,12 @@ qmckl_get_jastrow_factor_ee(qmckl_context context,
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze=ctx->electron.walk_num;
int64_t sze=ctx->electron.walker.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_ee",
"Array too small. Expected walk_num");
"Array too small. Expected walker.num");
}
memcpy(factor_ee, ctx->jastrow.factor_ee, sze*sizeof(double));
@ -1678,11 +1633,16 @@ qmckl_exit_code qmckl_provide_factor_ee(qmckl_context context)
/* Compute if necessary */
if (ctx->date > ctx->jastrow.factor_ee_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->jastrow.factor_ee);
ctx->jastrow.factor_ee = NULL;
}
/* Allocate array */
if (ctx->jastrow.factor_ee == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walk_num * sizeof(double);
mem_info.size = ctx->electron.walker.num * sizeof(double);
double* factor_ee = (double*) qmckl_malloc(context, mem_info);
if (factor_ee == NULL) {
@ -1695,7 +1655,7 @@ qmckl_exit_code qmckl_provide_factor_ee(qmckl_context context)
}
rc = qmckl_compute_factor_ee(context,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->electron.num,
ctx->electron.up_num,
ctx->jastrow.bord_num,
@ -1975,7 +1935,7 @@ qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context,
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num;
int64_t sze = ctx->electron.walker.num * 4 * ctx->electron.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
@ -2018,11 +1978,16 @@ qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context)
/* Compute if necessary */
if (ctx->date > ctx->jastrow.factor_ee_deriv_e_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->jastrow.factor_ee_deriv_e);
ctx->jastrow.factor_ee_deriv_e = NULL;
}
/* Allocate array */
if (ctx->jastrow.factor_ee_deriv_e == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walk_num * 4 * ctx->electron.num * sizeof(double);
mem_info.size = ctx->electron.walker.num * 4 * ctx->electron.num * sizeof(double);
double* factor_ee_deriv_e = (double*) qmckl_malloc(context, mem_info);
if (factor_ee_deriv_e == NULL) {
@ -2035,7 +2000,7 @@ qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context)
}
rc = qmckl_compute_factor_ee_deriv_e(context,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->electron.num,
ctx->electron.up_num,
ctx->jastrow.bord_num,
@ -2558,12 +2523,12 @@ qmckl_get_jastrow_factor_en(qmckl_context context,
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze=ctx->electron.walk_num;
int64_t sze=ctx->electron.walker.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_en",
"Array too small. Expected walk_num");
"Array too small. Expected walker.num");
}
memcpy(factor_en, ctx->jastrow.factor_en, sze*sizeof(double));
@ -2596,11 +2561,16 @@ qmckl_exit_code qmckl_provide_factor_en(qmckl_context context)
/* Compute if necessary */
if (ctx->date > ctx->jastrow.factor_en_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->jastrow.factor_en);
ctx->jastrow.factor_en = NULL;
}
/* Allocate array */
if (ctx->jastrow.factor_en == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walk_num * sizeof(double);
mem_info.size = ctx->electron.walker.num * sizeof(double);
double* factor_en = (double*) qmckl_malloc(context, mem_info);
if (factor_en == NULL) {
@ -2613,7 +2583,7 @@ qmckl_exit_code qmckl_provide_factor_en(qmckl_context context)
}
rc = qmckl_compute_factor_en(context,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow.type_nucl_num,
@ -2900,12 +2870,12 @@ qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context,
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num;
int64_t sze = ctx->electron.walker.num * 4 * ctx->electron.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_en_deriv_e",
"Array too small. Expected 4*walk_num*elec_num");
"Array too small. Expected 4*walker.num*elec_num");
}
memcpy(factor_en_deriv_e, ctx->jastrow.factor_en_deriv_e, sze*sizeof(double));
@ -2942,11 +2912,16 @@ qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context)
/* Compute if necessary */
if (ctx->date > ctx->jastrow.factor_en_deriv_e_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->jastrow.factor_en_deriv_e);
ctx->jastrow.factor_en_deriv_e = NULL;
}
/* Allocate array */
if (ctx->jastrow.factor_en_deriv_e == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walk_num * 4 * ctx->electron.num * sizeof(double);
mem_info.size = ctx->electron.walker.num * 4 * ctx->electron.num * sizeof(double);
double* factor_en_deriv_e = (double*) qmckl_malloc(context, mem_info);
if (factor_en_deriv_e == NULL) {
@ -2959,7 +2934,7 @@ qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context)
}
rc = qmckl_compute_factor_en_deriv_e(context,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow.type_nucl_num,
@ -3321,12 +3296,12 @@ qmckl_get_jastrow_een_rescaled_e(qmckl_context context,
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
int64_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walker.num * (ctx->jastrow.cord_num + 1);
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_een_rescaled_e",
"Array too small. Expected ctx->electron.num * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1)");
"Array too small. Expected ctx->electron.num * ctx->electron.num * ctx->electron.walker.num * (ctx->jastrow.cord_num + 1)");
}
memcpy(distance_rescaled, ctx->jastrow.een_rescaled_e, sze * sizeof(double));
@ -3358,12 +3333,17 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context)
/* Compute if necessary */
if (ctx->date > ctx->jastrow.een_rescaled_e_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->jastrow.een_rescaled_e);
ctx->jastrow.een_rescaled_e = NULL;
}
/* Allocate array */
if (ctx->jastrow.een_rescaled_e == 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 * (ctx->jastrow.cord_num + 1) * sizeof(double);
ctx->electron.walker.num * (ctx->jastrow.cord_num + 1) * sizeof(double);
double* een_rescaled_e = (double*) qmckl_malloc(context, mem_info);
if (een_rescaled_e == NULL) {
@ -3376,7 +3356,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context)
}
rc = qmckl_compute_een_rescaled_e(context,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->electron.num,
ctx->jastrow.cord_num,
ctx->electron.rescale_factor_kappa_ee,
@ -3827,12 +3807,12 @@ qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context,
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
int64_t sze = ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walker.num * (ctx->jastrow.cord_num + 1);
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_een_deriv_e",
"Array too small. Expected ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1)");
"Array too small. Expected ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walker.num * (ctx->jastrow.cord_num + 1)");
}
memcpy(distance_rescaled, ctx->jastrow.een_rescaled_e_deriv_e, sze * sizeof(double));
@ -3864,12 +3844,17 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context)
/* Compute if necessary */
if (ctx->date > ctx->jastrow.een_rescaled_e_deriv_e_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->jastrow.een_rescaled_e_deriv_e);
ctx->jastrow.een_rescaled_e_deriv_e = NULL;
}
/* Allocate array */
if (ctx->jastrow.een_rescaled_e_deriv_e == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.num * 4 * ctx->electron.num *
ctx->electron.walk_num * (ctx->jastrow.cord_num + 1) * sizeof(double);
ctx->electron.walker.num * (ctx->jastrow.cord_num + 1) * sizeof(double);
double* een_rescaled_e_deriv_e = (double*) qmckl_malloc(context, mem_info);
if (een_rescaled_e_deriv_e == NULL) {
@ -3882,11 +3867,11 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context)
}
rc = qmckl_compute_factor_een_rescaled_e_deriv_e(context,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->electron.num,
ctx->jastrow.cord_num,
ctx->electron.rescale_factor_kappa_ee,
ctx->electron.coord_new.data,
ctx->electron.walker.point.coord.data,
ctx->electron.ee_distance,
ctx->jastrow.een_rescaled_e,
ctx->jastrow.een_rescaled_e_deriv_e);
@ -3916,7 +3901,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context)
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
| ~cord_num~ | ~int64_t~ | in | Order of polynomials |
| ~rescale_factor_kappa_ee~ | ~double~ | in | Factor to rescale ee distances |
| ~coord_new~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates |
| ~coord_ee~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates |
| ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | in | Electron-electron distances |
| ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | in | Electron-electron distances |
| ~een_rescaled_e_deriv_e~ | ~double[walk_num][0:cord_num][elec_num][4][elec_num]~ | out | Electron-electron rescaled distances |
@ -3924,7 +3909,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context)
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f( &
context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, &
coord_new, ee_distance, een_rescaled_e, een_rescaled_e_deriv_e) &
coord_ee, ee_distance, een_rescaled_e, een_rescaled_e_deriv_e) &
result(info)
use qmckl
implicit none
@ -3933,7 +3918,7 @@ integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f( &
integer*8 , intent(in) :: elec_num
integer*8 , intent(in) :: cord_num
double precision , intent(in) :: rescale_factor_kappa_ee
double precision , intent(in) :: coord_new(elec_num,3,walk_num)
double precision , intent(in) :: coord_ee(elec_num,3,walk_num)
double precision , intent(in) :: ee_distance(elec_num,elec_num,walk_num)
double precision , intent(in) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
double precision , intent(out) :: een_rescaled_e_deriv_e(elec_num,4,elec_num,0:cord_num,walk_num)
@ -3972,7 +3957,7 @@ integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f( &
do i = 1, elec_num
rij_inv = 1.0d0 / ee_distance(i, j, nw)
do ii = 1, 3
elec_dist_deriv_e(ii, i, j) = (coord_new(i, ii, nw) - coord_new(j, ii, nw)) * rij_inv
elec_dist_deriv_e(ii, i, j) = (coord_ee(i, ii, nw) - coord_ee(j, ii, nw)) * rij_inv
end do
elec_dist_deriv_e(4, i, j) = 2.0d0 * rij_inv
end do
@ -4019,7 +4004,7 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f
const int64_t elec_num,
const int64_t cord_num,
const double rescale_factor_kappa_ee,
const double* coord_new,
const double* coord_ee,
const double* ee_distance,
const double* een_rescaled_e,
double* const een_rescaled_e_deriv_e );
@ -4036,7 +4021,7 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f
elec_num, &
cord_num, &
rescale_factor_kappa_ee, &
coord_new, &
coord_ee, &
ee_distance, &
een_rescaled_e, &
een_rescaled_e_deriv_e) &
@ -4050,7 +4035,7 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f
integer (c_int64_t) , intent(in) , value :: elec_num
integer (c_int64_t) , intent(in) , value :: cord_num
real (c_double ) , intent(in) , value :: rescale_factor_kappa_ee
real (c_double ) , intent(in) :: coord_new(elec_num,3,walk_num)
real (c_double ) , intent(in) :: coord_ee(elec_num,3,walk_num)
real (c_double ) , intent(in) :: ee_distance(elec_num,elec_num,walk_num)
real (c_double ) , intent(in) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
real (c_double ) , intent(out) :: een_rescaled_e_deriv_e(elec_num,4,elec_num,0:cord_num,walk_num)
@ -4062,7 +4047,7 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f
elec_num, &
cord_num, &
rescale_factor_kappa_ee, &
coord_new, &
coord_ee, &
ee_distance, &
een_rescaled_e, &
een_rescaled_e_deriv_e)
@ -4205,12 +4190,12 @@ qmckl_get_jastrow_een_rescaled_n(qmckl_context context,
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
int64_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow.cord_num + 1);
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_een_deriv_e",
"Array too small. Expected ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1)");
"Array too small. Expected ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow.cord_num + 1)");
}
memcpy(distance_rescaled, ctx->jastrow.een_rescaled_n, sze * sizeof(double));
@ -4242,12 +4227,17 @@ qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context)
/* Compute if necessary */
if (ctx->date > ctx->jastrow.een_rescaled_n_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->jastrow.een_rescaled_n);
ctx->jastrow.een_rescaled_n = NULL;
}
/* Allocate array */
if (ctx->jastrow.een_rescaled_n == 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 * (ctx->jastrow.cord_num + 1) * sizeof(double);
ctx->electron.walker.num * (ctx->jastrow.cord_num + 1) * sizeof(double);
double* een_rescaled_n = (double*) qmckl_malloc(context, mem_info);
if (een_rescaled_n == NULL) {
@ -4260,7 +4250,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context)
}
rc = qmckl_compute_een_rescaled_n(context,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow.cord_num,
@ -4538,12 +4528,12 @@ qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context,
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
int64_t sze = ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow.cord_num + 1);
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_jastrow_factor_een_deriv_e",
"Array too small. Expected ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1)");
"Array too small. Expected ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow.cord_num + 1)");
}
memcpy(distance_rescaled, ctx->jastrow.een_rescaled_n_deriv_e, sze * sizeof(double));
@ -4579,12 +4569,17 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context)
/* Compute if necessary */
if (ctx->date > ctx->jastrow.een_rescaled_n_deriv_e_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->jastrow.een_rescaled_n_deriv_e);
ctx->jastrow.een_rescaled_n_deriv_e = NULL;
}
/* Allocate array */
if (ctx->jastrow.een_rescaled_n_deriv_e == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.num * 4 * ctx->nucleus.num *
ctx->electron.walk_num * (ctx->jastrow.cord_num + 1) * sizeof(double);
ctx->electron.walker.num * (ctx->jastrow.cord_num + 1) * sizeof(double);
double* een_rescaled_n_deriv_e = (double*) qmckl_malloc(context, mem_info);
if (een_rescaled_n_deriv_e == NULL) {
@ -4597,12 +4592,12 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context)
}
rc = qmckl_compute_factor_een_rescaled_n_deriv_e(context,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow.cord_num,
ctx->electron.rescale_factor_kappa_en,
ctx->electron.coord_new.data,
ctx->electron.walker.point.coord.data,
ctx->nucleus.coord.data,
ctx->electron.en_distance,
ctx->jastrow.een_rescaled_n,
@ -4634,8 +4629,8 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context)
| ~nucl_num~ | ~int64_t~ | in | Number of atoms |
| ~cord_num~ | ~int64_t~ | in | Order of polynomials |
| ~rescale_factor_kappa_en~ | ~double~ | in | Factor to rescale ee distances |
| ~coord_new~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates |
| ~coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates |
| ~coord_ee~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates |
| ~coord_en~ | ~double[3][nucl_num]~ | in | Nuclear coordinates |
| ~en_distance~ | ~double[walk_num][elec_num][nucl_num]~ | in | Electron-nucleus distances |
| ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus distances |
| ~een_rescaled_n_deriv_e~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | out | Electron-nucleus rescaled distances |
@ -4644,7 +4639,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context)
integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f( &
context, walk_num, elec_num, nucl_num, &
cord_num, rescale_factor_kappa_en, &
coord_new, coord, en_distance, een_rescaled_n, een_rescaled_n_deriv_e) &
coord_ee, coord_en, en_distance, een_rescaled_n, een_rescaled_n_deriv_e) &
result(info)
use qmckl
implicit none
@ -4654,8 +4649,8 @@ integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f( &
integer*8 , intent(in) :: nucl_num
integer*8 , intent(in) :: cord_num
double precision , intent(in) :: rescale_factor_kappa_en
double precision , intent(in) :: coord_new(elec_num,3,walk_num)
double precision , intent(in) :: coord(nucl_num,3)
double precision , intent(in) :: coord_ee(elec_num,3,walk_num)
double precision , intent(in) :: coord_en(nucl_num,3)
double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num)
double precision , intent(in) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
double precision , intent(out) :: een_rescaled_n_deriv_e(elec_num,4,nucl_num,0:cord_num,walk_num)
@ -4701,7 +4696,7 @@ integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f( &
do i = 1, elec_num
ria_inv = 1.0d0 / en_distance(i, a, nw)
do ii = 1, 3
elnuc_dist_deriv_e(ii, i, a) = (coord_new(i, ii, nw) - coord(a, ii)) * ria_inv
elnuc_dist_deriv_e(ii, i, a) = (coord_ee(i, ii, nw) - coord_en(a, ii)) * ria_inv
end do
elnuc_dist_deriv_e(4, i, a) = 2.0d0 * ria_inv
end do
@ -4747,8 +4742,8 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f
const int64_t nucl_num,
const int64_t cord_num,
const double rescale_factor_kappa_en,
const double* coord_new,
const double* coord,
const double* coord_ee,
const double* coord_en,
const double* en_distance,
const double* een_rescaled_n,
double* const een_rescaled_n_deriv_e );
@ -4765,8 +4760,8 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f
nucl_num, &
cord_num, &
rescale_factor_kappa_en, &
coord_new, &
coord, &
coord_ee, &
coord_en, &
en_distance, &
een_rescaled_n, &
een_rescaled_n_deriv_e) &
@ -4781,8 +4776,8 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f
integer (c_int64_t) , intent(in) , value :: nucl_num
integer (c_int64_t) , intent(in) , value :: cord_num
real (c_double ) , intent(in) , value :: rescale_factor_kappa_en
real (c_double ) , intent(in) :: coord_new(elec_num,3,walk_num)
real (c_double ) , intent(in) :: coord(nucl_num,3)
real (c_double ) , intent(in) :: coord_ee(elec_num,3,walk_num)
real (c_double ) , intent(in) :: coord_en(nucl_num,3)
real (c_double ) , intent(in) :: en_distance(nucl_num,elec_num,walk_num)
real (c_double ) , intent(in) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
real (c_double ) , intent(out) :: een_rescaled_n_deriv_e(elec_num,4,nucl_num,0:cord_num,walk_num)
@ -4795,8 +4790,8 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f
nucl_num, &
cord_num, &
rescale_factor_kappa_en, &
coord_new, &
coord, &
coord_ee, &
coord_en, &
en_distance, &
een_rescaled_n, &
een_rescaled_n_deriv_e)
@ -4991,7 +4986,7 @@ qmckl_exit_code qmckl_get_jastrow_tmp_c(qmckl_context context, double* const tmp
assert (ctx != NULL);
size_t sze = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1)
* ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num;
* ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num;
memcpy(tmp_c, ctx->jastrow.tmp_c, sze * sizeof(double));
return QMCKL_SUCCESS;
@ -5018,7 +5013,7 @@ qmckl_exit_code qmckl_get_jastrow_dtmp_c(qmckl_context context, double* const dt
assert (ctx != NULL);
size_t sze = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1)
*4* ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num;
*4* ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num;
memcpy(dtmp_c, ctx->jastrow.dtmp_c, sze * sizeof(double));
return QMCKL_SUCCESS;
@ -5176,12 +5171,17 @@ qmckl_exit_code qmckl_provide_tmp_c(qmckl_context context)
/* Compute if necessary */
if (ctx->date > ctx->jastrow.tmp_c_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->jastrow.tmp_c);
ctx->jastrow.tmp_c = NULL;
}
/* Allocate array */
if (ctx->jastrow.tmp_c == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1)
,* ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * sizeof(double);
,* ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num * sizeof(double);
double* tmp_c = (double*) qmckl_malloc(context, mem_info);
if (tmp_c == NULL) {
@ -5207,7 +5207,7 @@ qmckl_exit_code qmckl_provide_tmp_c(qmckl_context context)
ctx->jastrow.cord_num,
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->jastrow.een_rescaled_e,
ctx->jastrow.een_rescaled_n,
ctx->jastrow.tmp_c);
@ -5216,7 +5216,7 @@ qmckl_exit_code qmckl_provide_tmp_c(qmckl_context context)
ctx->jastrow.cord_num,
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->jastrow.een_rescaled_e,
ctx->jastrow.een_rescaled_n,
ctx->jastrow.tmp_c);
@ -5225,7 +5225,7 @@ qmckl_exit_code qmckl_provide_tmp_c(qmckl_context context)
ctx->jastrow.cord_num,
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->jastrow.een_rescaled_e,
ctx->jastrow.een_rescaled_n,
ctx->jastrow.tmp_c);
@ -5237,7 +5237,7 @@ qmckl_exit_code qmckl_provide_tmp_c(qmckl_context context)
ctx->jastrow.cord_num,
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->jastrow.een_rescaled_e,
ctx->jastrow.een_rescaled_n,
ctx->jastrow.tmp_c);
@ -5266,12 +5266,17 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context)
/* Compute if necessary */
if (ctx->date > ctx->jastrow.dtmp_c_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->jastrow.dtmp_c);
ctx->jastrow.dtmp_c = NULL;
}
/* Allocate array */
if (ctx->jastrow.dtmp_c == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1)
,* 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * sizeof(double);
,* 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num * sizeof(double);
double* dtmp_c = (double*) qmckl_malloc(context, mem_info);
if (dtmp_c == NULL) {
@ -5296,7 +5301,7 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context)
ctx->jastrow.cord_num,
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->jastrow.een_rescaled_e_deriv_e,
ctx->jastrow.een_rescaled_n,
ctx->jastrow.dtmp_c);
@ -5305,7 +5310,7 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context)
ctx->jastrow.cord_num,
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->jastrow.een_rescaled_e_deriv_e,
ctx->jastrow.een_rescaled_n,
ctx->jastrow.dtmp_c);
@ -5314,7 +5319,7 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context)
ctx->jastrow.cord_num,
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->jastrow.een_rescaled_e_deriv_e,
ctx->jastrow.een_rescaled_n,
ctx->jastrow.dtmp_c);
@ -5326,7 +5331,7 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context)
ctx->jastrow.cord_num,
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->jastrow.een_rescaled_e_deriv_e,
ctx->jastrow.een_rescaled_n,
ctx->jastrow.dtmp_c);
@ -6942,7 +6947,6 @@ assert(fabs(tmp_c[0][0][1][0][0] - 2.7083473948352403) < 1e-12);
printf("%e\n%e\n", tmp_c[0][1][0][0][0],0.237440520852232);
assert(fabs(dtmp_c[0][1][0][0][0][0] - 0.237440520852232) < 1e-12);
return QMCKL_SUCCESS;
#+end_src
** Electron-electron-nucleus Jastrow \(f_{een}\)
@ -6978,7 +6982,7 @@ qmckl_get_jastrow_factor_een(qmckl_context context,
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = ctx->electron.walk_num;
int64_t sze = ctx->electron.walker.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
@ -7032,11 +7036,16 @@ qmckl_exit_code qmckl_provide_factor_een(qmckl_context context)
/* Compute if necessary */
if (ctx->date > ctx->jastrow.factor_een_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->jastrow.factor_een);
ctx->jastrow.factor_een = NULL;
}
/* Allocate array */
if (ctx->jastrow.factor_een == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walk_num * sizeof(double);
mem_info.size = ctx->electron.walker.num * sizeof(double);
double* factor_een = (double*) qmckl_malloc(context, mem_info);
if (factor_een == NULL) {
@ -7049,7 +7058,7 @@ qmckl_exit_code qmckl_provide_factor_een(qmckl_context context)
}
rc = qmckl_compute_factor_een(context,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow.cord_num,
@ -7441,7 +7450,6 @@ double factor_een[walk_num];
rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0]),walk_num);
assert(fabs(factor_een[0] + 0.37407972141304213) < 1e-12);
return QMCKL_SUCCESS;
#+end_src
** Electron-electron-nucleus Jastrow \(f_{een}\) derivative
@ -7477,7 +7485,7 @@ qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context,
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num;
int64_t sze = ctx->electron.walker.num * 4 * ctx->electron.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
@ -7543,11 +7551,16 @@ qmckl_exit_code qmckl_provide_factor_een_deriv_e(qmckl_context context)
/* Compute if necessary */
if (ctx->date > ctx->jastrow.factor_een_deriv_e_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->jastrow.factor_een_deriv_e);
ctx->jastrow.factor_een_deriv_e = NULL;
}
/* Allocate array */
if (ctx->jastrow.factor_een_deriv_e == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 4 * ctx->electron.num * ctx->electron.walk_num * sizeof(double);
mem_info.size = 4 * ctx->electron.num * ctx->electron.walker.num * sizeof(double);
double* factor_een_deriv_e = (double*) qmckl_malloc(context, mem_info);
if (factor_een_deriv_e == NULL) {
@ -7560,7 +7573,7 @@ qmckl_exit_code qmckl_provide_factor_een_deriv_e(qmckl_context context)
}
rc = qmckl_compute_factor_een_deriv_e(context,
ctx->electron.walk_num,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow.cord_num,

139
org/qmckl_local_energy.org
View File

@ -158,12 +158,12 @@ typedef struct qmckl_local_energy_struct {
double * accep_prob;
double * r_drift;
double * y_move;
int64_t e_kin_date;
int64_t e_pot_date;
int64_t e_local_date;
int64_t accep_prob_date;
int64_t r_drift_date;
int64_t y_move_date;
uint64_t e_kin_date;
uint64_t e_pot_date;
uint64_t e_local_date;
uint64_t accep_prob_date;
uint64_t r_drift_date;
uint64_t y_move_date;
int32_t uninitialized;
bool provided;
@ -217,10 +217,10 @@ qmckl_exit_code qmckl_get_kinetic_energy(qmckl_context context, double * const k
if(!qmckl_nucleus_provided(context)) return QMCKL_NOT_PROVIDED;
rc = qmckl_provide_ao_vgl(context);
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_vgl(context);
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_kinetic_energy(context);
@ -229,7 +229,7 @@ qmckl_exit_code qmckl_get_kinetic_energy(qmckl_context context, double * const k
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.walk_num * sizeof(double);
size_t sze = ctx->electron.walker.num * sizeof(double);
memcpy(kinetic_energy, ctx->local_energy.e_kin, sze);
return QMCKL_SUCCESS;
@ -304,13 +304,18 @@ qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context) {
}
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->local_energy.e_kin_date) {
if (ctx->electron.walker.point.date > ctx->local_energy.e_kin_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->local_energy.e_kin);
ctx->local_energy.e_kin = NULL;
}
/* Allocate array */
if (ctx->local_energy.e_kin == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walk_num * sizeof(double);
mem_info.size = ctx->electron.walker.num * sizeof(double);
double* e_kin = (double*) qmckl_malloc(context, mem_info);
if (e_kin == NULL) {
@ -324,7 +329,7 @@ qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context) {
if (ctx->det.type == 'G') {
rc = qmckl_compute_kinetic_energy(context,
ctx->det.walk_num,
ctx->electron.walker.num,
ctx->det.det_num_alpha,
ctx->det.det_num_beta,
ctx->electron.up_num,
@ -556,12 +561,9 @@ const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
rc = qmckl_set_electron_num (context, chbrclf_elec_up_num, chbrclf_elec_dn_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_electron_walk_num (context, chbrclf_walk_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord, chbrclf_walk_num*chbrclf_elec_num*3);
rc = qmckl_set_electron_coord (context, 'N', chbrclf_walk_num, elec_coord, chbrclf_walk_num*chbrclf_elec_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num);
@ -688,9 +690,6 @@ for(k = 0; k < det_num_beta; ++k)
rc = qmckl_set_determinant_type (context, typ);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_determinant_walk_num (context, chbrclf_walk_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_determinant_det_num_alpha (context, det_num_alpha);
assert (rc == QMCKL_SUCCESS);
@ -781,10 +780,10 @@ qmckl_exit_code qmckl_get_potential_energy(qmckl_context context, double * const
if(!qmckl_nucleus_provided(context)) return QMCKL_NOT_PROVIDED;
rc = qmckl_provide_ao_vgl(context);
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_vgl(context);
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_potential_energy(context);
@ -793,7 +792,7 @@ qmckl_exit_code qmckl_get_potential_energy(qmckl_context context, double * const
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.walk_num * sizeof(double);
size_t sze = ctx->electron.walker.num * sizeof(double);
memcpy(potential_energy, ctx->local_energy.e_pot, sze);
return QMCKL_SUCCESS;
@ -877,19 +876,24 @@ qmckl_exit_code qmckl_provide_potential_energy(qmckl_context context) {
}
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->local_energy.e_pot_date) {
if (ctx->electron.walker.point.date > ctx->local_energy.e_pot_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->local_energy.e_pot);
ctx->local_energy.e_pot = NULL;
}
/* Allocate array */
if (ctx->local_energy.e_pot == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walk_num * sizeof(double);
mem_info.size = ctx->electron.walker.num * sizeof(double);
double* e_pot = (double*) qmckl_malloc(context, mem_info);
if (e_pot == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_e_pot",
"qmckl_provide_potential_energy",
NULL);
}
ctx->local_energy.e_pot = e_pot;
@ -897,11 +901,11 @@ qmckl_exit_code qmckl_provide_potential_energy(qmckl_context context) {
if (ctx->det.type == 'G') {
rc = qmckl_compute_potential_energy(context,
ctx->det.walk_num,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.ee_pot,
ctx->electron.en_pot,
ctx->electron.ee_potential,
ctx->electron.en_potential,
ctx->nucleus.repulsion,
ctx->local_energy.e_pot);
} else {
@ -929,18 +933,18 @@ qmckl_exit_code qmckl_provide_potential_energy(qmckl_context context) {
:END:
#+NAME: qmckl_compute_potential_energy_args
| ~qmckl_context~ | ~context~ | in | Global state |
| ~int64_t~ | ~walk_num~ | in | Number of walkers |
| ~int64_t~ | ~elec_num~ | in | Number of electrons |
| ~int64_t~ | ~nucl_num~ | in | Number of MOs |
| ~double~ | ~ee_pot[walk_num]~ | in | ee potential |
| ~double~ | ~en_pot[walk_num]~ | in | en potential |
| ~double~ | ~repulsion~ | in | en potential |
| ~double~ | ~e_pot[walk_num]~ | out | Potential energy |
| ~qmckl_context~ | ~context~ | in | Global state |
| ~int64_t~ | ~walk_num~ | in | Number of walkers |
| ~int64_t~ | ~elec_num~ | in | Number of electrons |
| ~int64_t~ | ~nucl_num~ | in | Number of MOs |
| ~double~ | ~ee_potential[walk_num]~ | in | ee potential |
| ~double~ | ~en_potential[walk_num]~ | in | en potential |
| ~double~ | ~repulsion~ | in | en potential |
| ~double~ | ~e_pot[walk_num]~ | out | Potential energy |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_potential_energy_f(context, walk_num, &
elec_num, nucl_num, ee_pot, en_pot, repulsion, e_pot) &
elec_num, nucl_num, ee_potential, en_potential, repulsion, e_pot) &
result(info)
use qmckl
implicit none
@ -948,8 +952,8 @@ integer function qmckl_compute_potential_energy_f(context, walk_num, &
integer*8, intent(in) :: walk_num
integer*8, intent(in) :: elec_num
integer*8, intent(in) :: nucl_num
double precision, intent(in) :: ee_pot(walk_num)
double precision, intent(in) :: en_pot(walk_num)
double precision, intent(in) :: ee_potential(walk_num)
double precision, intent(in) :: en_potential(walk_num)
double precision, intent(in) :: repulsion
double precision, intent(inout) :: e_pot(walk_num)
integer*8 :: idet, iwalk, ielec, mo_id, imo
@ -971,9 +975,8 @@ integer function qmckl_compute_potential_energy_f(context, walk_num, &
return
endif
e_pot = 0.0d0 + repulsion
do iwalk = 1, walk_num
e_pot(iwalk) = e_pot(iwalk) + ee_pot(iwalk) + en_pot(iwalk)
e_pot(iwalk) = ee_potential(iwalk) + en_potential(iwalk) + repulsion
end do
end function qmckl_compute_potential_energy_f
@ -988,8 +991,8 @@ end function qmckl_compute_potential_energy_f
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const double* ee_pot,
const double* en_pot,
const double* ee_potential,
const double* en_potential,
const double repulsion,
double* const e_pot );
#+end_src
@ -999,7 +1002,7 @@ end function qmckl_compute_potential_energy_f
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_potential_energy &
(context, walk_num, elec_num, nucl_num, ee_pot, en_pot, repulsion, e_pot) &
(context, walk_num, elec_num, nucl_num, ee_potential, en_potential, repulsion, e_pot) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
@ -1009,14 +1012,14 @@ end function qmckl_compute_potential_energy_f
integer (c_int64_t) , intent(in) , value :: walk_num
integer (c_int64_t) , intent(in) , value :: elec_num
integer (c_int64_t) , intent(in) , value :: nucl_num
real (c_double ) , intent(in) :: ee_pot(walk_num)
real (c_double ) , intent(in) :: en_pot(walk_num)
real (c_double ) , intent(in) :: ee_potential(walk_num)
real (c_double ) , intent(in) :: en_potential(walk_num)
real (c_double ) , intent(in) , value :: repulsion
real (c_double ) , intent(out) :: e_pot(walk_num)
integer(c_int32_t), external :: qmckl_compute_potential_energy_f
info = qmckl_compute_potential_energy_f &
(context, walk_num, elec_num, nucl_num, ee_pot, en_pot, repulsion, e_pot)
(context, walk_num, elec_num, nucl_num, ee_potential, en_potential, repulsion, e_pot)
end function qmckl_compute_potential_energy
#+end_src
@ -1065,10 +1068,10 @@ qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double * const loc
if(!qmckl_nucleus_provided(context)) return QMCKL_NOT_PROVIDED;
rc = qmckl_provide_ao_vgl(context);
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_vgl(context);
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_local_energy(context);
@ -1077,12 +1080,12 @@ qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double * const loc
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->electron.walk_num;
const int64_t sze = ctx->electron.walker.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_local_energy",
"input array too small");
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_local_energy",
"input array too small");
}
memcpy(local_energy, ctx->local_energy.e_local, sze * sizeof(double));
@ -1159,13 +1162,18 @@ qmckl_exit_code qmckl_provide_local_energy(qmckl_context context) {
}
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->local_energy.e_local_date) {
if (ctx->electron.walker.point.date > ctx->local_energy.e_local_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->local_energy.e_local);
ctx->local_energy.e_local = NULL;
}
/* Allocate array */
if (ctx->local_energy.e_local == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walk_num * sizeof(double);
mem_info.size = ctx->electron.walker.num * sizeof(double);
double* local_energy = (double*) qmckl_malloc(context, mem_info);
if (local_energy == NULL) {
@ -1179,7 +1187,7 @@ qmckl_exit_code qmckl_provide_local_energy(qmckl_context context) {
if (ctx->det.type == 'G') {
rc = qmckl_compute_local_energy(context,
ctx->det.walk_num,
ctx->electron.walker.num,
ctx->local_energy.e_kin,
ctx->local_energy.e_pot,
ctx->local_energy.e_local);
@ -1327,10 +1335,10 @@ qmckl_exit_code qmckl_get_drift_vector(qmckl_context context, double * const dri
if(!qmckl_nucleus_provided(context)) return QMCKL_NOT_PROVIDED;
rc = qmckl_provide_ao_vgl(context);
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_vgl(context);
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_drift_vector(context);
@ -1339,7 +1347,7 @@ qmckl_exit_code qmckl_get_drift_vector(qmckl_context context, double * const dri
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.walk_num * ctx->electron.num * 3 * sizeof(double);
size_t sze = ctx->electron.walker.num * ctx->electron.num * 3 * sizeof(double);
memcpy(drift_vector, ctx->local_energy.r_drift, sze);
return QMCKL_SUCCESS;
@ -1398,13 +1406,18 @@ qmckl_exit_code qmckl_provide_drift_vector(qmckl_context context) {
}
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->local_energy.r_drift_date) {
if (ctx->electron.walker.point.date > ctx->local_energy.r_drift_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->local_energy.r_drift);
ctx->local_energy.r_drift = NULL;
}
/* Allocate array */
if (ctx->local_energy.r_drift == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walk_num * ctx->electron.num * 3 * sizeof(double);
mem_info.size = ctx->electron.walker.num * ctx->electron.num * 3 * sizeof(double);
double* r_drift = (double*) qmckl_malloc(context, mem_info);
if (r_drift == NULL) {
@ -1419,7 +1432,7 @@ qmckl_exit_code qmckl_provide_drift_vector(qmckl_context context) {
qmckl_exit_code rc;
if (ctx->det.type == 'G') {
rc = qmckl_compute_drift_vector(context,
ctx->det.walk_num,
ctx->electron.walker.num,
ctx->det.det_num_alpha,
ctx->det.det_num_beta,
ctx->electron.up_num,

97
org/qmckl_memory.org
View File

@ -156,7 +156,7 @@ void* qmckl_malloc(qmckl_context context, const qmckl_memory_info_struct info) {
assert (ctx->memory.element[pos].size == (size_t) 0);
/* Copy info at the new location */
ctx->memory.element[pos].size = info.size;
memcpy(&(ctx->memory.element[pos]), &info, sizeof(qmckl_memory_info_struct));
ctx->memory.element[pos].pointer = pointer;
ctx->memory.n_allocated += (size_t) 1;
}
@ -173,7 +173,7 @@ void* qmckl_malloc(qmckl_context context, const qmckl_memory_info_struct info) {
qmckl_context context = qmckl_context_create();
qmckl_memory_info_struct info = qmckl_memory_info_struct_zero;
info.size = (size_t) 3;
info.size = (size_t) 3*sizeof(int);
/* Allocate an array of ints */
int *a = (int*) qmckl_malloc(context, info);
@ -235,7 +235,7 @@ qmckl_exit_code qmckl_free(qmckl_context context, void * const ptr) {
/* Not found */
qmckl_unlock(context);
return qmckl_failwith(context,
QMCKL_FAILURE,
QMCKL_INVALID_ARG_2,
"qmckl_free",
"Pointer not found in context");
}
@ -272,7 +272,7 @@ assert(rc == QMCKL_SUCCESS);
/* Free again */
rc = qmckl_free(context, a);
assert(rc == QMCKL_FAILURE);
assert(rc == QMCKL_INVALID_ARG_2);
/* Clean up */
rc = qmckl_context_destroy(context);
@ -280,6 +280,95 @@ assert(rc == QMCKL_SUCCESS);
#+end_src
* Get the size of a memory block
All the blocks allocated with ~qmckl_malloc~ keep track of how many
bytes were allocated. Using ~qmckl_malloc_size~ allows to get this information.
# Header
#+begin_src c :tangle (eval h_private_func) :noexport
qmckl_exit_code
qmckl_get_malloc_info(qmckl_context context,
const void* pointer,
qmckl_memory_info_struct* info);
#+end_src
# Source
#+begin_src c :tangle (eval c)
qmckl_exit_code
qmckl_get_malloc_info(qmckl_context context,
const void* ptr,
qmckl_memory_info_struct* info)
{
assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT);
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
if (ptr == NULL) {
return qmckl_failwith(context,
QMCKL_INVALID_ARG_2,
"qmckl_get_malloc_info",
"Null pointer");
}
if (info == NULL) {
return qmckl_failwith(context,
QMCKL_INVALID_ARG_3,
"qmckl_get_malloc_info",
"Null pointer");
}
qmckl_lock(context);
{
/* Find the pointer entry */
size_t pos = (size_t) 0;
while ( pos < ctx->memory.array_size && ctx->memory.element[pos].pointer != ptr) {
pos += (size_t) 1;
}
if (pos >= ctx->memory.array_size) {
/* Not found */
qmckl_unlock(context);
return qmckl_failwith(context,
QMCKL_INVALID_ARG_2,
"qmckl_get_malloc_info",
"Pointer not found in context");
}
/* Copy info */
memcpy(info, &(ctx->memory.element[pos]), sizeof(qmckl_memory_info_struct));
}
qmckl_unlock(context);
return QMCKL_SUCCESS;
}
#+end_src
# Test :noexport:
#+begin_src c :tangle (eval c_test)
/* Create a context */
context = qmckl_context_create();
info = qmckl_memory_info_struct_zero;
info.size = (size_t) 3*sizeof(int);
/* Allocate an array of ints */
a = (int*) qmckl_malloc(context, info);
/* Check that the size of a is 3*sizeof(int) */
info = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, NULL, &info);
assert (rc == QMCKL_INVALID_ARG_2);
rc = qmckl_get_malloc_info(context, &rc, &info);
assert (rc == QMCKL_INVALID_ARG_2);
rc = qmckl_get_malloc_info(context, a, &info);
assert (rc == QMCKL_SUCCESS);
assert (info.size == 3*sizeof(int));
rc = qmckl_context_destroy(context);
#+end_src
* End of files :noexport:
#+begin_src c :comments org :tangle (eval h_private_func)

697
org/qmckl_mo.org
View File

@ -3,16 +3,16 @@
#+INCLUDE: ../tools/lib.org
The molecular orbitals (MOs) are defined in the basis of AOs along with a AO to MO
coefficient matrix \[C\]. Using these coefficients (e.g. from Hartree Fock SCF method)
coefficient matrix $C$. Using these coefficients (e.g. from Hartree Fock method)
the MOs are defined as follows:
\[
\phi_i(\mathbf{r}) = C_i * \chi_i (\mathbf{r})
\]
In this section we demonstrate how to use the QMCkl specific DGEMM
function to calculate the MOs.
By default, all the MOs are computed. A subset of MOs can be selected
for computation, for example to remove computation of the useless
virtual MOs present in a MO coefficient matrix.
* Headers :noexport:
@ -84,11 +84,11 @@ int main() {
The following arrays are stored in the context:
|-----------------+--------------------+----------------------------------------|
| ~mo_num~ | | Number of MOs |
| ~coefficient~ | ~[mo_num][ao_num]~ | Orbital coefficients |
| ~coefficient_t~ | ~[ao_num][mo_num]~ | Transposed of the Orbital coefficients |
|-----------------+--------------------+----------------------------------------|
|-------------------+--------------------+----------------------------------------|
| ~mo_num~ | | Number of MOs |
| ~coefficient~ | ~[mo_num][ao_num]~ | MO coefficients |
| ~coefficient_t~ | ~[ao_num][mo_num]~ | Transposed of the Orbital coefficients |
|-------------------+--------------------+----------------------------------------|
Computed data:
@ -103,7 +103,7 @@ int main() {
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_mo_basis_struct {
int64_t mo_num;
int64_t mo_num;
double * restrict coefficient;
double * restrict coefficient_t;
@ -144,150 +144,6 @@ qmckl_exit_code qmckl_init_mo_basis(qmckl_context context) {
}
#+end_src
** Access functions
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_mo_basis_mo_num (const qmckl_context context,
int64_t* mo_num);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_mo_basis_mo_num (const qmckl_context context,
int64_t* mo_num)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_mo_basis_mo_num",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1;
if ( (ctx->mo_basis.uninitialized & mask) != 0) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_get_mo_basis_mo_num",
NULL);
}
assert (ctx->mo_basis.mo_num > (int64_t) 0);
,*mo_num = ctx->mo_basis.mo_num;
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_mo_basis_coefficient (const qmckl_context context,
double* const coefficient,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code
qmckl_get_mo_basis_coefficient (const qmckl_context context,
double* const coefficient,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_mo_basis_coefficient",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 1;
if ( (ctx->ao_basis.uninitialized & mask) != 0) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_get_mo_basis_coefficient",
NULL);
}
if (coefficient == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_get_mo_basis_coefficient",
"NULL pointer");
}
if (size_max < ctx->ao_basis.ao_num * ctx->mo_basis.mo_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_mo_basis_coefficient",
"Array too small. Expected mo_num * ao_num");
}
assert (ctx->mo_basis.coefficient != NULL);
memcpy(coefficient, ctx->mo_basis.coefficient,
ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
When all the data for the AOs have been provided, the following
function returns ~true~.
#+begin_src c :comments org :tangle (eval h_func)
bool qmckl_mo_basis_provided (const qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :exports none
bool qmckl_mo_basis_provided(const qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
return ctx->mo_basis.provided;
}
#+end_src
*** Fortran interfaces
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_mo_basis_mo_num (context, &
mo_num) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(out) :: mo_num
end function qmckl_get_mo_basis_mo_num
end interface
interface
integer(c_int32_t) function qmckl_get_mo_basis_coefficient(context, &
coefficient, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
double precision, intent(out) :: coefficient(*)
integer (c_int64_t) , intent(in), value :: size_max
end function qmckl_get_mo_basis_coefficient
end interface
#+end_src
** Initialization functions
To set the basis set, all the following functions need to be
@ -398,6 +254,15 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (ctx->mo_basis.coefficient_t != NULL) {
qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient_t);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_finalize_mo_basis",
NULL);
}
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(double);
double* new_array = (double*) qmckl_malloc(context, mem_info);
@ -410,15 +275,6 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
assert (ctx->mo_basis.coefficient != NULL);
if (ctx->mo_basis.coefficient_t != NULL) {
qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_finalize_mo_basis",
NULL);
}
}
for (int64_t i=0 ; i<ctx->ao_basis.ao_num ; ++i) {
for (int64_t j=0 ; j<ctx->mo_basis.mo_num ; ++j) {
new_array[i*ctx->mo_basis.mo_num + j] = ctx->mo_basis.coefficient[j*ctx->ao_basis.ao_num + i];
@ -426,11 +282,269 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
}
ctx->mo_basis.coefficient_t = new_array;
qmckl_exit_code rc = QMCKL_SUCCESS;
return rc;
qmckl_exit_code rc;
if (ctx->mo_basis.mo_vgl != NULL) {
rc = qmckl_free(context, ctx->mo_basis.mo_vgl);
if (rc != QMCKL_SUCCESS) return rc;
ctx->mo_basis.mo_vgl = NULL;
ctx->mo_basis.mo_vgl_date = 0;
}
if (ctx->mo_basis.mo_value != NULL) {
rc = qmckl_free(context, ctx->mo_basis.mo_value);
if (rc != QMCKL_SUCCESS) return rc;
ctx->mo_basis.mo_value = NULL;
ctx->mo_basis.mo_value_date = 0;
}
return QMCKL_SUCCESS;
}
#+end_src
** Access functions
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_mo_basis_mo_num (const qmckl_context context,
int64_t* mo_num);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_mo_basis_mo_num (const qmckl_context context,
int64_t* mo_num)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_mo_basis_mo_num",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1;
if ( (ctx->mo_basis.uninitialized & mask) != 0) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_get_mo_basis_mo_num",
NULL);
}
assert (ctx->mo_basis.mo_num > (int64_t) 0);
,*mo_num = ctx->mo_basis.mo_num;
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_mo_basis_coefficient (const qmckl_context context,
double* const coefficient,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code
qmckl_get_mo_basis_coefficient (const qmckl_context context,
double* const coefficient,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_mo_basis_coefficient",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 1;
if ( (ctx->mo_basis.uninitialized & mask) != 0) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_get_mo_basis_coefficient",
NULL);
}
if (coefficient == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_get_mo_basis_coefficient",
"NULL pointer");
}
if (size_max < ctx->ao_basis.ao_num * ctx->mo_basis.mo_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_mo_basis_coefficient",
"Array too small: expected mo_num * ao_num.");
}
assert (ctx->mo_basis.coefficient != NULL);
memcpy(coefficient, ctx->mo_basis.coefficient,
ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
When all the data for the AOs have been provided, the following
function returns ~true~.
#+begin_src c :comments org :tangle (eval h_func)
bool qmckl_mo_basis_provided (const qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :exports none
bool qmckl_mo_basis_provided(const qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
return ctx->mo_basis.provided;
}
#+end_src
*** Fortran interfaces
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_mo_basis_mo_num (context, &
mo_num) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(out) :: mo_num
end function qmckl_get_mo_basis_mo_num
end interface
interface
integer(c_int32_t) function qmckl_get_mo_basis_coefficient(context, &
coefficient, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
double precision, intent(out) :: coefficient(*)
integer (c_int64_t) , intent(in), value :: size_max
end function qmckl_get_mo_basis_coefficient
end interface
#+end_src
** Update
Useless MOs can be removed, for instance virtual MOs in a single
determinant calculation.
To select a subset of MOs that will be kept, create an array of
integers of size =mo_num=. If the integer is zero, the MO is dropped,
otherwise it is kept.
#+begin_src c :comments org :tangle (eval h_func)
bool qmckl_mo_basis_select_mo (const qmckl_context context,
const int32_t* keep,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :exports none
bool qmckl_mo_basis_select_mo (const qmckl_context context,
const int32_t* keep,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_mo_basis_select_mo",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if ( !(qmckl_mo_basis_provided(context)) ) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_get_mo_basis_select_mo",
NULL);
}
if (keep == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_get_mo_basis_select_mo",
"NULL pointer");
}
const int64_t mo_num = ctx->mo_basis.mo_num;
const int64_t ao_num = ctx->ao_basis.ao_num;
if (size_max < mo_num) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_mo_basis_select_mo",
"Array too small: expected mo_num.");
}
int64_t mo_num_new = 0;
for (int64_t i=0 ; i<mo_num ; ++i) {
if (keep[i] != 0) ++mo_num_new;
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ao_num * mo_num_new * sizeof(double);
double* restrict coefficient = (double*) qmckl_malloc(context, mem_info);
int64_t k = 0;
for (int64_t i=0 ; i<mo_num ; ++i) {
if (keep[i] != 0) {
memcpy(&(coefficient[k*ao_num]), &(ctx->mo_basis.coefficient[i*ao_num]), ao_num*sizeof(double));
++k;
}
}
qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient);
if (rc != QMCKL_SUCCESS) return rc;
ctx->mo_basis.coefficient = coefficient;
ctx->mo_basis.mo_num = mo_num_new;
rc = qmckl_finalize_mo_basis(context);
return rc;
}
#+end_src
*** Fortran interface
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_mo_basis_select_mo (context, &
keep, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in), value :: context
integer (c_int32_t) , intent(in) :: keep(*)
integer (c_int64_t) , intent(in), value :: size_max
end function qmckl_mo_basis_select_mo
end interface
#+end_src
* Computation
** Computation of MOs: values only
@ -457,10 +571,7 @@ qmckl_get_mo_basis_mo_value(qmckl_context context,
qmckl_exit_code rc;
rc = qmckl_provide_ao_basis_ao_value(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_value(context);
rc = qmckl_provide_mo_basis_mo_value(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
@ -537,7 +648,7 @@ qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
ctx->mo_basis.mo_value = mo_value;
rc = qmckl_provide_mo_value(context);
rc = qmckl_provide_mo_basis_mo_value(context);
if (rc != QMCKL_SUCCESS) return rc;
ctx->mo_basis.mo_value = old_array;
@ -562,52 +673,62 @@ qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
*** Provide
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_mo_value(qmckl_context context);
#+end_src
#+CALL: write_provider_header( group="mo_basis", data="mo_value" )
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_mo_value(qmckl_context context)
#+RESULTS:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :export none
qmckl_exit_code qmckl_provide_mo_basis_mo_value(qmckl_context context);
#+end_src
#+CALL: write_provider_pre( group="mo_basis", data="mo_value", dimension="ctx->mo_basis.mo_num * ctx->point.num")
#+RESULTS:
#+begin_src c :comments org :tangle (eval c) :noweb yes :export none
qmckl_exit_code qmckl_provide_mo_basis_mo_value(qmckl_context context)
{
qmckl_exit_code rc;
qmckl_exit_code rc = QMCKL_SUCCESS;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_provide_mo_basis_mo_value",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_basis",
NULL);
}
rc = qmckl_provide_ao_basis_ao_value(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_value",
NULL);
}
if (!ctx->mo_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_mo_basis",
"qmckl_provide_mo_basis_mo_value",
NULL);
}
/* Compute if necessary */
if (ctx->point.date > ctx->mo_basis.mo_value_date) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->mo_basis.mo_num * ctx->point.num * sizeof(double);
if (ctx->mo_basis.mo_value != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->mo_basis.mo_value, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->mo_basis.mo_value);
assert (rc == QMCKL_SUCCESS);
ctx->mo_basis.mo_value = NULL;
}
}
/* Allocate array */
if (ctx->mo_basis.mo_value == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->point.num * ctx->mo_basis.mo_num * sizeof(double);
double* mo_value = (double*) qmckl_malloc(context, mem_info);
if (mo_value == NULL) {
@ -619,6 +740,10 @@ qmckl_exit_code qmckl_provide_mo_value(qmckl_context context)
ctx->mo_basis.mo_value = mo_value;
}
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
if (ctx->mo_basis.mo_vgl_date == ctx->point.date) {
// mo_vgl has been computed at this step: Just copy the data.
@ -635,6 +760,14 @@ qmckl_exit_code qmckl_provide_mo_value(qmckl_context context)
} else {
rc = qmckl_provide_ao_basis_ao_value(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_value",
NULL);
}
rc = qmckl_compute_mo_basis_mo_value(context,
ctx->ao_basis.ao_num,
ctx->mo_basis.mo_num,
@ -643,10 +776,17 @@ qmckl_exit_code qmckl_provide_mo_value(qmckl_context context)
ctx->ao_basis.ao_value,
ctx->mo_basis.mo_value);
if (rc != QMCKL_SUCCESS) {
return rc;
}
}
#+end_src
#+CALL: write_provider_post( group="mo_basis", data="mo_value" )
#+RESULTS:
#+begin_src c :comments org :tangle (eval c) :noweb yes :export none
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->mo_basis.mo_value_date = ctx->date;
@ -654,7 +794,7 @@ qmckl_exit_code qmckl_provide_mo_value(qmckl_context context)
return QMCKL_SUCCESS;
}
#+end_src
#+end_src
*** Compute
:PROPERTIES:
@ -704,8 +844,8 @@ integer function qmckl_compute_mo_basis_mo_value_doc_f(context, &
do j=1,point_num
mo_value(:,j) = 0.d0
do k=1,ao_num
if (ao_value(k,j) /= 0.d0) then
c1 = ao_value(k,j)
c1 = ao_value(k,j)
if (c1 /= 0.d0) then
do i=1,mo_num
mo_value(i,j) = mo_value(i,j) + coefficient_t(i,k) * c1
end do
@ -713,7 +853,7 @@ integer function qmckl_compute_mo_basis_mo_value_doc_f(context, &
end do
end do
else ! dgemm
else ! dgemm for checking
LDA = size(coefficient_t,1)
LDB = size(ao_value,1)
@ -851,7 +991,7 @@ qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
}
}
int64_t n;
int64_t n=0;
for (n=0 ; n < nidx-4 ; n+=4) {
const double* restrict ck1 = coefficient_t + idx[n ]*mo_num;
@ -872,8 +1012,7 @@ qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
}
}
const int64_t n0 = n < 0 ? 0 : n;
for (int64_t m=n0 ; m < nidx ; m+=1) {
for (int64_t m=n ; m < nidx ; m+=1) {
const double* restrict ck = coefficient_t + idx[m]*mo_num;
const double a1 = av1[m];
@ -914,10 +1053,7 @@ qmckl_get_mo_basis_mo_vgl(qmckl_context context,
qmckl_exit_code rc;
rc = qmckl_provide_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_vgl(context);
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
@ -994,7 +1130,7 @@ qmckl_get_mo_basis_mo_vgl_inplace (qmckl_context context,
ctx->mo_basis.mo_vgl = mo_vgl;
rc = qmckl_provide_mo_vgl(context);
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
ctx->mo_basis.mo_vgl = old_array;
@ -1019,52 +1155,62 @@ qmckl_get_mo_basis_mo_vgl_inplace (qmckl_context context,
*** Provide
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context);
#+end_src
#+CALL: write_provider_header( group="mo_basis", data="mo_vgl" )
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
#+RESULTS:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :export none
qmckl_exit_code qmckl_provide_mo_basis_mo_vgl(qmckl_context context);
#+end_src
#+CALL: write_provider_pre( group="mo_basis", data="mo_vgl", dimension="5 * ctx->mo_basis.mo_num * ctx->point.num")
#+RESULTS:
#+begin_src c :comments org :tangle (eval c) :noweb yes :export none
qmckl_exit_code qmckl_provide_mo_basis_mo_vgl(qmckl_context context)
{
qmckl_exit_code rc;
qmckl_exit_code rc = QMCKL_SUCCESS;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_provide_mo_basis_mo_vgl",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_basis",
NULL);
}
rc = qmckl_provide_ao_vgl(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_basis",
NULL);
}
if (!ctx->mo_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_mo_basis",
"qmckl_provide_mo_basis_mo_vgl",
NULL);
}
/* Compute if necessary */
if (ctx->point.date > ctx->mo_basis.mo_vgl_date) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 5 * ctx->mo_basis.mo_num * ctx->point.num * sizeof(double);
if (ctx->mo_basis.mo_vgl != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->mo_basis.mo_vgl, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->mo_basis.mo_vgl);
assert (rc == QMCKL_SUCCESS);
ctx->mo_basis.mo_vgl = NULL;
}
}
/* Allocate array */
if (ctx->mo_basis.mo_vgl == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 5 * ctx->point.num * ctx->mo_basis.mo_num * sizeof(double);
double* mo_vgl = (double*) qmckl_malloc(context, mem_info);
if (mo_vgl == NULL) {
@ -1076,6 +1222,17 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
ctx->mo_basis.mo_vgl = mo_vgl;
}
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_basis",
NULL);
}
rc = qmckl_compute_mo_basis_mo_vgl(context,
ctx->ao_basis.ao_num,
ctx->mo_basis.mo_num,
@ -1083,6 +1240,12 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
ctx->mo_basis.coefficient_t,
ctx->ao_basis.ao_vgl,
ctx->mo_basis.mo_vgl);
#+end_src
#+CALL: write_provider_post( group="mo_basis", data="mo_vgl" )
#+RESULTS:
#+begin_src c :comments org :tangle (eval c) :noweb yes :export none
if (rc != QMCKL_SUCCESS) {
return rc;
}
@ -1092,7 +1255,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
return QMCKL_SUCCESS;
}
#+end_src
#+end_src
*** Compute
:PROPERTIES:
@ -1308,7 +1471,7 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
}
}
int64_t n;
int64_t n=0;
for (n=0 ; n < nidx-4 ; n+=4) {
const double* restrict ck1 = coefficient_t + idx[n ]*mo_num;
@ -1353,8 +1516,7 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
}
}
const int64_t n0 = n < 0 ? 0 : n;
for (int64_t m=n0 ; m < nidx ; m+=1) {
for (int64_t m=n ; m < nidx ; m+=1) {
const double* restrict ck = coefficient_t + idx[m]*mo_num;
const double a1 = av1[m];
const double a2 = av2[m];
@ -1409,7 +1571,7 @@ elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] )
nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] )
nucl_2 = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] )
#double prim_vgl[prim_num][5][walk_num][elec_num];
#double prim_vgl[prim_num][5][point_num];
x = elec_26_w1 ; y = nucl_1
a = [( 8.236000E+03, -1.130000E-04 * 6.1616545431994848e+02 ),
( 1.235000E+03, -8.780000E-04 * 1.4847738511079908e+02 ),
@ -1457,15 +1619,14 @@ const int64_t nucl_num = chbrclf_nucl_num;
const double* nucl_charge = chbrclf_charge;
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
const int64_t point_num = walk_num*elec_num;
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_electron_walk_num (context, walk_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
rc = qmckl_set_point(context, 'N', point_num, elec_coord, point_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, nucl_num);
@ -1550,14 +1711,14 @@ assert(rc == QMCKL_SUCCESS);
assert(qmckl_ao_basis_provided(context));
double ao_vgl[walk_num*elec_num][5][chbrclf_ao_num];
double ao_vgl[point_num][5][chbrclf_ao_num];
rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]),
(int64_t) 5*walk_num*elec_num*chbrclf_ao_num);
(int64_t) 5*point_num*chbrclf_ao_num);
assert (rc == QMCKL_SUCCESS);
/* Set up MO data */
const int64_t mo_num = chbrclf_mo_num;
int64_t mo_num = chbrclf_mo_num;
rc = qmckl_set_mo_basis_mo_num(context, mo_num);
assert (rc == QMCKL_SUCCESS);
@ -1571,15 +1732,15 @@ assert(qmckl_mo_basis_provided(context));
rc = qmckl_context_touch(context);
assert (rc == QMCKL_SUCCESS);
double mo_value[walk_num*elec_num][chbrclf_mo_num];
rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), walk_num*elec_num*chbrclf_mo_num);
double mo_value[point_num][chbrclf_mo_num];
rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), point_num*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
double mo_vgl[walk_num*elec_num][5][chbrclf_mo_num];
rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), walk_num*elec_num*5*chbrclf_mo_num);
double mo_vgl[point_num][5][chbrclf_mo_num];
rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), point_num*5*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
for (int i=0 ; i< walk_num*elec_num; ++i) {
for (int i=0 ; i< point_num; ++i) {
for (int k=0 ; k< chbrclf_mo_num ; ++k) {
assert(fabs(mo_vgl[i][0][k] - mo_value[i][k]) < 1.e-12) ;
}
@ -1588,10 +1749,10 @@ for (int i=0 ; i< walk_num*elec_num; ++i) {
rc = qmckl_context_touch(context);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), walk_num*elec_num*chbrclf_mo_num);
rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), point_num*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
for (int i=0 ; i< walk_num*elec_num; ++i) {
for (int i=0 ; i< point_num; ++i) {
for (int k=0 ; k< chbrclf_mo_num ; ++k) {
assert(fabs(mo_vgl[i][0][k] - mo_value[i][k]) < 1.e-12) ;
}
@ -1629,7 +1790,7 @@ for (int i=0 ; i< walk_num*elec_num; ++i) {
// assert(rc == QMCKL_SUCCESS);
//
// // Calculate value of MO (1st electron)
// double mo_vgl[5][walk_num][elec_num][chbrclf_mo_num];
// double mo_vgl[5][point_num][chbrclf_mo_num];
// rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0][0]));
// assert (rc == QMCKL_SUCCESS);
// ovlmo1 += mo_vgl[0][0][0][0]*mo_vgl[0][0][0][0]*dr3;
@ -1653,6 +1814,32 @@ printf(" mo_vgl mo_vgl[1][26][223] %25.15e\n", mo_vgl[2][1][3]);
printf(" mo_vgl mo_vgl[0][26][224] %25.15e\n", mo_vgl[2][0][3]);
printf(" mo_vgl mo_vgl[1][26][224] %25.15e\n", mo_vgl[2][1][3]);
printf("\n");
/* Check selection of MOs */
int32_t keep[mo_num];
for (int i=0 ; i<mo_num ; ++i) {
keep[i] = 0;
}
keep[2] = 1;
keep[5] = 1;
rc = qmckl_mo_basis_select_mo(context, &(keep[0]), mo_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_get_mo_basis_mo_num(context, &mo_num);
printf(" mo_num: %ld\n", mo_num);
assert(mo_num == 2);
double mo_coefficient_new[mo_num][ao_num];
rc = qmckl_get_mo_basis_coefficient (context, &(mo_coefficient_new[0][0]), mo_num*ao_num);
for (int i=0 ; i<ao_num ; ++i) {
assert(mo_coefficient_new[0][i] == mo_coefficient[i + ao_num*2]);
assert(mo_coefficient_new[1][i] == mo_coefficient[i + ao_num*5]);
}
}
#+end_src

21
org/qmckl_point.org
View File

@ -2,10 +2,10 @@
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
This data structure contains cartesian coordinates where the functions
will be evaluated. For DFT codes these may be the integration grid
points. For QMC codes, these are the electron coordinates of all the
walkers.
This data structure contains cartesian coordinates where the
3-dimensional functions will be evaluated. For DFT codes these may be
the integration grid points. For QMC codes, these are the electron
coordinates of all the walkers.
* Headers :noexport:
#+begin_src elisp :noexport :results none
@ -19,7 +19,7 @@ walkers.
#include <stdbool.h>
#include "qmckl_blas_private_type.h"
#+end_src
#+begin_src c :tangle (eval h_private_func)
#ifndef QMCKL_POINT_HPF
#define QMCKL_POINT_HPF
@ -80,9 +80,7 @@ int main() {
| Variable | Type | Description |
|--------------+----------------+-------------------------------------------|
| ~num~ | ~int64_t~ | Total number of points |
| ~alloc_num~ | ~int64_t~ | Numer of allocated number of points |
| ~date~ | ~uint64_t~ | Last modification date of the coordinates |
| ~alloc_date~ | ~uint64_t~ | Last modification date of the allocation |
| ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix |
We consider that the matrix is stored 'transposed' and 'normal'
@ -93,9 +91,7 @@ int main() {
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_point_struct {
int64_t num;
int64_t alloc_num;
uint64_t date;
uint64_t alloc_date;
qmckl_matrix coord;
} qmckl_point_struct;
@ -312,7 +308,7 @@ qmckl_set_point (qmckl_context context,
assert (ctx != NULL);
qmckl_exit_code rc;
if (num > ctx->point.alloc_num) {
if (num != ctx->point.num) {
if (ctx->point.coord.data != NULL) {
rc = qmckl_matrix_free(context, &(ctx->point.coord));
@ -353,11 +349,6 @@ qmckl_set_point (qmckl_context context,
rc = qmckl_context_touch(context);
assert (rc == QMCKL_SUCCESS);
if (num > ctx->point.alloc_num) {
ctx->point.alloc_num = num;
ctx->point.alloc_date = ctx->point.date;
};
return QMCKL_SUCCESS;
}

4
org/qmckl_trexio.org
View File

@ -1137,7 +1137,7 @@ qmckl_trexio_read(const qmckl_context context, const char* file_name, const int6
rc = qmckl_failwith( context,
QMCKL_FAILURE,
"qmckl_trexio_read",
"QMCkl was not compiled without TREXIO");
"QMCkl was compiled without TREXIO");
#endif
return rc;
}
@ -1147,7 +1147,6 @@ qmckl_trexio_read(const qmckl_context context, const char* file_name, const int6
#+begin_src c :tangle (eval c_test)
#ifdef HAVE_TREXIO
#define walk_num 2
qmckl_exit_code rc;
char fname[256];
@ -1161,7 +1160,6 @@ strncpy(fname, QMCKL_TEST_DIR,255);
strncat(fname, "/chbrclf", 255);
printf("Test file: %s\n", fname);
rc = qmckl_set_electron_walk_num(context, walk_num);
rc = qmckl_trexio_read(context, fname, 255);
if (rc != QMCKL_SUCCESS) {

2
python/setup.py
View File

@ -46,7 +46,7 @@ qmckl_module = Extension(name = "._" + MODULE_NAME,
setup(name = MODULE_NAME,
version = "0.2.0",
version = "0.3.1",
author = "TREX-CoE",
author_email = "posenitskiy@irsamc.ups-tlse.fr",
description = """Python API of the QMCkl library""",

4
python/test/test_api.py
View File

@ -28,12 +28,10 @@ fname = join('data', 'Alz_small.h5')
pq.trexio_read(ctx, fname)
print('trexio_read: passed')
pq.set_electron_walk_num(ctx, walk_num)
mo_num = pq.get_mo_basis_mo_num(ctx)
assert mo_num == 404
pq.set_electron_coord(ctx, 'T', coord)
pq.set_electron_coord(ctx, 'T', walk_num, coord)
ao_type = pq.get_ao_basis_type(ctx)
assert 'G' in ao_type

1
share/doc/qmckl/html/emacs-htmlize Submodule

@ -0,0 +1 @@
Subproject commit dd27bc3f26efd728f2b1f01f9e4ac4f61f2ffbf9

22
tools/build_makefile.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/env python2
#!/usr/bin/env python
#
# Creates all the dependencies from the org-mode files
@ -248,10 +248,22 @@ def main():
"" ]
tmp = "EXTRA_DIST += "
r = subprocess.check_output("git ls-tree --name-only -r HEAD".split())
for line in r.splitlines():
if b"share/qmckl/test_data/" in line:
tmp += " \\\n " + line.decode('utf8')
# Manually get the benchmark data list
for item in glob(os.path.join("share/qmckl/test_data/", "*")):
if os.path.isfile(item):
#print(f'Depth 1, file item: {item}')
tmp += " \\\n " + item
elif os.path.isdir(item):
for item2 in glob(os.path.join(item, "*")):
if os.path.isfile(item2):
tmp += " \\\n " + item2
elif os.path.isdir(item2):
for item3 in glob(os.path.join(item2, "*")):
if os.path.isfile(item3):
tmp += " \\\n " + item3
elif os.path.isdir(item3):
print("Stopping at the depth level 3, not processing")
tmp += "\n"
output += tmp.split("\n")

23
tools/install_htmlize.sh
View File

@ -2,8 +2,25 @@
#
# Installs the htmlize Emacs plugin
${abs_srcdir}/tools/missing git clone "https://github.com/TREX-CoE/emacs-htmlize"
mv emacs-htmlize/htmlize.el $1
rm -rf emacs-htmlize
readonly HTMLIZE_EL=$1
# If the file already present - exit the script
if [ -f ${HTMLIZE_EL} ]
then
exit 0
fi
DOC_ROOT=${HTMLIZE_EL%/*}
readonly EMACS_HTMLIZE=${DOC_ROOT}/emacs-htmlize/htmlize.el
# Case 1: submodule cloned but htmlize.el is absent in DOC_ROOT
if [ -f ${EMACS_HTMLIZE} ]
then
cp ${EMACS_HTMLIZE} ${HTMLIZE_EL}
else
# Case 2: submodule has not been cloned
${abs_srcdir}/tools/missing git clone "https://github.com/TREX-CoE/emacs-htmlize"
mv emacs-htmlize/htmlize.el ${HTMLIZE_EL}
rm -rf emacs-htmlize
fi

17
tools/lib.org
View File

@ -278,7 +278,7 @@ return msg
template = """qmckl_exit_code qmckl_provide_{{ group }}_{{ data }}(qmckl_context context)
{
qmckl_exit_code rc;
qmckl_exit_code rc = QMCKL_SUCCESS;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
@ -300,8 +300,17 @@ template = """qmckl_exit_code qmckl_provide_{{ group }}_{{ data }}(qmckl_context
/* Compute if necessary */
if (ctx->point.date > ctx->{{ group }}.{{ data }}_date) {
if (ctx->point.alloc_date > ctx->{{ group }}.{{ data }}_date) {
if (ctx->{{ group }}.{{ data }} != NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = {{ dimension }} * sizeof(double);
if (ctx->{{ group }}.{{ data }} != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->{{ group }}.{{ data }}, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->{{ group }}.{{ data }});
assert (rc == QMCKL_SUCCESS);
ctx->{{ group }}.{{ data }} = NULL;
@ -311,8 +320,6 @@ template = """qmckl_exit_code qmckl_provide_{{ group }}_{{ data }}(qmckl_context
/* Allocate array */
if (ctx->{{ group }}.{{ data }} == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = {{ dimension }} * sizeof(double);
double* {{ data }} = (double*) qmckl_malloc(context, mem_info);
if ({{ data }} == NULL) {

118
tools/qmckl.scm Normal file
View File

@ -0,0 +1,118 @@
(define-module (qmckl)
#:use-module (guix packages)
#:use-module (gnu packages autotools)
#:use-module (gnu packages pkg-config)
#:use-module (gnu packages gcc)
#:use-module (gnu packages maths) ;; contains blas, lapack, hdf5
#:use-module (trexio) ;; custom trexio module, has to be provided via the -L command line option
#:use-module (guix download)
#:use-module (guix build-system gnu)
#:use-module (guix licenses)
;; the modules below are additional requirements for building the dev version
#:use-module (guix git-download)
#:use-module (gnu packages python)
#:use-module (gnu packages emacs)
#:use-module (gnu packages swig)
#:use-module (gnu packages m4))
(define-public qmckl-hpc-0.2.1
(package
(name "qmckl-hpc")
(version "0.2.1")
(source (origin
(method url-fetch)
(uri (string-append "https://github.com/TREX-CoE/qmckl/releases/download/v" version
"/qmckl-" version
".tar.gz"))
(sha256
(base32
;; the hash below is produced by guix download <url>
"18100fd4vp41saxiji734mq5lckjplbnmm1nz29da4azhxzbzki9"
))))
(build-system gnu-build-system)
(arguments
'(#:configure-flags
'("--enable-silent-rules"
"--enable-hpc"
"--with-openmp")))
(inputs
`(("trexio", trexio)
("gfortran", gfortran)
("openblas", openblas)
("lapack", lapack)
))
(synopsis "QMCkl: Quantum Monte Carlo Kernel Library")
(description "The main QMC algorithms are exposed in a simple language and provided a standard API
and tests to enable the development of high-performance QMCkl implementations taking
advantage of modern hardware.")
(home-page "https://trex-coe.github.io/qmckl/index.html")
(license bsd-3)))
;; Guix does not copy the .git folder so relying on it's presence is a bad practice !
(define-public qmckl-dev
(let ((commit "26f8a1b906c329fa92adc2480e1769b8a90347de")
(revision "1"))
(package
(name "qmckl-dev")
(version (git-version "0.2.2" revision commit))
(source (origin
(method git-fetch)
(uri (git-reference
(url "https://github.com/TREX-CoE/qmckl")
(commit commit)
(recursive? #t))) ;; tried recursive to get htmlize - fails !
(file-name (git-file-name name version))
(sha256
(base32
;; the hash below is produced by `guix hash -rx .`
"0px3880bnciky8mwivll56108j9ncxri3ic2bhavcwn1z12z7lcb"
))))
(build-system gnu-build-system)
(arguments
'(#:configure-flags
'("--enable-hpc"
"--with-openmp")
;; ignoring make errors is a hack for the build to pass
;; #:make-flags '("-i")))
#:phases
;; this is a workaround to activate QMCKL_DEVEL mode
(modify-phases %standard-phases
(add-after 'unpack 'set_devel
(lambda _
(mkdir-p ".git")))
)))
(inputs
`(("trexio", trexio)
("gfortran", gfortran)
("openblas", openblas)
("lapack", lapack)
("python", python-wrapper)
("emacs", emacs)
))
;; these inputs are required for autogen.sh to work properly
(native-inputs
`(("autoconf", autoconf)
("automake", automake)
("libtool", libtool)
("pkg-config", pkg-config)
("swig", swig)
("m4", m4)
))
(synopsis "QMCkl: Quantum Monte Carlo Kernel Library")
(description "The main QMC algorithms are exposed in a simple language and provided a standard API
and tests to enable the development of high-performance QMCkl implementations taking
advantage of modern hardware.")
(home-page "https://trex-coe.github.io/qmckl/index.html")
(license bsd-3))))
(define-public qmckl
;; Default version of QMCkl - change this to benchmark installation from Git
;; qmckl-dev)
qmckl-hpc-0.2.1)
;; return qmckl variable so that `quix package -f qmckl.scm` works
qmckl