1
0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-06-13 16:55:35 +02:00

Updated tests for Jastrow with kappa=0.6

This commit is contained in:
Anthony Scemama 2023-05-25 01:12:05 +02:00
parent 1e4bf9631f
commit edbe33f40f
2 changed files with 231 additions and 298 deletions

View File

@ -157,12 +157,12 @@ integer function qmckl_distance_sq_f(context, transa, transb, m, n, &
return
endif
if (iand(transab,2) == 0 .and. LDA < 3) then
if (iand(transab,2) == 0 .and. LDB < 3) then
info = QMCKL_INVALID_ARG_7
return
endif
if (iand(transab,2) == 2 .and. LDA < m) then
if (iand(transab,2) == 2 .and. LDB < n) then
info = QMCKL_INVALID_ARG_7
return
endif
@ -1328,37 +1328,29 @@ end function test_qmckl_dist_rescaled
\[
C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa
C(r_{ij}) = \left( 1 - \exp(-\kappa\, r_{ij})\right)/\kappa
\]
Here the gradient is defined as follows:
\[
\nabla (C_{ij}(\mathbf{r}_{ee})) = \left(\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta x},\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta y},\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta z} \right)
\nabla_i C(r_{ij}) = \left(\frac{\partial C(r_{ij})}{\partial x_i},\frac{\partial C(r_{ij})}{\partial y_i},\frac{\partial C(r_{ij})}{\partial z_i} \right)
\]
and the Laplacian is defined as follows:
\[
\triangle (C_{ij}(r_{ee})) = \frac{\delta^2}{\delta x^2} + \frac{\delta^2}{\delta y^2} + \frac{\delta^2}{\delta z^2}
\Delta_i C(r_{ij}) = \frac{\partial^2}{\partial x_i^2} + \frac{\partial^2}{\partial y_i^2} + \frac{\partial^2}{\partial z_i^2}
\]
Using the above three formulae, the expression for the gradient and Laplacian is
as follows:
Using the above three formulas, the expression for the gradient and Laplacian is:
\[
\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta x} = \frac{|(x_i - x_j)|}{r_{ij}} (1 - \kappa R_{ij})
\frac{\partial C(r_{ij})}{\partial x_i} = \frac{|(x_i -
x_j)|}{r_{ij}} \exp (- \kappa \, r_{ij})
\]
\[
\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta y} = \frac{|(y_i - y_j)|}{r_{ij}} (1 - \kappa R_{ij})
\]
\[
\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta z} = \frac{|(z_i - z_j)|}{r_{ij}} (1 - \kappa R_{ij})
\]
\[
\Delta(C_{ij}(r_{ee}) = \left[ \frac{2}{r_{ij}} - \kappa \right] (1-\kappa R_{ij})
\Delta C_{ij}(r_{ij}) = \left[ \frac{2}{r_{ij}} - \kappa \right] \exp (- \kappa \, r_{ij})
\]
If the input array is normal (~'N'~), the xyz coordinates are in
@ -1432,11 +1424,9 @@ integer function qmckl_distance_rescaled_gl_f(context, transa, transb, m, n, &
integer*8 :: i,j
real*8 :: x, y, z, dist, dist_inv
real*8 :: rescale_factor_kappa_inv, rij
real*8 :: rij
integer :: transab
rescale_factor_kappa_inv = 1.0d0/rescale_factor_kappa;
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
@ -1486,27 +1476,7 @@ integer function qmckl_distance_rescaled_gl_f(context, transa, transb, m, n, &
return
endif
if (iand(transab,2) == 0 .and. LDA < 3) then
info = QMCKL_INVALID_ARG_7
return
endif
if (iand(transab,2) == 2 .and. LDA < m) then
info = QMCKL_INVALID_ARG_7
return
endif
! check for LDB
if (iand(transab,1) == 0 .and. LDB < 3) then
info = QMCKL_INVALID_ARG_9
return
endif
if (iand(transab,1) == 1 .and. LDB < n) then
info = QMCKL_INVALID_ARG_9
return
endif
if (iand(transab,2) == 0 .and. LDB < 3) then
info = QMCKL_INVALID_ARG_9
return
@ -1535,11 +1505,11 @@ integer function qmckl_distance_rescaled_gl_f(context, transa, transb, m, n, &
z = A(3,i) - B(3,j)
dist = dsqrt(x*x + y*y + z*z)
dist_inv = 1.0d0/dist
rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij)
rij = dexp(-rescale_factor_kappa * dist)
C(1,i,j) = x * dist_inv * rij
C(2,i,j) = y * dist_inv * rij
C(3,i,j) = z * dist_inv * rij
C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa) * rij
end do
end do
@ -1552,11 +1522,11 @@ integer function qmckl_distance_rescaled_gl_f(context, transa, transb, m, n, &
z = A(i,3) - B(3,j)
dist = dsqrt(x*x + y*y + z*z)
dist_inv = 1.0d0/dist
rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij)
rij = dexp(-rescale_factor_kappa * dist)
C(1,i,j) = x * dist_inv * rij
C(2,i,j) = y * dist_inv * rij
C(3,i,j) = z * dist_inv * rij
C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa) * rij
end do
end do
@ -1569,11 +1539,11 @@ integer function qmckl_distance_rescaled_gl_f(context, transa, transb, m, n, &
z = A(3,i) - B(j,3)
dist = dsqrt(x*x + y*y + z*z)
dist_inv = 1.0d0/dist
rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij)
rij = dexp(-rescale_factor_kappa * dist)
C(1,i,j) = x * dist_inv * rij
C(2,i,j) = y * dist_inv * rij
C(3,i,j) = z * dist_inv * rij
C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa) * rij
end do
end do
@ -1586,11 +1556,11 @@ integer function qmckl_distance_rescaled_gl_f(context, transa, transb, m, n, &
z = A(i,3) - B(j,3)
dist = dsqrt(x*x + y*y + z*z)
dist_inv = 1.0d0/dist
rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij)
C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij)
rij = dexp(-rescale_factor_kappa * dist)
C(1,i,j) = x * dist_inv * rij
C(2,i,j) = y * dist_inv * rij
C(3,i,j) = z * dist_inv * rij
C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa) * rij
end do
end do

