TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2025-01-11 13:38:28 +01:00
Code Issues Projects Releases Wiki Activity

Merge branch 'master' of github.com:TREX-CoE/qmckl

Browse Source
This commit is contained in:
Anthony Scemama 2021-09-28 22:36:44 +02:00
commit 92b425fe85
6 changed files with 1629 additions and 360 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
org/qmckl_ao.org
View File

@ -14,10 +14,10 @@ Gaussian ($p=2$):
\exp \left( - \gamma_{ks} | \mathbf{r}-\mathbf{R}_A | ^p \right). \exp \left( - \gamma_{ks} | \mathbf{r}-\mathbf{R}_A | ^p \right).
\] \]
In the case of Gaussian functions, $n_s$ is always zero. In the case of Gaussian functions, $n_s$ is always zero. The
The normalization factor $\mathcal{N}_s$ ensures that all the functions normalization factor $\mathcal{N}_s$ ensures that all the functions of
of the shell are normalized to unity. Usually, basis sets are given the shell are normalized (integrate) to unity. Usually, basis sets are
a combination of normalized primitives, so the normalization given a combination of normalized primitives, so the normalization
coefficients of the primitives, $f_{ks}$, need also to be provided. coefficients of the primitives, $f_{ks}$, need also to be provided.
Atomic orbitals (AOs) are defined as Atomic orbitals (AOs) are defined as
@ -26,10 +26,11 @@ Atomic orbitals (AOs) are defined as
\chi_i (\mathbf{r}) = \mathcal{M}_i\, P_{\eta(i)}(\mathbf{r})\, R_{\theta(i)} (\mathbf{r}) \chi_i (\mathbf{r}) = \mathcal{M}_i\, P_{\eta(i)}(\mathbf{r})\, R_{\theta(i)} (\mathbf{r})
\] \]
where $\theta(i)$ returns the shell on which the AO is expanded, where $\theta(i)$ returns the shell on which the AO is expanded, and
and $\eta(i)$ denotes which angular function is chosen. $\eta(i)$ denotes which angular function is chosen and $P$ are the
Here, the parameter $\mathcal{M}_i$ is an extra parameter which allows generating functions of the given angular momentum $\eta(i)$. Here,
the normalization of the different functions of the same shell to be the parameter $\mathcal{M}_i$ is an extra parameter which allows the
normalization of the different functions of the same shell to be
different, as in GAMESS for example. different, as in GAMESS for example.
In this section we describe first how the basis set is stored in the In this section we describe first how the basis set is stored in the

606
org/qmckl_jastrow.org
View File

