1
0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-12-22 12:23:56 +01:00

Added spin-independent Jastrow

This commit is contained in:
Anthony Scemama 2023-09-26 17:32:51 +02:00
parent ad378103a5
commit 89a4a57c32

View File

@ -86,6 +86,7 @@
#include "config.h"
#endif
#include <stdbool.h>
#include <stdio.h>
#include "n2.h"
#include "qmckl_jastrow_champ_private_func.h"
@ -135,7 +136,7 @@ int main() {
#+NAME: qmckl_jastrow_args
| Variable | Type | Description |
|---------------------------+---------------------------------------+-------------------------------------------------------------------|
|---------------------+---------------------------------------+-------------------------------------------------------------------|
| ~uninitialized~ | ~int32_t~ | Keeps bits set for uninitialized data |
| ~rescale_factor_ee~ | ~double~ | The distance scaling factor |
| ~rescale_factor_en~ | ~double[type_nucl_num]~ | The distance scaling factor |
@ -147,6 +148,8 @@ int main() {
| ~a_vector~ | ~double[aord_num + 1][type_nucl_num]~ | a polynomial coefficients |
| ~b_vector~ | ~double[bord_num + 1]~ | b polynomial coefficients |
| ~c_vector~ | ~double[dim_c_vector][type_nucl_num]~ | c polynomial coefficients |
| ~c_vector~ | ~double[dim_c_vector][type_nucl_num]~ | c polynomial coefficients |
| ~spin_independent~ | ~bool~ | If true, use same parameters for parallel and anti-parallel pairs |
Computed data:
@ -368,6 +371,7 @@ typedef struct qmckl_jastrow_champ_struct{
uint64_t gl_date;
double rescale_factor_ee;
int32_t uninitialized;
bool spin_independent;
bool provided;
} qmckl_jastrow_champ_struct;
@ -394,17 +398,15 @@ qmckl_exit_code qmckl_init_jastrow_champ(qmckl_context context) {
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
ctx->jastrow_champ.uninitialized = (1 << 10) - 1;
ctx->jastrow_champ.uninitialized = (1 << 11) - 1;
/* Default values */
ctx->jastrow_champ.aord_num = -1;
ctx->jastrow_champ.bord_num = -1;
ctx->jastrow_champ.dim_c_vector = -1;
ctx->jastrow_champ.cord_num = -1;
ctx->jastrow_champ.dim_c_vector = -1;
ctx->jastrow_champ.type_nucl_num = -1;
ctx->jastrow_champ.spin_independent = false;
return QMCKL_SUCCESS;
}
@ -426,12 +428,13 @@ qmckl_exit_code qmckl_set_jastrow_champ_type_nucl_vector (qmckl_context contex
qmckl_exit_code qmckl_set_jastrow_champ_a_vector (qmckl_context context, const double * a_vector, const int64_t size_max);
qmckl_exit_code qmckl_set_jastrow_champ_b_vector (qmckl_context context, const double * b_vector, const int64_t size_max);
qmckl_exit_code qmckl_set_jastrow_champ_c_vector (qmckl_context context, const double * c_vector, const int64_t size_max);
qmckl_exit_code qmckl_set_jastrow_champ_spin_independent (qmckl_context context, const bool spin_independent);
#+end_src
#+NAME:pre2
#+begin_src c :exports none
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
return QMCKL_INVALID_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
@ -463,7 +466,7 @@ qmckl_set_jastrow_champ_aord_num(qmckl_context context, const int64_t aord_num)
int32_t mask = 1 << 0;
<<pre2>>
<<pre2>>
if (aord_num < 0) {
return qmckl_failwith( context,
@ -476,7 +479,7 @@ qmckl_set_jastrow_champ_aord_num(qmckl_context context, const int64_t aord_num)
ctx->jastrow_champ.uninitialized |= (1 << 5);
<<post2>>
}
}
qmckl_exit_code
qmckl_set_jastrow_champ_bord_num(qmckl_context context, const int64_t bord_num)
@ -484,7 +487,7 @@ qmckl_set_jastrow_champ_bord_num(qmckl_context context, const int64_t bord_num)
int32_t mask = 1 << 1;
<<pre2>>
<<pre2>>
if (bord_num < 0) {
return qmckl_failwith( context,
@ -497,7 +500,7 @@ qmckl_set_jastrow_champ_bord_num(qmckl_context context, const int64_t bord_num)
ctx->jastrow_champ.uninitialized |= (1 << 6);
<<post2>>
}
}
qmckl_exit_code
qmckl_set_jastrow_champ_cord_num(qmckl_context context, const int64_t cord_num)
@ -505,7 +508,7 @@ qmckl_set_jastrow_champ_cord_num(qmckl_context context, const int64_t cord_num)
int32_t mask = 1 << 2;
<<pre2>>
<<pre2>>
if (cord_num < 0) {
return qmckl_failwith( context,
@ -529,7 +532,7 @@ qmckl_set_jastrow_champ_cord_num(qmckl_context context, const int64_t cord_num)
}
<<post2>>
}
}
qmckl_exit_code
@ -537,7 +540,7 @@ qmckl_set_jastrow_champ_type_nucl_num(qmckl_context context, const int64_t type_
{
int32_t mask = 1 << 3;
<<pre2>>
<<pre2>>
if (type_nucl_num <= 0) {
return qmckl_failwith( context,
@ -549,7 +552,7 @@ qmckl_set_jastrow_champ_type_nucl_num(qmckl_context context, const int64_t type_
ctx->jastrow_champ.type_nucl_num = type_nucl_num;
<<post2>>
}
}
qmckl_exit_code
@ -560,7 +563,7 @@ qmckl_set_jastrow_champ_type_nucl_vector(qmckl_context context,
int32_t mask = 1 << 4;
<<pre2>>
<<pre2>>
int64_t type_nucl_num = ctx->jastrow_champ.type_nucl_num;
@ -616,7 +619,7 @@ qmckl_set_jastrow_champ_type_nucl_vector(qmckl_context context,
ctx->jastrow_champ.type_nucl_vector = new_array;
<<post2>>
}
}
qmckl_exit_code
@ -626,7 +629,7 @@ qmckl_set_jastrow_champ_a_vector(qmckl_context context,
{
int32_t mask = 1 << 5;
<<pre2>>
<<pre2>>
int64_t aord_num = ctx->jastrow_champ.aord_num;
if (aord_num < 0) {
@ -685,7 +688,7 @@ qmckl_set_jastrow_champ_a_vector(qmckl_context context,
ctx->jastrow_champ.a_vector = new_array;
<<post2>>
}
}
qmckl_exit_code
@ -695,7 +698,7 @@ qmckl_set_jastrow_champ_b_vector(qmckl_context context,
{
int32_t mask = 1 << 6;
<<pre2>>
<<pre2>>
int64_t bord_num = ctx->jastrow_champ.bord_num;
if (bord_num < 0) {
@ -745,7 +748,7 @@ qmckl_set_jastrow_champ_b_vector(qmckl_context context,
ctx->jastrow_champ.b_vector = new_array;
<<post2>>
}
}
qmckl_exit_code
@ -755,7 +758,7 @@ qmckl_set_jastrow_champ_c_vector(qmckl_context context,
{
int32_t mask = 1 << 7;
<<pre2>>
<<pre2>>
int64_t type_nucl_num = ctx->jastrow_champ.type_nucl_num;
if (type_nucl_num <= 0) {
@ -816,7 +819,7 @@ qmckl_set_jastrow_champ_c_vector(qmckl_context context,
ctx->jastrow_champ.c_vector = new_array;
<<post2>>
}
}
qmckl_exit_code
qmckl_set_jastrow_champ_rescale_factor_ee(qmckl_context context,
@ -836,7 +839,7 @@ qmckl_set_jastrow_champ_rescale_factor_ee(qmckl_context context,
ctx->jastrow_champ.rescale_factor_ee = rescale_factor_ee;
<<post2>>
}
}
qmckl_exit_code
@ -894,6 +897,20 @@ qmckl_set_jastrow_champ_rescale_factor_en(qmckl_context context,
<<post2>>
}
qmckl_exit_code
qmckl_set_jastrow_champ_spin_independent(qmckl_context context, const bool spin_independent)
{
int32_t mask = 1 << 10;
<<pre2>>
ctx->jastrow_champ.spin_independent = spin_independent;
<<post2>>
}
#+end_src
When the required information is completely entered, other data structures are
@ -908,7 +925,10 @@ qmckl_exit_code qmckl_finalize_jastrow_champ(qmckl_context context);
qmckl_exit_code qmckl_finalize_jastrow_champ(qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_finalize_jastrow_champ",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
@ -950,7 +970,7 @@ qmckl_exit_code qmckl_finalize_jastrow_champ(qmckl_context context) {
**** Fortran interface
#+begin_src f90 :tangle (eval fh_func) :comments org
interface
interface
integer(qmckl_exit_code) function qmckl_set_jastrow_champ_rescale_factor_ee (context, &
kappa_ee) bind(C)
use, intrinsic :: iso_c_binding
@ -1046,7 +1066,16 @@ interface
real(c_double), intent(in) :: c_vector(size_max)
end function qmckl_set_jastrow_champ_c_vector
end interface
integer(qmckl_exit_code) function qmckl_set_jastrow_champ_spin_independent(context, &
spin_independent) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (qmckl_context) , intent(in) , value :: context
logical(c_bool), intent(in), value :: spin_independent
end function qmckl_set_jastrow_champ_spin_independent
end interface
#+end_src
** Access functions
@ -1063,6 +1092,7 @@ qmckl_exit_code qmckl_get_jastrow_champ_c_vector (qmckl_context contex
qmckl_exit_code qmckl_get_jastrow_champ_rescale_factor_ee (const qmckl_context context, double* const rescale_factor_ee);
qmckl_exit_code qmckl_get_jastrow_champ_rescale_factor_en (const qmckl_context context, double* const rescale_factor_en, const int64_t size_max);
qmckl_exit_code qmckl_get_jastrow_champ_dim_c_vector (qmckl_context context, int64_t* const dim_c_vector);
qmckl_exit_code qmckl_get_jastrow_champ_spin_independent (qmckl_context context, bool* const spin_independent);
#+end_src
@ -1095,7 +1125,10 @@ bool qmckl_jastrow_champ_provided(const qmckl_context context) {
qmckl_exit_code qmckl_get_jastrow_champ_aord_num (const qmckl_context context, int64_t* const aord_num) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (char) 0;
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_jastrow_champ_aord_num",
NULL);
}
if (aord_num == NULL) {
@ -1122,7 +1155,10 @@ qmckl_exit_code qmckl_get_jastrow_champ_aord_num (const qmckl_context context, i
qmckl_exit_code qmckl_get_jastrow_champ_bord_num (const qmckl_context context, int64_t* const bord_num) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (char) 0;
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_jastrow_champ_bord_num",
NULL);
}
if (bord_num == NULL) {
@ -1149,7 +1185,10 @@ qmckl_exit_code qmckl_get_jastrow_champ_bord_num (const qmckl_context context, i
qmckl_exit_code qmckl_get_jastrow_champ_cord_num (const qmckl_context context, int64_t* const cord_num) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (char) 0;
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_jastrow_champ_cord_num",
NULL);
}
if (cord_num == NULL) {
@ -1176,7 +1215,10 @@ qmckl_exit_code qmckl_get_jastrow_champ_cord_num (const qmckl_context context, i
qmckl_exit_code qmckl_get_jastrow_champ_type_nucl_num (const qmckl_context context, int64_t* const type_nucl_num) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (char) 0;
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_jastrow_champ_type_nucl_num",
NULL);
}
if (type_nucl_num == NULL) {
@ -1207,7 +1249,10 @@ qmckl_get_jastrow_champ_type_nucl_vector (const qmckl_context context,
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (char) 0;
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_jastrow_champ_type_nucl_vector",
NULL);
}
if (type_nucl_vector == NULL) {
@ -1244,7 +1289,10 @@ qmckl_get_jastrow_champ_a_vector (const qmckl_context context,
const int64_t size_max) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (char) 0;
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_jastrow_champ_a_vector",
NULL);
}
if (a_vector == NULL) {
@ -1281,7 +1329,10 @@ qmckl_get_jastrow_champ_b_vector (const qmckl_context context,
const int64_t size_max) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (char) 0;
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_jastrow_champ_b_vector",
NULL);
}
if (b_vector == NULL) {
@ -1318,7 +1369,10 @@ qmckl_get_jastrow_champ_c_vector (const qmckl_context context,
const int64_t size_max) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (char) 0;
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_jastrow_champ_c_vector",
NULL);
}
if (c_vector == NULL) {
@ -1359,7 +1413,10 @@ qmckl_get_jastrow_champ_rescale_factor_ee (const qmckl_context context,
double* const rescale_factor_ee) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_jastrow_champ_rescale_factor_ee",
NULL);
}
if (rescale_factor_ee == NULL) {
@ -1389,7 +1446,10 @@ qmckl_get_jastrow_champ_rescale_factor_en (const qmckl_context context,
double* const rescale_factor_en,
const int64_t size_max) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_jastrow_champ_rescale_factor_en",
NULL);
}
if (rescale_factor_en == NULL) {
@ -1426,7 +1486,10 @@ qmckl_get_jastrow_champ_rescale_factor_en (const qmckl_context context,
qmckl_exit_code qmckl_get_jastrow_champ_dim_c_vector(qmckl_context context, int64_t* const dim_c_vector)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_jastrow_champ_dim_c_vector",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
@ -1434,6 +1497,35 @@ qmckl_exit_code qmckl_get_jastrow_champ_dim_c_vector(qmckl_context context, int6
,*dim_c_vector = ctx->jastrow_champ.dim_c_vector;
return QMCKL_SUCCESS;
}
qmckl_exit_code qmckl_get_jastrow_champ_spin_independent(const qmckl_context context, bool* const spin_independent) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_get_jastrow_champ_spin_independent",
NULL);
}
if (spin_independent == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_get_jastrow_champ_spin_independent",
"spin_independent is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 10;
if ( (ctx->jastrow_champ.uninitialized & mask) != 0) {
return QMCKL_NOT_PROVIDED;
}
,*spin_independent = ctx->jastrow_champ.spin_independent ;
return QMCKL_SUCCESS;
}
@ -1442,7 +1534,7 @@ qmckl_exit_code qmckl_get_jastrow_champ_dim_c_vector(qmckl_context context, int6
**** Fortran interface
#+begin_src f90 :tangle (eval fh_func) :comments org
interface
interface
integer(qmckl_exit_code) function qmckl_get_jastrow_champ_rescale_factor_ee (context, &
kappa_ee) bind(C)
use, intrinsic :: iso_c_binding
@ -1538,7 +1630,16 @@ interface
real(c_double), intent(out) :: c_vector(size_max)
end function qmckl_get_jastrow_champ_c_vector
end interface
integer(qmckl_exit_code) function qmckl_get_jastrow_champ_spin_independent(context, &
spin_independent) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (qmckl_context) , intent(in) , value :: context
logical(c_bool), intent(out) :: spin_independent
end function qmckl_get_jastrow_champ_spin_independent
end interface
#+end_src
** Test
@ -1660,10 +1761,12 @@ assert(qmckl_nucleus_provided(context));
** Electron-electron component
*** Asymptotic component
Calculate the asymptotic component ~asymp_jasb~ to be substracted from the
Calculate the asymptotic component ~asymp_jasb~ to be subtracted from the
electron-electron jastrow factor \(J_{\text{ee}}\). Two values are
computed. The first one is for parallel spin pairs, and the
second one for antiparallel spin pairs.
If the ~spin_independent~ variable is set to ~true~, then
$\delta^{\uparrow \downarrow}$ is always equal to one.
\[
J_{\text{ee}}^{\infty} = \frac{\frac{1}{2}(1+\delta^{\uparrow \downarrow})\,b_1 \kappa_\text{ee}^{-1}}{1 + b_2\,
@ -1766,6 +1869,7 @@ qmckl_exit_code qmckl_provide_jastrow_champ_asymp_jasb(qmckl_context context)
ctx->jastrow_champ.bord_num,
ctx->jastrow_champ.b_vector,
ctx->jastrow_champ.rescale_factor_ee,
ctx->jastrow_champ.spin_independent,
ctx->jastrow_champ.asymp_jasb);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -1787,16 +1891,17 @@ qmckl_exit_code qmckl_provide_jastrow_champ_asymp_jasb(qmckl_context context)
#+NAME: qmckl_asymp_jasb_args
| Variable | Type | In/Out | Description |
|---------------------+----------------------+--------+-------------------------|
|---------------------+----------------------+--------+---------------------------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~bord_num~ | ~int64_t~ | in | Order of the polynomial |
| ~b_vector~ | ~double[bord_num+1]~ | in | Values of b |
| ~rescale_factor_ee~ | ~double~ | in | Electron coordinates |
| ~spin_independent~ | ~bool~ | in | If true, same parameters for parallel and anti-parallel pairs |
| ~asymp_jasb~ | ~double[2]~ | out | Asymptotic value |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
function qmckl_compute_jastrow_champ_asymp_jasb_doc(context, &
bord_num, b_vector, rescale_factor_ee, asymp_jasb) &
bord_num, b_vector, rescale_factor_ee, spin_independent, asymp_jasb) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
@ -1807,6 +1912,7 @@ function qmckl_compute_jastrow_champ_asymp_jasb_doc(context, &
integer (c_int64_t) , intent(in) , value :: bord_num
real (c_double ) , intent(in) :: b_vector(bord_num+1)
real (c_double ) , intent(in) , value :: rescale_factor_ee
logical (c_bool ) , intent(in) , value :: spin_independent
real (c_double ) , intent(out) :: asymp_jasb(2)
integer(qmckl_exit_code) :: info
@ -1827,7 +1933,11 @@ function qmckl_compute_jastrow_champ_asymp_jasb_doc(context, &
endif
asym_one = b_vector(1) * kappa_inv / (1.0d0 + b_vector(2) * kappa_inv)
if (spin_independent) then
asymp_jasb(:) = (/asym_one, asym_one/)
else
asymp_jasb(:) = (/0.5d0*asym_one, asym_one/)
end if
x = kappa_inv
do p = 2, bord_num
@ -1845,7 +1955,8 @@ qmckl_exit_code qmckl_compute_jastrow_champ_asymp_jasb_doc (const qmckl_context
const int64_t bord_num,
const double* b_vector,
const double rescale_factor_ee,
double* const asymp_jasb );
const bool spin_independent,
double* const asymp_jasb);
#+end_src
@ -1854,6 +1965,7 @@ qmckl_exit_code qmckl_compute_jastrow_champ_asymp_jasb_hpc (const qmckl_context
const int64_t bord_num,
const double* b_vector,
const double rescale_factor_ee,
const bool spin_independent,
double* const asymp_jasb );
#+end_src
@ -1863,6 +1975,7 @@ qmckl_compute_jastrow_champ_asymp_jasb_hpc (const qmckl_context context,
const int64_t bord_num,
const double* b_vector,
const double rescale_factor_ee,
const bool spin_independent,
double* const asymp_jasb )
{
@ -1884,7 +1997,7 @@ qmckl_compute_jastrow_champ_asymp_jasb_hpc (const qmckl_context context,
f += b_vector[k]*x;
}
asymp_jasb[0] = 0.5 * asym_one + f;
asymp_jasb[0] = spin_independent ? asym_one + f : 0.5 * asym_one + f;
asymp_jasb[1] = asym_one + f;
return QMCKL_SUCCESS;
@ -1897,15 +2010,16 @@ qmckl_exit_code qmckl_compute_jastrow_champ_asymp_jasb (const qmckl_context cont
const int64_t bord_num,
const double* b_vector,
const double rescale_factor_ee,
const bool spin_independent,
double* const asymp_jasb );
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes
qmckl_exit_code qmckl_compute_jastrow_champ_asymp_jasb (
const qmckl_context context,
qmckl_exit_code qmckl_compute_jastrow_champ_asymp_jasb (const qmckl_context context,
const int64_t bord_num,
const double* b_vector,
const double rescale_factor_ee,
const bool spin_independent,
double* const asymp_jasb )
{
#ifdef HAVE_HPC
@ -1913,7 +2027,7 @@ qmckl_exit_code qmckl_compute_jastrow_champ_asymp_jasb (
#else
return qmckl_compute_jastrow_champ_asymp_jasb_doc
#endif
(context, bord_num, b_vector, rescale_factor_ee, asymp_jasb);
(context, bord_num, b_vector, rescale_factor_ee, spin_independent, asymp_jasb);
}
#+end_src
**** Test
@ -1958,6 +2072,10 @@ int64_t dim_c_vector=0;
assert(!qmckl_jastrow_champ_provided(context));
/* Set the data */
rc = qmckl_check(context,
qmckl_set_jastrow_champ_spin_independent(context, false)
);
rc = qmckl_check(context,
qmckl_set_jastrow_champ_aord_num(context, aord_num)
);
@ -2023,7 +2141,9 @@ for (int i=0 ; i<type_nucl_num ; ++i) {
assert(qmckl_jastrow_champ_provided(context));
double asymp_jasb[2];
rc = qmckl_get_jastrow_champ_asymp_jasb(context, asymp_jasb,2);
rc = qmckl_check(context,
qmckl_get_jastrow_champ_asymp_jasb(context, asymp_jasb,2)
);
// calculate asymp_jasb
assert(fabs(asymp_jasb[0]-0.7115733522582638) < 1.e-12);
@ -2035,6 +2155,8 @@ assert(fabs(asymp_jasb[1]-1.043287918508297 ) < 1.e-12);
Calculate the electron-electron jastrow component ~factor_ee~ using the ~asymp_jasb~
component and the electron-electron rescaled distances ~ee_distance_rescaled~.
If the ~spin_independent~ variable is set to ~true~, then
$\delta^{\uparrow \downarrow}$ is always equal to one.
\[
f_\text{ee} = \sum_{i,j<i} \left[
@ -2185,6 +2307,7 @@ qmckl_exit_code qmckl_provide_jastrow_champ_factor_ee(qmckl_context context)
ctx->jastrow_champ.b_vector,
ctx->jastrow_champ.ee_distance_rescaled,
ctx->jastrow_champ.asymp_jasb,
ctx->jastrow_champ.spin_independent,
ctx->jastrow_champ.factor_ee);
if (rc != QMCKL_SUCCESS) {
return rc;
@ -2221,8 +2344,8 @@ qmckl_exit_code qmckl_provide_jastrow_champ_factor_ee(qmckl_context context)
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
function qmckl_compute_jastrow_champ_factor_ee_doc(context, &
walk_num, elec_num, up_num, bord_num, &
b_vector, ee_distance_rescaled, asymp_jasb, factor_ee) &
walk_num, elec_num, up_num, bord_num, b_vector, &
ee_distance_rescaled, asymp_jasb, spin_independent, factor_ee) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
@ -2237,6 +2360,7 @@ function qmckl_compute_jastrow_champ_factor_ee_doc(context, &
real (c_double ) , intent(in) :: b_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)
logical (c_bool ) , intent(in), value :: spin_independent
real (c_double ) , intent(out) :: factor_ee(walk_num)
integer(qmckl_exit_code) :: info
@ -2266,17 +2390,22 @@ function qmckl_compute_jastrow_champ_factor_ee_doc(context, &
endif
do nw =1, walk_num
factor_ee(nw) = 0.0d0
do j=1,elec_num
do i=1,j-1
x = ee_distance_rescaled(i,j,nw)
if (spin_independent) then
factor_ee(nw) = factor_ee(nw) + b_vector(1) * x / (1.d0 + b_vector(2) * x) - asymp_jasb(2)
else
if ( (j <= up_num).or.(i > up_num) ) then
factor_ee(nw) = factor_ee(nw) + 0.5d0 * b_vector(1) * x / (1.d0 + b_vector(2) * x) - asymp_jasb(1)
else
factor_ee(nw) = factor_ee(nw) + b_vector(1) * x / (1.d0 + b_vector(2) * x) - asymp_jasb(2)
endif
endif
xk = x
do k=2,bord_num
@ -2292,8 +2421,8 @@ end function qmckl_compute_jastrow_champ_factor_ee_doc
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_doc (
const qmckl_context context,
qmckl_exit_code
qmckl_compute_jastrow_champ_factor_ee_doc (const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t up_num,
@ -2301,12 +2430,13 @@ qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_doc (
const double* b_vector,
const double* ee_distance_rescaled,
const double* asymp_jasb,
const bool spin_independent,
double* const factor_ee );
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_hpc (
const qmckl_context context,
qmckl_exit_code
qmckl_compute_jastrow_champ_factor_ee_hpc (const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t up_num,
@ -2314,12 +2444,13 @@ qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_hpc (
const double* b_vector,
const double* ee_distance_rescaled,
const double* asymp_jasb,
const bool spin_independent,
double* const factor_ee );
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes
qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_hpc (
const qmckl_context context,
qmckl_exit_code
qmckl_compute_jastrow_champ_factor_ee_hpc (const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t up_num,
@ -2327,7 +2458,9 @@ qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_hpc (
const double* b_vector,
const double* ee_distance_rescaled,
const double* asymp_jasb,
double* const factor_ee ) {
const bool spin_independent,
double* const factor_ee )
{
if (context == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
@ -2352,6 +2485,18 @@ qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_hpc (
factor_ee[nw] = 0.;
size_t ishift = nw * elec_num * elec_num;
if (spin_independent) {
for (int j = 0; j < elec_num; ++j ) {
const double* xj = &(ee_distance_rescaled[j * elec_num + ishift]);
for (int i = 0; i < j ; ++i) {
factor_ee[nw] += b_vector[0]*xj[i] / (1. + b_vector[1]*xj[i]);
}
}
} else {
for (int j = 0; j < up_num; ++j ) {
const double* xj = &(ee_distance_rescaled[j * elec_num + ishift]);
for (int i = 0; i < j ; ++i) {
@ -2369,6 +2514,7 @@ qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_hpc (
}
}
}
factor_ee[nw] -= fshift;
@ -2393,8 +2539,8 @@ qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_hpc (
# #+CALL: generate_c_header(table=qmckl_factor_ee_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee (
const qmckl_context context,
qmckl_exit_code
qmckl_compute_jastrow_champ_factor_ee (const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t up_num,
@ -2402,6 +2548,7 @@ qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee (
const double* b_vector,
const double* ee_distance_rescaled,
const double* asymp_jasb,
const bool spin_independent,
double* const factor_ee );
#+end_src
@ -2415,6 +2562,7 @@ qmckl_compute_jastrow_champ_factor_ee (const qmckl_context context,
const double* b_vector,
const double* ee_distance_rescaled,
const double* asymp_jasb,
const bool spin_independent,
double* const factor_ee )
{
@ -2424,7 +2572,7 @@ qmckl_compute_jastrow_champ_factor_ee (const qmckl_context context,
return qmckl_compute_jastrow_champ_factor_ee_doc
#endif
(context, walk_num, elec_num, up_num, bord_num, b_vector,
ee_distance_rescaled, asymp_jasb, factor_ee);
ee_distance_rescaled, asymp_jasb, spin_independent, factor_ee);
}
#+end_src
**** Test
@ -3686,7 +3834,7 @@ rc = qmckl_get_jastrow_champ_ee_distance_rescaled_gl(context, ee_distance_rescal
** Electron-nucleus component
*** Asymptotic component for
Calculate the asymptotic component ~asymp_jasa~ to be substracted from the final
Calculate the asymptotic component ~asymp_jasa~ to be subtracted from the final
electron-nucleus jastrow factor \(J_{\text{eN}}\). The asymptotic component is calculated
via the ~a_vector~ and the electron-nucleus rescale factors ~rescale_factor_en~.