View File

@ -199,6 +199,9 @@ int main() {
#+BEGIN_SRC python :results none :exports none
import numpy as np
kappa = 0.6
kappa_inv = 1./kappa
# For H2O we have the following data:
elec_num = 10
nucl_num = 2
@ -219,62 +222,20 @@ elec_coord = np.array( [[[-0.250655104764153 , 0.503070975550133 ,
[-0.108090166832043 , 0.189161729653261 , 2.15398313919894],
[ 0.397978144318712 , -0.254277292595981 , 2.54553335476344]]])
ee_distance_rescaled = np.array([
[ 0.000000000000000, 0.000000000000000, 0.000000000000000,
0.000000000000000, 0.000000000000000, 0.000000000000000,
0.000000000000000, 0.000000000000000, 0.000000000000000,
0.000000000000000 ],
[ 0.550227800352402, 0.000000000000000, 0.000000000000000,
0.000000000000000, 0.000000000000000, 0.000000000000000,
0.000000000000000, 0.000000000000000, 0.000000000000000,
0.000000000000000 ],
[ 0.919155060185168, 0.937695909123175, 0.000000000000000,
0.000000000000000, 0.000000000000000, 0.000000000000000,
0.000000000000000, 0.000000000000000, 0.000000000000000,
0.000000000000000 ],
[ 0.893325429242815, 0.851181978173561, 0.978501685226877,
0.000000000000000, 0.000000000000000, 0.000000000000000,
0.000000000000000, 0.000000000000000, 0.000000000000000,
0.000000000000000 ],
[ 0.982457268305353, 0.976125002619471, 0.994349933143149,
0.844077311588328, 0.000000000000000, 0.000000000000000,
0.000000000000000, 0.000000000000000, 0.000000000000000,
0.000000000000000 ],
[ 0.482407528408731, 0.414816073699124, 0.894716035479343,
0.876540187084407, 0.978921170036895, 0.000000000000000,
0.000000000000000, 0.000000000000000, 0.000000000000000,
0.000000000000000 ],
[ 0.459541909660400, 0.545007215761510, 0.883752955884551,
0.918958134888791, 0.986386936267237, 0.362209822236419,
0.000000000000000, 0.000000000000000, 0.000000000000000,
0.000000000000000 ],
[ 0.763732576854455, 0.817282762358449, 0.801802919535959,
0.900089095449775, 0.975704636491453, 0.707836537586060,
0.755705808346586, 0.000000000000000, 0.000000000000000,
0.000000000000000 ],
[ 0.904249454052971, 0.871097965261373, 0.982717262706270,
0.239901207363622, 0.836519456769083, 0.896135326270534,
0.930694340243023, 0.917708540815567, 0.000000000000000,
0.000000000000000 ],
[ 0.944400908070716, 0.922589018494961, 0.984615718580670,
0.514328661540623, 0.692362267147064, 0.931894098453677,
0.956034127544344, 0.931221472309472, 0.540903688625053,
0.000000000000000 ]])
ee_distance_rescaled = np.array(\
[ [(1.-np.exp(-kappa*np.linalg.norm(elec_coord[0,j,:]-elec_coord[0,i,:])))/kappa \
for i in range(elec_num) ]
for j in range(elec_num) ])
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 ]]))
en_distance_rescaled = \
np.array([ [(1.-np.exp(-kappa*np.linalg.norm(elec_coord[0,j,:]-nucl_coord[:,i])))/kappa \
for j in range(elec_num) ]
for i in range(nucl_num) ])
# symmetrize it
for i in range(elec_num):
for j in range(elec_num):
ee_distance_rescaled[i][j] = ee_distance_rescaled[j][i]
#for i in range(elec_num):
# for j in range(elec_num):
# ee_distance_rescaled[i][j] = ee_distance_rescaled[j][i]
# For N2, we have the following data:
type_nucl_num = 1
@ -283,6 +244,7 @@ bord_num = 5
cord_num = 5
dim_c_vector= 23
type_nucl_vector = [ 0, 0]
a_vector = np.array([
[0.000000000000000E+000],
[0.000000000000000E+000],
@ -343,8 +305,6 @@ lkpm_combined_index = [[1 , 1 , 2 , 0],
[3 , 0 , 5 , 1],
[1 , 0 , 5 , 2]]
kappa = 1.0
kappa_inv = 1.0/kappa
#+END_SRC
** Data structure
@ -1573,8 +1533,8 @@ 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;
int64_t nucl_num = n2_nucl_num;
double rescale_factor_ee = 1.0;
double rescale_factor_en[2] = { 1.0, 1.0 };
double rescale_factor_ee = 0.6;
double rescale_factor_en[2] = { 0.6, 0.6 };
double* elec_coord = &(n2_elec_coord[0][0][0]);
const double* nucl_charge = n2_charge;
@ -1980,11 +1940,11 @@ 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.43340325572525706
: asymp_jasb[0] : 0.31567342786262853
: asymp_jasb[1] : 0.5323750557252571
: asym_one : 0.6634291325000664
: asymp_jasb[0] : 0.7115733522582638
: asymp_jasb[1] : 1.043287918508297
#+begin_src c :tangle (eval c_test)
assert(qmckl_electron_provided(context));
@ -2070,8 +2030,8 @@ double asymp_jasb[2];
rc = qmckl_get_jastrow_champ_asymp_jasb(context, asymp_jasb,2);
// calculate asymp_jasb
assert(fabs(asymp_jasb[0]-0.31567342786262853) < 1.e-12);
assert(fabs(asymp_jasb[1]-0.5323750557252571) < 1.e-12);
assert(fabs(asymp_jasb[0]-0.7115733522582638) < 1.e-12);
assert(fabs(asymp_jasb[1]-1.043287918508297 ) < 1.e-12);
#+end_src
@ -2502,31 +2462,32 @@ 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
pow_ser = 0.0
for p in range(1,bord_num):
x = x * ee_distance_rescaled[i][j]
pow_ser = pow_ser + b_vector[p + 1] * x
pow_ser += b_vector[p+1] * x
if(i < up_num or j >= up_num):
spin_fact = 0.5
ipar = 0
else:
ipar = 1
spin_fact = 1.0
factor_ee = factor_ee + spin_fact * b_vector[0] * ee_distance_rescaled[i][j] \
factor_ee += spin_fact * b_vector[0] * ee_distance_rescaled[i][j] \
/ (1.0 + b_vector[1] * ee_distance_rescaled[i][j]) \
- asymp_jasb[ipar] + pow_ser
- 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
: asym_one : 0.6634291325000664
: asymp_jasb[0] : 0.7115733522582638
: asymp_jasb[1] : 1.043287918508297
: factor_ee : -16.83886184243964
#+begin_src c :tangle (eval c_test)
/* Check if Jastrow is properly initialized */
@ -2538,8 +2499,8 @@ rc = qmckl_check(context,
);
// calculate factor_ee
printf("%e\n%e\n\n",factor_ee[0],-4.282760865958113 );
assert(fabs(factor_ee[0]+4.282760865958113) < 1.e-12);
printf("%e\n%e\n\n",factor_ee[0],-16.83886184243964);
assert(fabs(factor_ee[0]+16.83886184243964) < 1.e-12);
#+end_src
@ -3025,7 +2986,7 @@ import numpy as np
<<asymp_jasb>>
kappa = 1.0
kappa = 0.6
dx = 1.e-3
elec_coord = np.array(elec_coord)[0]
@ -3039,6 +3000,8 @@ def make_dist(elec_coord):
return elec_dist
def make_dist_deriv(elec_coord):
elec_dist_d = np.zeros(shape=(4, elec_num, elec_num),dtype=float)
@ -3105,7 +3068,6 @@ for j in range(elec_num):
invden2 = invden * invden
invden3 = invden2 * invden
xinv = 1.0 / x
ipar = 1
dx[:] = ee_distance_rescaled_gl[:,j,i]
@ -3131,13 +3093,13 @@ print("factor_ee_gl[3][0]:",factor_ee_gl[3][0])
#+end_src
#+RESULTS:
: asym_one : 0.43340325572525706
: asymp_jasb[0] : 0.31567342786262853
: asymp_jasb[1] : 0.5323750557252571
: factor_ee_gl[0][0]: -0.16364894652107934
: factor_ee_gl[1][0]: 0.6927548119830084
: factor_ee_gl[2][0]: -0.073267755223968
: factor_ee_gl[3][0]: 1.5111672803213185
: asym_one : 0.6634291325000664
: asymp_jasb[0] : 0.7115733522582638
: asymp_jasb[1] : 1.043287918508297
: factor_ee_gl[0][0]:
: factor_ee_gl[1][0]:
: factor_ee_gl[2][0]:
: factor_ee_gl[3][0]:
#+begin_src c :tangle (eval c_test)
@ -3149,17 +3111,17 @@ double factor_ee_gl[walk_num][4][elec_num];
rc = qmckl_get_jastrow_champ_factor_ee_gl(context, &(factor_ee_gl[0][0][0]),walk_num*4*elec_num);
// check factor_ee_gl
printf("%f %f\n", factor_ee_gl[0][0][0], -0.16364894652107934);
assert(fabs(factor_ee_gl[0][0][0]+0.16364894652107934) < 1.e-12);
printf("%f %f\n", factor_ee_gl[0][0][0], -0.39319353942687446);
assert(fabs(factor_ee_gl[0][0][0]+0.39319353942687446) < 1.e-12);
printf("%f %f\n", factor_ee_gl[0][1][0], 0.6927548119830084 );
assert(fabs(factor_ee_gl[0][1][0]-0.6927548119830084 ) < 1.e-12);
printf("%f %f\n", factor_ee_gl[0][1][0], 1.0535615450668214);
assert(fabs(factor_ee_gl[0][1][0]-1.0535615450668214) < 1.e-12);
printf("%f %f\n", factor_ee_gl[0][2][0],-0.073267755223968 );
assert(fabs(factor_ee_gl[0][2][0]+0.073267755223968 ) < 1.e-12);
printf("%f %f\n", factor_ee_gl[0][2][0],-0.39098406960784515);
assert(fabs(factor_ee_gl[0][2][0]+0.39098406960784515) < 1.e-12);
printf("%f %f\n", factor_ee_gl[0][3][0], 1.5111672803213185 );
assert(fabs(factor_ee_gl[0][3][0]-1.5111672803213185 ) < 1.e-12);
printf("%f %f\n", factor_ee_gl[0][3][0],2.8650469630854483);
assert(fabs(factor_ee_gl[0][3][0]-2.8650469630854483) < 1.e-12);
#+end_src
*** Electron-electron rescaled distances
@ -3371,7 +3333,7 @@ qmckl_exit_code qmckl_compute_ee_distance_rescaled (
#+begin_src python :results output :exports none
import numpy as np
kappa = 1.0
kappa = 0.6
elec_1_w1 = np.array( [-0.250655104764153, 0.503070975550133 , -0.166554344502303])
elec_2_w1 = np.array( [-0.587812193472177, -0.128751981129274 , 0.187773606533075])
@ -3388,11 +3350,11 @@ print ( "[6][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_6_w1-elec_5_w1))
#+RESULTS:
: [0][0] : 0.0
: [0][1] : 0.5502278003524018
: [1][0] : 0.5502278003524018
: [0][1] :
: [1][0] : 0.6347507420688708
: [5][5] : 0.0
: [5][6] : 0.3622098222364193
: [6][5] : 0.3622098222364193
: [5][6] : 0.3941735387855409
: [6][5] :
#+begin_src c :tangle (eval c_test)
assert(qmckl_electron_provided(context));
@ -3409,7 +3371,7 @@ assert(ee_distance_rescaled[0] == 0.);
assert(ee_distance_rescaled[1] == ee_distance_rescaled[elec_num]);
// value of (1,0,0)
assert(fabs(ee_distance_rescaled[1]-0.5502278003524018) < 1.e-12);
assert(fabs(ee_distance_rescaled[1]-0.6347507420688708) < 1.e-12);
// (0,0,1) == 0.
assert(ee_distance_rescaled[5*elec_num + 5] == 0.);
@ -3418,7 +3380,7 @@ assert(ee_distance_rescaled[5*elec_num + 5] == 0.);
assert(ee_distance_rescaled[5*elec_num+6] == ee_distance_rescaled[6*elec_num+5]);
// value of (1,0,1)
assert(fabs(ee_distance_rescaled[5*elec_num+6]-0.3622098222364193) < 1.e-12);
assert(fabs(ee_distance_rescaled[5*elec_num+6]-0.3941735387855409) < 1.e-12);
#+end_src
@ -3914,15 +3876,15 @@ print("asymp_jasa[i] : ", asymp_jasa)
#+end_src
#+RESULTS: asymp_jasa
: asymp_jasa[i] : [0.548554]
: asymp_jasa[i] : [1.75529774]
#+begin_src c :tangle (eval c_test)
double asymp_jasa[2];
rc = qmckl_get_jastrow_champ_asymp_jasa(context, asymp_jasa, type_nucl_num);
// calculate asymp_jasb
printf("%e %e\n", asymp_jasa[0], 0.548554);
assert(fabs(0.548554 - asymp_jasa[0]) < 1.e-12);
printf("%e %e\n", asymp_jasa[0], 1.75529774);
assert(fabs(1.75529774 - asymp_jasa[0]) < 1.e-8);
#+end_src
@ -4296,11 +4258,11 @@ import numpy as np
factor_en = 0.0
for a in range(0,nucl_num):
for i in range(0,elec_num):
x = en_distance_rescaled[i][a]
x = en_distance_rescaled[a][i]
pow_ser = 0.0
for p in range(2,aord_num+1):
x = x * en_distance_rescaled[i][a]
x = x * en_distance_rescaled[a][i]
pow_ser = pow_ser + a_vector[(p-1) + 1][type_nucl_vector[a]] * x
factor_en = factor_en - a_vector[0][type_nucl_vector[a]] * x \
@ -4312,8 +4274,8 @@ print("factor_en :",factor_en)
#+end_src
#+RESULTS:
: asymp_jasa[i] : [0.548554]
: factor_en : -5.1052574308112755
: asymp_jasa[i] : [1.75529774]
: factor_en :
#+begin_src c :tangle (eval c_test)
@ -4324,8 +4286,8 @@ double factor_en[walk_num];
rc = qmckl_get_jastrow_champ_factor_en(context, factor_en,walk_num);
// calculate factor_en
printf("%f %f\n", factor_en[0], -5.1052574308112755);
assert(fabs(-5.1052574308112755 - factor_en[0]) < 1.e-12);
printf("%f %f\n", factor_en[0], -22.781375792083587);
assert(fabs(-22.781375792083587 - factor_en[0]) < 1.e-12);
#+end_src
@ -4375,6 +4337,22 @@ qmckl_get_jastrow_champ_factor_en_gl(qmckl_context context,
}
#+end_src
***** Fortran interface
#+begin_src f90 :tangle (eval fh_func) :comments org
interface
integer(qmckl_exit_code) function qmckl_get_jastrow_champ_factor_en_gl (context, &
factor_en_gl, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (qmckl_context) , intent(in), value :: context
integer(c_int64_t), intent(in), value :: size_max
double precision, intent(out) :: factor_en_gl(size_max)
end function qmckl_get_jastrow_champ_factor_en_gl
end interface
#+end_src
**** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_jastrow_champ_factor_en_gl(qmckl_context context);
@ -4472,17 +4450,17 @@ qmckl_exit_code qmckl_provide_jastrow_champ_factor_en_gl(qmckl_context context)
:END:
#+NAME: qmckl_factor_en_gl_args
| Variable | Type | In/Out | Description |
|--------------------------------+-------------------------------------------+--------+---------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~type_nucl_num~ | ~int64_t~ | in | Number of unique nuclei |
| ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | IDs of unique nuclei |
| ~aord_num~ | ~int64_t~ | in | Number of coefficients |
| ~a_vector~ | ~double[aord_num+1][type_nucl_num]~ | in | List of coefficients |
| ~en_distance_rescaled~ | ~double[walk_num][nucl_num][elec_num]~ | in | Electron-nucleus distances |
| Variable | Type | In/Out | Description |
|---------------------------+-------------------------------------------+--------+---------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~type_nucl_num~ | ~int64_t~ | in | Number of unique nuclei |
| ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | IDs of unique nuclei |
| ~aord_num~ | ~int64_t~ | in | Number of coefficients |
| ~a_vector~ | ~double[aord_num+1][type_nucl_num]~ | in | List of coefficients |
| ~en_distance_rescaled~ | ~double[walk_num][nucl_num][elec_num]~ | in | Electron-nucleus distances |
| ~en_distance_rescaled_gl~ | ~double[walk_num][4][nucl_num][elec_num]~ | in | Electron-nucleus distance derivatives |
| ~factor_en_gl~ | ~double[walk_num][4][elec_num]~ | out | Electron-nucleus jastrow |
@ -4662,7 +4640,7 @@ import numpy as np
<<jastrow_data>>
kappa = 1.0
kappa = 0.6
elec_coord = np.array(elec_coord)[0]
nucl_coord = np.array(nucl_coord)
@ -4682,7 +4660,7 @@ for a in range(nucl_num):
en_distance_rescaled_gl = 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]
f = 1.0 - kappa * en_distance_rescaled[a][i]
for ii in range(4):
en_distance_rescaled_gl[ii][i][a] = elnuc_dist_gl[ii][i][a]
en_distance_rescaled_gl[3][i][a] = en_distance_rescaled_gl[3][i][a] + \
@ -4698,7 +4676,7 @@ 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]
x = en_distance_rescaled[a][i]
if abs(x) < 1e-18:
continue
pow_ser_g = np.zeros(shape=(3),dtype=float)
@ -4715,7 +4693,7 @@ for a in range(nucl_num):
lap2 = 0.0
lap3 = 0.0
for ii in range(3):
x = en_distance_rescaled[i][a]
x = en_distance_rescaled[a][i]
if x < 1e-18:
continue
for p in range(2,aord_num+1):
@ -4723,7 +4701,7 @@ for a in range(nucl_num):
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]
x = x * en_distance_rescaled[a][i]
lap3 = lap3 - 2.0 * a_vector[1][type_nucl_vector[a]] * dx[ii] * dx[ii]
@ -4744,16 +4722,11 @@ print("factor_en_gl[3][0]:",factor_en_gl[3][0])
#+end_src
#+RESULTS:
: factor_en_gl[0][0]: -0.11609919541763383
: factor_en_gl[1][0]: 0.23301394780804574
: factor_en_gl[2][0]: -0.17548337641865783
: factor_en_gl[3][0]: 0.9667363412285741
: factor_en_gl[0][0]: 0.11609919541763383
: factor_en_gl[1][0]: -0.23301394780804574
: factor_en_gl[2][0]: 0.17548337641865783
: factor_en_gl[3][0]: -0.9667363412285741
#+begin_src c :tangle (eval c_test)
@ -4765,10 +4738,10 @@ double factor_en_gl[walk_num][4][elec_num];
rc = qmckl_get_jastrow_champ_factor_en_gl(context, &(factor_en_gl[0][0][0]),walk_num*4*elec_num);
// check factor_en_gl
assert(fabs(factor_en_gl[0][0][0]+0.11609919541763383) < 1.e-12);
assert(fabs(factor_en_gl[0][1][0]-0.23301394780804574) < 1.e-12);
assert(fabs(factor_en_gl[0][2][0]+0.17548337641865783) < 1.e-12);
assert(fabs(factor_en_gl[0][3][0]-0.9667363412285741 ) < 1.e-12);
assert(fabs( -0.19656663796630847 - factor_en_gl[0][0][0]) < 1.e-12);
assert(fabs( 0.3945140890522283 - factor_en_gl[0][1][0]) < 1.e-12);
assert(fabs( -0.5082964671286118 - factor_en_gl[0][2][0]) < 1.e-12);
assert(fabs( 1.8409460670666289 - factor_en_gl[0][3][0]) < 1.e-12);
#+end_src
@ -5034,7 +5007,7 @@ qmckl_exit_code qmckl_compute_en_distance_rescaled (
#+begin_src python :results output :exports none
import numpy as np
kappa = 1.0
kappa = 0.6
elec_1_w1 = np.array( [-0.250655104764153, 0.503070975550133 , -0.166554344502303])
elec_2_w1 = np.array( [-0.587812193472177, -0.128751981129274 , 0.187773606533075])
@ -5053,12 +5026,13 @@ print ( "[0][6] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_6_w1-nucl_1)) )/
#+end_src
#+RESULTS:
: [0][0] : 0.4435709484118112
: [1][0] : 0.8993601506374442
: [0][1] : 0.46760219699910477
: [0][5] : 0.1875631834682101
: [1][5] : 0.8840716589810682
: [0][6] : 0.42640469987268914
: [0][0] : 0.4942158656729477
: [1][0] : 1.2464137498005765
: [0][1] : 0.5248654474756858
: [0][5] : 0.19529459944794733
: [1][5] : 1.2091967687767369
: [0][6] : 0.4726452953409436
#+begin_src c :tangle (eval c_test)
@ -5074,22 +5048,18 @@ assert (rc == QMCKL_SUCCESS);
// (e,n,w) in Fortran notation
// (1,1,1)
assert(fabs(en_distance_rescaled[0][0][0] - 0.4435709484118112) < 1.e-12);
assert(fabs(en_distance_rescaled[0][0][0] - 0.4942158656729477) < 1.e-12);
// (1,2,1)
assert(fabs(en_distance_rescaled[0][1][0] - 0.8993601506374442) < 1.e-12);
assert(fabs(en_distance_rescaled[0][1][0] - 1.2464137498005765) < 1.e-12);
// (2,1,1)
assert(fabs(en_distance_rescaled[0][0][1] - 0.46760219699910477) < 1.e-12);
assert(fabs(en_distance_rescaled[0][0][1] - 0.5248654474756858) < 1.e-12);
// (1,1,2)
assert(fabs(en_distance_rescaled[0][0][5] - 0.1875631834682101) < 1.e-12);
assert(fabs(en_distance_rescaled[0][0][5] - 0.19529459944794733) < 1.e-12);
// (1,2,2)
assert(fabs(en_distance_rescaled[0][1][5] - 0.8840716589810682) < 1.e-12);
assert(fabs(en_distance_rescaled[0][1][5] - 1.2091967687767369) < 1.e-12);
// (2,1,2)
assert(fabs(en_distance_rescaled[0][0][6] - 0.42640469987268914) < 1.e-12);
assert(fabs(en_distance_rescaled[0][0][6] - 0.4726452953409436) < 1.e-12);
#+end_src
@ -5847,7 +5817,7 @@ 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
kappa = 0.6
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
@ -5887,12 +5857,12 @@ 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.017542731694647366
: een_rescaled_e[1, 3, 2] = 0.02214680362033448
: een_rescaled_e[1, 4, 2] = 0.0005700154999202759
: een_rescaled_e[1, 5, 2] = 0.3424402276009091
: een_rescaled_e[0, 2, 1] = 0.2211015082992776
: een_rescaled_e[0, 3, 1] = 0.2611178387068169
: een_rescaled_e[0, 4, 1] = 0.08840123507637472
: een_rescaled_e[1, 3, 2] = 0.10166855073546568
: een_rescaled_e[1, 4, 2] = 0.011311807324686948
: een_rescaled_e[1, 5, 2] = 0.5257156022077619
#+begin_src c :tangle (eval c_test)
assert(qmckl_electron_provided(context));
@ -5902,12 +5872,12 @@ double een_rescaled_e[walk_num][(cord_num + 1)][elec_num][elec_num];
rc = qmckl_get_jastrow_champ_een_rescaled_e(context, &(een_rescaled_e[0][0][0][0]),elec_num*elec_num*(cord_num+1)*walk_num);
// value of (0,2,1)
assert(fabs(een_rescaled_e[0][1][0][2]-0.08084493981483197) < 1.e-12);
assert(fabs(een_rescaled_e[0][1][0][3]-0.1066745707571846) < 1.e-12);
assert(fabs(een_rescaled_e[0][1][0][4]-0.01754273169464735) < 1.e-12);
assert(fabs(een_rescaled_e[0][2][1][3]-0.02214680362033448) < 1.e-12);
assert(fabs(een_rescaled_e[0][2][1][4]-0.0005700154999202759) < 1.e-12);
assert(fabs(een_rescaled_e[0][2][1][5]-0.3424402276009091) < 1.e-12);
assert(fabs(een_rescaled_e[0][1][0][2]- 0.2211015082992776 ) < 1.e-12);
assert(fabs(een_rescaled_e[0][1][0][3]- 0.2611178387068169 ) < 1.e-12);
assert(fabs(een_rescaled_e[0][1][0][4]- 0.0884012350763747 ) < 1.e-12);
assert(fabs(een_rescaled_e[0][2][1][3]- 0.1016685507354656 ) < 1.e-12);
assert(fabs(een_rescaled_e[0][2][1][4]- 0.0113118073246869 ) < 1.e-12);
assert(fabs(een_rescaled_e[0][2][1][5]- 0.5257156022077619 ) < 1.e-12);
#+end_src
*** Electron-electron rescaled distances derivatives in $J_\text{eeN}$
@ -6226,7 +6196,7 @@ for j in range(elec_num):
elec_dist_gl[:, j, j] = 0.0
kappa = 1.0
kappa = 0.6
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
@ -6268,22 +6238,14 @@ for l in range(0,cord_num+1):
for ii in range(0,4):
een_rescaled_e_gl[i,ii,j,l] = een_rescaled_e_gl[i,ii,j,l] * een_rescaled_e[i,j,l]
#print(" een_rescaled_e_gl[1, 1, 3, 1] = ",een_rescaled_e_gl[0, 0, 2, 1])
#print(" een_rescaled_e_gl[1, 1, 4, 1] = ",een_rescaled_e_gl[0, 0, 3, 1])
#print(" een_rescaled_e_gl[1, 1, 5, 1] = ",een_rescaled_e_gl[0, 0, 4, 1])
#print(" een_rescaled_e_gl[2, 1, 4, 2] = ",een_rescaled_e_gl[1, 0, 3, 2])
#print(" een_rescaled_e_gl[2, 1, 5, 2] = ",een_rescaled_e_gl[1, 0, 4, 2])
#print(" een_rescaled_e_gl[2, 1, 6, 2] = ",een_rescaled_e_gl[1, 0, 5, 2])
print(" een_rescaled_e_gl[1, 1, 3, 1] = ",een_rescaled_e_gl[0, 0, 2, 1])
print(" een_rescaled_e_gl[1, 1, 4, 1] = ",een_rescaled_e_gl[0, 0, 3, 1])
print(" een_rescaled_e_gl[1, 1, 5, 1] = ",een_rescaled_e_gl[0, 0, 4, 1])
print(" een_rescaled_e_gl[2, 1, 4, 2] = ",een_rescaled_e_gl[1, 0, 3, 2])
print(" een_rescaled_e_gl[2, 1, 5, 2] = ",een_rescaled_e_gl[1, 0, 4, 2])
print(" een_rescaled_e_gl[2, 1, 6, 2] = ",een_rescaled_e_gl[1, 0, 5, 2])
#+end_src
#+RESULTS: een_e_gl
: een_rescaled_e_gl[1, 1, 3, 1] = 0.05991352796887283
: een_rescaled_e_gl[1, 1, 4, 1] = 0.011714035071545248
: een_rescaled_e_gl[1, 1, 5, 1] = 0.00441398875758468
: een_rescaled_e_gl[2, 1, 4, 2] = 0.013553180060167595
: een_rescaled_e_gl[2, 1, 5, 2] = 0.00041342909359870457
: een_rescaled_e_gl[2, 1, 6, 2] = 0.5880599146214673
#+begin_src c :tangle (eval c_test)
double een_rescaled_e_gl[walk_num][(cord_num + 1)][elec_num][4][elec_num];
size_max=walk_num*(cord_num + 1)*elec_num*4*elec_num;
@ -6291,12 +6253,12 @@ rc = qmckl_get_jastrow_champ_een_rescaled_e_gl(context,
&(een_rescaled_e_gl[0][0][0][0][0]),size_max);
// value of (0,0,0,2,1)
assert(fabs(een_rescaled_e_gl[0][1][0][0][2] + 0.05991352796887283 ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][1][0][0][3] + 0.011714035071545248 ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][1][0][0][4] + 0.00441398875758468 ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][2][1][0][3] + 0.013553180060167595 ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][2][1][0][4] + 0.00041342909359870457) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][2][1][0][5] + 0.5880599146214673 ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][1][0][0][2] + 0.09831391870751387 ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][1][0][0][3] + 0.017204157459682526 ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][1][0][0][4] + 0.013345768421098641 ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][2][1][0][3] + 0.03733086358273962 ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][2][1][0][4] + 0.004922634822943517 ) < 1.e-12);
assert(fabs(een_rescaled_e_gl[0][2][1][0][5] + 0.5416751547830984 ) < 1.e-12);
#+end_src
*** Electron-nucleus rescaled distances in $J_\text{eeN}$
@ -6662,7 +6624,7 @@ 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
kappa = 0.6
een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float)
een_rescaled_n[:,:,0] = 1.0
@ -6685,12 +6647,12 @@ 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
: een_rescaled_n[0, 2, 1] =
: een_rescaled_n[0, 3, 1] =
: een_rescaled_n[0, 4, 1] =
: een_rescaled_n[1, 3, 2] =
: een_rescaled_n[1, 4, 2] =
: een_rescaled_n[1, 5, 2] =
#+begin_src c :tangle (eval c_test)
assert(qmckl_electron_provided(context));
@ -6700,12 +6662,12 @@ size_max=walk_num*(cord_num + 1)*nucl_num*elec_num;
rc = qmckl_get_jastrow_champ_een_rescaled_n(context, &(een_rescaled_n[0][0][0][0]),size_max);
// value of (0,2,1)
assert(fabs(een_rescaled_n[0][1][0][2]-0.10612983920006765) < 1.e-12);
assert(fabs(een_rescaled_n[0][1][0][3]-0.135652809635553) < 1.e-12);
assert(fabs(een_rescaled_n[0][1][0][4]-0.023391817607642338) < 1.e-12);
assert(fabs(een_rescaled_n[0][2][1][3]-0.880957224822116) < 1.e-12);
assert(fabs(een_rescaled_n[0][2][1][4]-0.027185942659395074) < 1.e-12);
assert(fabs(een_rescaled_n[0][2][1][5]-0.01343938025140174) < 1.e-12);
assert(fabs(een_rescaled_n[0][1][0][2]-0.2603169838750542 )< 1.e-12);
assert(fabs(een_rescaled_n[0][1][0][3]-0.3016180139679065 )< 1.e-12);
assert(fabs(een_rescaled_n[0][1][0][4]-0.10506023826192266)< 1.e-12);
assert(fabs(een_rescaled_n[0][2][1][3]-0.9267719759374164 )< 1.e-12);
assert(fabs(een_rescaled_n[0][2][1][4]-0.11497585238132658)< 1.e-12);
assert(fabs(een_rescaled_n[0][2][1][5]-0.07534033469115217)< 1.e-12);
#+end_src
*** Electron-nucleus rescaled distances derivatives in $J_\text{eeN}$
@ -7065,7 +7027,7 @@ for a in range(nucl_num):
elnuc_dist_gl[ii, i, a] = (elec_coord[i][ii] - nucl_coord[ii][a]) * rij_inv
elnuc_dist_gl[3, i, a] = 2.0 * rij_inv
kappa = 1.0
kappa = 0.6
een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float)
een_rescaled_n[:,:,0] = 1.0
@ -7103,12 +7065,12 @@ print(" een_rescaled_n_gl[2, 1, 6, 2] = ",een_rescaled_n_gl[5, 0, 1, 2])
#+end_src
#+RESULTS:
: een_rescaled_n_gl[1, 1, 3, 1] = -0.07633444246999128
: een_rescaled_n_gl[1, 1, 4, 1] = 0.00033282346259738276
: een_rescaled_n_gl[1, 1, 5, 1] = -0.004775370547333061
: een_rescaled_n_gl[2, 1, 4, 2] = 0.1362654644223866
: een_rescaled_n_gl[2, 1, 5, 2] = -0.0231253431662794
: een_rescaled_n_gl[2, 1, 6, 2] = 0.001593334817691633
: een_rescaled_n_gl[1, 1, 3, 1] =
: een_rescaled_n_gl[1, 1, 4, 1] =
: een_rescaled_n_gl[1, 1, 5, 1] =
: een_rescaled_n_gl[2, 1, 4, 2] =
: een_rescaled_n_gl[2, 1, 5, 2] =
: een_rescaled_n_gl[2, 1, 6, 2] =
#+begin_src c :tangle (eval c_test)
assert(qmckl_electron_provided(context));
@ -7118,12 +7080,12 @@ size_max=walk_num*(cord_num + 1)*nucl_num*4*elec_num;
rc = qmckl_get_jastrow_champ_een_rescaled_n_gl(context, &(een_rescaled_n_gl[0][0][0][0][0]),size_max);
// value of (0,2,1)
assert(fabs(een_rescaled_n_gl[0][1][0][0][2]+0.07633444246999128 ) < 1.e-12);
assert(fabs(een_rescaled_n_gl[0][1][0][0][3]-0.00033282346259738276) < 1.e-12);
assert(fabs(een_rescaled_n_gl[0][1][0][0][4]+0.004775370547333061 ) < 1.e-12);
assert(fabs(een_rescaled_n_gl[0][2][1][0][3]-0.1362654644223866 ) < 1.e-12);
assert(fabs(een_rescaled_n_gl[0][2][1][0][4]+0.0231253431662794 ) < 1.e-12);
assert(fabs(een_rescaled_n_gl[0][2][1][0][5]-0.001593334817691633 ) < 1.e-12);
assert(fabs( -0.11234061209936878 - een_rescaled_n_gl[0][1][0][0][2]) < 1.e-12);
assert(fabs( 0.0004440109367151707 - een_rescaled_n_gl[0][1][0][0][3]) < 1.e-12);
assert(fabs( -0.012868642597346566 - een_rescaled_n_gl[0][1][0][0][4]) < 1.e-12);
assert(fabs( 0.08601122289922644 - een_rescaled_n_gl[0][2][1][0][3]) < 1.e-12);
assert(fabs( -0.058681563677207206 - een_rescaled_n_gl[0][2][1][0][4]) < 1.e-12);
assert(fabs( 0.005359281880312882 - een_rescaled_n_gl[0][2][1][0][5]) < 1.e-12);
#+end_src
@ -8390,7 +8352,7 @@ 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
kappa = 0.6
een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float)
een_rescaled_n[:,:,0] = 1.0
@ -8409,7 +8371,7 @@ 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
kappa = 0.6
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
@ -8454,11 +8416,11 @@ rc = qmckl_get_jastrow_champ_tmp_c(context, &(tmp_c[0][0][0][0][0]));
double dtmp_c[walk_num][cord_num][cord_num+1][nucl_num][4][elec_num];
rc = qmckl_get_jastrow_champ_dtmp_c(context, &(dtmp_c[0][0][0][0][0][0]));
printf("%e\n%e\n", tmp_c[0][0][1][0][0], 2.7083473948352403);
assert(fabs(tmp_c[0][0][1][0][0] - 2.7083473948352403) < 1e-12);
printf("%e\n%e\n", tmp_c[0][0][1][0][0], 3.954384);
assert(fabs(tmp_c[0][0][1][0][0] - 3.954384) < 1e-6);
printf("%e\n%e\n", dtmp_c[0][1][0][0][0][0],0.237440520852232);
assert(fabs(dtmp_c[0][1][0][0][0][0] - 0.237440520852232) < 1e-12);
printf("%e\n%e\n", dtmp_c[0][1][0][0][0][0],3.278657e-01);
assert(fabs(dtmp_c[0][1][0][0][0][0] - 3.278657e-01 ) < 1e-6);
#+end_src
*** Electron-electron-nucleus Jastrow $f_{een}$
@ -8993,7 +8955,7 @@ import numpy as np
<<helper_funcs>>
kappa = 1.0
kappa = 0.6
factor_een = 0.0
@ -9018,7 +8980,7 @@ print("factor_een:",factor_een)
#+end_src
#+RESULTS:
: factor_een: -0.37407972141304213
: factor_een: -0.382580260174321
#+begin_src c :tangle (eval c_test)
@ -9028,7 +8990,7 @@ assert(qmckl_jastrow_champ_provided(context));
double factor_een[walk_num];
rc = qmckl_get_jastrow_champ_factor_een(context, &(factor_een[0]),walk_num);
assert(fabs(factor_een[0] + 0.37407972141304213) < 1e-12);
assert(fabs(factor_een[0] + 0.382580260174321) < 1e-12);
#+end_src
*** Electron-electron-nucleus Jastrow $f_{een}$ derivative
@ -9798,7 +9760,7 @@ import numpy as np
<<helper_funcs>>
kappa = 1.0
kappa = 0.6
factor_een = 0.0
@ -9816,27 +9778,28 @@ for n in range(0, dim_c_vector):
for j in range(0, elec_num):
accu = 0.0
accu2 = 0.0
daccu = 0.0
daccu2 = 0.0
for i in range(0, elec_num):
accu = accu + een_rescaled_e[i,j,k] * \
een_rescaled_n[a,i,m]
accu2 = accu2 + een_rescaled_e[i,j,k] * \
een_rescaled_n[a,i,m+l]
# daccu[0:4] = daccu[0:4] + een_rescaled_e_gl_t[k,j,0:4,i,k] * \
# een_rescaled_n[a,i,m]
# daccu[0:4] = daccu[0:4] + een_rescaled_e_gl_t[k,j,0:4,i,k] * \
# een_rescaled_n[a,i,m]
accu2 = accu2 + accu * een_rescaled_n[a,j,m+l]
# factor_een = factor_een + accu2 * cn
accu2 = accu2 + accu * een_rescaled_n[a,j,m+l]
factor_een = factor_een + accu2 * cn
print("factor_een:",factor_een)
#+end_src
#+RESULTS:
: een_rescaled_e_gl[1, 1, 3, 1] = 0.09831391870751387
: een_rescaled_e_gl[1, 1, 4, 1] = 0.017204157459682526
: een_rescaled_e_gl[1, 1, 5, 1] = 0.013345768421098641
: een_rescaled_e_gl[2, 1, 4, 2] = 0.03733086358273962
: een_rescaled_e_gl[2, 1, 5, 2] = 0.004922634822943517
: een_rescaled_e_gl[2, 1, 6, 2] = 0.5416751547830984
: (6, 10, 4, 10)
: factor_een: 0.0
: factor_een: -14.956095654486404
#+begin_src c :tangle (eval c_test)
@ -9847,16 +9810,16 @@ double factor_een_gl[4][walk_num][elec_num];
rc = qmckl_get_jastrow_champ_factor_een_gl(context, &(factor_een_gl[0][0][0]),4*walk_num*elec_num);
printf("%20.15e\n", factor_een_gl[0][0][0]);
assert(fabs(factor_een_gl[0][0][0] - (-5.481671107220383e-04)) < 1e-12);
assert(fabs(8.967809309100624e-02 - factor_een_gl[0][0][0]) < 1e-12);
printf("%20.15e\n", factor_een_gl[1][0][1]);
assert(fabs(factor_een_gl[1][0][1] - (-5.402107832095666e-02)) < 1e-12);
assert(fabs(-3.401545636077585e-02 - factor_een_gl[1][0][1]) < 1e-12);
printf("%20.15e\n", factor_een_gl[2][0][2]);
assert(fabs(factor_een_gl[2][0][2] - (-1.648945927082279e-01)) < 1e-12);
assert(fabs(-2.631321052321952e-01 - factor_een_gl[2][0][2]) < 1e-12);
printf("%20.15e\n", factor_een_gl[3][0][3]);
assert(fabs(factor_een_gl[3][0][3] - (-1.269746119491287e+00)) < 1e-12);
assert(fabs(-1.016785559040419e+00 - factor_een_gl[3][0][3]) < 1e-12);
#+end_src
** Total Jastrow