TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2025-01-08 20:33:40 +01:00
Code Issues Projects Releases Wiki Activity

Finished factor_en. #22

Browse Source
This commit is contained in:
vijay gopal chilkuri 2021-07-06 15:51:51 +05:30
parent 0df816c0ba
commit 7c226d0a99
3 changed files with 444 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
org/qmckl_distance.org
View File

@ -923,7 +923,7 @@ integer function qmckl_distance_rescaled_f(context, transa, transb, m, n, &
if (transb == 'N' .or. transb == 'n') then
continue
else if (transa == 'T' .or. transa == 't') then
else if (transb == 'T' .or. transb == 't') then
transab = transab + 2
else
transab = -100

482
org/qmckl_jastrow.org
View File

@ -72,27 +72,28 @@ int main() {
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 |
| ~uint64_t~ | ~type_nucl_num~ | in | Number of Nucleii types |
| ~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 |
| ~double~ | ~factor_ee_date~ | out | Jastrow factor: electron-electron part |
| ~double~ | ~factor_en[walk_num]~ | out | Jastrow factor: electron-nucleus part |
| ~double~ | ~factor_en_date~ | out | Jastrow factor: electron-nucleus part |
| ~double~ | ~factor_een[walk_num]~ | out | Jastrow factor: electron-electron-nucleus part |
| ~double~ | ~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 |
| ~double~ | ~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 |
| ~double~ | ~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 |
| ~double~ | ~factor_een_deriv_e_date~ | out | Keep track of the date for the een derivative |
|-----------+--------------------------------------------+-----+-------------------------------------------------------------------|
| ~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 |
| ~double~ | ~factor_ee_date~ | out | Jastrow factor: electron-electron part |
| ~double~ | ~factor_en[walk_num]~ | out | Jastrow factor: electron-nucleus part |
| ~double~ | ~factor_en_date~ | out | Jastrow factor: electron-nucleus part |
| ~double~ | ~factor_een[walk_num]~ | out | Jastrow factor: electron-electron-nucleus part |
| ~double~ | ~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 |
| ~double~ | ~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 |
| ~double~ | ~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 |
| ~double~ | ~factor_een_deriv_e_date~ | out | Keep track of the date for the een derivative |
computed data:
@ -172,6 +173,16 @@ ee_distance_rescaled = [
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):
@ -182,8 +193,15 @@ aord_num = 5
bord_num = 5
cord_num = 23
dim_cord_vec = 23
aord_vector = [ 0.000000000000000E+000, 0.000000000000000E+000, -0.380512000000000E+000,
-0.157996000000000E+000, -3.155800000000000E-002, 2.151200000000000E-002]
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,
@ -260,6 +278,7 @@ typedef struct qmckl_jastrow_struct{
int64_t factor_ee_deriv_e_date;
int64_t factor_en_deriv_e_date;
int64_t factor_een_deriv_e_date;
int64_t* type_nucl_vector;
double * aord_vector;
double * bord_vector;
double * cord_vector;
@ -301,7 +320,7 @@ qmckl_exit_code qmckl_init_jastrow(qmckl_context context) {
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
ctx->jastrow.uninitialized = (1 << 5) - 1;
ctx->jastrow.uninitialized = (1 << 6) - 1;
/* Default values */
@ -312,13 +331,14 @@ qmckl_exit_code qmckl_init_jastrow(qmckl_context context) {
** 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_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);
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
@ -462,6 +482,33 @@ qmckl_exit_code qmckl_get_jastrow_type_nucl_num (const qmckl_context context, in
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) {
@ -478,7 +525,7 @@ qmckl_exit_code qmckl_get_jastrow_aord_vector (const qmckl_context context, doub
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
int32_t mask = 1 << 2;
int32_t mask = 1 << 3;
if ( (ctx->jastrow.uninitialized & mask) != 0) {
return QMCKL_NOT_PROVIDED;
@ -505,7 +552,7 @@ qmckl_exit_code qmckl_get_jastrow_bord_vector (const qmckl_context context, doub
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
int32_t mask = 1 << 3;
int32_t mask = 1 << 4;
if ( (ctx->jastrow.uninitialized & mask) != 0) {
return QMCKL_NOT_PROVIDED;
@ -532,7 +579,7 @@ qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, doub
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
int32_t mask = 1 << 4;
int32_t mask = 1 << 5;
if ( (ctx->jastrow.uninitialized & mask) != 0) {
return QMCKL_NOT_PROVIDED;
@ -551,12 +598,13 @@ qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, doub
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_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);
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
@ -630,11 +678,61 @@ qmckl_exit_code qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int
<<post2>>
}
qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double const * aord_vector) {
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;
@ -687,7 +785,7 @@ qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double cons
qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double const * bord_vector) {
<<pre2>>
int32_t mask = 1 << 3;
int32_t mask = 1 << 4;
int64_t bord_num;
qmckl_exit_code rc = qmckl_get_jastrow_bord_num(context, &bord_num);
@ -737,7 +835,7 @@ qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double cons
qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double const * cord_vector) {
<<pre2>>
int32_t mask = 1 << 4;
int32_t mask = 1 << 5;
int64_t cord_num;
qmckl_exit_code rc = qmckl_get_jastrow_cord_num(context, &cord_num);
@ -807,7 +905,7 @@ qmckl_exit_code qmckl_set_jastrow_dependencies(qmckl_context context) {
NULL);
}
int32_t mask = 1 << 5;
int32_t mask = 1 << 6;
<<post2>>
}
@ -1272,6 +1370,7 @@ print("asymp_jasb[1] : ", asymp_jasb[1])
#+begin_src c :tangle (eval c_test)
assert(qmckl_electron_provided(context));
int64_t* type_nucl_vector = &(n2_type_nucl_vector[0]);
double* aord_vector = &(n2_aord_vector[0][0]);
double* bord_vector = &(n2_bord_vector[0]);
double* cord_vector = &(n2_cord_vector[0][0]);
@ -1285,6 +1384,8 @@ rc = qmckl_set_jastrow_ord_num(context, n2_aord_num, n2_bord_num, n2_cord_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_jastrow_type_nucl_num(context, n2_type_nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_jastrow_type_nucl_vector(context, n2_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);
@ -2029,6 +2130,299 @@ 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[type_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(type_nucl_num)
double precision , intent(in) :: aord_vector(aord_num, nucl_num)
double precision , intent(in) :: en_distance_rescaled(walk_num, elec_num, nucl_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, i, a)
power_ser = 0.0d0
do p = 2, aord_num
x = x * en_distance_rescaled(nw, i, a)
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, i, a) / &
(1.0d0 + aord_vector(2, type_nucl_vector(a)) * &
en_distance_rescaled(nw, i, a)) &
+ 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(type_nucl_num)
integer (c_int64_t) , intent(in) , value :: aord_num
real (c_double ) , intent(in) :: aord_vector(type_nucl_num,aord_num + 1)
real (c_double ) , intent(in) :: en_distance_rescaled(elec_num,nucl_num,walk_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
* End of files :noexport:
#+begin_src c :tangle (eval h_private_type)

6
tests/n2.h
View File

@ -5,7 +5,7 @@ double n2_charge[n2_nucl_num] = { 5., 5.};
double n2_nucl_coord[3][n2_nucl_num] =
{ {0.000000, 0.000000 },
{0.000000, 0.000000 },
{2.059801, 2.059801 } };
{0.000000, 2.059801 } };
#define n2_elec_up_num ((int64_t) 5)
#define n2_elec_dn_num ((int64_t) 5)
@ -32,6 +32,10 @@ double n2_elec_coord[n2_walk_num][n2_elec_num][3] = { {
#define n2_cord_num ((int64_t) 23)
#define n2_dim_cord_vec ((int64_t) 23)
int64_t n2_type_nucl_vector[n2_nucl_num] = {
1,
1};
double n2_aord_vector[n2_aord_num + 1][n2_type_nucl_num] = {
{ 0. },
{ 0. },