mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-08-17 10:51:44 +02:00
5530 lines
194 KiB
Org Mode
5530 lines
194 KiB
Org Mode
#+TITLE: Jastrow Factor
|
|
#+SETUPFILE: ../tools/theme.setup
|
|
#+INCLUDE: ../tools/lib.org
|
|
|
|
Functions for the calculation of the Jastrow factor \(f_{ee}, f_{en}, f_{een}\).
|
|
These are stored in the ~factor_ee~, ~factor_en~, and ~factor_een~ variables.
|
|
The ~jastrow~ structure contains all the information required to build
|
|
these factors along with their derivatives.
|
|
|
|
* Headers :noexport:
|
|
#+begin_src elisp :noexport :results none
|
|
(org-babel-lob-ingest "../tools/lib.org")
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :tangle (eval h_private_type)
|
|
#ifndef QMCKL_JASTROW_HPT
|
|
#define QMCKL_JASTROW_HPT
|
|
#include <stdbool.h>
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval c_test) :noweb yes
|
|
#include "qmckl.h"
|
|
#include <assert.h>
|
|
#include <math.h>
|
|
#ifdef HAVE_CONFIG_H
|
|
#include "config.h"
|
|
#endif
|
|
|
|
#include <stdio.h>
|
|
#include "n2.h"
|
|
|
|
int main() {
|
|
qmckl_context context;
|
|
context = qmckl_context_create();
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval c)
|
|
#ifdef HAVE_CONFIG_H
|
|
#include "config.h"
|
|
#endif
|
|
|
|
#ifdef HAVE_STDINT_H
|
|
#include <stdint.h>
|
|
#elif HAVE_INTTYPES_H
|
|
#include <inttypes.h>
|
|
#endif
|
|
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
#include <stdbool.h>
|
|
#include <assert.h>
|
|
#include <math.h>
|
|
|
|
#include <stdio.h>
|
|
|
|
#include "qmckl.h"
|
|
#include "qmckl_context_private_type.h"
|
|
#include "qmckl_memory_private_type.h"
|
|
#include "qmckl_memory_private_func.h"
|
|
#include "qmckl_jastrow_private_func.h"
|
|
#include "qmckl_jastrow_private_type.h"
|
|
#+end_src
|
|
|
|
* Context
|
|
:PROPERTIES:
|
|
:Name: qmckl_jastrow
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
The following data stored in the context:
|
|
|
|
#+NAME: qmckl_jastrow_args
|
|
|------------+--------------------------------------------+-----+-------------------------------------------------------------------|
|
|
| ~int32_t~ | ~uninitialized~ | in | Keeps bit set for uninitialized data |
|
|
| ~int64_t~ | ~aord_num~ | in | The number of a coeffecients |
|
|
| ~int64_t~ | ~bord_num~ | in | The number of b coeffecients |
|
|
| ~int64_t~ | ~cord_num~ | in | The number of c coeffecients |
|
|
| ~int64_t~ | ~type_nucl_num~ | in | Number of Nucleii types |
|
|
| ~int64_t~ | ~type_nucl_vector[nucl_num]~ | in | IDs of types of Nucleii |
|
|
| ~double~ | ~aord_vector[aord_num + 1][type_nucl_num]~ | in | Order of a polynomial coefficients |
|
|
| ~double~ | ~bord_vector[bord_num + 1]~ | in | Order of b polynomial coefficients |
|
|
| ~double~ | ~cord_vector[cord_num][type_nucl_num]~ | in | Order of c polynomial coefficients |
|
|
| ~double~ | ~factor_ee[walk_num]~ | out | Jastrow factor: electron-electron part |
|
|
| ~uint64_t~ | ~factor_ee_date~ | out | Jastrow factor: electron-electron part |
|
|
| ~double~ | ~factor_en[walk_num]~ | out | Jastrow factor: electron-nucleus part |
|
|
| ~uint64_t~ | ~factor_en_date~ | out | Jastrow factor: electron-nucleus part |
|
|
| ~double~ | ~factor_een[walk_num]~ | out | Jastrow factor: electron-electron-nucleus part |
|
|
| ~uint64_t~ | ~factor_een_date~ | out | Jastrow factor: electron-electron-nucleus part |
|
|
| ~double~ | ~factor_ee_deriv_e[4][nelec][walk_num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part |
|
|
| ~uint64_t~ | ~factor_ee_deriv_e_date~ | out | Keep track of the date for the derivative |
|
|
| ~double~ | ~factor_en_deriv_e[4][nelec][walk_num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part |
|
|
| ~uint64_t~ | ~factor_en_deriv_e_date~ | out | Keep track of the date for the en derivative |
|
|
| ~double~ | ~factor_een_deriv_e[4][nelec][walk_num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part |
|
|
| ~uint64_t~ | ~factor_een_deriv_e_date~ | out | Keep track of the date for the een derivative |
|
|
|
|
computed data:
|
|
|
|
|------------+-----------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------|
|
|
| ~int64_t~ | ~dim_cord_vect~ | Number of unique C coefficients |
|
|
| ~uint64_t~ | ~dim_cord_vect_date~ | Number of unique C coefficients |
|
|
| ~double~ | ~asymp_jasb[2]~ | Asymptotic component |
|
|
| ~uint64_t~ | ~asymp_jasb_date~ | Asymptotic component |
|
|
| ~double~ | ~cord_vect_full[dim_cord_vect][nucl_num]~ | vector of non-zero coefficients |
|
|
| ~uint64_t~ | ~cord_vect_full_date~ | Keep track of changes here |
|
|
| ~int64_t~ | ~lkpm_combined_index[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices |
|
|
| ~uint64_t~ | ~lkpm_combined_index_date~ | Transform l,k,p, and m into consecutive indices |
|
|
| ~double~ | ~tmp_c[elec_num][nucl_num][ncord + 1][ncord][walk_num]~ | vector of non-zero coefficients |
|
|
| ~double~ | ~dtmp_c[elec_num][4][nucl_num][ncord + 1][ncord][walk_num]~ | vector of non-zero coefficients |
|
|
| ~double~ | ~een_rescaled_e[walk_num][elec_num][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord |
|
|
| ~uint64_t~ | ~een_rescaled_e_date~ | Keep track of the date of creation |
|
|
| ~double~ | ~een_rescaled_n[walk_num][elec_num][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord |
|
|
| ~uint64_t~ | ~een_rescaled_n_date~ | Keep track of the date of creation |
|
|
| ~double~ | ~een_rescaled_e_deriv_e[walk_num][elec_num][4][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons |
|
|
| ~uint64_t~ | ~een_rescaled_e_deriv_e_date~ | Keep track of the date of creation |
|
|
| ~double~ | ~een_rescaled_n_deriv_e[walk_num][elec_num][4][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons |
|
|
| ~uint64_t~ | ~een_rescaled_n_deriv_e_date~ | Keep track of the date of creation |
|
|
|
|
For H2O we have the following data:
|
|
|
|
#+NAME: jastrow_data
|
|
#+BEGIN_SRC python :results output
|
|
import numpy as np
|
|
|
|
elec_num = 10
|
|
nucl_num = 2
|
|
up_num = 5
|
|
down_num = 5
|
|
nucl_coord = np.array([ [0.000000, 0.000000 ],
|
|
[0.000000, 0.000000 ],
|
|
[0.000000, 2.059801 ] ])
|
|
|
|
elec_coord = [[[-0.250655104764153 , 0.503070975550133 , -0.166554344502303],
|
|
[-0.587812193472177 , -0.128751981129274 , 0.187773606533075],
|
|
[ 1.61335569047166 , -0.615556732874863 , -1.43165470979934 ],
|
|
[-4.901239896295210E-003 , -1.120440036458986E-002 , 1.99761909330422 ],
|
|
[ 0.766647499681200 , -0.293515395797937 , 3.66454589201239 ],
|
|
[-0.127732483187947 , -0.138975497694196 , -8.669850480215846E-002],
|
|
[-0.232271834949124 , -1.059321673434182E-002 , -0.504862241464867],
|
|
[ 1.09360863531826 , -2.036103063808752E-003 , -2.702796910818986E-002],
|
|
[-0.108090166832043 , 0.189161729653261 , 2.15398313919894],
|
|
[ 0.397978144318712 , -0.254277292595981 , 2.54553335476344]]];
|
|
|
|
ee_distance_rescaled = [
|
|
[ 0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
|
|
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
|
|
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
|
|
0.000000000000000E+000],
|
|
[ 0.550227800352402 ,0.000000000000000E+000 ,0.000000000000000E+000,
|
|
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
|
|
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
|
|
0.000000000000000E+000],
|
|
[ 0.919155060185168 ,0.937695909123175 ,0.000000000000000E+000,
|
|
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
|
|
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
|
|
0.000000000000000E+000],
|
|
[ 0.893325429242815 ,0.851181978173561 ,0.978501685226877 ,
|
|
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
|
|
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
|
|
0.000000000000000E+000],
|
|
[ 0.982457268305353 ,0.976125002619471 ,0.994349933143149 ,
|
|
0.844077311588328 ,0.000000000000000E+000 ,0.000000000000000E+000,
|
|
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
|
|
0.000000000000000E+000],
|
|
[ 0.482407528408731 ,0.414816073699124 ,0.894716035479343 ,
|
|
0.876540187084407 ,0.978921170036895 ,0.000000000000000E+000,
|
|
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
|
|
0.000000000000000E+000],
|
|
[ 0.459541909660400 ,0.545007215761510 ,0.883752955884551 ,
|
|
0.918958134888791 ,0.986386936267237 ,0.362209822236419 ,
|
|
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
|
|
0.000000000000000E+000],
|
|
[ 0.763732576854455 ,0.817282762358449 ,0.801802919535959 ,
|
|
0.900089095449775 ,0.975704636491453 ,0.707836537586060 ,
|
|
0.755705808346586 ,0.000000000000000E+000 ,0.000000000000000E+000,
|
|
0.000000000000000E+000],
|
|
[ 0.904249454052971 ,0.871097965261373 ,0.982717262706270 ,
|
|
0.239901207363622 ,0.836519456769083 ,0.896135326270534 ,
|
|
0.930694340243023 ,0.917708540815567 ,0.000000000000000E+000,
|
|
0.000000000000000E+000],
|
|
[ 0.944400908070716 ,0.922589018494961 ,0.984615718580670 ,
|
|
0.514328661540623 ,0.692362267147064 ,0.931894098453677 ,
|
|
0.956034127544344 ,0.931221472309472 ,0.540903688625053 ,
|
|
0.000000000000000E+000]]
|
|
|
|
en_distance_rescaled = np.transpose(np.array([
|
|
[ 0.443570948411811 , 0.467602196999105 , 0.893870160799932 ,
|
|
0.864347190364447 , 0.976608182392358 , 0.187563183468210 ,
|
|
0.426404699872689 , 0.665107090128166 , 0.885246991424583 ,
|
|
0.924902909715270 ],
|
|
[ 0.899360150637444 , 0.860035135365386 , 0.979659405613798 ,
|
|
6.140678415933776E-002, 0.835118398056681 , 0.884071658981068 ,
|
|
0.923860000907362 , 0.905203414522289 , 0.211286300932359 ,
|
|
0.492104840907350 ]]))
|
|
|
|
# symmetrize it
|
|
for i in range(elec_num):
|
|
for j in range(elec_num):
|
|
ee_distance_rescaled[i][j] = ee_distance_rescaled[j][i]
|
|
|
|
type_nucl_num = 1
|
|
aord_num = 5
|
|
bord_num = 5
|
|
cord_num = 23
|
|
dim_cord_vect= 23
|
|
type_nucl_vector = [ 1, 1]
|
|
aord_vector = [
|
|
[0.000000000000000E+000],
|
|
[0.000000000000000E+000],
|
|
[-0.380512000000000E+000],
|
|
[-0.157996000000000E+000],
|
|
[-3.155800000000000E-002],
|
|
[2.151200000000000E-002]]
|
|
|
|
bord_vector = [ 0.500000000000000E-000, 0.153660000000000E-000, 6.722620000000000E-002,
|
|
2.157000000000000E-002, 7.309600000000000E-003, 2.866000000000000E-003]
|
|
cord_vector = [ 0.571702000000000E-000, -0.514253000000000E-000, -0.513043000000000E-000,
|
|
9.486000000000000E-003, -4.205000000000000E-003, 0.426325800000000E-000,
|
|
8.288150000000000E-002, 5.118600000000000E-003, -2.997800000000000E-003,
|
|
-5.270400000000000E-003, -7.499999999999999E-005, -8.301649999999999E-002,
|
|
1.454340000000000E-002, 5.143510000000000E-002, 9.250000000000000E-004,
|
|
-4.099100000000000E-003, 4.327600000000000E-003, -1.654470000000000E-003,
|
|
2.614000000000000E-003, -1.477000000000000E-003, -1.137000000000000E-003,
|
|
-4.010475000000000E-002, 6.106710000000000E-003 ]
|
|
cord_vector_full = [
|
|
[ 0.571702000000000E-000, -0.514253000000000E-000, -0.513043000000000E-000,
|
|
9.486000000000000E-003, -4.205000000000000E-003, 0.426325800000000E-000,
|
|
8.288150000000000E-002, 5.118600000000000E-003, -2.997800000000000E-003,
|
|
-5.270400000000000E-003, -7.499999999999999E-005, -8.301649999999999E-002,
|
|
1.454340000000000E-002, 5.143510000000000E-002, 9.250000000000000E-004,
|
|
-4.099100000000000E-003, 4.327600000000000E-003, -1.654470000000000E-003,
|
|
2.614000000000000E-003, -1.477000000000000E-003, -1.137000000000000E-003,
|
|
-4.010475000000000E-002, 6.106710000000000E-003 ],
|
|
[ 0.571702000000000E-000, -0.514253000000000E-000, -0.513043000000000E-000,
|
|
9.486000000000000E-003, -4.205000000000000E-003, 0.426325800000000E-000,
|
|
8.288150000000000E-002, 5.118600000000000E-003, -2.997800000000000E-003,
|
|
-5.270400000000000E-003, -7.499999999999999E-005, -8.301649999999999E-002,
|
|
1.454340000000000E-002, 5.143510000000000E-002, 9.250000000000000E-004,
|
|
-4.099100000000000E-003, 4.327600000000000E-003, -1.654470000000000E-003,
|
|
2.614000000000000E-003, -1.477000000000000E-003, -1.137000000000000E-003,
|
|
-4.010475000000000E-002, 6.106710000000000E-003 ],
|
|
]
|
|
lkpm_combined_index = [[1 , 1 , 2 , 0],
|
|
[0 , 0 , 2 , 1],
|
|
[1 , 2 , 3 , 0],
|
|
[2 , 1 , 3 , 0],
|
|
[0 , 1 , 3 , 1],
|
|
[1 , 0 , 3 , 1],
|
|
[1 , 3 , 4 , 0],
|
|
[2 , 2 , 4 , 0],
|
|
[0 , 2 , 4 , 1],
|
|
[3 , 1 , 4 , 0],
|
|
[1 , 1 , 4 , 1],
|
|
[2 , 0 , 4 , 1],
|
|
[0 , 0 , 4 , 2],
|
|
[1 , 4 , 5 , 0],
|
|
[2 , 3 , 5 , 0],
|
|
[0 , 3 , 5 , 1],
|
|
[3 , 2 , 5 , 0],
|
|
[1 , 2 , 5 , 1],
|
|
[4 , 1 , 5 , 0],
|
|
[2 , 1 , 5 , 1],
|
|
[0 , 1 , 5 , 2],
|
|
[3 , 0 , 5 , 1],
|
|
[1 , 0 , 5 , 2]]
|
|
|
|
kappa = 1.0
|
|
kappa_inv = 1.0/kappa
|
|
#+END_SRC
|
|
|
|
#+RESULTS: jastrow_data
|
|
|
|
** Data structure
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_type)
|
|
typedef struct qmckl_jastrow_struct{
|
|
int32_t uninitialized;
|
|
int64_t aord_num;
|
|
int64_t bord_num;
|
|
int64_t cord_num;
|
|
int64_t type_nucl_num;
|
|
uint64_t asymp_jasb_date;
|
|
uint64_t tmp_c_date;
|
|
uint64_t dtmp_c_date;
|
|
uint64_t factor_ee_date;
|
|
uint64_t factor_en_date;
|
|
uint64_t factor_een_date;
|
|
uint64_t factor_ee_deriv_e_date;
|
|
uint64_t factor_en_deriv_e_date;
|
|
uint64_t factor_een_deriv_e_date;
|
|
int64_t* type_nucl_vector;
|
|
double * aord_vector;
|
|
double * bord_vector;
|
|
double * cord_vector;
|
|
double * asymp_jasb;
|
|
double * factor_ee;
|
|
double * factor_en;
|
|
double * factor_een;
|
|
double * factor_ee_deriv_e;
|
|
double * factor_en_deriv_e;
|
|
double * factor_een_deriv_e;
|
|
int64_t dim_cord_vect;
|
|
uint64_t dim_cord_vect_date;
|
|
double * cord_vect_full;
|
|
uint64_t cord_vect_full_date;
|
|
int64_t* lkpm_combined_index;
|
|
uint64_t lkpm_combined_index_date;
|
|
double * tmp_c;
|
|
double * dtmp_c;
|
|
double * een_rescaled_e;
|
|
double * een_rescaled_n;
|
|
uint64_t een_rescaled_e_date;
|
|
uint64_t een_rescaled_n_date;
|
|
double * een_rescaled_e_deriv_e;
|
|
double * een_rescaled_n_deriv_e;
|
|
uint64_t een_rescaled_e_deriv_e_date;
|
|
uint64_t een_rescaled_n_deriv_e_date;
|
|
bool provided;
|
|
char * type;
|
|
} qmckl_jastrow_struct;
|
|
#+end_src
|
|
|
|
|
|
The ~uninitialized~ integer contains one bit set to one for each
|
|
initialization function which has not been called. It becomes equal
|
|
to zero after all initialization functions have been called. The
|
|
struct is then initialized and ~provided == true~.
|
|
Some values are initialized by default, and are not concerned by
|
|
this mechanism.
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code qmckl_init_jastrow(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c)
|
|
qmckl_exit_code qmckl_init_jastrow(qmckl_context context) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return false;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
ctx->jastrow.uninitialized = (1 << 6) - 1;
|
|
|
|
/* Default values */
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
** Access functions
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :exports none
|
|
qmckl_exit_code qmckl_get_jastrow_aord_num (qmckl_context context, int64_t* const aord_num);
|
|
qmckl_exit_code qmckl_get_jastrow_bord_num (qmckl_context context, int64_t* const bord_num);
|
|
qmckl_exit_code qmckl_get_jastrow_cord_num (qmckl_context context, int64_t* const bord_num);
|
|
qmckl_exit_code qmckl_get_jastrow_type_nucl_num (qmckl_context context, int64_t* const type_nucl_num);
|
|
qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (qmckl_context context, int64_t* const type_nucl_num);
|
|
qmckl_exit_code qmckl_get_jastrow_aord_vector (qmckl_context context, double * const aord_vector);
|
|
qmckl_exit_code qmckl_get_jastrow_bord_vector (qmckl_context context, double * const bord_vector);
|
|
qmckl_exit_code qmckl_get_jastrow_cord_vector (qmckl_context context, double * const cord_vector);
|
|
#+end_src
|
|
|
|
Along with these core functions, calculation of the jastrow factor
|
|
requires the following additional information to be set:
|
|
|
|
|
|
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_jastrow_provided (const qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
bool qmckl_jastrow_provided(const qmckl_context context) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return false;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
return ctx->jastrow.provided;
|
|
}
|
|
#+end_src
|
|
|
|
#+NAME:post
|
|
#+begin_src c :exports none
|
|
if ( (ctx->jastrow.uninitialized & mask) != 0) {
|
|
return NULL;
|
|
}
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_jastrow_aord_num (const qmckl_context context, int64_t* const aord_num) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return (char) 0;
|
|
}
|
|
|
|
if (aord_num == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_jastrow_aord_num",
|
|
"aord_num is a null pointer");
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 0;
|
|
|
|
if ( (ctx->jastrow.uninitialized & mask) != 0) {
|
|
return QMCKL_NOT_PROVIDED;
|
|
}
|
|
|
|
assert (ctx->jastrow.aord_num > 0);
|
|
*aord_num = ctx->jastrow.aord_num;
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
qmckl_exit_code qmckl_get_jastrow_bord_num (const qmckl_context context, int64_t* const bord_num) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return (char) 0;
|
|
}
|
|
|
|
if (bord_num == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_jastrow_bord_num",
|
|
"aord_num is a null pointer");
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 0;
|
|
|
|
if ( (ctx->jastrow.uninitialized & mask) != 0) {
|
|
return QMCKL_NOT_PROVIDED;
|
|
}
|
|
|
|
assert (ctx->jastrow.bord_num > 0);
|
|
*bord_num = ctx->jastrow.bord_num;
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
qmckl_exit_code qmckl_get_jastrow_cord_num (const qmckl_context context, int64_t* const cord_num) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return (char) 0;
|
|
}
|
|
|
|
if (cord_num == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_jastrow_cord_num",
|
|
"aord_num is a null pointer");
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 0;
|
|
|
|
if ( (ctx->jastrow.uninitialized & mask) != 0) {
|
|
return QMCKL_NOT_PROVIDED;
|
|
}
|
|
|
|
assert (ctx->jastrow.cord_num > 0);
|
|
*cord_num = ctx->jastrow.cord_num;
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
qmckl_exit_code qmckl_get_jastrow_type_nucl_num (const qmckl_context context, int64_t* const type_nucl_num) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return (char) 0;
|
|
}
|
|
|
|
if (type_nucl_num == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_jastrow_type_nucl_num",
|
|
"type_nucl_num is a null pointer");
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 1;
|
|
|
|
if ( (ctx->jastrow.uninitialized & mask) != 0) {
|
|
return QMCKL_NOT_PROVIDED;
|
|
}
|
|
|
|
assert (ctx->jastrow.type_nucl_num > 0);
|
|
*type_nucl_num = ctx->jastrow.type_nucl_num;
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (const qmckl_context context, int64_t * const type_nucl_vector) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return (char) 0;
|
|
}
|
|
|
|
if (type_nucl_vector == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_jastrow_type_nucl_vector",
|
|
"type_nucl_vector is a null pointer");
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 2;
|
|
|
|
if ( (ctx->jastrow.uninitialized & mask) != 0) {
|
|
return QMCKL_NOT_PROVIDED;
|
|
}
|
|
|
|
assert (ctx->jastrow.type_nucl_vector != NULL);
|
|
memcpy(type_nucl_vector, ctx->jastrow.type_nucl_vector, ctx->jastrow.type_nucl_num*sizeof(int64_t));
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
qmckl_exit_code qmckl_get_jastrow_aord_vector (const qmckl_context context, double * const aord_vector) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return (char) 0;
|
|
}
|
|
|
|
if (aord_vector == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_jastrow_aord_vector",
|
|
"aord_vector is a null pointer");
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 3;
|
|
|
|
if ( (ctx->jastrow.uninitialized & mask) != 0) {
|
|
return QMCKL_NOT_PROVIDED;
|
|
}
|
|
|
|
assert (ctx->jastrow.aord_vector != NULL);
|
|
memcpy(aord_vector, ctx->jastrow.aord_vector, ctx->jastrow.aord_num*sizeof(double));
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
qmckl_exit_code qmckl_get_jastrow_bord_vector (const qmckl_context context, double * const bord_vector) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return (char) 0;
|
|
}
|
|
|
|
if (bord_vector == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_jastrow_bord_vector",
|
|
"bord_vector is a null pointer");
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 4;
|
|
|
|
if ( (ctx->jastrow.uninitialized & mask) != 0) {
|
|
return QMCKL_NOT_PROVIDED;
|
|
}
|
|
|
|
assert (ctx->jastrow.bord_vector != NULL);
|
|
memcpy(bord_vector, ctx->jastrow.bord_vector, ctx->jastrow.bord_num*sizeof(double));
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, double * const cord_vector) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return (char) 0;
|
|
}
|
|
|
|
if (cord_vector == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_get_jastrow_cord_vector",
|
|
"cord_vector is a null pointer");
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int32_t mask = 1 << 5;
|
|
|
|
if ( (ctx->jastrow.uninitialized & mask) != 0) {
|
|
return QMCKL_NOT_PROVIDED;
|
|
}
|
|
|
|
assert (ctx->jastrow.cord_vector != NULL);
|
|
memcpy(cord_vector, ctx->jastrow.cord_vector, ctx->jastrow.cord_num*sizeof(double));
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
#+end_src
|
|
|
|
** Initialization functions
|
|
|
|
To prepare for the Jastrow and its derivative, all the following functions need to be
|
|
called.
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func)
|
|
qmckl_exit_code qmckl_set_jastrow_ord_num (qmckl_context context, const int64_t aord_num, const int64_t bord_num, const int64_t cord_num);
|
|
qmckl_exit_code qmckl_set_jastrow_type_nucl_num (qmckl_context context, const int64_t type_nucl_num);
|
|
qmckl_exit_code qmckl_set_jastrow_type_nucl_vector (qmckl_context context, const int64_t* type_nucl_vector, const int64_t nucl_num);
|
|
qmckl_exit_code qmckl_set_jastrow_aord_vector (qmckl_context context, const double * aord_vector);
|
|
qmckl_exit_code qmckl_set_jastrow_bord_vector (qmckl_context context, const double * bord_vector);
|
|
qmckl_exit_code qmckl_set_jastrow_cord_vector (qmckl_context context, const double * cord_vector);
|
|
qmckl_exit_code qmckl_set_jastrow_dependencies (qmckl_context context);
|
|
#+end_src
|
|
|
|
#+NAME:pre2
|
|
#+begin_src c :exports none
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
#+end_src
|
|
|
|
#+NAME:post2
|
|
#+begin_src c :exports none
|
|
ctx->jastrow.uninitialized &= ~mask;
|
|
ctx->jastrow.provided = (ctx->jastrow.uninitialized == 0);
|
|
if (ctx->jastrow.provided) {
|
|
//qmckl_exit_code rc_ = qmckl_set_jastrow_dependencies(context);
|
|
//if (rc_ != QMCKL_SUCCESS) return rc_;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
#+end_src
|
|
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_set_jastrow_ord_num(qmckl_context context, const int64_t aord_num, const int64_t bord_num, const int64_t cord_num) {
|
|
<<pre2>>
|
|
|
|
if (aord_num <= 0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_jastrow_ord_num",
|
|
"aord_num <= 0");
|
|
}
|
|
|
|
if (bord_num <= 0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_jastrow_ord_num",
|
|
"bord_num <= 0");
|
|
}
|
|
|
|
if (cord_num <= 0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_jastrow_ord_num",
|
|
"cord_num <= 0");
|
|
}
|
|
|
|
int32_t mask = 1 << 0;
|
|
ctx->jastrow.aord_num = aord_num;
|
|
ctx->jastrow.bord_num = bord_num;
|
|
ctx->jastrow.cord_num = cord_num;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
qmckl_exit_code qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int64_t type_nucl_num) {
|
|
<<pre2>>
|
|
|
|
if (type_nucl_num <= 0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_jastrow_type_nucl_num",
|
|
"type_nucl_num < 0");
|
|
}
|
|
|
|
int32_t mask = 1 << 1;
|
|
ctx->jastrow.type_nucl_num = type_nucl_num;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
qmckl_exit_code qmckl_set_jastrow_type_nucl_vector(qmckl_context context, int64_t const * type_nucl_vector, const int64_t nucl_num) {
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 2;
|
|
|
|
int64_t type_nucl_num;
|
|
qmckl_exit_code rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
if (type_nucl_num == 0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_jastrow_type_nucl_vector",
|
|
"type_nucl_num is not set");
|
|
}
|
|
|
|
if (type_nucl_vector == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_jastrow_type_nucl_vector",
|
|
"type_nucl_vector = NULL");
|
|
}
|
|
|
|
if (ctx->jastrow.type_nucl_vector != NULL) {
|
|
qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.type_nucl_vector);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc,
|
|
"qmckl_set_type_nucl_vector",
|
|
NULL);
|
|
}
|
|
}
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = nucl_num * sizeof(int64_t);
|
|
int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);
|
|
|
|
if(new_array == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_set_jastrow_type_nucl_vector",
|
|
NULL);
|
|
}
|
|
|
|
memcpy(new_array, type_nucl_vector, mem_info.size);
|
|
|
|
ctx->jastrow.type_nucl_vector = new_array;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double const * aord_vector) {
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 3;
|
|
|
|
int64_t aord_num;
|
|
qmckl_exit_code rc = qmckl_get_jastrow_aord_num(context, &aord_num);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
int64_t type_nucl_num;
|
|
rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
if (aord_num == 0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_jastrow_coefficient",
|
|
"aord_num is not set");
|
|
}
|
|
|
|
if (aord_vector == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_jastrow_aord_vector",
|
|
"aord_vector = NULL");
|
|
}
|
|
|
|
if (ctx->jastrow.aord_vector != NULL) {
|
|
qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.aord_vector);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc,
|
|
"qmckl_set_ord_vector",
|
|
NULL);
|
|
}
|
|
}
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = (aord_num + 1) * type_nucl_num * sizeof(double);
|
|
double* new_array = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if(new_array == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_set_jastrow_coefficient",
|
|
NULL);
|
|
}
|
|
|
|
memcpy(new_array, aord_vector, mem_info.size);
|
|
|
|
ctx->jastrow.aord_vector = new_array;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double const * bord_vector) {
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 4;
|
|
|
|
int64_t bord_num;
|
|
qmckl_exit_code rc = qmckl_get_jastrow_bord_num(context, &bord_num);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
if (bord_num == 0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_jastrow_coefficient",
|
|
"bord_num is not set");
|
|
}
|
|
|
|
if (bord_vector == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_jastrow_bord_vector",
|
|
"bord_vector = NULL");
|
|
}
|
|
|
|
if (ctx->jastrow.bord_vector != NULL) {
|
|
qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.bord_vector);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc,
|
|
"qmckl_set_ord_vector",
|
|
NULL);
|
|
}
|
|
}
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = (bord_num + 1) * sizeof(double);
|
|
double* new_array = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if(new_array == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_set_jastrow_coefficient",
|
|
NULL);
|
|
}
|
|
|
|
memcpy(new_array, bord_vector, mem_info.size);
|
|
|
|
ctx->jastrow.bord_vector = new_array;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double const * cord_vector) {
|
|
<<pre2>>
|
|
|
|
int32_t mask = 1 << 5;
|
|
|
|
int64_t cord_num;
|
|
qmckl_exit_code rc = qmckl_get_jastrow_cord_num(context, &cord_num);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
int64_t type_nucl_num;
|
|
rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
if (cord_num == 0) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_FAILURE,
|
|
"qmckl_set_jastrow_coefficient",
|
|
"cord_num is not set");
|
|
}
|
|
|
|
if (cord_vector == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_jastrow_cord_vector",
|
|
"cord_vector = NULL");
|
|
}
|
|
|
|
if (ctx->jastrow.cord_vector != NULL) {
|
|
qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.cord_vector);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return qmckl_failwith( context, rc,
|
|
"qmckl_set_ord_vector",
|
|
NULL);
|
|
}
|
|
}
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = cord_num * type_nucl_num * sizeof(double);
|
|
double* new_array = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if(new_array == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_set_jastrow_coefficient",
|
|
NULL);
|
|
}
|
|
|
|
memcpy(new_array, cord_vector, mem_info.size);
|
|
|
|
ctx->jastrow.cord_vector = new_array;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
qmckl_exit_code qmckl_set_jastrow_dependencies(qmckl_context context) {
|
|
<<pre2>>
|
|
|
|
/* Check for electron data */
|
|
if (!(ctx->electron.provided)) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_provide_ee_distance",
|
|
NULL);
|
|
}
|
|
|
|
/* Check for nucleus data */
|
|
if (!(ctx->nucleus.provided)) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_provide_en_distance",
|
|
NULL);
|
|
}
|
|
|
|
int32_t mask = 1 << 6;
|
|
|
|
<<post2>>
|
|
}
|
|
|
|
#+end_src
|
|
|
|
When the required information is completely entered, other data structures are
|
|
computed to accelerate the calculations. The intermediates factors
|
|
are precontracted using BLAS LEVEL 3 operations for an optimal FLOP count.
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_INVALID_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
/* ----------------------------------- */
|
|
/* Check for the necessary information */
|
|
/* ----------------------------------- */
|
|
|
|
/* Check for the electron data
|
|
1. elec_num
|
|
2. ee_distances_rescaled
|
|
*/
|
|
if (!(ctx->electron.provided)) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_electron",
|
|
NULL);
|
|
}
|
|
|
|
/* Check for the nucleus data
|
|
1. nucl_num
|
|
2. en_distances_rescaled
|
|
*/
|
|
if (!(ctx->nucleus.provided)) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_NOT_PROVIDED,
|
|
"qmckl_nucleus",
|
|
NULL);
|
|
}
|
|
|
|
qmckl_exit_code rc = QMCKL_FAILURE;
|
|
|
|
|
|
/* ----------------------------------- */
|
|
/* Start calculation of data */
|
|
/* ----------------------------------- */
|
|
|
|
|
|
}
|
|
#+end_src
|
|
|
|
** Test
|
|
#+begin_src c :tangle (eval c_test)
|
|
/* Reference input data */
|
|
int64_t walk_num = n2_walk_num;
|
|
int64_t elec_num = n2_elec_num;
|
|
int64_t elec_up_num = n2_elec_up_num;
|
|
int64_t elec_dn_num = n2_elec_dn_num;
|
|
double rescale_factor_kappa_ee = 1.0;
|
|
double rescale_factor_kappa_en = 1.0;
|
|
double nucl_rescale_factor_kappa = 1.0;
|
|
double* elec_coord = &(n2_elec_coord[0][0][0]);
|
|
|
|
const double* nucl_charge = n2_charge;
|
|
int64_t nucl_num = n2_nucl_num;
|
|
double* charge = n2_charge;
|
|
double* nucl_coord = &(n2_nucl_coord[0][0]);
|
|
|
|
/* Provide Electron data */
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
assert(!qmckl_electron_provided(context));
|
|
|
|
int64_t n;
|
|
rc = qmckl_get_electron_num (context, &n);
|
|
assert(rc == QMCKL_NOT_PROVIDED);
|
|
|
|
rc = qmckl_get_electron_up_num (context, &n);
|
|
assert(rc == QMCKL_NOT_PROVIDED);
|
|
|
|
rc = qmckl_get_electron_down_num (context, &n);
|
|
assert(rc == QMCKL_NOT_PROVIDED);
|
|
|
|
|
|
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(!qmckl_electron_provided(context));
|
|
|
|
rc = qmckl_get_electron_up_num (context, &n);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(n == elec_up_num);
|
|
|
|
rc = qmckl_get_electron_down_num (context, &n);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(n == elec_dn_num);
|
|
|
|
rc = qmckl_get_electron_num (context, &n);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(n == elec_num);
|
|
|
|
double k_ee = 0.;
|
|
double k_en = 0.;
|
|
rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(k_ee == 1.0);
|
|
|
|
rc = qmckl_get_electron_rescale_factor_en (context, &k_en);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(k_en == 1.0);
|
|
|
|
rc = qmckl_set_electron_rescale_factor_en(context, rescale_factor_kappa_en);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_set_electron_rescale_factor_ee(context, rescale_factor_kappa_ee);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(k_ee == rescale_factor_kappa_ee);
|
|
|
|
rc = qmckl_get_electron_rescale_factor_en (context, &k_en);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(k_en == rescale_factor_kappa_en);
|
|
|
|
|
|
int64_t w;
|
|
rc = qmckl_get_electron_walk_num (context, &w);
|
|
assert(rc == QMCKL_NOT_PROVIDED);
|
|
|
|
|
|
rc = qmckl_set_electron_walk_num (context, walk_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_get_electron_walk_num (context, &w);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(w == walk_num);
|
|
|
|
assert(qmckl_electron_provided(context));
|
|
|
|
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
double elec_coord2[walk_num*3*elec_num];
|
|
|
|
rc = qmckl_get_electron_coord (context, 'N', elec_coord2);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
for (int64_t i=0 ; i<3*elec_num ; ++i) {
|
|
assert( elec_coord[i] == elec_coord2[i] );
|
|
}
|
|
|
|
|
|
/* Provide Nucleus data */
|
|
|
|
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);
|
|
assert(k == 1.0);
|
|
|
|
|
|
rc = qmckl_set_nucleus_rescale_factor (context, nucl_rescale_factor_kappa);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_get_nucleus_rescale_factor (context, &k);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
assert(k == nucl_rescale_factor_kappa);
|
|
|
|
double nucl_coord2[3*nucl_num];
|
|
|
|
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2);
|
|
assert(rc == QMCKL_NOT_PROVIDED);
|
|
|
|
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
assert(!qmckl_nucleus_provided(context));
|
|
|
|
rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
for (int64_t k=0 ; k<3 ; ++k) {
|
|
for (int64_t i=0 ; i<nucl_num ; ++i) {
|
|
assert( nucl_coord[nucl_num*k+i] == nucl_coord2[3*i+k] );
|
|
}
|
|
}
|
|
|
|
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
for (int64_t i=0 ; i<3*nucl_num ; ++i) {
|
|
assert( nucl_coord[i] == nucl_coord2[i] );
|
|
}
|
|
|
|
double nucl_charge2[nucl_num];
|
|
|
|
rc = qmckl_get_nucleus_charge(context, nucl_charge2);
|
|
assert(rc == QMCKL_NOT_PROVIDED);
|
|
|
|
rc = qmckl_set_nucleus_charge(context, nucl_charge);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
rc = qmckl_get_nucleus_charge(context, nucl_charge2);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
for (int64_t i=0 ; i<nucl_num ; ++i) {
|
|
assert( nucl_charge[i] == nucl_charge2[i] );
|
|
}
|
|
assert(qmckl_nucleus_provided(context));
|
|
|
|
#+end_src
|
|
|
|
* Computation
|
|
|
|
The computed data is stored in the context so that it can be reused
|
|
by different kernels. To ensure that the data is valid, for each
|
|
computed data the date of the context is stored when it is computed.
|
|
To know if some data needs to be recomputed, we check if the date of
|
|
the dependencies are more recent than the date of the data to
|
|
compute. If it is the case, then the data is recomputed and the
|
|
current date is stored.
|
|
|
|
** Asymptotic component for \(f_{ee}\)
|
|
|
|
Calculate the asymptotic component ~asymp_jasb~ to be substracted from the final
|
|
electron-electron jastrow factor \(f_{ee}\). The asymptotic componenet is calculated
|
|
via the ~bord_vector~ and the electron-electron rescale factor ~rescale_factor_kappa~.
|
|
|
|
\[
|
|
J_{asymp} = \frac{b_1 \kappa^-1}{1 + b_2 \kappa^-1}
|
|
\]
|
|
|
|
*** Get
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_jastrow_asymp_jasb(qmckl_context context, double* const asymp_jasb);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_jastrow_asymp_jasb(qmckl_context context, double* const asymp_jasb)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_asymp_jasb(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = 2;
|
|
memcpy(asymp_jasb, ctx->jastrow.asymp_jasb, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_asymp_jasb(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_asymp_jasb(qmckl_context context)
|
|
{
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
/* Check if ee kappa is provided */
|
|
double rescale_factor_kappa_ee;
|
|
rc = qmckl_get_electron_rescale_factor_ee(context, &rescale_factor_kappa_ee);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->date > ctx->jastrow.asymp_jasb_date) {
|
|
|
|
/* Allocate array */
|
|
if (ctx->jastrow.asymp_jasb == NULL) {
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = 2 * sizeof(double);
|
|
double* asymp_jasb = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (asymp_jasb == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_asymp_jasb",
|
|
NULL);
|
|
}
|
|
ctx->jastrow.asymp_jasb = asymp_jasb;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_asymp_jasb(context,
|
|
ctx->jastrow.bord_num,
|
|
ctx->jastrow.bord_vector,
|
|
rescale_factor_kappa_ee,
|
|
ctx->jastrow.asymp_jasb);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->jastrow.asymp_jasb_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_asymp_jasb
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_asymp_jasb_args
|
|
| qmckl_context | context | in | Global state |
|
|
| int64_t | bord_num | in | Number of electrons |
|
|
| double | bord_vector[bord_num + 1] | in | Number of walkers |
|
|
| double | rescale_factor_kappa_ee | in | Electron coordinates |
|
|
| double | asymp_jasb[2] | out | Electron-electron distances |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_asymp_jasb_f(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: bord_num
|
|
double precision , intent(in) :: bord_vector(bord_num + 1)
|
|
double precision , intent(in) :: rescale_factor_kappa_ee
|
|
double precision , intent(out) :: asymp_jasb(2)
|
|
|
|
integer*8 :: i, p
|
|
double precision :: kappa_inv, x, asym_one
|
|
kappa_inv = 1.0d0 / rescale_factor_kappa_ee
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (bord_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
asym_one = bord_vector(1) * kappa_inv / (1.0d0 + bord_vector(2) * kappa_inv)
|
|
asymp_jasb(:) = (/asym_one, 0.5d0 * asym_one/)
|
|
|
|
do i = 1, 2
|
|
x = kappa_inv
|
|
do p = 2, bord_num
|
|
x = x * kappa_inv
|
|
asymp_jasb(i) = asymp_jasb(i) + bord_vector(p + 1) * x
|
|
end do
|
|
end do
|
|
|
|
end function qmckl_compute_asymp_jasb_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_header(table=qmckl_asymp_jasb_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_private_func) :comments org
|
|
qmckl_exit_code qmckl_compute_asymp_jasb (
|
|
const qmckl_context context,
|
|
const int64_t bord_num,
|
|
const double* bord_vector,
|
|
const double rescale_factor_kappa_ee,
|
|
double* const asymp_jasb );
|
|
#+end_src
|
|
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_asymp_jasb_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_asymp_jasb &
|
|
(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: bord_num
|
|
real (c_double ) , intent(in) :: bord_vector(bord_num + 1)
|
|
real (c_double ) , intent(in) , value :: rescale_factor_kappa_ee
|
|
real (c_double ) , intent(out) :: asymp_jasb(2)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_asymp_jasb_f
|
|
info = qmckl_compute_asymp_jasb_f &
|
|
(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb)
|
|
|
|
end function qmckl_compute_asymp_jasb
|
|
#+end_src
|
|
|
|
*** Test
|
|
#+name: asymp_jasb
|
|
#+begin_src python :results output :exports none :noweb yes
|
|
import numpy as np
|
|
|
|
<<jastrow_data>>
|
|
|
|
asym_one = bord_vector[0] * kappa_inv / (1.0 + bord_vector[1]*kappa_inv)
|
|
asymp_jasb = np.array([asym_one, 0.5 * asym_one])
|
|
|
|
for i in range(2):
|
|
x = kappa_inv
|
|
for p in range(1,bord_num):
|
|
x = x * kappa_inv
|
|
asymp_jasb[i] += bord_vector[p + 1] * x
|
|
|
|
print("asym_one : ", asym_one)
|
|
print("asymp_jasb[0] : ", asymp_jasb[0])
|
|
print("asymp_jasb[1] : ", asymp_jasb[1])
|
|
#+end_src
|
|
|
|
#+RESULTS: asymp_jasb
|
|
: asym_one : 0.6634291325000664
|
|
: asymp_jasb[0] : 1.043287918508297
|
|
: asymp_jasb[1] : 0.7115733522582638
|
|
|
|
#+RESULTS:
|
|
: asym_one : 0.43340325572525706
|
|
: asymp_jasb[0] : 0.5323750557252571
|
|
: asymp_jasb[1] : 0.31567342786262853
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
assert(qmckl_electron_provided(context));
|
|
|
|
int64_t type_nucl_num = n2_type_nucl_num;
|
|
int64_t* type_nucl_vector = &(n2_type_nucl_vector[0]);
|
|
int64_t aord_num = n2_aord_num;
|
|
int64_t bord_num = n2_bord_num;
|
|
int64_t cord_num = n2_cord_num;
|
|
double* aord_vector = &(n2_aord_vector[0][0]);
|
|
double* bord_vector = &(n2_bord_vector[0]);
|
|
double* cord_vector = &(n2_cord_vector[0][0]);
|
|
|
|
/* Initialize the Jastrow data */
|
|
rc = qmckl_init_jastrow(context);
|
|
assert(!qmckl_jastrow_provided(context));
|
|
|
|
/* Set the data */
|
|
rc = qmckl_set_jastrow_ord_num(context, aord_num, bord_num, cord_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
rc = qmckl_set_jastrow_type_nucl_num(context, type_nucl_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
rc = qmckl_set_jastrow_type_nucl_vector(context, type_nucl_vector, nucl_num);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
rc = qmckl_set_jastrow_aord_vector(context, aord_vector);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
rc = qmckl_set_jastrow_bord_vector(context, bord_vector);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
rc = qmckl_set_jastrow_cord_vector(context, cord_vector);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
rc = qmckl_set_jastrow_dependencies(context);
|
|
assert(rc == QMCKL_SUCCESS);
|
|
|
|
/* Check if Jastrow is properly initialized */
|
|
assert(qmckl_jastrow_provided(context));
|
|
|
|
double asymp_jasb[2];
|
|
rc = qmckl_get_jastrow_asymp_jasb(context, asymp_jasb);
|
|
|
|
// calculate asymp_jasb
|
|
assert(fabs(asymp_jasb[0]-0.5323750557252571) < 1.e-12);
|
|
assert(fabs(asymp_jasb[1]-0.31567342786262853) < 1.e-12);
|
|
|
|
#+end_src
|
|
|
|
** Electron-electron component \(f_{ee}\)
|
|
|
|
Calculate the electron-electron jastrow component ~factor_ee~ using the ~asymp_jasb~
|
|
componenet and the electron-electron rescaled distances ~ee_distance_rescaled~.
|
|
|
|
\[
|
|
f_{ee} = \sum_{i,j<i} \left\{ \frac{ \eta B_0 C_{ij}}{1 - B_1 C_{ij}} - J_{asymp} + \sum^{nord}_{k}B_k C_{ij}^k \right\}
|
|
\]
|
|
|
|
|
|
*** Get
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_jastrow_factor_ee(qmckl_context context, double* const factor_ee);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_jastrow_factor_ee(qmckl_context context, double* const factor_ee)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_factor_ee(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
memcpy(factor_ee, ctx->jastrow.factor_ee, ctx->electron.walk_num*sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_factor_ee(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_factor_ee(qmckl_context context)
|
|
{
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
/* Check if ee rescaled distance is provided */
|
|
rc = qmckl_provide_ee_distance_rescaled(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->date > ctx->jastrow.factor_ee_date) {
|
|
|
|
/* 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);
|
|
double* factor_ee = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (factor_ee == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_provide_factor_ee",
|
|
NULL);
|
|
}
|
|
ctx->jastrow.factor_ee = factor_ee;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_factor_ee(context,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.num,
|
|
ctx->electron.up_num,
|
|
ctx->jastrow.bord_num,
|
|
ctx->jastrow.bord_vector,
|
|
ctx->electron.ee_distance_rescaled,
|
|
ctx->jastrow.asymp_jasb,
|
|
ctx->jastrow.factor_ee);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->jastrow.factor_ee_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_factor_ee
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_factor_ee_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 | up_num | in | Number of alpha electrons |
|
|
| int64_t | bord_num | in | Number of coefficients |
|
|
| double | bord_vector[bord_num + 1] | in | List of coefficients |
|
|
| double | ee_distance_rescaled[walk_num][elec_num][elec_num] | in | Electron-electron distances |
|
|
| double | asymp_jasb[2] | in | Electron-electron distances |
|
|
| double | factor_ee[walk_num] | out | Electron-electron distances |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num, bord_num, &
|
|
bord_vector, ee_distance_rescaled, asymp_jasb, factor_ee) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: walk_num, elec_num, bord_num, up_num
|
|
double precision , intent(in) :: bord_vector(bord_num + 1)
|
|
double precision , intent(in) :: ee_distance_rescaled(walk_num, elec_num, elec_num)
|
|
double precision , intent(in) :: asymp_jasb(2)
|
|
double precision , intent(out) :: factor_ee(walk_num)
|
|
|
|
integer*8 :: i, j, p, ipar, nw
|
|
double precision :: pow_ser, x, spin_fact, power_ser
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
if (bord_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_4
|
|
return
|
|
endif
|
|
|
|
factor_ee = 0.0d0
|
|
|
|
do nw =1, walk_num
|
|
do j = 1, elec_num
|
|
do i = 1, j - 1
|
|
x = ee_distance_rescaled(nw,i,j)
|
|
power_ser = 0.0d0
|
|
spin_fact = 1.0d0
|
|
ipar = 1
|
|
|
|
do p = 2, bord_num
|
|
x = x * ee_distance_rescaled(nw,i,j)
|
|
power_ser = power_ser + bord_vector(p + 1) * x
|
|
end do
|
|
|
|
if(j .LE. up_num .OR. i .GT. up_num) then
|
|
spin_fact = 0.5d0
|
|
ipar = 2
|
|
endif
|
|
|
|
factor_ee(nw) = factor_ee(nw) + spin_fact * bord_vector(1) * &
|
|
ee_distance_rescaled(nw,i,j) / &
|
|
(1.0d0 + bord_vector(2) * &
|
|
ee_distance_rescaled(nw,i,j)) &
|
|
-asymp_jasb(ipar) + power_ser
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end function qmckl_compute_factor_ee_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_header(table=qmckl_factor_ee_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_compute_factor_ee (
|
|
const qmckl_context context,
|
|
const int64_t walk_num,
|
|
const int64_t elec_num,
|
|
const int64_t up_num,
|
|
const int64_t bord_num,
|
|
const double* bord_vector,
|
|
const double* ee_distance_rescaled,
|
|
const double* asymp_jasb,
|
|
double* const factor_ee );
|
|
#+end_src
|
|
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_factor_ee_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_factor_ee &
|
|
(context, &
|
|
walk_num, &
|
|
elec_num, &
|
|
up_num, &
|
|
bord_num, &
|
|
bord_vector, &
|
|
ee_distance_rescaled, &
|
|
asymp_jasb, &
|
|
factor_ee) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
integer (c_int64_t) , intent(in) , value :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: up_num
|
|
integer (c_int64_t) , intent(in) , value :: bord_num
|
|
real (c_double ) , intent(in) :: bord_vector(bord_num + 1)
|
|
real (c_double ) , intent(in) :: ee_distance_rescaled(elec_num,elec_num,walk_num)
|
|
real (c_double ) , intent(in) :: asymp_jasb(2)
|
|
real (c_double ) , intent(out) :: factor_ee(walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_factor_ee_f
|
|
info = qmckl_compute_factor_ee_f &
|
|
(context, &
|
|
walk_num, &
|
|
elec_num, &
|
|
up_num, &
|
|
bord_num, &
|
|
bord_vector, &
|
|
ee_distance_rescaled, &
|
|
asymp_jasb, &
|
|
factor_ee)
|
|
|
|
end function qmckl_compute_factor_ee
|
|
#+end_src
|
|
|
|
*** Test
|
|
#+begin_src python :results output :exports none :noweb yes
|
|
import numpy as np
|
|
|
|
<<jastrow_data>>
|
|
|
|
<<asymp_jasb>>
|
|
|
|
factor_ee = 0.0
|
|
for i in range(0,elec_num):
|
|
for j in range(0,i):
|
|
x = ee_distance_rescaled[i][j]
|
|
pow_ser = 0.0
|
|
spin_fact = 1.0
|
|
ipar = 0
|
|
|
|
for p in range(1,bord_num):
|
|
x = x * ee_distance_rescaled[i][j]
|
|
pow_ser = pow_ser + bord_vector[p + 1] * x
|
|
|
|
if(i < up_num or j >= up_num):
|
|
spin_fact = 0.5
|
|
ipar = 1
|
|
|
|
factor_ee = factor_ee + spin_fact * bord_vector[0] * ee_distance_rescaled[i][j] \
|
|
/ (1.0 + bord_vector[1] * ee_distance_rescaled[i][j]) \
|
|
- asymp_jasb[ipar] + pow_ser
|
|
print("factor_ee :",factor_ee)
|
|
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: asym_one : 0.43340325572525706
|
|
: asymp_jasb[0] : 0.5323750557252571
|
|
: asymp_jasb[1] : 0.31567342786262853
|
|
: factor_ee : -4.282760865958113
|
|
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
/* Check if Jastrow is properly initialized */
|
|
assert(qmckl_jastrow_provided(context));
|
|
|
|
double factor_ee[walk_num];
|
|
rc = qmckl_get_jastrow_factor_ee(context, factor_ee);
|
|
|
|
// calculate factor_ee
|
|
assert(fabs(factor_ee[0]+4.282760865958113) < 1.e-12);
|
|
|
|
#+end_src
|
|
|
|
** Electron-electron component derivative \(f'_{ee}\)
|
|
|
|
Calculate the derivative of the ~factor_ee~ using the ~ee_distance_rescaled~ and
|
|
the electron-electron rescaled distances derivatives ~ee_distance_rescaled_deriv_e~.
|
|
There are four components, the gradient which has 3 components in the \(x, y, z\)
|
|
directions and the laplacian as the last component.
|
|
|
|
TODO: Add equation
|
|
|
|
|
|
*** Get
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, double* const factor_ee_deriv_e);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, double* const factor_ee_deriv_e)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_factor_ee_deriv_e(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num;
|
|
memcpy(factor_ee_deriv_e, ctx->jastrow.factor_ee_deriv_e, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context)
|
|
{
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
/* Check if ee rescaled distance is provided */
|
|
rc = qmckl_provide_ee_distance_rescaled(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Check if ee rescaled distance deriv e is provided */
|
|
rc = qmckl_provide_ee_distance_rescaled_deriv_e(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->date > ctx->jastrow.factor_ee_deriv_e_date) {
|
|
|
|
/* 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);
|
|
double* factor_ee_deriv_e = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (factor_ee_deriv_e == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_provide_factor_ee_deriv_e",
|
|
NULL);
|
|
}
|
|
ctx->jastrow.factor_ee_deriv_e = factor_ee_deriv_e;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_factor_ee_deriv_e(context,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.num,
|
|
ctx->electron.up_num,
|
|
ctx->jastrow.bord_num,
|
|
ctx->jastrow.bord_vector,
|
|
ctx->electron.ee_distance_rescaled,
|
|
ctx->electron.ee_distance_rescaled_deriv_e,
|
|
ctx->jastrow.asymp_jasb,
|
|
ctx->jastrow.factor_ee_deriv_e);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->jastrow.factor_ee_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_factor_ee_deriv_e
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_factor_ee_deriv_e_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 | up_num | in | Number of alpha electrons |
|
|
| int64_t | bord_num | in | Number of coefficients |
|
|
| double | bord_vector[bord_num + 1] | in | List of coefficients |
|
|
| double | ee_distance_rescaled[walk_num][elec_num][elec_num] | in | Electron-electron distances |
|
|
| double | ee_distance_rescaled_deriv_e[walk_num][4][elec_num][elec_num] | in | Electron-electron distances |
|
|
| double | asymp_jasb[2] | in | Electron-electron distances |
|
|
| double | factor_ee_deriv_e[walk_num][4][elec_num] | out | Electron-electron distances |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, up_num, bord_num, &
|
|
bord_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e, &
|
|
asymp_jasb, factor_ee_deriv_e) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: walk_num, elec_num, bord_num, up_num
|
|
double precision , intent(in) :: bord_vector(bord_num + 1)
|
|
double precision , intent(in) :: ee_distance_rescaled(walk_num, elec_num, elec_num)
|
|
double precision , intent(in) :: ee_distance_rescaled_deriv_e(walk_num, 4, elec_num, elec_num)
|
|
double precision , intent(in) :: asymp_jasb(2)
|
|
double precision , intent(out) :: factor_ee_deriv_e(elec_num,4,walk_num)
|
|
|
|
integer*8 :: i, j, p, ipar, nw, ii
|
|
double precision :: x, spin_fact, y
|
|
double precision :: den, invden, invden2, invden3, xinv
|
|
double precision :: lap1, lap2, lap3, third
|
|
double precision, dimension(3) :: pow_ser_g
|
|
double precision, dimension(4) :: dx
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
if (bord_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_4
|
|
return
|
|
endif
|
|
|
|
factor_ee_deriv_e = 0.0d0
|
|
third = 1.0d0 / 3.0d0
|
|
|
|
do nw =1, walk_num
|
|
do j = 1, elec_num
|
|
do i = 1, elec_num
|
|
x = ee_distance_rescaled(nw, i, j)
|
|
if(abs(x) < 1.0d-18) cycle
|
|
pow_ser_g = 0.0d0
|
|
spin_fact = 1.0d0
|
|
den = 1.0d0 + bord_vector(2) * x
|
|
invden = 1.0d0 / den
|
|
invden2 = invden * invden
|
|
invden3 = invden2 * invden
|
|
xinv = 1.0d0 / (x + 1.0d-18)
|
|
ipar = 1
|
|
|
|
do ii = 1, 4
|
|
dx(ii) = ee_distance_rescaled_deriv_e(nw, ii, i, j)
|
|
end do
|
|
|
|
if((i .LE. up_num .AND. j .LE. up_num ) .OR. &
|
|
(i .GT. up_num .AND. j .GT. up_num)) then
|
|
spin_fact = 0.5d0
|
|
endif
|
|
|
|
lap1 = 0.0d0
|
|
lap2 = 0.0d0
|
|
lap3 = 0.0d0
|
|
do ii = 1, 3
|
|
x = ee_distance_rescaled(nw, i, j)
|
|
if(abs(x) < 1.0d-18) cycle
|
|
do p = 2, bord_num
|
|
y = p * bord_vector(p + 1) * x
|
|
pow_ser_g(ii) = pow_ser_g(ii) + y * dx(ii)
|
|
lap1 = lap1 + (p - 1) * y * xinv * dx(ii) * dx(ii)
|
|
lap2 = lap2 + y
|
|
x = x * ee_distance_rescaled(nw, i, j)
|
|
end do
|
|
|
|
lap3 = lap3 - 2.0d0 * bord_vector(2) * dx(ii) * dx(ii)
|
|
|
|
factor_ee_deriv_e( j, ii, nw) = factor_ee_deriv_e( j, ii, nw) + spin_fact * bord_vector(1) * &
|
|
dx(ii) * invden2 + pow_ser_g(ii)
|
|
end do
|
|
|
|
ii = 4
|
|
lap2 = lap2 * dx(ii) * third
|
|
lap3 = lap3 + den * dx(ii)
|
|
lap3 = lap3 * (spin_fact * bord_vector(1) * invden3)
|
|
factor_ee_deriv_e( j, ii, nw) = factor_ee_deriv_e( j, ii, nw) + lap1 + lap2 + lap3
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end function qmckl_compute_factor_ee_deriv_e_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_header(table=qmckl_factor_ee_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_compute_factor_ee_deriv_e (
|
|
const qmckl_context context,
|
|
const int64_t walk_num,
|
|
const int64_t elec_num,
|
|
const int64_t up_num,
|
|
const int64_t bord_num,
|
|
const double* bord_vector,
|
|
const double* ee_distance_rescaled,
|
|
const double* ee_distance_rescaled_deriv_e,
|
|
const double* asymp_jasb,
|
|
double* const factor_ee_deriv_e );
|
|
#+end_src
|
|
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_factor_ee_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_factor_ee_deriv_e &
|
|
(context, &
|
|
walk_num, &
|
|
elec_num, &
|
|
up_num, &
|
|
bord_num, &
|
|
bord_vector, &
|
|
ee_distance_rescaled, &
|
|
ee_distance_rescaled_deriv_e, &
|
|
asymp_jasb, &
|
|
factor_ee_deriv_e) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
integer (c_int64_t) , intent(in) , value :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: up_num
|
|
integer (c_int64_t) , intent(in) , value :: bord_num
|
|
real (c_double ) , intent(in) :: bord_vector(bord_num + 1)
|
|
real (c_double ) , intent(in) :: ee_distance_rescaled(elec_num,elec_num,walk_num)
|
|
real (c_double ) , intent(in) :: ee_distance_rescaled_deriv_e(elec_num,elec_num,4,walk_num)
|
|
real (c_double ) , intent(in) :: asymp_jasb(2)
|
|
real (c_double ) , intent(out) :: factor_ee_deriv_e(elec_num,4,walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_factor_ee_deriv_e_f
|
|
info = qmckl_compute_factor_ee_deriv_e_f &
|
|
(context, &
|
|
walk_num, &
|
|
elec_num, &
|
|
up_num, &
|
|
bord_num, &
|
|
bord_vector, &
|
|
ee_distance_rescaled, &
|
|
ee_distance_rescaled_deriv_e, &
|
|
asymp_jasb, &
|
|
factor_ee_deriv_e)
|
|
|
|
end function qmckl_compute_factor_ee_deriv_e
|
|
#+end_src
|
|
|
|
*** Test
|
|
#+begin_src python :results output :exports none :noweb yes
|
|
import numpy as np
|
|
|
|
<<jastrow_data>>
|
|
|
|
<<asymp_jasb>>
|
|
|
|
kappa = 1.0
|
|
|
|
elec_coord = np.array(elec_coord)[0]
|
|
elec_dist = np.zeros(shape=(elec_num, elec_num),dtype=float)
|
|
for i in range(elec_num):
|
|
for j in range(elec_num):
|
|
elec_dist[i, j] = np.linalg.norm(elec_coord[i] - elec_coord[j])
|
|
|
|
elec_dist_deriv_e = np.zeros(shape=(4,elec_num, elec_num),dtype=float)
|
|
for j in range(elec_num):
|
|
for i in range(elec_num):
|
|
rij_inv = 1.0 / elec_dist[i, j]
|
|
for ii in range(3):
|
|
elec_dist_deriv_e[ii, i, j] = (elec_coord[j][ii] - elec_coord[i][ii]) * rij_inv
|
|
elec_dist_deriv_e[3, i, j] = 2.0 * rij_inv
|
|
elec_dist_deriv_e[:, j, j] = 0.0
|
|
|
|
ee_distance_rescaled_deriv_e = np.zeros(shape=(4,elec_num,elec_num),dtype=float)
|
|
for j in range(elec_num):
|
|
for i in range(elec_num):
|
|
f = 1.0 - kappa * ee_distance_rescaled[i][j]
|
|
for ii in range(4):
|
|
ee_distance_rescaled_deriv_e[ii][i][j] = elec_dist_deriv_e[ii][i][j]
|
|
ee_distance_rescaled_deriv_e[3][i][j] = ee_distance_rescaled_deriv_e[3][i][j] + \
|
|
(-kappa * ee_distance_rescaled_deriv_e[0][i][j] * ee_distance_rescaled_deriv_e[0][i][j]) + \
|
|
(-kappa * ee_distance_rescaled_deriv_e[1][i][j] * ee_distance_rescaled_deriv_e[1][i][j]) + \
|
|
(-kappa * ee_distance_rescaled_deriv_e[2][i][j] * ee_distance_rescaled_deriv_e[2][i][j])
|
|
for ii in range(4):
|
|
ee_distance_rescaled_deriv_e[ii][i][j] = ee_distance_rescaled_deriv_e[ii][i][j] * f
|
|
|
|
third = 1.0 / 3.0
|
|
factor_ee_deriv_e = np.zeros(shape=(4,elec_num),dtype=float)
|
|
dx = np.zeros(shape=(4),dtype=float)
|
|
pow_ser_g = np.zeros(shape=(4),dtype=float)
|
|
for j in range(elec_num):
|
|
for i in range(elec_num):
|
|
x = ee_distance_rescaled[j][i]
|
|
if abs(x) < 1e-18:
|
|
continue
|
|
pow_ser_g = np.zeros(shape=(4),dtype=float)
|
|
spin_fact = 1.0
|
|
den = 1.0 + bord_vector[1] * ee_distance_rescaled[j][i]
|
|
invden = 1.0 / den
|
|
invden2 = invden * invden
|
|
invden3 = invden2 * invden
|
|
xinv = 1.0 / (ee_distance_rescaled[j][i] + 1.0E-18)
|
|
ipar = 1
|
|
|
|
for ii in range(4):
|
|
dx[ii] = ee_distance_rescaled_deriv_e[ii][j][i]
|
|
|
|
if((i <= (up_num-1) and j <= (up_num-1) ) or \
|
|
(i > (up_num-1) and j > (up_num-1))):
|
|
spin_fact = 0.5
|
|
|
|
lap1 = 0.0
|
|
lap2 = 0.0
|
|
lap3 = 0.0
|
|
for ii in range(3):
|
|
x = ee_distance_rescaled[j][i]
|
|
if x < 1e-18:
|
|
continue
|
|
for p in range(2,bord_num+1):
|
|
y = p * bord_vector[(p-1) + 1] * x
|
|
pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii]
|
|
lap1 = lap1 + (p - 1) * y * xinv * dx[ii] * dx[ii]
|
|
lap2 = lap2 + y
|
|
x = x * ee_distance_rescaled[j][i]
|
|
|
|
lap3 = lap3 - 2.0 * bord_vector[1] * dx[ii] * dx[ii]
|
|
|
|
factor_ee_deriv_e[ii][j] = factor_ee_deriv_e[ii][j] + spin_fact * bord_vector[0] * \
|
|
dx[ii] * invden2 + pow_ser_g[ii]
|
|
|
|
ii = 3
|
|
lap2 = lap2 * dx[ii] * third
|
|
lap3 = lap3 + den * dx[ii]
|
|
lap3 = lap3 * (spin_fact * bord_vector[0] * invden3)
|
|
factor_ee_deriv_e[ii][j] = factor_ee_deriv_e[ii][j] + lap1 + lap2 + lap3
|
|
|
|
print("factor_ee_deriv_e[0][0]:",factor_ee_deriv_e[0][0])
|
|
print("factor_ee_deriv_e[1][0]:",factor_ee_deriv_e[1][0])
|
|
print("factor_ee_deriv_e[2][0]:",factor_ee_deriv_e[2][0])
|
|
print("factor_ee_deriv_e[3][0]:",factor_ee_deriv_e[3][0])
|
|
print(factor_ee_deriv_e)
|
|
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
#+begin_example
|
|
asym_one : 0.43340325572525706
|
|
asymp_jasb[0] : 0.5323750557252571
|
|
asymp_jasb[1] : 0.31567342786262853
|
|
factor_ee_deriv_e[0][0]: 0.16364894652107934
|
|
factor_ee_deriv_e[1][0]: -0.6927548119830084
|
|
factor_ee_deriv_e[2][0]: 0.073267755223968
|
|
factor_ee_deriv_e[3][0]: 1.5111672803213185
|
|
[[ 0.16364895 0.60354957 -0.19825547 0.02359797 -0.13123153 -0.18789233
|
|
0.07762515 -0.42459184 0.27920265 -0.2056531 ]
|
|
[-0.69275481 0.15690393 0.09831069 0.18490587 0.04361723 0.3250686
|
|
0.12657961 -0.01736522 -0.40149005 0.17622416]
|
|
[ 0.07326776 -0.27532276 0.22396943 0.18771633 -0.34506246 0.07298062
|
|
0.63302352 -0.00910198 -0.30238713 -0.25908332]
|
|
[ 1.51116728 1.5033247 0.00325003 2.89377255 0.1338393 2.15893795
|
|
1.74732003 0.23561147 2.67455607 0.82810434]]
|
|
#+end_example
|
|
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
/* Check if Jastrow is properly initialized */
|
|
assert(qmckl_jastrow_provided(context));
|
|
|
|
// calculate factor_ee_deriv_e
|
|
double factor_ee_deriv_e[walk_num][4][elec_num];
|
|
rc = qmckl_get_jastrow_factor_ee_deriv_e(context, &(factor_ee_deriv_e[0][0][0]));
|
|
|
|
// check factor_ee_deriv_e
|
|
assert(fabs(factor_ee_deriv_e[0][0][0]-0.16364894652107934) < 1.e-12);
|
|
assert(fabs(factor_ee_deriv_e[0][1][0]+0.6927548119830084 ) < 1.e-12);
|
|
assert(fabs(factor_ee_deriv_e[0][2][0]-0.073267755223968 ) < 1.e-12);
|
|
assert(fabs(factor_ee_deriv_e[0][3][0]-1.5111672803213185 ) < 1.e-12);
|
|
|
|
#+end_src
|
|
|
|
** Electron-nucleus component \(f_{en}\)
|
|
|
|
Calculate the electron-electron jastrow component ~factor_en~ using the ~aord_vector~
|
|
coeffecients and the electron-nucleus rescaled distances ~en_distance_rescaled~.
|
|
|
|
\[
|
|
f_{en} = \sum_{i,j<i} \left\{ \frac{ A_0 C_{ij}}{1 - A_1 C_{ij}} + \sum^{nord}_{k}A_k C_{ij}^k \right\}
|
|
\]
|
|
|
|
|
|
*** Get
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_jastrow_factor_en(qmckl_context context, double* const factor_en);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_jastrow_factor_en(qmckl_context context, double* const factor_en)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_factor_en(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
memcpy(factor_en, ctx->jastrow.factor_en, ctx->electron.walk_num*sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_factor_en(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_factor_en(qmckl_context context)
|
|
{
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
/* Check if en rescaled distance is provided */
|
|
rc = qmckl_provide_en_distance_rescaled(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->date > ctx->jastrow.factor_en_date) {
|
|
|
|
/* 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);
|
|
double* factor_en = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (factor_en == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_provide_factor_en",
|
|
NULL);
|
|
}
|
|
ctx->jastrow.factor_en = factor_en;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_factor_en(context,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.num,
|
|
ctx->nucleus.num,
|
|
ctx->jastrow.type_nucl_num,
|
|
ctx->jastrow.type_nucl_vector,
|
|
ctx->jastrow.aord_num,
|
|
ctx->jastrow.aord_vector,
|
|
ctx->electron.en_distance_rescaled,
|
|
ctx->jastrow.factor_en);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->jastrow.factor_en_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_factor_en
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_factor_en_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 nucleii |
|
|
| int64_t | type_nucl_num | in | Number of unique nuclei |
|
|
| int64_t | type_nucl_vector[nucl_num] | in | IDs of unique nucleii |
|
|
| int64_t | aord_num | in | Number of coefficients |
|
|
| double | aord_vector[aord_num + 1][type_nucl_num] | in | List of coefficients |
|
|
| double | en_distance_rescaled[walk_num][nucl_num][elec_num] | in | Electron-nucleus distances |
|
|
| double | factor_en[walk_num] | out | Electron-nucleus jastrow |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_factor_en_f(context, walk_num, elec_num, nucl_num, type_nucl_num, &
|
|
type_nucl_vector, aord_num, aord_vector, &
|
|
en_distance_rescaled, factor_en) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: walk_num, elec_num, aord_num, nucl_num, type_nucl_num
|
|
integer*8 , intent(in) :: type_nucl_vector(nucl_num)
|
|
double precision , intent(in) :: aord_vector(aord_num + 1, type_nucl_num)
|
|
double precision , intent(in) :: en_distance_rescaled(walk_num, nucl_num, elec_num)
|
|
double precision , intent(out) :: factor_en(walk_num)
|
|
|
|
integer*8 :: i, a, p, ipar, nw
|
|
double precision :: x, spin_fact, power_ser
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
if (nucl_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_4
|
|
return
|
|
endif
|
|
|
|
if (aord_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_7
|
|
return
|
|
endif
|
|
|
|
factor_en = 0.0d0
|
|
|
|
do nw =1, walk_num
|
|
do a = 1, nucl_num
|
|
do i = 1, elec_num
|
|
x = en_distance_rescaled(nw, a, i)
|
|
power_ser = 0.0d0
|
|
|
|
do p = 2, aord_num
|
|
x = x * en_distance_rescaled(nw, a, i)
|
|
power_ser = power_ser + aord_vector(p + 1, type_nucl_vector(a)) * x
|
|
end do
|
|
|
|
factor_en(nw) = factor_en(nw) + aord_vector(1, type_nucl_vector(a)) * &
|
|
en_distance_rescaled(nw, a, i) / &
|
|
(1.0d0 + aord_vector(2, type_nucl_vector(a)) * &
|
|
en_distance_rescaled(nw, a, i)) &
|
|
+ power_ser
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end function qmckl_compute_factor_en_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_header(table=qmckl_factor_en_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_compute_factor_en (
|
|
const qmckl_context context,
|
|
const int64_t walk_num,
|
|
const int64_t elec_num,
|
|
const int64_t nucl_num,
|
|
const int64_t type_nucl_num,
|
|
const int64_t* type_nucl_vector,
|
|
const int64_t aord_num,
|
|
const double* aord_vector,
|
|
const double* en_distance_rescaled,
|
|
double* const factor_en );
|
|
#+end_src
|
|
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_factor_en_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_factor_en &
|
|
(context, &
|
|
walk_num, &
|
|
elec_num, &
|
|
nucl_num, &
|
|
type_nucl_num, &
|
|
type_nucl_vector, &
|
|
aord_num, &
|
|
aord_vector, &
|
|
en_distance_rescaled, &
|
|
factor_en) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
integer (c_int64_t) , intent(in) , value :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
integer (c_int64_t) , intent(in) , value :: type_nucl_num
|
|
integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num)
|
|
integer (c_int64_t) , intent(in) , value :: aord_num
|
|
real (c_double ) , intent(in) :: aord_vector(aord_num + 1, type_nucl_num)
|
|
real (c_double ) , intent(in) :: en_distance_rescaled(walk_num, nucl_num, elec_num)
|
|
real (c_double ) , intent(out) :: factor_en(walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_factor_en_f
|
|
info = qmckl_compute_factor_en_f &
|
|
(context, &
|
|
walk_num, &
|
|
elec_num, &
|
|
nucl_num, &
|
|
type_nucl_num, &
|
|
type_nucl_vector, &
|
|
aord_num, &
|
|
aord_vector, &
|
|
en_distance_rescaled, &
|
|
factor_en)
|
|
|
|
end function qmckl_compute_factor_en
|
|
#+end_src
|
|
|
|
*** Test
|
|
#+begin_src python :results output :exports none :noweb yes
|
|
import numpy as np
|
|
|
|
<<jastrow_data>>
|
|
|
|
factor_en = 0.0
|
|
for a in range(0,nucl_num):
|
|
for i in range(0,elec_num):
|
|
x = en_distance_rescaled[i][a]
|
|
pow_ser = 0.0
|
|
|
|
for p in range(2,aord_num+1):
|
|
x = x * en_distance_rescaled[i][a]
|
|
pow_ser = pow_ser + aord_vector[(p-1) + 1][type_nucl_vector[a]-1] * x
|
|
|
|
factor_en = factor_en + aord_vector[0][type_nucl_vector[a]-1] * en_distance_rescaled[i][a] \
|
|
/ (1.0 + aord_vector[1][type_nucl_vector[a]-1] * en_distance_rescaled[i][a]) \
|
|
+ pow_ser
|
|
print("factor_en :",factor_en)
|
|
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: factor_en : -5.865822569188727
|
|
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
/* Check if Jastrow is properly initialized */
|
|
assert(qmckl_jastrow_provided(context));
|
|
|
|
double factor_en[walk_num];
|
|
rc = qmckl_get_jastrow_factor_en(context, factor_en);
|
|
|
|
// calculate factor_en
|
|
assert(fabs(factor_en[0]+5.865822569188727) < 1.e-12);
|
|
|
|
#+end_src
|
|
|
|
** Electron-nucleus component derivative \(f'_{en}\)
|
|
Calculate the electron-electron jastrow component ~factor_en_deriv_e~ derivative
|
|
with respect to the electron coordinates using the ~en_distance_rescaled~ and ~en_distance_rescaled_deriv_e~ which are already calculated previously.
|
|
|
|
TODO: write equations.
|
|
|
|
*** Get
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, double* const factor_en_deriv_e);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, double* const factor_en_deriv_e)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_factor_en_deriv_e(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num;
|
|
memcpy(factor_en_deriv_e, ctx->jastrow.factor_en_deriv_e, sze*sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context)
|
|
{
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
/* Check if en rescaled distance is provided */
|
|
rc = qmckl_provide_en_distance_rescaled(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Check if en rescaled distance derivatives is provided */
|
|
rc = qmckl_provide_en_distance_rescaled_deriv_e(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->date > ctx->jastrow.factor_en_deriv_e_date) {
|
|
|
|
/* 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);
|
|
double* factor_en_deriv_e = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (factor_en_deriv_e == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_provide_factor_en_deriv_e",
|
|
NULL);
|
|
}
|
|
ctx->jastrow.factor_en_deriv_e = factor_en_deriv_e;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_factor_en_deriv_e(context,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.num,
|
|
ctx->nucleus.num,
|
|
ctx->jastrow.type_nucl_num,
|
|
ctx->jastrow.type_nucl_vector,
|
|
ctx->jastrow.aord_num,
|
|
ctx->jastrow.aord_vector,
|
|
ctx->electron.en_distance_rescaled,
|
|
ctx->electron.en_distance_rescaled_deriv_e,
|
|
ctx->jastrow.factor_en_deriv_e);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->jastrow.factor_en_deriv_e_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_factor_en_deriv_e
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_factor_en_deriv_e_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 nucleii |
|
|
| int64_t | type_nucl_num | in | Number of unique nuclei |
|
|
| int64_t | type_nucl_vector[nucl_num] | in | IDs of unique nucleii |
|
|
| int64_t | aord_num | in | Number of coefficients |
|
|
| double | aord_vector[aord_num + 1][type_nucl_num] | in | List of coefficients |
|
|
| double | en_distance_rescaled[walk_num][nucl_num][elec_num] | in | Electron-nucleus distances |
|
|
| double | en_distance_rescaled_deriv_e[walk_num][4][nucl_num][elec_num] | in | Electron-nucleus distance derivatives |
|
|
| double | factor_en_deriv_e[walk_num][4][elec_num] | out | Electron-nucleus jastrow |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num, nucl_num, type_nucl_num, &
|
|
type_nucl_vector, aord_num, aord_vector, &
|
|
en_distance_rescaled, en_distance_rescaled_deriv_e, factor_en_deriv_e) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: walk_num, elec_num, aord_num, nucl_num, type_nucl_num
|
|
integer*8 , intent(in) :: type_nucl_vector(nucl_num)
|
|
double precision , intent(in) :: aord_vector(aord_num + 1, type_nucl_num)
|
|
double precision , intent(in) :: en_distance_rescaled(walk_num, elec_num, nucl_num)
|
|
double precision , intent(in) :: en_distance_rescaled_deriv_e(walk_num, 4, elec_num, nucl_num)
|
|
double precision , intent(out) :: factor_en_deriv_e(elec_num,4,walk_num)
|
|
|
|
integer*8 :: i, a, p, ipar, nw, ii
|
|
double precision :: x, spin_fact, den, invden, invden2, invden3, xinv
|
|
double precision :: y, lap1, lap2, lap3, third
|
|
double precision, dimension(3) :: power_ser_g
|
|
double precision, dimension(4) :: dx
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
if (nucl_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_4
|
|
return
|
|
endif
|
|
|
|
if (aord_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_7
|
|
return
|
|
endif
|
|
|
|
factor_en_deriv_e = 0.0d0
|
|
third = 1.0d0 / 3.0d0
|
|
|
|
do nw =1, walk_num
|
|
do a = 1, nucl_num
|
|
do i = 1, elec_num
|
|
x = en_distance_rescaled(nw, i, a)
|
|
if(abs(x) < 1.0d-18) continue
|
|
power_ser_g = 0.0d0
|
|
den = 1.0d0 + aord_vector(2, type_nucl_vector(a)) * x
|
|
invden = 1.0d0 / den
|
|
invden2 = invden * invden
|
|
invden3 = invden2 * invden
|
|
xinv = 1.0d0 / x
|
|
|
|
do ii = 1, 4
|
|
dx(ii) = en_distance_rescaled_deriv_e(nw, ii, i, a)
|
|
end do
|
|
|
|
lap1 = 0.0d0
|
|
lap2 = 0.0d0
|
|
lap3 = 0.0d0
|
|
do ii = 1, 3
|
|
x = en_distance_rescaled(nw, i, a)
|
|
do p = 2, aord_num
|
|
y = p * aord_vector(p + 1, type_nucl_vector(a)) * x
|
|
power_ser_g(ii) = power_ser_g(ii) + y * dx(ii)
|
|
lap1 = lap1 + (p - 1) * y * xinv * dx(ii) * dx(ii)
|
|
lap2 = lap2 + y
|
|
x = x * en_distance_rescaled(nw, i, a)
|
|
end do
|
|
|
|
lap3 = lap3 - 2.0d0 * aord_vector(2, type_nucl_vector(a)) * dx(ii) * dx(ii)
|
|
|
|
factor_en_deriv_e(i, ii, nw) = factor_en_deriv_e(i, ii, nw) + aord_vector(1, type_nucl_vector(a)) &
|
|
* dx(ii) * invden2 &
|
|
+ power_ser_g(ii)
|
|
|
|
end do
|
|
|
|
ii = 4
|
|
lap2 = lap2 * dx(ii) * third
|
|
lap3 = lap3 + den * dx(ii)
|
|
lap3 = lap3 * aord_vector(1, type_nucl_vector(a)) * invden3
|
|
factor_en_deriv_e(i, ii, nw) = factor_en_deriv_e(i, ii, nw) + lap1 + lap2 + lap3
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end function qmckl_compute_factor_en_deriv_e_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_header(table=qmckl_factor_en_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_compute_factor_en_deriv_e (
|
|
const qmckl_context context,
|
|
const int64_t walk_num,
|
|
const int64_t elec_num,
|
|
const int64_t nucl_num,
|
|
const int64_t type_nucl_num,
|
|
const int64_t* type_nucl_vector,
|
|
const int64_t aord_num,
|
|
const double* aord_vector,
|
|
const double* en_distance_rescaled,
|
|
const double* en_distance_rescaled_deriv_e,
|
|
double* const factor_en_deriv_e );
|
|
#+end_src
|
|
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_factor_en_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_factor_en_deriv_e &
|
|
(context, &
|
|
walk_num, &
|
|
elec_num, &
|
|
nucl_num, &
|
|
type_nucl_num, &
|
|
type_nucl_vector, &
|
|
aord_num, &
|
|
aord_vector, &
|
|
en_distance_rescaled, &
|
|
en_distance_rescaled_deriv_e, &
|
|
factor_en_deriv_e) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
integer (c_int64_t) , intent(in) , value :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
integer (c_int64_t) , intent(in) , value :: type_nucl_num
|
|
integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num)
|
|
integer (c_int64_t) , intent(in) , value :: aord_num
|
|
real (c_double ) , intent(in) :: aord_vector(aord_num + 1, type_nucl_num)
|
|
real (c_double ) , intent(in) :: en_distance_rescaled(elec_num,nucl_num,walk_num)
|
|
real (c_double ) , intent(in) :: en_distance_rescaled_deriv_e(elec_num,nucl_num,4,walk_num)
|
|
real (c_double ) , intent(out) :: factor_en_deriv_e(elec_num,4,walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_factor_en_deriv_e_f
|
|
info = qmckl_compute_factor_en_deriv_e_f &
|
|
(context, &
|
|
walk_num, &
|
|
elec_num, &
|
|
nucl_num, &
|
|
type_nucl_num, &
|
|
type_nucl_vector, &
|
|
aord_num, &
|
|
aord_vector, &
|
|
en_distance_rescaled, &
|
|
en_distance_rescaled_deriv_e, &
|
|
factor_en_deriv_e)
|
|
|
|
end function qmckl_compute_factor_en_deriv_e
|
|
#+end_src
|
|
|
|
*** Test
|
|
#+begin_src python :results output :exports none :noweb yes
|
|
import numpy as np
|
|
|
|
<<jastrow_data>>
|
|
|
|
kappa = 1.0
|
|
|
|
elec_coord = np.array(elec_coord)[0]
|
|
nucl_coord = np.array(nucl_coord)
|
|
elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float)
|
|
for i in range(elec_num):
|
|
for j in range(nucl_num):
|
|
elnuc_dist[i, j] = np.linalg.norm(elec_coord[i] - nucl_coord[:,j])
|
|
|
|
elnuc_dist_deriv_e = np.zeros(shape=(4, elec_num, nucl_num),dtype=float)
|
|
for a in range(nucl_num):
|
|
for i in range(elec_num):
|
|
rij_inv = 1.0 / elnuc_dist[i, a]
|
|
for ii in range(3):
|
|
elnuc_dist_deriv_e[ii, i, a] = (elec_coord[i][ii] - nucl_coord[ii][a]) * rij_inv
|
|
elnuc_dist_deriv_e[3, i, a] = 2.0 * rij_inv
|
|
|
|
en_distance_rescaled_deriv_e = np.zeros(shape=(4,elec_num,nucl_num),dtype=float)
|
|
for a in range(nucl_num):
|
|
for i in range(elec_num):
|
|
f = 1.0 - kappa * en_distance_rescaled[i][a]
|
|
for ii in range(4):
|
|
en_distance_rescaled_deriv_e[ii][i][a] = elnuc_dist_deriv_e[ii][i][a]
|
|
en_distance_rescaled_deriv_e[3][i][a] = en_distance_rescaled_deriv_e[3][i][a] + \
|
|
(-kappa * en_distance_rescaled_deriv_e[0][i][a] * en_distance_rescaled_deriv_e[0][i][a]) + \
|
|
(-kappa * en_distance_rescaled_deriv_e[1][i][a] * en_distance_rescaled_deriv_e[1][i][a]) + \
|
|
(-kappa * en_distance_rescaled_deriv_e[2][i][a] * en_distance_rescaled_deriv_e[2][i][a])
|
|
for ii in range(4):
|
|
en_distance_rescaled_deriv_e[ii][i][a] = en_distance_rescaled_deriv_e[ii][i][a] * f
|
|
|
|
third = 1.0 / 3.0
|
|
factor_en_deriv_e = np.zeros(shape=(4,elec_num),dtype=float)
|
|
dx = np.zeros(shape=(4),dtype=float)
|
|
pow_ser_g = np.zeros(shape=(3),dtype=float)
|
|
for a in range(nucl_num):
|
|
for i in range(elec_num):
|
|
x = en_distance_rescaled[i][a]
|
|
if abs(x) < 1e-18:
|
|
continue
|
|
pow_ser_g = np.zeros(shape=(3),dtype=float)
|
|
den = 1.0 + aord_vector[1][type_nucl_vector[a]-1] * x
|
|
invden = 1.0 / den
|
|
invden2 = invden * invden
|
|
invden3 = invden2 * invden
|
|
xinv = 1.0 / (x + 1.0E-18)
|
|
|
|
for ii in range(4):
|
|
dx[ii] = en_distance_rescaled_deriv_e[ii][i][a]
|
|
|
|
lap1 = 0.0
|
|
lap2 = 0.0
|
|
lap3 = 0.0
|
|
for ii in range(3):
|
|
x = en_distance_rescaled[i][a]
|
|
if x < 1e-18:
|
|
continue
|
|
for p in range(2,aord_num+1):
|
|
y = p * aord_vector[(p-1) + 1][type_nucl_vector[a]-1] * x
|
|
pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii]
|
|
lap1 = lap1 + (p - 1) * y * xinv * dx[ii] * dx[ii]
|
|
lap2 = lap2 + y
|
|
x = x * en_distance_rescaled[i][a]
|
|
|
|
lap3 = lap3 - 2.0 * aord_vector[1][type_nucl_vector[a]-1] * dx[ii] * dx[ii]
|
|
|
|
factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + aord_vector[0][type_nucl_vector[a]-1] * \
|
|
dx[ii] * invden2 + pow_ser_g[ii]
|
|
|
|
ii = 3
|
|
lap2 = lap2 * dx[ii] * third
|
|
lap3 = lap3 + den * dx[ii]
|
|
lap3 = lap3 * (aord_vector[0][type_nucl_vector[a]-1] * invden3)
|
|
factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + lap1 + lap2 + lap3
|
|
|
|
print("factor_en_deriv_e[0][0]:",factor_en_deriv_e[0][0])
|
|
print("factor_en_deriv_e[1][0]:",factor_en_deriv_e[1][0])
|
|
print("factor_en_deriv_e[2][0]:",factor_en_deriv_e[2][0])
|
|
print("factor_en_deriv_e[3][0]:",factor_en_deriv_e[3][0])
|
|
|
|
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: factor_en_deriv_e[0][0]: 0.11609919541763383
|
|
: factor_en_deriv_e[1][0]: -0.23301394780804574
|
|
: factor_en_deriv_e[2][0]: 0.17548337641865783
|
|
: factor_en_deriv_e[3][0]: -0.9667363412285741
|
|
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
/* Check if Jastrow is properly initialized */
|
|
assert(qmckl_jastrow_provided(context));
|
|
|
|
// calculate factor_en_deriv_e
|
|
double factor_en_deriv_e[walk_num][4][elec_num];
|
|
rc = qmckl_get_jastrow_factor_en_deriv_e(context, &(factor_en_deriv_e[0][0][0]));
|
|
|
|
// check factor_en_deriv_e
|
|
assert(fabs(factor_en_deriv_e[0][0][0]-0.11609919541763383) < 1.e-12);
|
|
assert(fabs(factor_en_deriv_e[0][1][0]+0.23301394780804574) < 1.e-12);
|
|
assert(fabs(factor_en_deriv_e[0][2][0]-0.17548337641865783) < 1.e-12);
|
|
assert(fabs(factor_en_deriv_e[0][3][0]+0.9667363412285741 ) < 1.e-12);
|
|
|
|
#+end_src
|
|
|
|
** Electron-electron rescaled distances for each order ~een_rescaled_e~ stores the table of the rescaled distances between all
|
|
pairs of electrons and raised to the power \(p\) defined by ~cord_num~:
|
|
|
|
\[
|
|
C_{ij,p} = \left( 1 - \exp{-\kappa C_{ij}} \right)^p
|
|
\]
|
|
|
|
where \(C_{ij}\) is the matrix of electron-electron distances.
|
|
|
|
*** Get
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double* const distance_rescaled);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double* const distance_rescaled)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_een_rescaled_e(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
|
|
memcpy(distance_rescaled, ctx->jastrow.een_rescaled_e, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
/* Check if ee distance is provided */
|
|
qmckl_exit_code rc = qmckl_provide_ee_distance(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->date > ctx->jastrow.een_rescaled_e_date) {
|
|
|
|
/* 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);
|
|
double* een_rescaled_e = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (een_rescaled_e == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_een_rescaled_e",
|
|
NULL);
|
|
}
|
|
ctx->jastrow.een_rescaled_e = een_rescaled_e;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_een_rescaled_e(context,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.num,
|
|
ctx->jastrow.cord_num,
|
|
ctx->electron.rescale_factor_kappa_ee,
|
|
ctx->electron.ee_distance,
|
|
ctx->jastrow.een_rescaled_e);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->jastrow.een_rescaled_e_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_een_rescaled_e
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_factor_een_rescaled_e_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 | cord_num | in | Order of polynomials |
|
|
| double | rescale_factor_kappa_ee | in | Factor to rescale ee distances |
|
|
| double | ee_distance[walk_num][elec_num][elec_num] | in | Electron-electron distances |
|
|
| double | een_rescaled_e[walk_num][elec_num][elec_num][0:cord_num] | out | Electron-electron rescaled distances |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, &
|
|
ee_distance, een_rescaled_e) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: walk_num
|
|
integer*8 , intent(in) :: elec_num
|
|
integer*8 , intent(in) :: cord_num
|
|
double precision , intent(in) :: rescale_factor_kappa_ee
|
|
double precision , intent(in) :: ee_distance(elec_num,elec_num,walk_num)
|
|
double precision , intent(out) :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num)
|
|
double precision,dimension(:,:),allocatable :: een_rescaled_e_ij
|
|
double precision :: x
|
|
integer*8 :: i, j, k, l, nw
|
|
|
|
allocate(een_rescaled_e_ij(elec_num * (elec_num - 1) / 2, cord_num + 1))
|
|
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
if (cord_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_4
|
|
return
|
|
endif
|
|
|
|
! Prepare table of exponentiated distances raised to appropriate power
|
|
een_rescaled_e = 0.0d0
|
|
do nw = 1, walk_num
|
|
een_rescaled_e_ij = 0.0d0
|
|
een_rescaled_e_ij(:, 1) = 1.0d0
|
|
|
|
k = 0
|
|
do j = 1, elec_num
|
|
do i = 1, j - 1
|
|
k = k + 1
|
|
een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_kappa_ee * ee_distance(i, j, nw))
|
|
end do
|
|
end do
|
|
|
|
do l = 2, cord_num
|
|
do k = 1, elec_num * (elec_num - 1)/2
|
|
een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 2)
|
|
end do
|
|
end do
|
|
|
|
! prepare the actual een table
|
|
een_rescaled_e(0, :, :, nw) = 1.0d0
|
|
do l = 1, cord_num
|
|
k = 0
|
|
do j = 1, elec_num
|
|
do i = 1, j - 1
|
|
k = k + 1
|
|
x = een_rescaled_e_ij(k, l + 1)
|
|
een_rescaled_e(l, i, j, nw) = x
|
|
een_rescaled_e(l, j, i, nw) = x
|
|
end do
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end function qmckl_compute_een_rescaled_e_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_header(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_compute_een_rescaled_e (
|
|
const qmckl_context context,
|
|
const int64_t walk_num,
|
|
const int64_t elec_num,
|
|
const int64_t cord_num,
|
|
const double rescale_factor_kappa_ee,
|
|
const double* ee_distance,
|
|
double* const een_rescaled_e );
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_een_rescaled_e &
|
|
(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
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) :: ee_distance(elec_num,elec_num,walk_num)
|
|
real (c_double ) , intent(out) :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_een_rescaled_e_f
|
|
info = qmckl_compute_een_rescaled_e_f &
|
|
(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e)
|
|
|
|
end function qmckl_compute_een_rescaled_e
|
|
#+end_src
|
|
|
|
*** Test
|
|
|
|
#+begin_src python :results output :exports none :noweb yes
|
|
import numpy as np
|
|
|
|
<<jastrow_data>>
|
|
|
|
elec_coord = np.array(elec_coord)[0]
|
|
elec_dist = np.zeros(shape=(elec_num, elec_num),dtype=float)
|
|
for i in range(elec_num):
|
|
for j in range(elec_num):
|
|
elec_dist[i, j] = np.linalg.norm(elec_coord[i] - elec_coord[j])
|
|
|
|
kappa = 1.0
|
|
|
|
een_rescaled_e_ij = np.zeros(shape=(elec_num * (elec_num - 1)//2, cord_num+1), dtype=float)
|
|
een_rescaled_e_ij[:,0] = 1.0
|
|
|
|
k = 0
|
|
for j in range(elec_num):
|
|
for i in range(j):
|
|
een_rescaled_e_ij[k, 1] = np.exp(-kappa * elec_dist[i, j])
|
|
k = k + 1
|
|
|
|
for l in range(2, cord_num + 1):
|
|
for k in range(elec_num * (elec_num - 1)//2):
|
|
een_rescaled_e_ij[k, l] = een_rescaled_e_ij[k, l - 1] * een_rescaled_e_ij[k, 1]
|
|
|
|
een_rescaled_e = np.zeros(shape=(elec_num, elec_num, cord_num + 1), dtype=float)
|
|
een_rescaled_e[:,:,0] = 1.0
|
|
|
|
for l in range(1,cord_num+1):
|
|
k = 0
|
|
for j in range(elec_num):
|
|
for i in range(j):
|
|
x = een_rescaled_e_ij[k, l]
|
|
een_rescaled_e[i, j, l] = x
|
|
een_rescaled_e[j, i, l] = x
|
|
k = k + 1
|
|
|
|
print(" een_rescaled_e[0, 2, 1] = ",een_rescaled_e[0, 2, 1])
|
|
print(" een_rescaled_e[0, 3, 1] = ",een_rescaled_e[0, 3, 1])
|
|
print(" een_rescaled_e[0, 4, 1] = ",een_rescaled_e[0, 4, 1])
|
|
print(" een_rescaled_e[1, 3, 2] = ",een_rescaled_e[1, 3, 2])
|
|
print(" een_rescaled_e[1, 4, 2] = ",een_rescaled_e[1, 4, 2])
|
|
print(" een_rescaled_e[1, 5, 2] = ",een_rescaled_e[1, 5, 2])
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: een_rescaled_e[0, 2, 1] = 0.08084493981483197
|
|
: een_rescaled_e[0, 3, 1] = 0.1066745707571846
|
|
: een_rescaled_e[0, 4, 1] = 0.01754273169464735
|
|
: een_rescaled_e[1, 3, 2] = 0.02214680362033448
|
|
: een_rescaled_e[1, 4, 2] = 0.0005700154999202759
|
|
: een_rescaled_e[1, 5, 2] = 0.3424402276009091
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
assert(qmckl_electron_provided(context));
|
|
|
|
|
|
double een_rescaled_e[walk_num][elec_num][elec_num][(cord_num + 1)];
|
|
rc = qmckl_get_jastrow_een_rescaled_e(context, &(een_rescaled_e[0][0][0][0]));
|
|
|
|
// value of (0,2,1)
|
|
assert(fabs(een_rescaled_e[0][0][2][1]-0.08084493981483197) < 1.e-12);
|
|
assert(fabs(een_rescaled_e[0][0][3][1]-0.1066745707571846) < 1.e-12);
|
|
assert(fabs(een_rescaled_e[0][0][4][1]-0.01754273169464735) < 1.e-12);
|
|
assert(fabs(een_rescaled_e[0][1][3][2]-0.02214680362033448) < 1.e-12);
|
|
assert(fabs(een_rescaled_e[0][1][4][2]-0.0005700154999202759) < 1.e-12);
|
|
assert(fabs(een_rescaled_e[0][1][5][2]-0.3424402276009091) < 1.e-12);
|
|
|
|
#+end_src
|
|
|
|
** Electron-electron rescaled distances for each order and derivatives ~een_rescaled_e~ stores the table of the rescaled distances between all
|
|
pairs of electrons and raised to the power \(p\) defined by ~cord_num~.
|
|
Here we take its derivatives required for the een jastrow.
|
|
|
|
TODO: write formulae
|
|
|
|
|
|
*** Get
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context, double* const distance_rescaled);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context, double* const distance_rescaled)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_een_rescaled_e_deriv_e(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
|
|
memcpy(distance_rescaled, ctx->jastrow.een_rescaled_e_deriv_e, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
/* Check if ee distance is provided */
|
|
qmckl_exit_code rc = qmckl_provide_een_rescaled_e(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->date > ctx->jastrow.een_rescaled_e_deriv_e_date) {
|
|
|
|
/* 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);
|
|
double* een_rescaled_e_deriv_e = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (een_rescaled_e_deriv_e == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_een_rescaled_e_deriv_e",
|
|
NULL);
|
|
}
|
|
ctx->jastrow.een_rescaled_e_deriv_e = een_rescaled_e_deriv_e;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_factor_een_rescaled_e_deriv_e(context,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.num,
|
|
ctx->jastrow.cord_num,
|
|
ctx->electron.rescale_factor_kappa_ee,
|
|
ctx->electron.coord_new,
|
|
ctx->electron.ee_distance,
|
|
ctx->jastrow.een_rescaled_e,
|
|
ctx->jastrow.een_rescaled_e_deriv_e);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->jastrow.een_rescaled_e_deriv_e_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_een_rescaled_e_deriv_e
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_factor_een_rescaled_e_deriv_e_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 | cord_num | in | Order of polynomials |
|
|
| double | rescale_factor_kappa_ee | in | Factor to rescale ee distances |
|
|
| double | coord_new[walk_num][3][elec_num] | in | Electron coordinates |
|
|
| double | ee_distance[walk_num][elec_num][elec_num] | in | Electron-electron distances |
|
|
| double | een_rescaled_e[walk_num][elec_num][elec_num][0:cord_num] | in | Electron-electron distances |
|
|
| double | een_rescaled_e_deriv_e[walk_num][elec_num][4][elec_num][0:cord_num] | out | Electron-electron rescaled distances |
|
|
|
|
#+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) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: walk_num
|
|
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) :: ee_distance(elec_num,elec_num,walk_num)
|
|
double precision , intent(in) :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num)
|
|
double precision , intent(out) :: een_rescaled_e_deriv_e(0:cord_num,elec_num,4,elec_num,walk_num)
|
|
double precision,dimension(:,:,:),allocatable :: elec_dist_deriv_e
|
|
double precision :: x, rij_inv, kappa_l
|
|
integer*8 :: i, j, k, l, nw, ii
|
|
|
|
allocate(elec_dist_deriv_e(4,elec_num,elec_num))
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
if (cord_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_4
|
|
return
|
|
endif
|
|
|
|
! Prepare table of exponentiated distances raised to appropriate power
|
|
een_rescaled_e_deriv_e = 0.0d0
|
|
do nw = 1, walk_num
|
|
do j = 1, elec_num
|
|
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
|
|
end do
|
|
elec_dist_deriv_e(4, i, j) = 2.0d0 * rij_inv
|
|
end do
|
|
elec_dist_deriv_e(:, j, j) = 0.0d0
|
|
end do
|
|
|
|
! prepare the actual een table
|
|
do l = 1, cord_num
|
|
kappa_l = - dble(l) * rescale_factor_kappa_ee
|
|
do j = 1, elec_num
|
|
do i = 1, elec_num
|
|
do ii = 1, 4
|
|
een_rescaled_e_deriv_e(l, i, ii, j, nw) = kappa_l * elec_dist_deriv_e(ii, i, j)
|
|
end do
|
|
|
|
een_rescaled_e_deriv_e(l, i, 4, j, nw) = een_rescaled_e_deriv_e(l, i, 4, j, nw) &
|
|
+ een_rescaled_e_deriv_e(l, i, 1, j, nw) * een_rescaled_e_deriv_e(l, i, 1, j, nw) &
|
|
+ een_rescaled_e_deriv_e(l, i, 2, j, nw) * een_rescaled_e_deriv_e(l, i, 2, j, nw) &
|
|
+ een_rescaled_e_deriv_e(l, i, 3, j, nw) * een_rescaled_e_deriv_e(l, i, 3, j, nw)
|
|
|
|
do ii = 1, 4
|
|
een_rescaled_e_deriv_e(l, i, ii, j, nw) = een_rescaled_e_deriv_e(l, i, ii, j, nw) * &
|
|
een_rescaled_e(l, i, j, nw)
|
|
end do
|
|
end do
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end function qmckl_compute_factor_een_rescaled_e_deriv_e_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_header(table=qmckl_factor_een_rescaled_e_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_compute_factor_een_rescaled_e_deriv_e (
|
|
const qmckl_context context,
|
|
const int64_t walk_num,
|
|
const int64_t elec_num,
|
|
const int64_t cord_num,
|
|
const double rescale_factor_kappa_ee,
|
|
const double* coord_new,
|
|
const double* ee_distance,
|
|
const double* een_rescaled_e,
|
|
double* const een_rescaled_e_deriv_e );
|
|
#+end_src
|
|
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_e_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_factor_een_rescaled_e_deriv_e &
|
|
(context, &
|
|
walk_num, &
|
|
elec_num, &
|
|
cord_num, &
|
|
rescale_factor_kappa_ee, &
|
|
coord_new, &
|
|
ee_distance, &
|
|
een_rescaled_e, &
|
|
een_rescaled_e_deriv_e) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
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) :: ee_distance(elec_num,elec_num,walk_num)
|
|
real (c_double ) , intent(in) :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num)
|
|
real (c_double ) , intent(out) :: een_rescaled_e_deriv_e(0:cord_num,elec_num,4,elec_num,walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_factor_een_rescaled_e_deriv_e_f
|
|
info = 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)
|
|
|
|
end function qmckl_compute_factor_een_rescaled_e_deriv_e
|
|
#+end_src
|
|
|
|
*** Test
|
|
|
|
#+begin_src python :results output :exports none :noweb yes
|
|
import numpy as np
|
|
|
|
<<jastrow_data>>
|
|
|
|
elec_coord = np.array(elec_coord)[0]
|
|
elec_dist = np.zeros(shape=(elec_num, elec_num),dtype=float)
|
|
for i in range(elec_num):
|
|
for j in range(elec_num):
|
|
elec_dist[i, j] = np.linalg.norm(elec_coord[i] - elec_coord[j])
|
|
|
|
kappa = 1.0
|
|
|
|
een_rescaled_e_ij = np.zeros(shape=(elec_num * (elec_num - 1)//2, cord_num+1), dtype=float)
|
|
een_rescaled_e_ij[:,0] = 1.0
|
|
|
|
k = 0
|
|
for j in range(elec_num):
|
|
for i in range(j):
|
|
een_rescaled_e_ij[k, 1] = np.exp(-kappa * elec_dist[i, j])
|
|
k = k + 1
|
|
|
|
for l in range(2, cord_num + 1):
|
|
for k in range(elec_num * (elec_num - 1)//2):
|
|
een_rescaled_e_ij[k, l] = een_rescaled_e_ij[k, l - 1] * een_rescaled_e_ij[k, 1]
|
|
|
|
een_rescaled_e = np.zeros(shape=(elec_num, elec_num, cord_num + 1), dtype=float)
|
|
een_rescaled_e[:,:,0] = 1.0
|
|
|
|
for l in range(1,cord_num+1):
|
|
k = 0
|
|
for j in range(elec_num):
|
|
for i in range(j):
|
|
x = een_rescaled_e_ij[k, l]
|
|
een_rescaled_e[i, j, l] = x
|
|
een_rescaled_e[j, i, l] = x
|
|
k = k + 1
|
|
|
|
print(" een_rescaled_e[0, 2, 1] = ",een_rescaled_e[0, 2, 1])
|
|
print(" een_rescaled_e[0, 3, 1] = ",een_rescaled_e[0, 3, 1])
|
|
print(" een_rescaled_e[0, 4, 1] = ",een_rescaled_e[0, 4, 1])
|
|
print(" een_rescaled_e[1, 3, 2] = ",een_rescaled_e[1, 3, 2])
|
|
print(" een_rescaled_e[1, 4, 2] = ",een_rescaled_e[1, 4, 2])
|
|
print(" een_rescaled_e[1, 5, 2] = ",een_rescaled_e[1, 5, 2])
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: een_rescaled_e[0, 2, 1] = 0.08084493981483197
|
|
: een_rescaled_e[0, 3, 1] = 0.1066745707571846
|
|
: een_rescaled_e[0, 4, 1] = 0.01754273169464735
|
|
: een_rescaled_e[1, 3, 2] = 0.02214680362033448
|
|
: een_rescaled_e[1, 4, 2] = 0.0005700154999202759
|
|
: een_rescaled_e[1, 5, 2] = 0.3424402276009091
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
//assert(qmckl_electron_provided(context));
|
|
|
|
#+end_src
|
|
|
|
** Electron-nucleus rescaled distances for each order ~een_rescaled_n~ stores the table of the rescaled distances between
|
|
electrons and nucleii raised to the power \(p\) defined by ~cord_num~:
|
|
|
|
\[
|
|
C_{ia,p} = \left( 1 - \exp{-\kappa C_{ia}} \right)^p
|
|
\]
|
|
|
|
where \(C_{ia}\) is the matrix of electron-nucleus distances.
|
|
|
|
*** Get
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double* const distance_rescaled);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double* const distance_rescaled)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_een_rescaled_n(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
|
|
memcpy(distance_rescaled, ctx->jastrow.een_rescaled_n, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
/* Check if ee distance is provided */
|
|
qmckl_exit_code rc = qmckl_provide_en_distance(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->date > ctx->jastrow.een_rescaled_n_date) {
|
|
|
|
/* 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);
|
|
double* een_rescaled_n = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (een_rescaled_n == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_een_rescaled_n",
|
|
NULL);
|
|
}
|
|
ctx->jastrow.een_rescaled_n = een_rescaled_n;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_een_rescaled_n(context,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.num,
|
|
ctx->nucleus.num,
|
|
ctx->jastrow.cord_num,
|
|
ctx->electron.rescale_factor_kappa_en,
|
|
ctx->electron.en_distance,
|
|
ctx->jastrow.een_rescaled_n);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->jastrow.een_rescaled_n_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_een_rescaled_n
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_factor_een_rescaled_n_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 atoms |
|
|
| int64_t | cord_num | in | Order of polynomials |
|
|
| double | rescale_factor_kappa_en | in | Factor to rescale ee distances |
|
|
| double | en_distance[walk_num][elec_num][nucl_num] | in | Electron-nucleus distances |
|
|
| double | een_rescaled_n[walk_num][elec_num][nucl_num][0:cord_num] | out | Electron-nucleus rescaled distances |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_een_rescaled_n_f(context, walk_num, elec_num, nucl_num, cord_num, rescale_factor_kappa_en, &
|
|
en_distance, een_rescaled_n) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: walk_num
|
|
integer*8 , intent(in) :: elec_num
|
|
integer*8 , intent(in) :: nucl_num
|
|
integer*8 , intent(in) :: cord_num
|
|
double precision , intent(in) :: rescale_factor_kappa_en
|
|
double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num)
|
|
double precision , intent(out) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
|
|
double precision :: x
|
|
integer*8 :: i, a, k, l, nw
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
if (nucl_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_4
|
|
return
|
|
endif
|
|
|
|
if (cord_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_5
|
|
return
|
|
endif
|
|
|
|
! Prepare table of exponentiated distances raised to appropriate power
|
|
een_rescaled_n = 0.0d0
|
|
do nw = 1, walk_num
|
|
|
|
! prepare the actual een table
|
|
een_rescaled_n(0, :, :, nw) = 1.0d0
|
|
|
|
do a = 1, nucl_num
|
|
do i = 1, elec_num
|
|
een_rescaled_n(1, a, i, nw) = dexp(-rescale_factor_kappa_en * en_distance(i, a, nw))
|
|
end do
|
|
end do
|
|
|
|
do l = 2, cord_num
|
|
do a = 1, nucl_num
|
|
do i = 1, elec_num
|
|
een_rescaled_n(l, a, i, nw) = een_rescaled_n(l - 1, a, i, nw) * een_rescaled_n(1, a, i, nw)
|
|
end do
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end function qmckl_compute_een_rescaled_n_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_header(table=qmckl_factor_een_rescaled_n_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_compute_een_rescaled_n (
|
|
const qmckl_context context,
|
|
const int64_t walk_num,
|
|
const int64_t elec_num,
|
|
const int64_t nucl_num,
|
|
const int64_t cord_num,
|
|
const double rescale_factor_kappa_en,
|
|
const double* en_distance,
|
|
double* const een_rescaled_n );
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_n_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_een_rescaled_n &
|
|
(context, &
|
|
walk_num, &
|
|
elec_num, &
|
|
nucl_num, &
|
|
cord_num, &
|
|
rescale_factor_kappa_en, &
|
|
en_distance, &
|
|
een_rescaled_n) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
integer (c_int64_t) , intent(in) , value :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
integer (c_int64_t) , intent(in) , value :: cord_num
|
|
real (c_double ) , intent(in) , value :: rescale_factor_kappa_en
|
|
real (c_double ) , intent(in) :: en_distance(nucl_num,elec_num,walk_num)
|
|
real (c_double ) , intent(out) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_een_rescaled_n_f
|
|
info = qmckl_compute_een_rescaled_n_f &
|
|
(context, &
|
|
walk_num, &
|
|
elec_num, &
|
|
nucl_num, &
|
|
cord_num, &
|
|
rescale_factor_kappa_en, &
|
|
en_distance, &
|
|
een_rescaled_n)
|
|
|
|
end function qmckl_compute_een_rescaled_n
|
|
#+end_src
|
|
|
|
*** Test
|
|
|
|
#+begin_src python :results output :exports none :noweb yes
|
|
import numpy as np
|
|
|
|
<<jastrow_data>>
|
|
|
|
elec_coord = np.array(elec_coord)[0]
|
|
nucl_coord = np.array(nucl_coord)
|
|
elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float)
|
|
for i in range(elec_num):
|
|
for a in range(nucl_num):
|
|
elnuc_dist[i, a] = np.linalg.norm(elec_coord[i] - nucl_coord[:,a])
|
|
|
|
kappa = 1.0
|
|
|
|
een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float)
|
|
een_rescaled_n[:,:,0] = 1.0
|
|
|
|
for a in range(nucl_num):
|
|
for i in range(elec_num):
|
|
een_rescaled_n[a, i, 1] = np.exp(-kappa * elnuc_dist[i, a])
|
|
|
|
for l in range(2,cord_num+1):
|
|
for a in range(nucl_num):
|
|
for i in range(elec_num):
|
|
een_rescaled_n[a, i, l] = een_rescaled_n[a, i, l - 1] * een_rescaled_n[a, i, 1]
|
|
|
|
print(" een_rescaled_n[0, 2, 1] = ",een_rescaled_n[0, 2, 1])
|
|
print(" een_rescaled_n[0, 3, 1] = ",een_rescaled_n[0, 3, 1])
|
|
print(" een_rescaled_n[0, 4, 1] = ",een_rescaled_n[0, 4, 1])
|
|
print(" een_rescaled_n[1, 3, 2] = ",een_rescaled_n[1, 3, 2])
|
|
print(" een_rescaled_n[1, 4, 2] = ",een_rescaled_n[1, 4, 2])
|
|
print(" een_rescaled_n[1, 5, 2] = ",een_rescaled_n[1, 5, 2])
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: een_rescaled_n[0, 2, 1] = 0.10612983920006765
|
|
: een_rescaled_n[0, 3, 1] = 0.135652809635553
|
|
: een_rescaled_n[0, 4, 1] = 0.023391817607642338
|
|
: een_rescaled_n[1, 3, 2] = 0.880957224822116
|
|
: een_rescaled_n[1, 4, 2] = 0.027185942659395074
|
|
: een_rescaled_n[1, 5, 2] = 0.01343938025140174
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
assert(qmckl_electron_provided(context));
|
|
|
|
double een_rescaled_n[walk_num][elec_num][nucl_num][(cord_num + 1)];
|
|
rc = qmckl_get_jastrow_een_rescaled_n(context, &(een_rescaled_n[0][0][0][0]));
|
|
|
|
// value of (0,2,1)
|
|
assert(fabs(een_rescaled_n[0][2][0][1]-0.10612983920006765) < 1.e-12);
|
|
assert(fabs(een_rescaled_n[0][3][0][1]-0.135652809635553) < 1.e-12);
|
|
assert(fabs(een_rescaled_n[0][4][0][1]-0.023391817607642338) < 1.e-12);
|
|
assert(fabs(een_rescaled_n[0][3][1][2]-0.880957224822116) < 1.e-12);
|
|
assert(fabs(een_rescaled_n[0][4][1][2]-0.027185942659395074) < 1.e-12);
|
|
assert(fabs(een_rescaled_n[0][5][1][2]-0.01343938025140174) < 1.e-12);
|
|
|
|
#+end_src
|
|
|
|
** Electron-nucleus rescaled distances for each order and derivatives ~een_rescaled_n_deriv_e~ stores the table of the rescaled distances between
|
|
electrons and nucleii raised to the power \(p\) defined by ~cord_num~:
|
|
|
|
|
|
*** Get
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context, double* const distance_rescaled);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context, double* const distance_rescaled)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_een_rescaled_n_deriv_e(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
|
|
memcpy(distance_rescaled, ctx->jastrow.een_rescaled_n_deriv_e, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
/* Check if ee distance is provided */
|
|
qmckl_exit_code rc = qmckl_provide_en_distance(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Check if ee distance is provided */
|
|
rc = qmckl_provide_een_rescaled_n(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->date > ctx->jastrow.een_rescaled_n_deriv_e_date) {
|
|
|
|
/* 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);
|
|
double* een_rescaled_n_deriv_e = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (een_rescaled_n_deriv_e == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_een_rescaled_n_deriv_e",
|
|
NULL);
|
|
}
|
|
ctx->jastrow.een_rescaled_n_deriv_e = een_rescaled_n_deriv_e;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_factor_een_rescaled_n_deriv_e(context,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.num,
|
|
ctx->nucleus.num,
|
|
ctx->jastrow.cord_num,
|
|
ctx->electron.rescale_factor_kappa_en,
|
|
ctx->electron.coord_new,
|
|
ctx->nucleus.coord,
|
|
ctx->electron.en_distance,
|
|
ctx->jastrow.een_rescaled_n,
|
|
ctx->jastrow.een_rescaled_n_deriv_e);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->jastrow.een_rescaled_n_deriv_e_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_factor_een_rescaled_n_deriv_e
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_compute_factor_een_rescaled_n_deriv_e_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 atoms |
|
|
| int64_t | cord_num | in | Order of polynomials |
|
|
| double | rescale_factor_kappa_en | in | Factor to rescale ee distances |
|
|
| double | coord_new[walk_num][3][elec_num] | in | Electron coordinates |
|
|
| double | coord[3][nucl_num] | in | Nuclear coordinates |
|
|
| double | en_distance[walk_num][elec_num][nucl_num] | in | Electron-nucleus distances |
|
|
| double | een_rescaled_n[walk_num][elec_num][nucl_num][0:cord_num] | in | Electron-nucleus distances |
|
|
| double | een_rescaled_n_deriv_e[walk_num][elec_num][4][nucl_num][0:cord_num] | out | Electron-nucleus rescaled distances |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
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) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: walk_num
|
|
integer*8 , intent(in) :: elec_num
|
|
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) :: en_distance(elec_num,nucl_num,walk_num)
|
|
double precision , intent(in) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
|
|
double precision , intent(out) :: een_rescaled_n_deriv_e(0:cord_num,nucl_num,4,elec_num,walk_num)
|
|
double precision,dimension(:,:,:),allocatable :: elnuc_dist_deriv_e
|
|
double precision :: x, ria_inv, kappa_l
|
|
integer*8 :: i, a, k, l, nw, ii
|
|
|
|
allocate(elnuc_dist_deriv_e(4, elec_num, nucl_num))
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
if (nucl_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_4
|
|
return
|
|
endif
|
|
|
|
if (cord_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_5
|
|
return
|
|
endif
|
|
|
|
! Prepare table of exponentiated distances raised to appropriate power
|
|
een_rescaled_n_deriv_e = 0.0d0
|
|
do nw = 1, walk_num
|
|
|
|
! prepare the actual een table
|
|
do a = 1, nucl_num
|
|
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
|
|
end do
|
|
elnuc_dist_deriv_e(4, i, a) = 2.0d0 * ria_inv
|
|
end do
|
|
end do
|
|
|
|
do l = 0, cord_num
|
|
kappa_l = - dble(l) * rescale_factor_kappa_en
|
|
do a = 1, nucl_num
|
|
do i = 1, elec_num
|
|
do ii = 1, 4
|
|
een_rescaled_n_deriv_e(l, a, ii, i, nw) = kappa_l * elnuc_dist_deriv_e(ii, i, a)
|
|
end do
|
|
|
|
een_rescaled_n_deriv_e(l, a, 4, i, nw) = een_rescaled_n_deriv_e(l, a, 4, i, nw) &
|
|
+ een_rescaled_n_deriv_e(l, a, 1, i, nw) * een_rescaled_n_deriv_e(l, a, 1, i, nw) &
|
|
+ een_rescaled_n_deriv_e(l, a, 2, i, nw) * een_rescaled_n_deriv_e(l, a, 2, i, nw) &
|
|
+ een_rescaled_n_deriv_e(l, a, 3, i, nw) * een_rescaled_n_deriv_e(l, a, 3, i, nw)
|
|
|
|
do ii = 1, 4
|
|
een_rescaled_n_deriv_e(l, a, ii, i, nw) = een_rescaled_n_deriv_e(l, a, ii, i, nw) * &
|
|
een_rescaled_n(l, a, i, nw)
|
|
end do
|
|
end do
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end function qmckl_compute_factor_een_rescaled_n_deriv_e_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_header(table=qmckl_compute_factor_een_rescaled_n_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_compute_factor_een_rescaled_n_deriv_e (
|
|
const qmckl_context context,
|
|
const int64_t walk_num,
|
|
const int64_t elec_num,
|
|
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* en_distance,
|
|
const double* een_rescaled_n,
|
|
double* const een_rescaled_n_deriv_e );
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_compute_factor_een_rescaled_n_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_factor_een_rescaled_n_deriv_e &
|
|
(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) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
integer (c_int64_t) , intent(in) , value :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
integer (c_int64_t) , intent(in) , value :: 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) :: en_distance(nucl_num,elec_num,walk_num)
|
|
real (c_double ) , intent(in) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
|
|
real (c_double ) , intent(out) :: een_rescaled_n_deriv_e(0:cord_num,nucl_num,4,elec_num,walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_factor_een_rescaled_n_deriv_e_f
|
|
info = 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)
|
|
|
|
end function qmckl_compute_factor_een_rescaled_n_deriv_e
|
|
#+end_src
|
|
|
|
*** Test
|
|
|
|
#+begin_src python :results output :exports none :noweb yes
|
|
import numpy as np
|
|
|
|
<<jastrow_data>>
|
|
|
|
elec_coord = np.array(elec_coord)[0]
|
|
nucl_coord = np.array(nucl_coord)
|
|
elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float)
|
|
for i in range(elec_num):
|
|
for a in range(nucl_num):
|
|
elnuc_dist[i, a] = np.linalg.norm(elec_coord[i] - nucl_coord[:,a])
|
|
|
|
kappa = 1.0
|
|
|
|
een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float)
|
|
een_rescaled_n[:,:,0] = 1.0
|
|
|
|
for a in range(nucl_num):
|
|
for i in range(elec_num):
|
|
een_rescaled_n[a, i, 1] = np.exp(-kappa * elnuc_dist[i, a])
|
|
|
|
for l in range(2,cord_num+1):
|
|
for a in range(nucl_num):
|
|
for i in range(elec_num):
|
|
een_rescaled_n[a, i, l] = een_rescaled_n[a, i, l - 1] * een_rescaled_n[a, i, 1]
|
|
|
|
print(" een_rescaled_n[0, 2, 1] = ",een_rescaled_n[0, 2, 1])
|
|
print(" een_rescaled_n[0, 3, 1] = ",een_rescaled_n[0, 3, 1])
|
|
print(" een_rescaled_n[0, 4, 1] = ",een_rescaled_n[0, 4, 1])
|
|
print(" een_rescaled_n[1, 3, 2] = ",een_rescaled_n[1, 3, 2])
|
|
print(" een_rescaled_n[1, 4, 2] = ",een_rescaled_n[1, 4, 2])
|
|
print(" een_rescaled_n[1, 5, 2] = ",een_rescaled_n[1, 5, 2])
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: een_rescaled_n[0, 2, 1] = 0.10612983920006765
|
|
: een_rescaled_n[0, 3, 1] = 0.135652809635553
|
|
: een_rescaled_n[0, 4, 1] = 0.023391817607642338
|
|
: een_rescaled_n[1, 3, 2] = 0.880957224822116
|
|
: een_rescaled_n[1, 4, 2] = 0.027185942659395074
|
|
: een_rescaled_n[1, 5, 2] = 0.01343938025140174
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
//assert(qmckl_electron_provided(context));
|
|
|
|
#+end_src
|
|
|
|
** Prepare for electron-electron-nucleus Jastrow \(f_{een}\)
|
|
|
|
Prepare ~cord_vect_full~ and ~lkpm_combined_index~ tables required for the
|
|
calculation of the three-body jastrow ~factor_een~ and its derivative ~factor_een_deriv_e~.
|
|
|
|
*** Get
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_jastrow_dim_cord_vect(qmckl_context context, int64_t* const dim_cord_vect);
|
|
qmckl_exit_code qmckl_get_jastrow_cord_vect_full(qmckl_context context, double* const cord_vect_full);
|
|
qmckl_exit_code qmckl_get_jastrow_lkpm_combined_index(qmckl_context context, int64_t* const lkpm_combined_index);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_jastrow_dim_cord_vect(qmckl_context context, int64_t* const dim_cord_vect)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_dim_cord_vect(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
*dim_cord_vect = ctx->jastrow.dim_cord_vect;
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
qmckl_exit_code qmckl_get_jastrow_cord_vect_full(qmckl_context context, double* const cord_vect_full)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_dim_cord_vect(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
rc = qmckl_provide_cord_vect_full(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = ctx->jastrow.dim_cord_vect * ctx->nucleus.num;
|
|
memcpy(cord_vect_full, ctx->jastrow.cord_vect_full, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
qmckl_exit_code qmckl_get_jastrow_lkpm_combined_index(qmckl_context context, int64_t* const lkpm_combined_index)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_dim_cord_vect(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
rc = qmckl_provide_cord_vect_full(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
size_t sze = ctx->jastrow.dim_cord_vect * 4;
|
|
memcpy(lkpm_combined_index, ctx->jastrow.lkpm_combined_index, sze * sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_dim_cord_vect(qmckl_context context);
|
|
qmckl_exit_code qmckl_provide_cord_vect_full(qmckl_context context);
|
|
qmckl_exit_code qmckl_provide_lkpm_combined_index(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_dim_cord_vect(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->date > ctx->jastrow.dim_cord_vect_date) {
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_dim_cord_vect(context,
|
|
ctx->jastrow.cord_num,
|
|
&(ctx->jastrow.dim_cord_vect));
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->jastrow.dim_cord_vect_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
qmckl_exit_code qmckl_provide_cord_vect_full(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
/* Check if dim_cord_vect is provided */
|
|
qmckl_exit_code rc = qmckl_provide_dim_cord_vect(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->date > ctx->jastrow.cord_vect_full_date) {
|
|
|
|
/* Allocate array */
|
|
if (ctx->jastrow.cord_vect_full == NULL) {
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = ctx->jastrow.dim_cord_vect * ctx->nucleus.num * sizeof(double);
|
|
double* cord_vect_full = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (cord_vect_full == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_provide_cord_vect_full",
|
|
NULL);
|
|
}
|
|
ctx->jastrow.cord_vect_full = cord_vect_full;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_cord_vect_full(context,
|
|
ctx->nucleus.num,
|
|
ctx->jastrow.cord_num,
|
|
ctx->jastrow.dim_cord_vect,
|
|
ctx->jastrow.type_nucl_num,
|
|
ctx->jastrow.type_nucl_vector,
|
|
ctx->jastrow.cord_vector,
|
|
ctx->jastrow.cord_vect_full);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->jastrow.cord_vect_full_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
|
|
qmckl_exit_code qmckl_provide_lkpm_combined_index(qmckl_context context)
|
|
{
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
/* Check if dim_cord_vect is provided */
|
|
qmckl_exit_code rc = qmckl_provide_dim_cord_vect(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->date > ctx->jastrow.lkpm_combined_index_date) {
|
|
|
|
/* Allocate array */
|
|
if (ctx->jastrow.lkpm_combined_index == NULL) {
|
|
|
|
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
|
mem_info.size = 4 * ctx->jastrow.dim_cord_vect * sizeof(int64_t);
|
|
int64_t* lkpm_combined_index = (int64_t*) qmckl_malloc(context, mem_info);
|
|
|
|
if (lkpm_combined_index == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_provide_lkpm_combined_index",
|
|
NULL);
|
|
}
|
|
ctx->jastrow.lkpm_combined_index = lkpm_combined_index;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_lkpm_combined_index(context,
|
|
ctx->jastrow.cord_num,
|
|
ctx->jastrow.dim_cord_vect,
|
|
ctx->jastrow.lkpm_combined_index);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->jastrow.lkpm_combined_index_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute dim_cord_vect
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_dim_cord_vect
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_factor_dim_cord_vect_args
|
|
| qmckl_context | context | in | Global state |
|
|
| int64_t | cord_num | in | Order of polynomials |
|
|
| int64_t | dim_cord_vect | out | dimension of cord_vect_full table |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_dim_cord_vect_f(context, cord_num, dim_cord_vect) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: cord_num
|
|
integer*8 , intent(out) :: dim_cord_vect
|
|
double precision :: x
|
|
integer*8 :: i, a, k, l, p, lmax
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (cord_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
dim_cord_vect = 0
|
|
|
|
do p = 2, cord_num
|
|
do k = p - 1, 0, -1
|
|
if (k .ne. 0) then
|
|
lmax = p - k
|
|
else
|
|
lmax = p - k - 2
|
|
endif
|
|
do l = lmax, 0, -1
|
|
if (iand(p - k - l, 1_8) == 1) cycle
|
|
dim_cord_vect = dim_cord_vect + 1
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end function qmckl_compute_dim_cord_vect_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_header(table=qmckl_factor_dim_cord_vect_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_compute_dim_cord_vect (
|
|
const qmckl_context context,
|
|
const int64_t cord_num,
|
|
int64_t* const dim_cord_vect );
|
|
#+end_src
|
|
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_factor_dim_cord_vect_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_dim_cord_vect &
|
|
(context, cord_num, dim_cord_vect) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: cord_num
|
|
integer (c_int64_t) , intent(out) :: dim_cord_vect
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_dim_cord_vect_f
|
|
info = qmckl_compute_dim_cord_vect_f &
|
|
(context, cord_num, dim_cord_vect)
|
|
|
|
end function qmckl_compute_dim_cord_vect
|
|
#+end_src
|
|
|
|
*** Compute cord_vect_full
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_cord_vect_full
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_factor_cord_vect_full_args
|
|
| qmckl_context | context | in | Global state |
|
|
| int64_t | nucl_num | in | Number of atoms |
|
|
| int64_t | cord_num | in | Order of polynomials |
|
|
| int64_t | dim_cord_vect | in | dimension of cord full table |
|
|
| int64_t | type_nucl_num | in | dimension of cord full table |
|
|
| int64_t | type_nucl_vector[nucl_num] | in | dimension of cord full table |
|
|
| double | cord_vector[cord_num][type_nucl_num] | in | dimension of cord full table |
|
|
| double | cord_vect_full[dim_cord_vect][nucl_num] | out | Full list of coefficients |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_cord_vect_full_f(context, nucl_num, cord_num, dim_cord_vect, type_nucl_num, &
|
|
type_nucl_vector, cord_vector, cord_vect_full) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: nucl_num
|
|
integer*8 , intent(in) :: cord_num
|
|
integer*8 , intent(in) :: dim_cord_vect
|
|
integer*8 , intent(in) :: type_nucl_num
|
|
integer*8 , intent(in) :: type_nucl_vector(nucl_num)
|
|
double precision , intent(in) :: cord_vector(cord_num, type_nucl_num)
|
|
double precision , intent(out) :: cord_vect_full(nucl_num,dim_cord_vect)
|
|
double precision :: x
|
|
integer*8 :: i, a, k, l, nw
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (nucl_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (cord_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
if (type_nucl_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_4
|
|
return
|
|
endif
|
|
|
|
if (dim_cord_vect <= 0) then
|
|
info = QMCKL_INVALID_ARG_5
|
|
return
|
|
endif
|
|
|
|
|
|
do a = 1, nucl_num
|
|
cord_vect_full(1:dim_cord_vect,a) = cord_vector(1:dim_cord_vect,type_nucl_vector(a))
|
|
end do
|
|
|
|
end function qmckl_compute_cord_vect_full_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_header(table=qmckl_factor_cord_vect_full_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_compute_cord_vect_full (
|
|
const qmckl_context context,
|
|
const int64_t nucl_num,
|
|
const int64_t cord_num,
|
|
const int64_t dim_cord_vect,
|
|
const int64_t type_nucl_num,
|
|
const int64_t* type_nucl_vector,
|
|
const double* cord_vector,
|
|
double* const cord_vect_full );
|
|
#+end_src
|
|
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_factor_cord_vect_full_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_cord_vect_full &
|
|
(context, &
|
|
nucl_num, &
|
|
cord_num, &
|
|
dim_cord_vect, &
|
|
type_nucl_num, &
|
|
type_nucl_vector, &
|
|
cord_vector, &
|
|
cord_vect_full) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
integer (c_int64_t) , intent(in) , value :: cord_num
|
|
integer (c_int64_t) , intent(in) , value :: dim_cord_vect
|
|
integer (c_int64_t) , intent(in) , value :: type_nucl_num
|
|
integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num)
|
|
real (c_double ) , intent(in) :: cord_vector(type_nucl_num,cord_num)
|
|
real (c_double ) , intent(out) :: cord_vect_full(nucl_num,dim_cord_vect)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_cord_vect_full_f
|
|
info = qmckl_compute_cord_vect_full_f &
|
|
(context, &
|
|
nucl_num, &
|
|
cord_num, &
|
|
dim_cord_vect, &
|
|
type_nucl_num, &
|
|
type_nucl_vector, &
|
|
cord_vector, &
|
|
cord_vect_full)
|
|
|
|
end function qmckl_compute_cord_vect_full
|
|
#+end_src
|
|
|
|
*** Compute lkpm_combined_index
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_lkpm_combined_index
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_factor_lkpm_combined_index_args
|
|
| qmckl_context | context | in | Global state |
|
|
| int64_t | cord_num | in | Order of polynomials |
|
|
| int64_t | dim_cord_vect | in | dimension of cord full table |
|
|
| int64_t | lpkm_combined_index[4][dim_cord_vect] | out | Full list of combined indices |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_lkpm_combined_index_f(context, cord_num, dim_cord_vect, &
|
|
lkpm_combined_index) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: cord_num
|
|
integer*8 , intent(in) :: dim_cord_vect
|
|
integer*8 , intent(out) :: lkpm_combined_index(dim_cord_vect, 4)
|
|
double precision :: x
|
|
integer*8 :: i, a, k, l, kk, p, lmax, m
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (cord_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (dim_cord_vect <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
|
|
kk = 0
|
|
do p = 2, cord_num
|
|
do k = p - 1, 0, -1
|
|
if (k .ne. 0) then
|
|
lmax = p - k
|
|
else
|
|
lmax = p - k - 2
|
|
end if
|
|
do l = lmax, 0, -1
|
|
if (iand(p - k - l, 1_8) .eq. 1) cycle
|
|
m = (p - k - l)/2
|
|
kk = kk + 1
|
|
lkpm_combined_index(kk, 1) = l
|
|
lkpm_combined_index(kk, 2) = k
|
|
lkpm_combined_index(kk, 3) = p
|
|
lkpm_combined_index(kk, 4) = m
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end function qmckl_compute_lkpm_combined_index_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_header(table=qmckl_factor_lkpm_combined_index_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_compute_lkpm_combined_index (
|
|
const qmckl_context context,
|
|
const int64_t cord_num,
|
|
const int64_t dim_cord_vect,
|
|
int64_t* const lpkm_combined_index );
|
|
#+end_src
|
|
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_factor_lkpm_combined_index_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_lkpm_combined_index &
|
|
(context, cord_num, dim_cord_vect, lpkm_combined_index) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: cord_num
|
|
integer (c_int64_t) , intent(in) , value :: dim_cord_vect
|
|
integer (c_int64_t) , intent(out) :: lpkm_combined_index(dim_cord_vect,4)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_lkpm_combined_index_f
|
|
info = qmckl_compute_lkpm_combined_index_f &
|
|
(context, cord_num, dim_cord_vect, lpkm_combined_index)
|
|
|
|
end function qmckl_compute_lkpm_combined_index
|
|
#+end_src
|
|
|
|
|
|
*** Test
|
|
|
|
#+begin_src python :results output :exports none :noweb yes
|
|
import numpy as np
|
|
|
|
<<jastrow_data>>
|
|
|
|
elec_coord = np.array(elec_coord)[0]
|
|
nucl_coord = np.array(nucl_coord)
|
|
elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float)
|
|
for i in range(elec_num):
|
|
for a in range(nucl_num):
|
|
elnuc_dist[i, a] = np.linalg.norm(elec_coord[i] - nucl_coord[:,a])
|
|
|
|
kappa = 1.0
|
|
|
|
een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float)
|
|
een_rescaled_n[:,:,0] = 1.0
|
|
|
|
for a in range(nucl_num):
|
|
for i in range(elec_num):
|
|
een_rescaled_n[a, i, 1] = np.exp(-kappa * elnuc_dist[i, a])
|
|
|
|
for l in range(2,cord_num+1):
|
|
for a in range(nucl_num):
|
|
for i in range(elec_num):
|
|
een_rescaled_n[a, i, l] = een_rescaled_n[a, i, l - 1] * een_rescaled_n[a, i, 1]
|
|
|
|
print(" een_rescaled_n[0, 2, 1] = ",een_rescaled_n[0, 2, 1])
|
|
print(" een_rescaled_n[0, 3, 1] = ",een_rescaled_n[0, 3, 1])
|
|
print(" een_rescaled_n[0, 4, 1] = ",een_rescaled_n[0, 4, 1])
|
|
print(" een_rescaled_n[1, 3, 2] = ",een_rescaled_n[1, 3, 2])
|
|
print(" een_rescaled_n[1, 4, 2] = ",een_rescaled_n[1, 4, 2])
|
|
print(" een_rescaled_n[1, 5, 2] = ",een_rescaled_n[1, 5, 2])
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: een_rescaled_n[0, 2, 1] = 0.10612983920006765
|
|
: een_rescaled_n[0, 3, 1] = 0.135652809635553
|
|
: een_rescaled_n[0, 4, 1] = 0.023391817607642338
|
|
: een_rescaled_n[1, 3, 2] = 0.880957224822116
|
|
: een_rescaled_n[1, 4, 2] = 0.027185942659395074
|
|
: een_rescaled_n[1, 5, 2] = 0.01343938025140174
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
//assert(qmckl_electron_provided(context));
|
|
//
|
|
|
|
#+end_src
|
|
|
|
** Electron-electron-nucleus Jastrow \(f_{een}\)
|
|
|
|
Calculate the electron-electron-nuclear three-body jastrow component ~factor_een~
|
|
using the above prepared tables.
|
|
|
|
TODO: write equations.
|
|
|
|
*** Get
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* const factor_een);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* const factor_een)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_factor_een(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int64_t sze = ctx->electron.walk_num * ctx->electron.num;
|
|
memcpy(factor_een, ctx->jastrow.factor_een, sze*sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_factor_een(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_factor_een(qmckl_context context)
|
|
{
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
/* Check if en rescaled distance is provided */
|
|
rc = qmckl_provide_een_rescaled_e(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Check if en rescaled distance derivatives is provided */
|
|
rc = qmckl_provide_een_rescaled_n(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Check if en rescaled distance derivatives is provided */
|
|
rc = qmckl_provide_cord_vect_full(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Check if en rescaled distance derivatives is provided */
|
|
rc = qmckl_provide_lkpm_combined_index(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->date > ctx->jastrow.factor_een_date) {
|
|
|
|
/* 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);
|
|
double* factor_een = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (factor_een == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_provide_factor_een",
|
|
NULL);
|
|
}
|
|
ctx->jastrow.factor_een = factor_een;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_factor_een(context,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.num,
|
|
ctx->nucleus.num,
|
|
ctx->jastrow.cord_num,
|
|
ctx->jastrow.dim_cord_vect,
|
|
ctx->jastrow.cord_vect_full,
|
|
ctx->jastrow.lkpm_combined_index,
|
|
ctx->jastrow.een_rescaled_e,
|
|
ctx->jastrow.een_rescaled_n,
|
|
ctx->jastrow.factor_een);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->jastrow.factor_een_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_factor_een
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_factor_een_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 nucleii |
|
|
| int64_t | cord_num | in | order of polynomials |
|
|
| int64_t | dim_cord_vect | in | dimension of full coefficient vector |
|
|
| double | cord_vect_full[dim_cord_vect][nucl_num] | in | full coefficient vector |
|
|
| int64_t | lkpm_combined_index[4][dim_cord_vect] | in | combined indices |
|
|
| double | een_rescaled_e[walk_num][elec_num][elec_num][0:cord_num] | in | Electron-nucleus rescaled |
|
|
| double | een_rescaled_n[walk_num][elec_num][nucl_num][0:cord_num] | in | Electron-nucleus rescaled factor |
|
|
| double | factor_een[walk_num] | out | Electron-nucleus jastrow |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_factor_een_f(context, walk_num, elec_num, nucl_num, cord_num, dim_cord_vect, &
|
|
cord_vect_full, lkpm_combined_index, &
|
|
een_rescaled_e, een_rescaled_n, factor_een) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_cord_vect
|
|
integer*8 , intent(in) :: lkpm_combined_index(4,dim_cord_vect)
|
|
double precision , intent(in) :: cord_vect_full(dim_cord_vect, nucl_num)
|
|
double precision , intent(in) :: een_rescaled_e(walk_num, elec_num, elec_num, 0:cord_num)
|
|
double precision , intent(in) :: een_rescaled_n(walk_num, elec_num, nucl_num, 0:cord_num)
|
|
double precision , intent(out) :: factor_een(walk_num)
|
|
|
|
integer*8 :: i, a, j, l, k, p, m, n, nw
|
|
double precision :: accu, accu2, cn
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
if (nucl_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_4
|
|
return
|
|
endif
|
|
|
|
if (cord_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_5
|
|
return
|
|
endif
|
|
|
|
factor_een = 0.0d0
|
|
|
|
do nw =1, walk_num
|
|
do n = 1, dim_cord_vect
|
|
l = lkpm_combined_index(1, n)
|
|
k = lkpm_combined_index(2, n)
|
|
p = lkpm_combined_index(3, n)
|
|
m = lkpm_combined_index(4, n)
|
|
|
|
do a = 1, nucl_num
|
|
accu2 = 0.0d0
|
|
cn = cord_vect_full(n, a)
|
|
do j = 1, elec_num
|
|
accu = 0.0d0
|
|
do i = 1, elec_num
|
|
accu = accu + een_rescaled_e(nw, i, j, k) * &
|
|
een_rescaled_n(nw, i, a, m)
|
|
end do
|
|
accu2 = accu2 + accu * een_rescaled_n(nw, j, a, m + l)
|
|
end do
|
|
factor_een(nw) = factor_een(nw) + accu2 + cn
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end function qmckl_compute_factor_een_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_header(table=qmckl_factor_een_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_compute_factor_een (
|
|
const qmckl_context context,
|
|
const int64_t walk_num,
|
|
const int64_t elec_num,
|
|
const int64_t nucl_num,
|
|
const int64_t cord_num,
|
|
const int64_t dim_cord_vect,
|
|
const double* cord_vect_full,
|
|
const int64_t* lkpm_combined_index,
|
|
const double* een_rescaled_e,
|
|
const double* een_rescaled_n,
|
|
double* const factor_een );
|
|
#+end_src
|
|
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_factor_een_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_factor_een &
|
|
(context, &
|
|
walk_num, &
|
|
elec_num, &
|
|
nucl_num, &
|
|
cord_num, &
|
|
dim_cord_vect, &
|
|
cord_vect_full, &
|
|
lkpm_combined_index, &
|
|
een_rescaled_e, &
|
|
een_rescaled_n, &
|
|
factor_een) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
integer (c_int64_t) , intent(in) , value :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
integer (c_int64_t) , intent(in) , value :: cord_num
|
|
integer (c_int64_t) , intent(in) , value :: dim_cord_vect
|
|
real (c_double ) , intent(in) :: cord_vect_full(nucl_num,dim_cord_vect)
|
|
integer (c_int64_t) , intent(in) :: lkpm_combined_index(dim_cord_vect,4)
|
|
real (c_double ) , intent(in) :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num)
|
|
real (c_double ) , intent(in) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
|
|
real (c_double ) , intent(out) :: factor_een(walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_factor_een_f
|
|
info = qmckl_compute_factor_een_f &
|
|
(context, &
|
|
walk_num, &
|
|
elec_num, &
|
|
nucl_num, &
|
|
cord_num, &
|
|
dim_cord_vect, &
|
|
cord_vect_full, &
|
|
lkpm_combined_index, &
|
|
een_rescaled_e, &
|
|
een_rescaled_n, &
|
|
factor_een)
|
|
|
|
end function qmckl_compute_factor_een
|
|
#+end_src
|
|
|
|
*** Test
|
|
#+begin_src python :results output :exports none :noweb yes
|
|
import numpy as np
|
|
|
|
<<jastrow_data>>
|
|
|
|
kappa = 1.0
|
|
|
|
elec_coord = np.array(elec_coord)[0]
|
|
nucl_coord = np.array(nucl_coord)
|
|
elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float)
|
|
for i in range(elec_num):
|
|
for j in range(nucl_num):
|
|
elnuc_dist[i, j] = np.linalg.norm(elec_coord[i] - nucl_coord[:,j])
|
|
|
|
elnuc_dist_deriv_e = np.zeros(shape=(4, elec_num, nucl_num),dtype=float)
|
|
for a in range(nucl_num):
|
|
for i in range(elec_num):
|
|
rij_inv = 1.0 / elnuc_dist[i, a]
|
|
for ii in range(3):
|
|
elnuc_dist_deriv_e[ii, i, a] = (elec_coord[i][ii] - nucl_coord[ii][a]) * rij_inv
|
|
elnuc_dist_deriv_e[3, i, a] = 2.0 * rij_inv
|
|
|
|
en_distance_rescaled_deriv_e = np.zeros(shape=(4,elec_num,nucl_num),dtype=float)
|
|
for a in range(nucl_num):
|
|
for i in range(elec_num):
|
|
f = 1.0 - kappa * en_distance_rescaled[i][a]
|
|
for ii in range(4):
|
|
en_distance_rescaled_deriv_e[ii][i][a] = elnuc_dist_deriv_e[ii][i][a]
|
|
en_distance_rescaled_deriv_e[3][i][a] = en_distance_rescaled_deriv_e[3][i][a] + \
|
|
(-kappa * en_distance_rescaled_deriv_e[0][i][a] * en_distance_rescaled_deriv_e[0][i][a]) + \
|
|
(-kappa * en_distance_rescaled_deriv_e[1][i][a] * en_distance_rescaled_deriv_e[1][i][a]) + \
|
|
(-kappa * en_distance_rescaled_deriv_e[2][i][a] * en_distance_rescaled_deriv_e[2][i][a])
|
|
for ii in range(4):
|
|
en_distance_rescaled_deriv_e[ii][i][a] = en_distance_rescaled_deriv_e[ii][i][a] * f
|
|
|
|
third = 1.0 / 3.0
|
|
factor_en_deriv_e = np.zeros(shape=(4,elec_num),dtype=float)
|
|
dx = np.zeros(shape=(4),dtype=float)
|
|
pow_ser_g = np.zeros(shape=(3),dtype=float)
|
|
for a in range(nucl_num):
|
|
for i in range(elec_num):
|
|
x = en_distance_rescaled[i][a]
|
|
if abs(x) < 1e-18:
|
|
continue
|
|
pow_ser_g = np.zeros(shape=(3),dtype=float)
|
|
den = 1.0 + aord_vector[1][type_nucl_vector[a]-1] * x
|
|
invden = 1.0 / den
|
|
invden2 = invden * invden
|
|
invden3 = invden2 * invden
|
|
xinv = 1.0 / (x + 1.0E-18)
|
|
|
|
for ii in range(4):
|
|
dx[ii] = en_distance_rescaled_deriv_e[ii][i][a]
|
|
|
|
lap1 = 0.0
|
|
lap2 = 0.0
|
|
lap3 = 0.0
|
|
for ii in range(3):
|
|
x = en_distance_rescaled[i][a]
|
|
if x < 1e-18:
|
|
continue
|
|
for p in range(2,aord_num+1):
|
|
y = p * aord_vector[(p-1) + 1][type_nucl_vector[a]-1] * x
|
|
pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii]
|
|
lap1 = lap1 + (p - 1) * y * xinv * dx[ii] * dx[ii]
|
|
lap2 = lap2 + y
|
|
x = x * en_distance_rescaled[i][a]
|
|
|
|
lap3 = lap3 - 2.0 * aord_vector[1][type_nucl_vector[a]-1] * dx[ii] * dx[ii]
|
|
|
|
factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + aord_vector[0][type_nucl_vector[a]-1] * \
|
|
dx[ii] * invden2 + pow_ser_g[ii]
|
|
|
|
ii = 3
|
|
lap2 = lap2 * dx[ii] * third
|
|
lap3 = lap3 + den * dx[ii]
|
|
lap3 = lap3 * (aord_vector[0][type_nucl_vector[a]-1] * invden3)
|
|
factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + lap1 + lap2 + lap3
|
|
|
|
print("factor_en_deriv_e[0][0]:",factor_en_deriv_e[0][0])
|
|
print("factor_en_deriv_e[1][0]:",factor_en_deriv_e[1][0])
|
|
print("factor_en_deriv_e[2][0]:",factor_en_deriv_e[2][0])
|
|
print("factor_en_deriv_e[3][0]:",factor_en_deriv_e[3][0])
|
|
|
|
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: factor_en_deriv_e[0][0]: 0.11609919541763383
|
|
: factor_en_deriv_e[1][0]: -0.23301394780804574
|
|
: factor_en_deriv_e[2][0]: 0.17548337641865783
|
|
: factor_en_deriv_e[3][0]: -0.9667363412285741
|
|
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
/* Check if Jastrow is properly initialized */
|
|
//assert(qmckl_jastrow_provided(context));
|
|
//
|
|
//// calculate factor_en_deriv_e
|
|
//double factor_een[walk_num];
|
|
//rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0]));
|
|
|
|
//// check factor_en_deriv_e
|
|
//assert(fabs(factor_en_deriv_e[0][0][0]-0.11609919541763383) < 1.e-12);
|
|
//assert(fabs(factor_en_deriv_e[0][1][0]+0.23301394780804574) < 1.e-12);
|
|
//assert(fabs(factor_en_deriv_e[0][2][0]-0.17548337641865783) < 1.e-12);
|
|
//assert(fabs(factor_en_deriv_e[0][3][0]+0.9667363412285741 ) < 1.e-12);
|
|
|
|
#+end_src
|
|
|
|
** Electron-electron-nucleus Jastrow \(f_{een}\) derivative
|
|
|
|
Calculate the electron-electron-nuclear three-body jastrow component ~factor_een_deriv_e~
|
|
using the above prepared tables.
|
|
|
|
TODO: write equations.
|
|
|
|
*** Get
|
|
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
|
qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, double* const factor_een_deriv_e);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, double* const factor_een_deriv_e)
|
|
{
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
rc = qmckl_provide_factor_een_deriv_e(context);
|
|
if (rc != QMCKL_SUCCESS) return rc;
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
int64_t sze = ctx->electron.walk_num * ctx->electron.num;
|
|
memcpy(factor_een_deriv_e, ctx->jastrow.factor_een_deriv_e, sze*sizeof(double));
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Provide :noexport:
|
|
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_factor_een_deriv_e(qmckl_context context);
|
|
#+end_src
|
|
|
|
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
|
qmckl_exit_code qmckl_provide_factor_een_deriv_e(qmckl_context context)
|
|
{
|
|
|
|
qmckl_exit_code rc;
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return QMCKL_NULL_CONTEXT;
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
|
assert (ctx != NULL);
|
|
|
|
/* Check if en rescaled distance is provided */
|
|
rc = qmckl_provide_een_rescaled_e(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Check if en rescaled distance derivatives is provided */
|
|
rc = qmckl_provide_een_rescaled_n(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Check if en rescaled distance is provided */
|
|
rc = qmckl_provide_een_rescaled_e_deriv_e(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Check if en rescaled distance derivatives is provided */
|
|
rc = qmckl_provide_een_rescaled_n_deriv_e(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Check if en rescaled distance derivatives is provided */
|
|
rc = qmckl_provide_cord_vect_full(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Check if en rescaled distance derivatives is provided */
|
|
rc = qmckl_provide_lkpm_combined_index(context);
|
|
if(rc != QMCKL_SUCCESS) return rc;
|
|
|
|
/* Compute if necessary */
|
|
if (ctx->date > ctx->jastrow.factor_een_deriv_e_date) {
|
|
|
|
/* 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);
|
|
double* factor_een_deriv_e = (double*) qmckl_malloc(context, mem_info);
|
|
|
|
if (factor_een_deriv_e == NULL) {
|
|
return qmckl_failwith( context,
|
|
QMCKL_ALLOCATION_FAILED,
|
|
"qmckl_provide_factor_een_deriv_e",
|
|
NULL);
|
|
}
|
|
ctx->jastrow.factor_een_deriv_e = factor_een_deriv_e;
|
|
}
|
|
|
|
qmckl_exit_code rc =
|
|
qmckl_compute_factor_een_deriv_e(context,
|
|
ctx->electron.walk_num,
|
|
ctx->electron.num,
|
|
ctx->nucleus.num,
|
|
ctx->jastrow.cord_num,
|
|
ctx->jastrow.dim_cord_vect,
|
|
ctx->jastrow.cord_vect_full,
|
|
ctx->jastrow.lkpm_combined_index,
|
|
ctx->jastrow.een_rescaled_e,
|
|
ctx->jastrow.een_rescaled_n,
|
|
ctx->jastrow.een_rescaled_e_deriv_e,
|
|
ctx->jastrow.een_rescaled_n_deriv_e,
|
|
ctx->jastrow.factor_een_deriv_e);
|
|
if (rc != QMCKL_SUCCESS) {
|
|
return rc;
|
|
}
|
|
|
|
ctx->jastrow.factor_een_deriv_e_date = ctx->date;
|
|
}
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
*** Compute
|
|
:PROPERTIES:
|
|
:Name: qmckl_compute_factor_een_deriv_e
|
|
:CRetType: qmckl_exit_code
|
|
:FRetType: qmckl_exit_code
|
|
:END:
|
|
|
|
#+NAME: qmckl_factor_een_deriv_e_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 nucleii |
|
|
| int64_t | cord_num | in | order of polynomials |
|
|
| int64_t | dim_cord_vect | in | dimension of full coefficient vector |
|
|
| double | cord_vect_full[dim_cord_vect][nucl_num] | in | full coefficient vector |
|
|
| int64_t | lkpm_combined_index[4][dim_cord_vect] | in | combined indices |
|
|
| double | een_rescaled_e[walk_num][elec_num][elec_num][0:cord_num] | in | Electron-nucleus rescaled |
|
|
| double | een_rescaled_n[walk_num][elec_num][nucl_num][0:cord_num] | in | Electron-nucleus rescaled factor |
|
|
| double | een_rescaled_e_deriv_e[walk_num][elec_num][4][elec_num][0:cord_num] | in | Electron-nucleus rescaled |
|
|
| double | een_rescaled_n_deriv_e[walk_num][elec_num][4][nucl_num][0:cord_num] | in | Electron-nucleus rescaled factor |
|
|
| double | factor_een_deriv_e[walk_num][4][elec_num] | out | Electron-nucleus jastrow |
|
|
|
|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
|
integer function qmckl_compute_factor_een_deriv_e_f(context, walk_num, elec_num, nucl_num, cord_num, dim_cord_vect, &
|
|
cord_vect_full, lkpm_combined_index, &
|
|
een_rescaled_e, een_rescaled_n, &
|
|
een_rescaled_e_deriv_e, een_rescaled_n_deriv_e, factor_een_deriv_e) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in) :: context
|
|
integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_cord_vect
|
|
integer*8 , intent(in) :: lkpm_combined_index(4,dim_cord_vect)
|
|
double precision , intent(in) :: cord_vect_full(dim_cord_vect, nucl_num)
|
|
double precision , intent(in) :: een_rescaled_e(walk_num, elec_num, elec_num, 0:cord_num)
|
|
double precision , intent(in) :: een_rescaled_n(walk_num, elec_num, nucl_num, 0:cord_num)
|
|
double precision , intent(in) :: een_rescaled_e_deriv_e(walk_num, elec_num, 4, elec_num, 0:cord_num)
|
|
double precision , intent(in) :: een_rescaled_n_deriv_e(walk_num, elec_num, 4, nucl_num, 0:cord_num)
|
|
double precision , intent(out) :: factor_een_deriv_e(elec_num, 4, walk_num)
|
|
|
|
integer*8 :: i, a, j, l, k, p, m, n, nw
|
|
double precision :: accu, accu2, cn
|
|
double precision :: daccu(1:4), daccu2(1:4)
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (walk_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (elec_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_3
|
|
return
|
|
endif
|
|
|
|
if (nucl_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_4
|
|
return
|
|
endif
|
|
|
|
if (cord_num <= 0) then
|
|
info = QMCKL_INVALID_ARG_5
|
|
return
|
|
endif
|
|
|
|
factor_een_deriv_e = 0.0d0
|
|
|
|
do nw =1, walk_num
|
|
do n = 1, dim_cord_vect
|
|
l = lkpm_combined_index(1, n)
|
|
k = lkpm_combined_index(2, n)
|
|
p = lkpm_combined_index(3, n)
|
|
m = lkpm_combined_index(4, n)
|
|
|
|
do a = 1, nucl_num
|
|
cn = cord_vect_full(n, a)
|
|
do j = 1, elec_num
|
|
accu = 0.0d0
|
|
accu2 = 0.0d0
|
|
daccu = 0.0d0
|
|
daccu2 = 0.0d0
|
|
do i = 1, elec_num
|
|
accu = accu + een_rescaled_e(nw, i, j, k) * &
|
|
een_rescaled_n(nw, i, a, m)
|
|
accu2 = accu2 + een_rescaled_e(nw, i, j, k) * &
|
|
een_rescaled_n(nw, i, a, m + l)
|
|
daccu(1:4) = daccu(1:4) + een_rescaled_e_deriv_e(nw, j, 1:4, i, k) * &
|
|
een_rescaled_n(nw, i, a, m)
|
|
daccu2(1:4) = daccu2(1:4) + een_rescaled_e_deriv_e(nw, j, 1:4, i, k) * &
|
|
een_rescaled_n(nw, i, a, m + l)
|
|
end do
|
|
factor_een_deriv_e(j, 1:4, nw) = factor_een_deriv_e(j, 1:4, nw) + &
|
|
(accu * een_rescaled_n_deriv_e(nw, j, 1:4, a, m + l) &
|
|
+ daccu(1:4) * een_rescaled_n(nw, j, a, m + l) &
|
|
+ daccu2(1:4) * een_rescaled_n(nw, j, a, m) &
|
|
+ accu2 * een_rescaled_n_deriv_e(nw, j, 1:4, a, m)) * cn
|
|
|
|
factor_een_deriv_e(j, 4, nw) = factor_een_deriv_e(j, 4, nw) + 2.0d0 * ( &
|
|
daccu (1) * een_rescaled_n_deriv_e(nw, j, 1, a, m + l) + &
|
|
daccu (2) * een_rescaled_n_deriv_e(nw, j, 2, a, m + l) + &
|
|
daccu (3) * een_rescaled_n_deriv_e(nw, j, 3, a, m + l) + &
|
|
daccu2(1) * een_rescaled_n_deriv_e(nw, j, 1, a, m ) + &
|
|
daccu2(2) * een_rescaled_n_deriv_e(nw, j, 2, a, m ) + &
|
|
daccu2(3) * een_rescaled_n_deriv_e(nw, j, 3, a, m ) ) * cn
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end function qmckl_compute_factor_een_deriv_e_f
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_header(table=qmckl_factor_een_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_compute_factor_een_deriv_e (
|
|
const qmckl_context context,
|
|
const int64_t walk_num,
|
|
const int64_t elec_num,
|
|
const int64_t nucl_num,
|
|
const int64_t cord_num,
|
|
const int64_t dim_cord_vect,
|
|
const double* cord_vect_full,
|
|
const int64_t* lkpm_combined_index,
|
|
const double* een_rescaled_e,
|
|
const double* een_rescaled_n,
|
|
const double* een_rescaled_e_deriv_e,
|
|
const double* een_rescaled_n_deriv_e,
|
|
double* const factor_een_deriv_e );
|
|
#+end_src
|
|
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_factor_een_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_compute_factor_een_deriv_e &
|
|
(context, &
|
|
walk_num, &
|
|
elec_num, &
|
|
nucl_num, &
|
|
cord_num, &
|
|
dim_cord_vect, &
|
|
cord_vect_full, &
|
|
lkpm_combined_index, &
|
|
een_rescaled_e, &
|
|
een_rescaled_n, &
|
|
een_rescaled_e_deriv_e, &
|
|
een_rescaled_n_deriv_e, &
|
|
factor_een_deriv_e) &
|
|
bind(C) result(info)
|
|
|
|
use, intrinsic :: iso_c_binding
|
|
implicit none
|
|
|
|
integer (c_int64_t) , intent(in) , value :: context
|
|
integer (c_int64_t) , intent(in) , value :: walk_num
|
|
integer (c_int64_t) , intent(in) , value :: elec_num
|
|
integer (c_int64_t) , intent(in) , value :: nucl_num
|
|
integer (c_int64_t) , intent(in) , value :: cord_num
|
|
integer (c_int64_t) , intent(in) , value :: dim_cord_vect
|
|
real (c_double ) , intent(in) :: cord_vect_full(nucl_num,dim_cord_vect)
|
|
integer (c_int64_t) , intent(in) :: lkpm_combined_index(dim_cord_vect,4)
|
|
real (c_double ) , intent(in) :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num)
|
|
real (c_double ) , intent(in) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
|
|
real (c_double ) , intent(in) :: een_rescaled_e_deriv_e(0:cord_num,elec_num,4,elec_num,walk_num)
|
|
real (c_double ) , intent(in) :: een_rescaled_n_deriv_e(0:cord_num,nucl_num,4,elec_num,walk_num)
|
|
real (c_double ) , intent(out) :: factor_een_deriv_e(elec_num,4,walk_num)
|
|
|
|
integer(c_int32_t), external :: qmckl_compute_factor_een_deriv_e_f
|
|
info = qmckl_compute_factor_een_deriv_e_f &
|
|
(context, &
|
|
walk_num, &
|
|
elec_num, &
|
|
nucl_num, &
|
|
cord_num, &
|
|
dim_cord_vect, &
|
|
cord_vect_full, &
|
|
lkpm_combined_index, &
|
|
een_rescaled_e, &
|
|
een_rescaled_n, &
|
|
een_rescaled_e_deriv_e, &
|
|
een_rescaled_n_deriv_e, &
|
|
factor_een_deriv_e)
|
|
|
|
end function qmckl_compute_factor_een_deriv_e
|
|
#+end_src
|
|
|
|
*** Test
|
|
#+begin_src python :results output :exports none :noweb yes
|
|
import numpy as np
|
|
|
|
<<jastrow_data>>
|
|
|
|
kappa = 1.0
|
|
|
|
elec_coord = np.array(elec_coord)[0]
|
|
nucl_coord = np.array(nucl_coord)
|
|
elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float)
|
|
for i in range(elec_num):
|
|
for j in range(nucl_num):
|
|
elnuc_dist[i, j] = np.linalg.norm(elec_coord[i] - nucl_coord[:,j])
|
|
|
|
elnuc_dist_deriv_e = np.zeros(shape=(4, elec_num, nucl_num),dtype=float)
|
|
for a in range(nucl_num):
|
|
for i in range(elec_num):
|
|
rij_inv = 1.0 / elnuc_dist[i, a]
|
|
for ii in range(3):
|
|
elnuc_dist_deriv_e[ii, i, a] = (elec_coord[i][ii] - nucl_coord[ii][a]) * rij_inv
|
|
elnuc_dist_deriv_e[3, i, a] = 2.0 * rij_inv
|
|
|
|
en_distance_rescaled_deriv_e = np.zeros(shape=(4,elec_num,nucl_num),dtype=float)
|
|
for a in range(nucl_num):
|
|
for i in range(elec_num):
|
|
f = 1.0 - kappa * en_distance_rescaled[i][a]
|
|
for ii in range(4):
|
|
en_distance_rescaled_deriv_e[ii][i][a] = elnuc_dist_deriv_e[ii][i][a]
|
|
en_distance_rescaled_deriv_e[3][i][a] = en_distance_rescaled_deriv_e[3][i][a] + \
|
|
(-kappa * en_distance_rescaled_deriv_e[0][i][a] * en_distance_rescaled_deriv_e[0][i][a]) + \
|
|
(-kappa * en_distance_rescaled_deriv_e[1][i][a] * en_distance_rescaled_deriv_e[1][i][a]) + \
|
|
(-kappa * en_distance_rescaled_deriv_e[2][i][a] * en_distance_rescaled_deriv_e[2][i][a])
|
|
for ii in range(4):
|
|
en_distance_rescaled_deriv_e[ii][i][a] = en_distance_rescaled_deriv_e[ii][i][a] * f
|
|
|
|
third = 1.0 / 3.0
|
|
factor_en_deriv_e = np.zeros(shape=(4,elec_num),dtype=float)
|
|
dx = np.zeros(shape=(4),dtype=float)
|
|
pow_ser_g = np.zeros(shape=(3),dtype=float)
|
|
for a in range(nucl_num):
|
|
for i in range(elec_num):
|
|
x = en_distance_rescaled[i][a]
|
|
if abs(x) < 1e-18:
|
|
continue
|
|
pow_ser_g = np.zeros(shape=(3),dtype=float)
|
|
den = 1.0 + aord_vector[1][type_nucl_vector[a]-1] * x
|
|
invden = 1.0 / den
|
|
invden2 = invden * invden
|
|
invden3 = invden2 * invden
|
|
xinv = 1.0 / (x + 1.0E-18)
|
|
|
|
for ii in range(4):
|
|
dx[ii] = en_distance_rescaled_deriv_e[ii][i][a]
|
|
|
|
lap1 = 0.0
|
|
lap2 = 0.0
|
|
lap3 = 0.0
|
|
for ii in range(3):
|
|
x = en_distance_rescaled[i][a]
|
|
if x < 1e-18:
|
|
continue
|
|
for p in range(2,aord_num+1):
|
|
y = p * aord_vector[(p-1) + 1][type_nucl_vector[a]-1] * x
|
|
pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii]
|
|
lap1 = lap1 + (p - 1) * y * xinv * dx[ii] * dx[ii]
|
|
lap2 = lap2 + y
|
|
x = x * en_distance_rescaled[i][a]
|
|
|
|
lap3 = lap3 - 2.0 * aord_vector[1][type_nucl_vector[a]-1] * dx[ii] * dx[ii]
|
|
|
|
factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + aord_vector[0][type_nucl_vector[a]-1] * \
|
|
dx[ii] * invden2 + pow_ser_g[ii]
|
|
|
|
ii = 3
|
|
lap2 = lap2 * dx[ii] * third
|
|
lap3 = lap3 + den * dx[ii]
|
|
lap3 = lap3 * (aord_vector[0][type_nucl_vector[a]-1] * invden3)
|
|
factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + lap1 + lap2 + lap3
|
|
|
|
print("factor_en_deriv_e[0][0]:",factor_en_deriv_e[0][0])
|
|
print("factor_en_deriv_e[1][0]:",factor_en_deriv_e[1][0])
|
|
print("factor_en_deriv_e[2][0]:",factor_en_deriv_e[2][0])
|
|
print("factor_en_deriv_e[3][0]:",factor_en_deriv_e[3][0])
|
|
|
|
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: factor_en_deriv_e[0][0]: 0.11609919541763383
|
|
: factor_en_deriv_e[1][0]: -0.23301394780804574
|
|
: factor_en_deriv_e[2][0]: 0.17548337641865783
|
|
: factor_en_deriv_e[3][0]: -0.9667363412285741
|
|
|
|
|
|
#+begin_src c :tangle (eval c_test)
|
|
///* Check if Jastrow is properly initialized */
|
|
//assert(qmckl_jastrow_provided(context));
|
|
//
|
|
//// calculate factor_en_deriv_e
|
|
//double factor_een[walk_num];
|
|
//rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0]));
|
|
//
|
|
//// check factor_en_deriv_e
|
|
//assert(fabs(factor_en_deriv_e[0][0][0]-0.11609919541763383) < 1.e-12);
|
|
//assert(fabs(factor_en_deriv_e[0][1][0]+0.23301394780804574) < 1.e-12);
|
|
//assert(fabs(factor_en_deriv_e[0][2][0]-0.17548337641865783) < 1.e-12);
|
|
//assert(fabs(factor_en_deriv_e[0][3][0]+0.9667363412285741 ) < 1.e-12);
|
|
|
|
#+end_src
|
|
|
|
* End of files :noexport:
|
|
|
|
#+begin_src c :tangle (eval h_private_type)
|
|
#endif
|
|
#+end_src
|
|
|
|
*** Test
|
|
#+begin_src c :tangle (eval c_test)
|
|
rc = qmckl_context_destroy(context);
|
|
assert (rc == QMCKL_SUCCESS);
|
|
|
|
return 0;
|
|
}
|
|
#+end_src
|
|
|
|
|
|
# -*- mode: org -*-
|
|
# vim: syntax=c
|
|
|
|
|
|
|