@ -202,7 +202,7 @@ for i in range(elec_num):
type_nucl_num = 1 type_nucl_num = 1
aord_num = 5 aord_num = 5
bord_num = 5 bord_num = 5
cord_num = 23 cord_num = 5
dim_cord_vect= 23 dim_cord_vect= 23
type_nucl_vector = [ 1, 1] type_nucl_vector = [ 1, 1]
aord_vector = [ aord_vector = [
@ -609,7 +609,7 @@ qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, doub
} }
assert (ctx->jastrow.cord_vector != NULL); assert (ctx->jastrow.cord_vector != NULL);
memcpy(cord_vector, ctx->jastrow.cord_vector, ctx->jastrow.cord_num*sizeof(double)); memcpy(cord_vector, ctx->jastrow.cord_vector, ctx->jastrow.dim_cord_vect*sizeof(double));
return QMCKL_SUCCESS; return QMCKL_SUCCESS;
} }
@ -860,19 +860,22 @@ qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double cons
int32_t mask = 1 << 5; int32_t mask = 1 << 5;
int64_t cord_num; qmckl_exit_code rc = qmckl_provide_dim_cord_vect(context);
qmckl_exit_code rc = qmckl_get_jastrow_cord_num(context, &cord_num); if (rc != QMCKL_SUCCESS) return rc;
int64_t dim_cord_vect;
rc = qmckl_get_jastrow_dim_cord_vect(context, &dim_cord_vect);
if (rc != QMCKL_SUCCESS) return rc; if (rc != QMCKL_SUCCESS) return rc;
int64_t type_nucl_num; int64_t type_nucl_num;
rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num); rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num);
if (rc != QMCKL_SUCCESS) return rc; if (rc != QMCKL_SUCCESS) return rc;
if (cord_num == 0) { if (dim_cord_vect == 0) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_FAILURE, QMCKL_FAILURE,
"qmckl_set_jastrow_coefficient", "qmckl_set_jastrow_coefficient",
"cord_num is not set"); "dim_cord_vect is not set");
} }
if (cord_vector == NULL) { if (cord_vector == NULL) {
@ -892,7 +895,7 @@ qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double cons
} }
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = cord_num * type_nucl_num * sizeof(double); mem_info.size = dim_cord_vect * type_nucl_num * sizeof(double);
double* new_array = (double*) qmckl_malloc(context, mem_info); double* new_array = (double*) qmckl_malloc(context, mem_info);
if(new_array == NULL) { if(new_array == NULL) {
@ -1324,20 +1327,20 @@ end function qmckl_compute_asymp_jasb_f
#+CALL: generate_c_header(table=qmckl_asymp_jasb_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+CALL: generate_c_header(table=qmckl_asymp_jasb_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+RESULTS: #+RESULTS:
#+begin_src c :tangle (eval h_private_func) :comments org #+BEGIN_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_asymp_jasb ( qmckl_exit_code qmckl_compute_asymp_jasb (
const qmckl_context context, const qmckl_context context,
const int64_t bord_num, const int64_t bord_num,
const double* bord_vector, const double* bord_vector,
const double rescale_factor_kappa_ee, const double rescale_factor_kappa_ee,
double* const asymp_jasb ); double* const asymp_jasb );
#+end_src #+END_src
#+CALL: generate_c_interface(table=qmckl_asymp_jasb_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+CALL: generate_c_interface(table=qmckl_asymp_jasb_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+RESULTS: #+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none #+BEGIN_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_asymp_jasb & integer(c_int32_t) function qmckl_compute_asymp_jasb &
(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) & (context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) &
bind(C) result(info) bind(C) result(info)
@ -1356,7 +1359,7 @@ end function qmckl_compute_asymp_jasb_f
(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) (context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb)
end function qmckl_compute_asymp_jasb end function qmckl_compute_asymp_jasb
#+end_src #+END_src
*** Test *** Test
#+name: asymp_jasb #+name: asymp_jasb
@ -1380,11 +1383,6 @@ print("asymp_jasb[1] : ", asymp_jasb[1])
#+end_src #+end_src
#+RESULTS: asymp_jasb #+RESULTS: asymp_jasb
: asym_one : 0.6634291325000664
: asymp_jasb[0] : 1.043287918508297
: asymp_jasb[1] : 0.7115733522582638
#+RESULTS:
: asym_one : 0.43340325572525706 : asym_one : 0.43340325572525706
: asymp_jasb[0] : 0.5323750557252571 : asymp_jasb[0] : 0.5323750557252571
: asymp_jasb[1] : 0.31567342786262853 : asymp_jasb[1] : 0.31567342786262853
@ -1416,6 +1414,8 @@ rc = qmckl_set_jastrow_aord_vector(context, aord_vector);
assert(rc == QMCKL_SUCCESS); assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_jastrow_bord_vector(context, bord_vector); rc = qmckl_set_jastrow_bord_vector(context, bord_vector);
assert(rc == QMCKL_SUCCESS); 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); rc = qmckl_set_jastrow_cord_vector(context, cord_vector);
assert(rc == QMCKL_SUCCESS); assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_jastrow_dependencies(context); rc = qmckl_set_jastrow_dependencies(context);
@ -2114,28 +2114,16 @@ 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[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[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[3][0]:",factor_ee_deriv_e[3][0])
print(factor_ee_deriv_e)
#+end_src #+end_src
#+RESULTS: #+RESULTS:
#+begin_example : asym_one : 0.43340325572525706
asym_one : 0.43340325572525706 : asymp_jasb[0] : 0.5323750557252571
asymp_jasb[0] : 0.5323750557252571 : asymp_jasb[1] : 0.31567342786262853
asymp_jasb[1] : 0.31567342786262853 : factor_ee_deriv_e[0][0]: 0.16364894652107934
factor_ee_deriv_e[0][0]: 0.16364894652107934 : factor_ee_deriv_e[1][0]: -0.6927548119830084
factor_ee_deriv_e[1][0]: -0.6927548119830084 : factor_ee_deriv_e[2][0]: 0.073267755223968
factor_ee_deriv_e[2][0]: 0.073267755223968 : factor_ee_deriv_e[3][0]: 1.5111672803213185
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) #+begin_src c :tangle (eval c_test)
@ -3039,6 +3027,13 @@ integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cor
end do end do
end do end do
end do end do
do l = 0, cord_num
do j = 1, elec_num
een_rescaled_e(l, j, j, nw) = 0.0d0
end do
end do
end do end do
end function qmckl_compute_een_rescaled_e_f end function qmckl_compute_een_rescaled_e_f
@ -3124,6 +3119,10 @@ for l in range(1,cord_num+1):
een_rescaled_e[j, i, l] = x een_rescaled_e[j, i, l] = x
k = k + 1 k = k + 1
for l in range(0,cord_num+1):
for j in range(0, elec_num):
een_rescaled_e[j,j,l] = 0.0
print(" een_rescaled_e[0, 2, 1] = ",een_rescaled_e[0, 2, 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, 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[0, 4, 1] = ",een_rescaled_e[0, 4, 1])
@ -3135,7 +3134,7 @@ print(" een_rescaled_e[1, 5, 2] = ",een_rescaled_e[1, 5, 2])
#+RESULTS: #+RESULTS:
: een_rescaled_e[0, 2, 1] = 0.08084493981483197 : een_rescaled_e[0, 2, 1] = 0.08084493981483197
: een_rescaled_e[0, 3, 1] = 0.1066745707571846 : een_rescaled_e[0, 3, 1] = 0.1066745707571846
: een_rescaled_e[0, 4, 1] = 0.01754273169464735 : een_rescaled_e[0, 4, 1] = 0.017542731694647366
: een_rescaled_e[1, 3, 2] = 0.02214680362033448 : een_rescaled_e[1, 3, 2] = 0.02214680362033448
: een_rescaled_e[1, 4, 2] = 0.0005700154999202759 : een_rescaled_e[1, 4, 2] = 0.0005700154999202759
: een_rescaled_e[1, 5, 2] = 0.3424402276009091 : een_rescaled_e[1, 5, 2] = 0.3424402276009091
@ -3159,9 +3158,10 @@ assert(fabs(een_rescaled_e[0][1][5][2]-0.3424402276009091) < 1.e-12);
** Electron-electron rescaled distances for each order and derivatives ** Electron-electron rescaled distances for each order and derivatives
~een_rescaled_e~ stores the table of the rescaled distances between all ~een_rescaled_e_deriv_e~ stores the table of the derivatives of the
pairs of electrons and raised to the power \(p\) defined by ~cord_num~. rescaled distances between all pairs of electrons and raised to the
Here we take its derivatives required for the een jastrow. power \(p\) defined by ~cord_num~. Here we take its derivatives
required for the een jastrow.
TODO: write formulae TODO: write formulae
@ -3419,7 +3419,7 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f
#+end_src #+end_src
*** Test *** Test
#+name: een_e_deriv_e
#+begin_src python :results output :exports none :noweb yes #+begin_src python :results output :exports none :noweb yes
import numpy as np import numpy as np
@ -3431,6 +3431,16 @@ for i in range(elec_num):
for j 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[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
kappa = 1.0 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 = np.zeros(shape=(elec_num * (elec_num - 1)//2, cord_num+1), dtype=float)
@ -3458,25 +3468,49 @@ for l in range(1,cord_num+1):
een_rescaled_e[j, i, l] = x een_rescaled_e[j, i, l] = x
k = k + 1 k = k + 1
print(" een_rescaled_e[0, 2, 1] = ",een_rescaled_e[0, 2, 1]) een_rescaled_e_deriv_e = np.zeros(shape=(elec_num,4,elec_num,cord_num+1),dtype=float)
print(" een_rescaled_e[0, 3, 1] = ",een_rescaled_e[0, 3, 1]) for l in range(0,cord_num+1):
print(" een_rescaled_e[0, 4, 1] = ",een_rescaled_e[0, 4, 1]) kappa_l = -1.0 * kappa * l
print(" een_rescaled_e[1, 3, 2] = ",een_rescaled_e[1, 3, 2]) for j in range(0,elec_num):
print(" een_rescaled_e[1, 4, 2] = ",een_rescaled_e[1, 4, 2]) for i in range(0,elec_num):
print(" een_rescaled_e[1, 5, 2] = ",een_rescaled_e[1, 5, 2]) for ii in range(0,4):
een_rescaled_e_deriv_e[i,ii,j,l] = kappa_l * elec_dist_deriv_e[ii,i,j]
een_rescaled_e_deriv_e[i,3,j,l] = een_rescaled_e_deriv_e[i,3,j,l] + \
een_rescaled_e_deriv_e[i,0,j,l] * een_rescaled_e_deriv_e[i,0,j,l] + \
een_rescaled_e_deriv_e[i,1,j,l] * een_rescaled_e_deriv_e[i,1,j,l] + \
een_rescaled_e_deriv_e[i,2,j,l] * een_rescaled_e_deriv_e[i,2,j,l]
for ii in range(0,4):
een_rescaled_e_deriv_e[i,ii,j,l] = een_rescaled_e_deriv_e[i,ii,j,l] * een_rescaled_e[i,j,l]
#print(" een_rescaled_e_deriv_e[1, 1, 3, 1] = ",een_rescaled_e_deriv_e[0, 0, 2, 1])
#print(" een_rescaled_e_deriv_e[1, 1, 4, 1] = ",een_rescaled_e_deriv_e[0, 0, 3, 1])
#print(" een_rescaled_e_deriv_e[1, 1, 5, 1] = ",een_rescaled_e_deriv_e[0, 0, 4, 1])
#print(" een_rescaled_e_deriv_e[2, 1, 4, 2] = ",een_rescaled_e_deriv_e[1, 0, 3, 2])
#print(" een_rescaled_e_deriv_e[2, 1, 5, 2] = ",een_rescaled_e_deriv_e[1, 0, 4, 2])
#print(" een_rescaled_e_deriv_e[2, 1, 6, 2] = ",een_rescaled_e_deriv_e[1, 0, 5, 2])
#+end_src #+end_src
#+RESULTS: #+RESULTS: een_e_deriv_e
: een_rescaled_e[0, 2, 1] = 0.08084493981483197 : een_rescaled_e_deriv_e[1, 1, 3, 1] = 0.05991352796887283
: een_rescaled_e[0, 3, 1] = 0.1066745707571846 : een_rescaled_e_deriv_e[1, 1, 4, 1] = 0.011714035071545248
: een_rescaled_e[0, 4, 1] = 0.01754273169464735 : een_rescaled_e_deriv_e[1, 1, 5, 1] = 0.00441398875758468
: een_rescaled_e[1, 3, 2] = 0.02214680362033448 : een_rescaled_e_deriv_e[2, 1, 4, 2] = 0.013553180060167595
: een_rescaled_e[1, 4, 2] = 0.0005700154999202759 : een_rescaled_e_deriv_e[2, 1, 5, 2] = 0.00041342909359870457
: een_rescaled_e[1, 5, 2] = 0.3424402276009091 : een_rescaled_e_deriv_e[2, 1, 6, 2] = 0.5880599146214673
#+begin_src c :tangle (eval c_test) #+begin_src c :tangle (eval c_test)
//assert(qmckl_electron_provided(context)); //assert(qmckl_electron_provided(context));
double een_rescaled_e_deriv_e[walk_num][elec_num][4][elec_num][(cord_num + 1)];
rc = qmckl_get_jastrow_een_rescaled_e_deriv_e(context, &(een_rescaled_e_deriv_e[0][0][0][0][0]));
// value of (0,0,0,2,1)
assert(fabs(een_rescaled_e_deriv_e[0][0][0][2][1] + 0.05991352796887283 ) < 1.e-12);
assert(fabs(een_rescaled_e_deriv_e[0][0][0][3][1] + 0.011714035071545248 ) < 1.e-12);
assert(fabs(een_rescaled_e_deriv_e[0][0][0][4][1] + 0.00441398875758468 ) < 1.e-12);
assert(fabs(een_rescaled_e_deriv_e[0][1][0][3][2] + 0.013553180060167595 ) < 1.e-12);
assert(fabs(een_rescaled_e_deriv_e[0][1][0][4][2] + 0.00041342909359870457) < 1.e-12);
assert(fabs(een_rescaled_e_deriv_e[0][1][0][5][2] + 0.5880599146214673 ) < 1.e-12);
#+end_src #+end_src
** Electron-nucleus rescaled distances for each order ** Electron-nucleus rescaled distances for each order
@ -4076,6 +4110,14 @@ for i in range(elec_num):
for a in range(nucl_num): for a in range(nucl_num):
elnuc_dist[i, a] = np.linalg.norm(elec_coord[i] - nucl_coord[:,a]) elnuc_dist[i, a] = np.linalg.norm(elec_coord[i] - nucl_coord[:,a])
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
kappa = 1.0 kappa = 1.0
een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float) een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float)
@ -4090,24 +4132,50 @@ for l in range(2,cord_num+1):
for i in range(elec_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] 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]) een_rescaled_n_deriv_e = np.zeros(shape=(elec_num,4,nucl_num,cord_num+1),dtype=float)
print(" een_rescaled_n[0, 3, 1] = ",een_rescaled_n[0, 3, 1]) for l in range(0,cord_num+1):
print(" een_rescaled_n[0, 4, 1] = ",een_rescaled_n[0, 4, 1]) kappa_l = -1.0 * kappa * l
print(" een_rescaled_n[1, 3, 2] = ",een_rescaled_n[1, 3, 2]) for j in range(0,elec_num):
print(" een_rescaled_n[1, 4, 2] = ",een_rescaled_n[1, 4, 2]) for a in range(0,nucl_num):
print(" een_rescaled_n[1, 5, 2] = ",een_rescaled_n[1, 5, 2]) for ii in range(0,4):
een_rescaled_n_deriv_e[j,ii,a,l] = kappa_l * elnuc_dist_deriv_e[ii,j,a]
een_rescaled_n_deriv_e[j,3,a,l] = een_rescaled_n_deriv_e[j,3,a,l] + \
een_rescaled_n_deriv_e[j,0,a,l] * een_rescaled_n_deriv_e[j,0,a,l] + \
een_rescaled_n_deriv_e[j,1,a,l] * een_rescaled_n_deriv_e[j,1,a,l] + \
een_rescaled_n_deriv_e[j,2,a,l] * een_rescaled_n_deriv_e[j,2,a,l]
for ii in range(0,4):
een_rescaled_n_deriv_e[j,ii,a,l] = een_rescaled_n_deriv_e[j,ii,a,l] * een_rescaled_n[a,j,l]
print(" een_rescaled_n_deriv_e[1, 1, 3, 1] = ",een_rescaled_n_deriv_e[2, 0, 0, 1])
print(" een_rescaled_n_deriv_e[1, 1, 4, 1] = ",een_rescaled_n_deriv_e[3, 0, 0, 1])
print(" een_rescaled_n_deriv_e[1, 1, 5, 1] = ",een_rescaled_n_deriv_e[4, 0, 0, 1])
print(" een_rescaled_n_deriv_e[2, 1, 4, 2] = ",een_rescaled_n_deriv_e[3, 0, 1, 2])
print(" een_rescaled_n_deriv_e[2, 1, 5, 2] = ",een_rescaled_n_deriv_e[4, 0, 1, 2])
print(" een_rescaled_n_deriv_e[2, 1, 6, 2] = ",een_rescaled_n_deriv_e[5, 0, 1, 2])
#+end_src #+end_src
#+RESULTS: #+RESULTS:
: een_rescaled_n[0, 2, 1] = 0.10612983920006765 : een_rescaled_n_deriv_e[1, 1, 3, 1] = -0.07633444246999128
: een_rescaled_n[0, 3, 1] = 0.135652809635553 : een_rescaled_n_deriv_e[1, 1, 4, 1] = 0.00033282346259738276
: een_rescaled_n[0, 4, 1] = 0.023391817607642338 : een_rescaled_n_deriv_e[1, 1, 5, 1] = -0.004775370547333061
: een_rescaled_n[1, 3, 2] = 0.880957224822116 : een_rescaled_n_deriv_e[2, 1, 4, 2] = 0.1362654644223866
: een_rescaled_n[1, 4, 2] = 0.027185942659395074 : een_rescaled_n_deriv_e[2, 1, 5, 2] = -0.0231253431662794
: een_rescaled_n[1, 5, 2] = 0.01343938025140174 : een_rescaled_n_deriv_e[2, 1, 6, 2] = 0.001593334817691633
#+begin_src c :tangle (eval c_test) #+begin_src c :tangle (eval c_test)
//assert(qmckl_electron_provided(context)); assert(qmckl_electron_provided(context));
double een_rescaled_n_deriv_e[walk_num][elec_num][4][nucl_num][(cord_num + 1)];
rc = qmckl_get_jastrow_een_rescaled_n_deriv_e(context, &(een_rescaled_n_deriv_e[0][0][0][0][0]));
// value of (0,2,1)
assert(fabs(een_rescaled_n_deriv_e[0][2][0][0][1]+0.07633444246999128 ) < 1.e-12);
assert(fabs(een_rescaled_n_deriv_e[0][3][0][0][1]-0.00033282346259738276) < 1.e-12);
assert(fabs(een_rescaled_n_deriv_e[0][4][0][0][1]+0.004775370547333061 ) < 1.e-12);
assert(fabs(een_rescaled_n_deriv_e[0][3][0][1][2]-0.1362654644223866 ) < 1.e-12);
assert(fabs(een_rescaled_n_deriv_e[0][4][0][1][2]+0.0231253431662794 ) < 1.e-12);
assert(fabs(een_rescaled_n_deriv_e[0][5][0][1][2]-0.001593334817691633 ) < 1.e-12);
#+end_src #+end_src
@ -4264,7 +4332,6 @@ qmckl_exit_code qmckl_provide_cord_vect_full(qmckl_context context)
qmckl_exit_code rc = qmckl_exit_code rc =
qmckl_compute_cord_vect_full(context, qmckl_compute_cord_vect_full(context,
ctx->nucleus.num, ctx->nucleus.num,
ctx->jastrow.cord_num,
ctx->jastrow.dim_cord_vect, ctx->jastrow.dim_cord_vect,
ctx->jastrow.type_nucl_num, ctx->jastrow.type_nucl_num,
ctx->jastrow.type_nucl_vector, ctx->jastrow.type_nucl_vector,
@ -4426,26 +4493,24 @@ end function qmckl_compute_dim_cord_vect_f
#+NAME: qmckl_factor_cord_vect_full_args #+NAME: qmckl_factor_cord_vect_full_args
| qmckl_context | context | in | Global state | | qmckl_context | context | in | Global state |
| int64_t | nucl_num | in | Number of atoms | | 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 | dim_cord_vect | in | dimension of cord full table |
| int64_t | type_nucl_num | 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 | | 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_vector[dim_cord_vect][type_nucl_num] | in | dimension of cord full table |
| double | cord_vect_full[dim_cord_vect][nucl_num] | out | Full list of coefficients | | double | cord_vect_full[dim_cord_vect][nucl_num] | out | Full list of coefficients |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes #+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, & integer function qmckl_compute_cord_vect_full_f(context, nucl_num, dim_cord_vect, type_nucl_num, &
type_nucl_vector, cord_vector, cord_vect_full) & type_nucl_vector, cord_vector, cord_vect_full) &
result(info) result(info)
use qmckl use qmckl
implicit none implicit none
integer(qmckl_context), intent(in) :: context integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: nucl_num
integer*8 , intent(in) :: cord_num
integer*8 , intent(in) :: dim_cord_vect integer*8 , intent(in) :: dim_cord_vect
integer*8 , intent(in) :: type_nucl_num integer*8 , intent(in) :: type_nucl_num
integer*8 , intent(in) :: type_nucl_vector(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(in) :: cord_vector(type_nucl_num, dim_cord_vect)
double precision , intent(out) :: cord_vect_full(nucl_num,dim_cord_vect) double precision , intent(out) :: cord_vect_full(nucl_num,dim_cord_vect)
double precision :: x double precision :: x
integer*8 :: i, a, k, l, nw integer*8 :: i, a, k, l, nw
@ -4462,11 +4527,6 @@ integer function qmckl_compute_cord_vect_full_f(context, nucl_num, cord_num, dim
return return
endif endif
if (cord_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif
if (type_nucl_num <= 0) then if (type_nucl_num <= 0) then
info = QMCKL_INVALID_ARG_4 info = QMCKL_INVALID_ARG_4
return return
@ -4479,7 +4539,7 @@ integer function qmckl_compute_cord_vect_full_f(context, nucl_num, cord_num, dim
do a = 1, nucl_num do a = 1, nucl_num
cord_vect_full(1:dim_cord_vect,a) = cord_vector(1:dim_cord_vect,type_nucl_vector(a)) cord_vect_full(a,1:dim_cord_vect) = cord_vector(type_nucl_vector(a),1:dim_cord_vect)
end do end do
end function qmckl_compute_cord_vect_full_f end function qmckl_compute_cord_vect_full_f
@ -4492,7 +4552,6 @@ end function qmckl_compute_cord_vect_full_f
qmckl_exit_code qmckl_compute_cord_vect_full ( qmckl_exit_code qmckl_compute_cord_vect_full (
const qmckl_context context, const qmckl_context context,
const int64_t nucl_num, const int64_t nucl_num,
const int64_t cord_num,
const int64_t dim_cord_vect, const int64_t dim_cord_vect,
const int64_t type_nucl_num, const int64_t type_nucl_num,
const int64_t* type_nucl_vector, const int64_t* type_nucl_vector,
@ -4506,14 +4565,7 @@ end function qmckl_compute_cord_vect_full_f
#+RESULTS: #+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none #+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_cord_vect_full & integer(c_int32_t) function qmckl_compute_cord_vect_full &
(context, & (context, nucl_num, dim_cord_vect, type_nucl_num, type_nucl_vector, cord_vector, cord_vect_full) &
nucl_num, &
cord_num, &
dim_cord_vect, &
type_nucl_num, &
type_nucl_vector, &
cord_vector, &
cord_vect_full) &
bind(C) result(info) bind(C) result(info)
use, intrinsic :: iso_c_binding use, intrinsic :: iso_c_binding
@ -4521,23 +4573,15 @@ end function qmckl_compute_cord_vect_full_f
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: nucl_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 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) , value :: type_nucl_num
integer (c_int64_t) , intent(in) :: type_nucl_vector(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(in) :: cord_vector(type_nucl_num,dim_cord_vect)
real (c_double ) , intent(out) :: cord_vect_full(nucl_num,dim_cord_vect) real (c_double ) , intent(out) :: cord_vect_full(nucl_num,dim_cord_vect)
integer(c_int32_t), external :: qmckl_compute_cord_vect_full_f integer(c_int32_t), external :: qmckl_compute_cord_vect_full_f
info = qmckl_compute_cord_vect_full_f & info = qmckl_compute_cord_vect_full_f &
(context, & (context, nucl_num, dim_cord_vect, type_nucl_num, type_nucl_vector, cord_vector, cord_vect_full)
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 function qmckl_compute_cord_vect_full
#+end_src #+end_src
@ -4644,9 +4688,9 @@ end function qmckl_compute_lkpm_combined_index_f
end function qmckl_compute_lkpm_combined_index end function qmckl_compute_lkpm_combined_index
#+end_src #+end_src
*** Test *** Test
#+name: helper_funcs
#+begin_src python :results output :exports none :noweb yes #+begin_src python :results output :exports none :noweb yes
import numpy as np import numpy as np
@ -4673,21 +4717,46 @@ for l in range(2,cord_num+1):
for i in range(elec_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] 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]) elec_dist = np.zeros(shape=(elec_num, elec_num),dtype=float)
print(" een_rescaled_n[0, 3, 1] = ",een_rescaled_n[0, 3, 1]) for i in range(elec_num):
print(" een_rescaled_n[0, 4, 1] = ",een_rescaled_n[0, 4, 1]) for j in range(elec_num):
print(" een_rescaled_n[1, 3, 2] = ",een_rescaled_n[1, 3, 2]) elec_dist[i, j] = np.linalg.norm(elec_coord[i] - elec_coord[j])
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]) 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
for l in range(0,cord_num+1):
for j in range(0, elec_num):
een_rescaled_e[j,j,l] = 0.0
lkpm_of_cindex = np.array(lkpm_combined_index).T
#+end_src #+end_src
#+RESULTS: #+RESULTS: helper_funcs
: 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) #+begin_src c :tangle (eval c_test)
//assert(qmckl_electron_provided(context)); //assert(qmckl_electron_provided(context));
@ -4834,10 +4903,10 @@ integer function qmckl_compute_factor_een_f(context, walk_num, elec_num, nucl_nu
implicit none implicit none
integer(qmckl_context), intent(in) :: context 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) :: walk_num, elec_num, cord_num, nucl_num, dim_cord_vect
integer*8 , intent(in) :: lkpm_combined_index(4,dim_cord_vect) integer*8 , intent(in) :: lkpm_combined_index(dim_cord_vect,4)
double precision , intent(in) :: cord_vect_full(dim_cord_vect, nucl_num) double precision , intent(in) :: cord_vect_full(nucl_num, dim_cord_vect)
double precision , intent(in) :: een_rescaled_e(walk_num, elec_num, elec_num, 0:cord_num) double precision , intent(in) :: een_rescaled_e(0:cord_num, elec_num, elec_num, walk_num)
double precision , intent(in) :: een_rescaled_n(walk_num, elec_num, nucl_num, 0:cord_num) double precision , intent(in) :: een_rescaled_n(0:cord_num, nucl_num, elec_num, walk_num)
double precision , intent(out) :: factor_een(walk_num) double precision , intent(out) :: factor_een(walk_num)
integer*8 :: i, a, j, l, k, p, m, n, nw integer*8 :: i, a, j, l, k, p, m, n, nw
@ -4874,23 +4943,27 @@ integer function qmckl_compute_factor_een_f(context, walk_num, elec_num, nucl_nu
do nw =1, walk_num do nw =1, walk_num
do n = 1, dim_cord_vect do n = 1, dim_cord_vect
l = lkpm_combined_index(1, n) l = lkpm_combined_index(n, 1)
k = lkpm_combined_index(2, n) k = lkpm_combined_index(n, 2)
p = lkpm_combined_index(3, n) p = lkpm_combined_index(n, 3)
m = lkpm_combined_index(4, n) m = lkpm_combined_index(n, 4)
do a = 1, nucl_num do a = 1, nucl_num
accu2 = 0.0d0 accu2 = 0.0d0
cn = cord_vect_full(n, a) cn = cord_vect_full(a, n)
do j = 1, elec_num do j = 1, elec_num
accu = 0.0d0 accu = 0.0d0
do i = 1, elec_num do i = 1, elec_num
accu = accu + een_rescaled_e(nw, i, j, k) * & accu = accu + een_rescaled_e(k,i,j,nw) * &
een_rescaled_n(nw, i, a, m) een_rescaled_n(m,a,i,nw)
!if(nw .eq. 1) then
! print *,l,k,p,m,j,i,een_rescaled_e(k,i,j,nw), een_rescaled_n(m,a,i,nw), accu
!endif
end do end do
accu2 = accu2 + accu * een_rescaled_n(nw, j, a, m + l) accu2 = accu2 + accu * een_rescaled_n(m + l,a,j,nw)
!print *, l,m,nw,accu, accu2, een_rescaled_n(m + l, a, j, nw), cn, factor_een(nw)
end do end do
factor_een(nw) = factor_een(nw) + accu2 + cn factor_een(nw) = factor_een(nw) + accu2 * cn
end do end do
end do end do
end do end do
@ -4916,7 +4989,6 @@ end function qmckl_compute_factor_een_f
double* const factor_een ); double* const factor_een );
#+end_src #+end_src
#+CALL: generate_c_interface(table=qmckl_factor_een_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+CALL: generate_c_interface(table=qmckl_factor_een_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+RESULTS: #+RESULTS:
@ -4973,100 +5045,44 @@ import numpy as np
<<jastrow_data>> <<jastrow_data>>
<<helper_funcs>>
kappa = 1.0 kappa = 1.0
elec_coord = np.array(elec_coord)[0] factor_een = 0.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 n in range(0, dim_cord_vect):
for a in range(nucl_num): l = lkpm_of_cindex[0,n]
for i in range(elec_num): k = lkpm_of_cindex[1,n]
rij_inv = 1.0 / elnuc_dist[i, a] p = lkpm_of_cindex[2,n]
for ii in range(3): m = lkpm_of_cindex[3,n]
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])
for a in range(0, nucl_num):
accu2 = 0.0
cn = cord_vector_full[a][n]
for j in range(0, elec_num):
accu = 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 + accu * een_rescaled_n[a,j,m+l]
factor_een = factor_een + accu2 * cn
print("factor_een:",factor_een)
#+end_src #+end_src
#+RESULTS: #+RESULTS:
: factor_en_deriv_e[0][0]: 0.11609919541763383 : factor_een: -0.37407972141304213
: 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) #+begin_src c :tangle (eval c_test)
/* Check if Jastrow is properly initialized */ /* Check if Jastrow is properly initialized */
//assert(qmckl_jastrow_provided(context)); assert(qmckl_jastrow_provided(context));
//
double factor_een[walk_num];
rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0]));
assert(fabs(factor_een[0] + 0.37407972141304213) < 1e-12);
#+end_src #+end_src
** Electron-electron-nucleus Jastrow \(f_{een}\) derivative ** Electron-electron-nucleus Jastrow \(f_{een}\) derivative
@ -5221,12 +5237,12 @@ integer function qmckl_compute_factor_een_deriv_e_f(context, walk_num, elec_num,
implicit none implicit none
integer(qmckl_context), intent(in) :: context 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) :: walk_num, elec_num, cord_num, nucl_num, dim_cord_vect
integer*8 , intent(in) :: lkpm_combined_index(4,dim_cord_vect) integer*8 , intent(in) :: lkpm_combined_index(dim_cord_vect, 4)
double precision , intent(in) :: cord_vect_full(dim_cord_vect, nucl_num) double precision , intent(in) :: cord_vect_full(nucl_num, dim_cord_vect)
double precision , intent(in) :: een_rescaled_e(walk_num, elec_num, elec_num, 0:cord_num) double precision , intent(in) :: een_rescaled_e(0:cord_num, elec_num, elec_num, walk_num)
double precision , intent(in) :: een_rescaled_n(walk_num, elec_num, nucl_num, 0:cord_num) double precision , intent(in) :: een_rescaled_n(0:cord_num, nucl_num, elec_num, walk_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_e_deriv_e(0:cord_num, elec_num, 4, elec_num, walk_num)
double precision , intent(in) :: een_rescaled_n_deriv_e(walk_num, elec_num, 4, nucl_num, 0:cord_num) double precision , intent(in) :: een_rescaled_n_deriv_e(0:cord_num, nucl_num, 4, elec_num, walk_num)
double precision , intent(out) :: factor_een_deriv_e(elec_num, 4, walk_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 integer*8 :: i, a, j, l, k, p, m, n, nw
@ -5264,41 +5280,41 @@ integer function qmckl_compute_factor_een_deriv_e_f(context, walk_num, elec_num,
do nw =1, walk_num do nw =1, walk_num
do n = 1, dim_cord_vect do n = 1, dim_cord_vect
l = lkpm_combined_index(1, n) l = lkpm_combined_index(n, 1)
k = lkpm_combined_index(2, n) k = lkpm_combined_index(n, 2)
p = lkpm_combined_index(3, n) p = lkpm_combined_index(n, 3)
m = lkpm_combined_index(4, n) m = lkpm_combined_index(n, 4)
do a = 1, nucl_num do a = 1, nucl_num
cn = cord_vect_full(n, a) cn = cord_vect_full(a, n)
do j = 1, elec_num do j = 1, elec_num
accu = 0.0d0 accu = 0.0d0
accu2 = 0.0d0 accu2 = 0.0d0
daccu = 0.0d0 daccu = 0.0d0
daccu2 = 0.0d0 daccu2 = 0.0d0
do i = 1, elec_num do i = 1, elec_num
accu = accu + een_rescaled_e(nw, i, j, k) * & accu = accu + een_rescaled_e(k, i, j, nw) * &
een_rescaled_n(nw, i, a, m) een_rescaled_n(m, a, i, nw)
accu2 = accu2 + een_rescaled_e(nw, i, j, k) * & accu2 = accu2 + een_rescaled_e(k, i, j, nw) * &
een_rescaled_n(nw, i, a, m + l) een_rescaled_n(m + l, a, i, nw)
daccu(1:4) = daccu(1:4) + een_rescaled_e_deriv_e(nw, j, 1:4, i, k) * & daccu(1:4) = daccu(1:4) + een_rescaled_e_deriv_e(k, j, 1:4, i, nw) * &
een_rescaled_n(nw, i, a, m) een_rescaled_n(m, a, i, nw)
daccu2(1:4) = daccu2(1:4) + een_rescaled_e_deriv_e(nw, j, 1:4, i, k) * & daccu2(1:4) = daccu2(1:4) + een_rescaled_e_deriv_e(k, j, 1:4, i, nw) * &
een_rescaled_n(nw, i, a, m + l) een_rescaled_n(m + l, a, i, nw)
end do end do
factor_een_deriv_e(j, 1:4, nw) = factor_een_deriv_e(j, 1:4, nw) + & 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) & (accu * een_rescaled_n_deriv_e(m + l, a, 1:4, j, nw) &
+ daccu(1:4) * een_rescaled_n(nw, j, a, m + l) & + daccu(1:4) * een_rescaled_n(m + l, a, j, nw) &
+ daccu2(1:4) * een_rescaled_n(nw, j, a, m) & + daccu2(1:4) * een_rescaled_n(m, a, j, nw) &
+ accu2 * een_rescaled_n_deriv_e(nw, j, 1:4, a, m)) * cn + accu2 * een_rescaled_n_deriv_e(m, a, 1:4, j, nw)) * cn
factor_een_deriv_e(j, 4, nw) = factor_een_deriv_e(j, 4, nw) + 2.0d0 * ( & 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 (1) * een_rescaled_n_deriv_e(m + l, a, 1, j, nw) + &
daccu (2) * een_rescaled_n_deriv_e(nw, j, 2, a, m + l) + & daccu (2) * een_rescaled_n_deriv_e(m + l, a, 2, j, nw) + &
daccu (3) * een_rescaled_n_deriv_e(nw, j, 3, a, m + l) + & daccu (3) * een_rescaled_n_deriv_e(m + l, a, 3, j, nw) + &
daccu2(1) * een_rescaled_n_deriv_e(nw, j, 1, a, m ) + & daccu2(1) * een_rescaled_n_deriv_e(m, a, 1, j, nw ) + &
daccu2(2) * een_rescaled_n_deriv_e(nw, j, 2, a, m ) + & daccu2(2) * een_rescaled_n_deriv_e(m, a, 2, j, nw ) + &
daccu2(3) * een_rescaled_n_deriv_e(nw, j, 3, a, m ) ) * cn daccu2(3) * een_rescaled_n_deriv_e(m, a, 3, j, nw ) ) * cn
end do end do
end do end do
@ -5391,98 +5407,60 @@ import numpy as np
<<jastrow_data>> <<jastrow_data>>
<<een_e_deriv_e>>
<<helper_funcs>>
kappa = 1.0 kappa = 1.0
elec_coord = np.array(elec_coord)[0] factor_een = 0.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) daccu = np.zeros(4, dtype=float)
for a in range(nucl_num): daccu2 = np.zeros(4, dtype=float)
for i in range(elec_num): een_rescaled_e_deriv_e_t = een_rescaled_e_deriv_e.T
rij_inv = 1.0 / elnuc_dist[i, a] print(een_rescaled_e_deriv_e_t.shape)
for ii in range(3): for n in range(0, dim_cord_vect):
elnuc_dist_deriv_e[ii, i, a] = (elec_coord[i][ii] - nucl_coord[ii][a]) * rij_inv l = lkpm_of_cindex[0,n]
elnuc_dist_deriv_e[3, i, a] = 2.0 * rij_inv k = lkpm_of_cindex[1,n]
p = lkpm_of_cindex[2,n]
m = lkpm_of_cindex[3,n]
en_distance_rescaled_deriv_e = np.zeros(shape=(4,elec_num,nucl_num),dtype=float) for a in range(0, nucl_num):
for a in range(nucl_num): cn = cord_vector_full[a][n]
for i in range(elec_num): for j in range(0, elec_num):
f = 1.0 - kappa * en_distance_rescaled[i][a] accu = 0.0
for ii in range(4): accu2 = 0.0
en_distance_rescaled_deriv_e[ii][i][a] = elnuc_dist_deriv_e[ii][i][a] daccu = 0.0
en_distance_rescaled_deriv_e[3][i][a] = en_distance_rescaled_deriv_e[3][i][a] + \ daccu2 = 0.0
(-kappa * en_distance_rescaled_deriv_e[0][i][a] * en_distance_rescaled_deriv_e[0][i][a]) + \ for i in range(0, elec_num):
(-kappa * en_distance_rescaled_deriv_e[1][i][a] * en_distance_rescaled_deriv_e[1][i][a]) + \ accu = accu + een_rescaled_e[i,j,k] * \
(-kappa * en_distance_rescaled_deriv_e[2][i][a] * en_distance_rescaled_deriv_e[2][i][a]) een_rescaled_n[a,i,m]
for ii in range(4): accu2 = accu2 + een_rescaled_e[i,j,k] * \
en_distance_rescaled_deriv_e[ii][i][a] = en_distance_rescaled_deriv_e[ii][i][a] * f een_rescaled_n[a,i,m+l]
# daccu[0:4] = daccu[0:4] + een_rescaled_e_deriv_e_t[k,j,0:4,i,k] * \
third = 1.0 / 3.0 # een_rescaled_n[a,i,m]
factor_en_deriv_e = np.zeros(shape=(4,elec_num),dtype=float) # daccu[0:4] = daccu[0:4] + een_rescaled_e_deriv_e_t[k,j,0:4,i,k] * \
dx = np.zeros(shape=(4),dtype=float) # een_rescaled_n[a,i,m]
pow_ser_g = np.zeros(shape=(3),dtype=float) accu2 = accu2 + accu * een_rescaled_n[a,j,m+l]
for a in range(nucl_num): # factor_een = factor_een + accu2 * cn
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])
print("factor_een:",factor_een)
#+end_src #+end_src
#+RESULTS: #+RESULTS:
: factor_en_deriv_e[0][0]: 0.11609919541763383 : (6, 10, 4, 10)
: factor_en_deriv_e[1][0]: -0.23301394780804574 : factor_een: 0.0
: factor_en_deriv_e[2][0]: 0.17548337641865783
: factor_en_deriv_e[3][0]: -0.9667363412285741
#+begin_src c :tangle (eval c_test) #+begin_src c :tangle (eval c_test)
///* Check if Jastrow is properly initialized */ /* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_provided(context));
double factor_een_deriv_e[walk_num][elec_num];
rc = qmckl_get_jastrow_factor_een_deriv_e(context, &(factor_een_deriv_e[0][0]));
assert(fabs(factor_een_deriv_e[0][0] + 0.0005481671107226865) < 1e-12);
#+end_src #+end_src
* End of files :noexport: * End of files :noexport:

1306
org/qmckl_sherman_morrison_woodbury.org Normal file
View File

File diff suppressed because one or more lines are too long

4
org/qmckl_tests.org
View File

@ -1191,7 +1191,7 @@ double n2_elec_coord[n2_walk_num][n2_elec_num][3] = { {
#define n2_type_nucl_num ((int64_t) 1) #define n2_type_nucl_num ((int64_t) 1)
#define n2_aord_num ((int64_t) 5) #define n2_aord_num ((int64_t) 5)
#define n2_bord_num ((int64_t) 5) #define n2_bord_num ((int64_t) 5)
#define n2_cord_num ((int64_t) 23) #define n2_cord_num ((int64_t) 5)
#define n2_dim_cord_vec ((int64_t) 23) #define n2_dim_cord_vec ((int64_t) 23)
int64_t n2_type_nucl_vector[n2_nucl_num] = { int64_t n2_type_nucl_vector[n2_nucl_num] = {
@ -1214,7 +1214,7 @@ double n2_bord_vector[n2_bord_num + 1] = {
0.0073096 , 0.0073096 ,
0.002866 }; 0.002866 };
double n2_cord_vector[n2_cord_num][n2_type_nucl_num] = { double n2_cord_vector[n2_dim_cord_vec][n2_type_nucl_num] = {
{ 5.717020e-01}, { 5.717020e-01},
{-5.142530e-01}, {-5.142530e-01},
{-5.130430e-01}, {-5.130430e-01},

1
org/table_of_contents
View File

@ -7,6 +7,7 @@ qmckl_nucleus.org
qmckl_electron.org qmckl_electron.org
qmckl_ao.org qmckl_ao.org
qmckl_jastrow.org qmckl_jastrow.org
qmckl_sherman_morrison_woodbury.org
qmckl_distance.org qmckl_distance.org
qmckl_utils.org qmckl_utils.org
qmckl_tests.org qmckl_tests.org

29
tools/lib.org
View File

@ -11,11 +11,10 @@
** Function to get the value of a property. ** Function to get the value of a property.
#+NAME: get_value #+NAME: get_value
#+begin_src elisp :var key="Type" #+begin_src elisp :var key="Type"
(setq x (org-property-values key)) (org-with-point-at org-babel-current-src-block-location
(pop x) (org-entry-get nil key t))
#+end_src #+end_src
#+RESULTS: get_value
** Table of function arguments ** Table of function arguments
@ -42,6 +41,8 @@ f_of_c_d = { '' : ''
, 'qmckl_exit_code' : 'integer (c_int32_t)' , 'qmckl_exit_code' : 'integer (c_int32_t)'
, 'int32_t' : 'integer (c_int32_t)' , 'int32_t' : 'integer (c_int32_t)'
, 'int64_t' : 'integer (c_int64_t)' , 'int64_t' : 'integer (c_int64_t)'
, 'uint32_t' : 'integer (c_int32_t)'
, 'uint64_t' : 'integer (c_int64_t)'
, 'float' : 'real (c_float )' , 'float' : 'real (c_float )'
, 'double' : 'real (c_double )' , 'double' : 'real (c_double )'
, 'char' : 'character' , 'char' : 'character'
@ -55,6 +56,8 @@ ctypeid_d = { '' : ''
, 'qmckl_exit_code' : 'integer(c_int32_t)' , 'qmckl_exit_code' : 'integer(c_int32_t)'
, 'integer' : 'integer(c_int32_t)' , 'integer' : 'integer(c_int32_t)'
, 'integer*8' : 'integer(c_int64_t)' , 'integer*8' : 'integer(c_int64_t)'
, 'integer' : 'integer(c_uint32_t)'
, 'integer*8' : 'integer(c_uint64_t)'
, 'real' : 'real(c_float)' , 'real' : 'real(c_float)'
, 'real*8' : 'real(c_double)' , 'real*8' : 'real(c_double)'
, 'character' : 'character(c_char)' , 'character' : 'character(c_char)'
@ -184,23 +187,6 @@ results='\n'.join(results)
return results return results
#+END_SRC #+END_SRC
#+RESULTS: generate_c_interface
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function [] &
() &
bind(C) result(info)
use, intrinsic :: iso_c_binding
implicit none
integer(c_int32_t), external :: []_f
info = []_f &
()
end function []
#+end_src
*** Generates a Fortran interface to the C function *** Generates a Fortran interface to the C function
#+NAME: generate_f_interface #+NAME: generate_f_interface
@ -252,8 +238,5 @@ results='\n'.join(results)
return results return results
#+END_SRC #+END_SRC
#+RESULTS: generate_f_interface
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
#+end_src