mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-12-22 12:23:56 +01:00
Fixed elec_coord bugs
This commit is contained in:
parent
e009322467
commit
63d831ae04
@ -68,10 +68,11 @@
|
||||
** Reformulation of the three-body part
|
||||
|
||||
To accelerate the computation of the three-body part, the Jastrow
|
||||
factor is re-expressed as follows:
|
||||
factor is re-expressed as follows, with $m=(p-k)/2 -l/2$:
|
||||
|
||||
|
||||
\begin{eqnarray*}
|
||||
& & \sum_{\alpha=1}^{N_{\text{nucl}}}
|
||||
J_{kpl} & = & \sum_{\alpha=1}^{N_{\text{nucl}}}
|
||||
\sum_{i=1}^{N_{\text{elec}}}
|
||||
\sum_{j=1}^{i-1}
|
||||
c_{lkp\alpha} \left[ g_\text{e}({r}_{ij}) \right]^k
|
||||
@ -93,8 +94,35 @@
|
||||
\left[ \left[ g_\alpha({R}_{i\alpha}) \right]^{l+m} \left[ g_\alpha({R}_{j\alpha}) \right]^{m} +
|
||||
\left[ g_\alpha({R}_{i\alpha}) \right]^{l} \left[
|
||||
g_\alpha({R}_{j\alpha}) \right]^{l+m} \right] \\
|
||||
& = &
|
||||
\sum_{\alpha=1}^{N_{\text{nucl}}} c_{lkp\alpha}
|
||||
\sum_{i=1}^{N_{\text{elec}}}
|
||||
\sum_{j=1}^{N_{\text{elec}}}
|
||||
\left[ g_\alpha({R}_{i\alpha}) \right]^{l+m}
|
||||
\left[ g_\text{e}({r}_{ij}) \right]^k
|
||||
\left[ g_\alpha({R}_{j\alpha}) \right]^{m} \\
|
||||
& = &
|
||||
\sum_{\alpha=1}^{N_{\text{nucl}}} c_{lkp\alpha}
|
||||
\sum_{i=1}^{N_{\text{elec}}}
|
||||
\left[ g_\alpha({R}_{i\alpha}) \right]^{l+m}
|
||||
P_{i\alpha}^{km}, \text{ with }
|
||||
P_{i\alpha}^{km} =
|
||||
\sum_{j=1}^{N_{\text{elec}}}
|
||||
\left[ g_\text{e}({r}_{ij}) \right]^k
|
||||
\left[ g_\alpha({R}_{j\alpha}) \right]^{m}. \\
|
||||
J & = &
|
||||
\sum_{p=2}^{N_{\text{ord}}}
|
||||
\sum_{k=0}^{p-1}
|
||||
\sum_{l=0}^{p-k-2\delta_{k,0}}
|
||||
\sum_{\alpha=1}^{N_{\text{nucl}}} c_{lkp\alpha}
|
||||
\sum_{i=1}^{N_{\text{elec}}}
|
||||
\left[ g_\alpha({R}_{i\alpha}) \right]^{(p-k+l)/2}
|
||||
P_{i\alpha}^{k, (p-k-l)/2}
|
||||
\end{eqnarray*}
|
||||
|
||||
The computation of $P$ scales as $\mathcal{O}(N_\text{elec}^2 N_\text{nucl}n^2)$, and
|
||||
the computation of $J$ scales as $\mathcal{O}(N_\text{elec}N_\text{nucl}n^2)$.
|
||||
|
||||
* Headers :noexport:
|
||||
#+begin_src elisp :noexport :results none
|
||||
(org-babel-lob-ingest "../tools/lib.org")
|
||||
@ -126,6 +154,11 @@
|
||||
#include "qmckl_context_private_type.h"
|
||||
#include "qmckl_jastrow_champ_private_func.h"
|
||||
|
||||
#ifdef __GNUC__
|
||||
#pragma GCC push_options
|
||||
#pragma GCC optimize ("O0")
|
||||
#endif
|
||||
|
||||
int main() {
|
||||
qmckl_context context;
|
||||
context = qmckl_context_create();
|
||||
@ -1798,7 +1831,6 @@ for (int64_t i=0 ; i<nucl_num ; ++i) {
|
||||
assert( nucl_charge[i] == nucl_charge2[i] );
|
||||
}
|
||||
assert(qmckl_nucleus_provided(context));
|
||||
|
||||
#+end_src
|
||||
|
||||
* Computation
|
||||
@ -2081,6 +2113,7 @@ qmckl_exit_code qmckl_compute_jastrow_champ_asymp_jasb (const qmckl_context cont
|
||||
(context, bord_num, b_vector, rescale_factor_ee, spin_independent, asymp_jasb);
|
||||
}
|
||||
#+end_src
|
||||
|
||||
**** Test :noexport:
|
||||
#+name: asymp_jasb
|
||||
#+begin_src python :results output :exports none :noweb yes
|
||||
@ -3486,6 +3519,7 @@ qmckl_compute_jastrow_champ_factor_ee (const qmckl_context context,
|
||||
ee_distance_rescaled, asymp_jasb, spin_independent, factor_ee);
|
||||
}
|
||||
#+end_src
|
||||
|
||||
**** Test :noexport:
|
||||
#+begin_src python :results output :exports none :noweb yes
|
||||
import numpy as np
|
||||
@ -3515,37 +3549,14 @@ for i in range(0,elec_num):
|
||||
/ (1.0 + b_vector[1] * ee_distance_rescaled[i][j]) \
|
||||
- asymp_jasb[ipar] + pow_ser
|
||||
print("factor_ee :",factor_ee)
|
||||
print("ee_distance_rescaled :",ee_distance_rescaled)
|
||||
|
||||
#+end_src
|
||||
|
||||
#+RESULTS:
|
||||
#+begin_example
|
||||
asym_one : 0.6634291325000664
|
||||
asymp_jasb[0] : 0.7115733522582638
|
||||
asymp_jasb[1] : 1.043287918508297
|
||||
factor_ee : -16.83886184243964
|
||||
ee_distance_rescaled : [[0. 0.63475074 1.29816415 1.23147027 1.51933127 0.54402406
|
||||
0.51452479 0.96538731 1.25878564 1.3722995 ]
|
||||
[0.63475074 0. 1.35148664 1.13524156 1.48940503 0.4582292
|
||||
0.62758076 1.06560856 1.179133 1.30763703]
|
||||
[1.29816415 1.35148664 0. 1.50021375 1.59200788 1.23488312
|
||||
1.20844259 1.0355537 1.52064535 1.53049239]
|
||||
[1.23147027 1.13524156 1.50021375 0. 1.12016142 1.19158954
|
||||
1.29762585 1.24824277 0.25292267 0.58609336]
|
||||
[1.51933127 1.48940503 1.59200788 1.12016142 0. 1.50217017
|
||||
1.54012828 1.48753895 1.10441805 0.84504381]
|
||||
[0.54402406 0.4582292 1.23488312 1.19158954 1.50217017 0.
|
||||
0.39417354 0.87009603 1.23838502 1.33419121]
|
||||
[0.51452479 0.62758076 1.20844259 1.29762585 1.54012828 0.39417354
|
||||
0. 0.95118809 1.33068934 1.41097406]
|
||||
[0.96538731 1.06560856 1.0355537 1.24824277 1.48753895 0.87009603
|
||||
0.95118809 0. 1.29422213 1.33222493]
|
||||
[1.25878564 1.179133 1.52064535 0.25292267 1.10441805 1.23838502
|
||||
1.33068934 1.29422213 0. 0.62196802]
|
||||
[1.3722995 1.30763703 1.53049239 0.58609336 0.84504381 1.33419121
|
||||
1.41097406 1.33222493 0.62196802 0. ]]
|
||||
#+end_example
|
||||
: 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 */
|
||||
@ -4218,10 +4229,10 @@ print("factor_ee_gl[3][0]:",factor_ee_gl[3][0])
|
||||
: 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]:
|
||||
: factor_ee_gl[0][0]: -0.39319353942687446
|
||||
: factor_ee_gl[1][0]: 1.0535615450668214
|
||||
: factor_ee_gl[2][0]: -0.39098406960784515
|
||||
: factor_ee_gl[3][0]: 2.8650469630854483
|
||||
|
||||
#+begin_src c :tangle (eval c_test)
|
||||
/* Check if Jastrow is properly initialized */
|
||||
@ -5555,7 +5566,7 @@ assert(qmckl_electron_provided(context));
|
||||
coeffecients and the electron-nucleus rescaled distances ~en_distance_rescaled~.
|
||||
|
||||
\[
|
||||
f_{\alpha}(R_{i\alpha}) = - \sum_{i,j<i} \left[ \frac{ A_0 C_{ij}}{1 - A_1 C_{ij}} + \sum^{N^\alpha_{\text{ord}}}_{k}A_k C_{ij}^k \right]
|
||||
f_{\alpha}(R_{i\alpha}) = - \sum_{i,j<i} \left[ \frac{ A_0 C_{ij}}{1 + A_1 C_{ij}} + \sum_{k=2}^{N^\alpha_{\text{ord}}}A_k C_{ij}^k \right]
|
||||
\]
|
||||
|
||||
|
||||
@ -6023,7 +6034,21 @@ assert(fabs(22.781375792083587 - factor_en[0]) < 1.e-12);
|
||||
with respect to the electron coordinates using the ~en_distance_rescaled~ and
|
||||
~en_distance_rescaled_gl~ which are already calculated previously.
|
||||
|
||||
TODO: write equations.
|
||||
The derivative is calculated in the function ~qmckl_compute_jastrow_champ_factor_en_gl~.
|
||||
The formula is given by:
|
||||
\[
|
||||
\nabla_i f_{\alpha}(R_{i\alpha}) = \sum_{j=1}^{N_\text{elec}} \left[
|
||||
\frac{ A_0\, \nabla_i C_{ij} }{(1 + A_1 C_{ij})^2} +
|
||||
\sum_{k=2}^{N^\alpha_{\text{ord}}} A_k\, k\, C_{ij}^{k-1}\,\nabla_i C_{ij} \right]
|
||||
\]
|
||||
|
||||
\[
|
||||
\Delta_i f_{\alpha}(R_{i\alpha}) = \sum_{j=1}^{N_\text{elec}} \left[
|
||||
\frac{ A_0\, \Delta_i C_{ij} }{(1 + A_1 C_{ij})^2} -
|
||||
\frac{ 2 A_0\, A_1 (\nabla_i C_{ij})^2}{(1 + A_1 C_{ij})^3} +
|
||||
\sum_{k=2}^{N^\alpha_{\text{ord}}} A_k\, k\, C_{ij}^{k-1} C_{ij}^{k-2}
|
||||
\left[ \Delta_i C_{ij} + (k-1)(\nabla_i C_{ij})^2 \right] \right]
|
||||
\]
|
||||
|
||||
**** Get
|
||||
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
||||
@ -6186,6 +6211,7 @@ 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 |
|
||||
@ -6199,6 +6225,7 @@ qmckl_exit_code qmckl_provide_jastrow_champ_factor_en_gl(qmckl_context context)
|
||||
| ~en_distance_rescaled~ | ~double[walk_num][nucl_num][elec_num]~ | in | Electron-nucleus distances |
|
||||
| ~en_distance_rescaled_gl~ | ~double[walk_num][nucl_num][elec_num][4]~ | in | Electron-nucleus distance derivatives |
|
||||
| ~factor_en_gl~ | ~double[walk_num][4][elec_num]~ | out | Electron-nucleus jastrow |
|
||||
|---------------------------+-------------------------------------------+--------+---------------------------------------|
|
||||
|
||||
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
||||
function qmckl_compute_jastrow_champ_factor_en_gl_doc( &
|
||||
@ -6714,13 +6741,34 @@ assert(qmckl_jastrow_champ_provided(context));
|
||||
#+end_src
|
||||
|
||||
** Electron-electron-nucleus component
|
||||
|
||||
#+name: en_dist
|
||||
#+begin_src python :results output :exports none :noweb yes
|
||||
elec_coord = np.array(elec_coord)[0]
|
||||
nucl_coord = np.array(nucl_coord)
|
||||
elnuc_dist = np.zeros(shape=(nucl_num, elec_num),dtype=float)
|
||||
for i in range(elec_num):
|
||||
for a in range(nucl_num):
|
||||
elnuc_dist[a,i] = np.linalg.norm(elec_coord[i] - nucl_coord[:,a])
|
||||
|
||||
kappa_n = 0.6
|
||||
|
||||
een_rescaled_n = np.zeros(shape=(cord_num+1, nucl_num, elec_num), dtype=float)
|
||||
|
||||
for p in range(0, cord_num+1):
|
||||
for a in range(nucl_num):
|
||||
for i in range(elec_num):
|
||||
een_rescaled_n[p,a,i] = np.exp(-kappa_n * p * elnuc_dist[a,i])
|
||||
|
||||
#+end_src
|
||||
|
||||
*** Electron-electron rescaled distances in $J_\text{eeN}$
|
||||
|
||||
~een_rescaled_e~ stores the table of the rescaled distances between all
|
||||
pairs of electrons and raised to the power \(p\) defined by ~cord_num~:
|
||||
|
||||
\[
|
||||
C_{ij,p} = \left[ \exp\left(-\kappa_\text{e}\, r_{ij}\right) \right]^p
|
||||
[g_e(r)_{ij}]^p = \exp\left(-p\,\kappa_\text{e}\, r_{ij}\right)
|
||||
\]
|
||||
|
||||
where \(r_{ij}\) is the matrix of electron-electron distances.
|
||||
@ -6880,7 +6928,6 @@ function qmckl_compute_een_rescaled_e_doc( &
|
||||
real (c_double ) , intent(out) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
|
||||
integer(qmckl_exit_code) :: info
|
||||
|
||||
double precision,dimension(:,:),allocatable :: een_rescaled_e_ij
|
||||
double precision :: x
|
||||
integer*8 :: i, j, k, l, nw
|
||||
|
||||
@ -6906,49 +6953,16 @@ function qmckl_compute_een_rescaled_e_doc( &
|
||||
return
|
||||
endif
|
||||
|
||||
allocate(een_rescaled_e_ij(elec_num * (elec_num - 1) / 2, cord_num + 1))
|
||||
|
||||
! Prepare table of exponentiated distances raised to appropriate power
|
||||
do nw = 1, walk_num
|
||||
een_rescaled_e_ij(:, 1) = 1.0d0
|
||||
|
||||
|
||||
k = 0
|
||||
do l = 0, cord_num
|
||||
do j = 1, elec_num
|
||||
do i = 1, j - 1
|
||||
k = k + 1
|
||||
een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_ee * ee_distance(i, j, nw))
|
||||
do i = 1, elec_num
|
||||
een_rescaled_e(i, j, l, nw) = dexp(-rescale_factor_ee * ee_distance(i, j, nw))**l
|
||||
end do
|
||||
end do
|
||||
|
||||
|
||||
do l = 2, cord_num
|
||||
do k = 1, elec_num * (elec_num - 1)/2
|
||||
een_rescaled_e_ij(k, l+1) = een_rescaled_e_ij(k, l) * een_rescaled_e_ij(k, 2)
|
||||
end do
|
||||
end do
|
||||
|
||||
! prepare the actual een table
|
||||
een_rescaled_e(:, :, 0, nw) = 1.0d0
|
||||
do j = 1, elec_num
|
||||
een_rescaled_e(j, j, 0, nw) = 0.0d0
|
||||
end do
|
||||
|
||||
do l = 1, cord_num
|
||||
k = 0
|
||||
do j = 1, elec_num
|
||||
do i = 1, j - 1
|
||||
k = k + 1
|
||||
x = een_rescaled_e_ij(k, l+1)
|
||||
een_rescaled_e(i, j, l, nw) = x
|
||||
een_rescaled_e(j, i, l, nw) = x
|
||||
end do
|
||||
een_rescaled_e(j, j, l, nw) = 0.0d0
|
||||
end do
|
||||
end do
|
||||
|
||||
end do
|
||||
|
||||
end function qmckl_compute_een_rescaled_e_doc
|
||||
#+end_src
|
||||
|
||||
@ -6974,6 +6988,13 @@ qmckl_exit_code qmckl_compute_een_rescaled_e_hpc (const qmckl_context context,
|
||||
const double* ee_distance,
|
||||
double* const een_rescaled_e ) {
|
||||
|
||||
return qmckl_compute_een_rescaled_e_doc (context,
|
||||
walk_num,
|
||||
elec_num,
|
||||
cord_num,
|
||||
rescale_factor_ee,
|
||||
ee_distance,
|
||||
een_rescaled_e ) ;
|
||||
if (context == QMCKL_NULL_CONTEXT) {
|
||||
return QMCKL_INVALID_CONTEXT;
|
||||
}
|
||||
@ -7153,78 +7174,46 @@ qmckl_exit_code qmckl_compute_een_rescaled_e_hpc (const qmckl_context context,
|
||||
|
||||
**** Test
|
||||
|
||||
#+name: ee_dist
|
||||
#+begin_src python :results output :exports none :noweb yes
|
||||
import numpy as np
|
||||
|
||||
<<jastrow_data>>
|
||||
|
||||
elec_coord = np.array(elec_coord)[0]
|
||||
elec_dist = np.zeros(shape=(elec_num, elec_num),dtype=float)
|
||||
|
||||
for i in range(elec_num):
|
||||
for j in range(elec_num):
|
||||
elec_dist[i, j] = np.linalg.norm(elec_coord[i] - elec_coord[j])
|
||||
elec_dist[i,j] = np.linalg.norm(elec_coord[i] - elec_coord[j])
|
||||
|
||||
kappa = 0.6
|
||||
kappa_e = 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
|
||||
|
||||
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
|
||||
een_rescaled_e = np.zeros(shape=(cord_num+1, elec_num, elec_num), dtype=float)
|
||||
for i in range(elec_num):
|
||||
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
|
||||
|
||||
print(" een_rescaled_e[0, 2, 1] = ",een_rescaled_e[0, 2, 1])
|
||||
print(" een_rescaled_e[0, 3, 1] = ",een_rescaled_e[0, 3, 1])
|
||||
print(" een_rescaled_e[0, 4, 1] = ",een_rescaled_e[0, 4, 1])
|
||||
print(" een_rescaled_e[1, 3, 2] = ",een_rescaled_e[1, 3, 2])
|
||||
print(" een_rescaled_e[1, 4, 2] = ",een_rescaled_e[1, 4, 2])
|
||||
print(" een_rescaled_e[1, 5, 2] = ",een_rescaled_e[1, 5, 2])
|
||||
for p in range(0, cord_num+1):
|
||||
een_rescaled_e[p,i,j] = np.exp(-kappa_e * p * elec_dist[i,j])
|
||||
#+end_src
|
||||
|
||||
#+RESULTS:
|
||||
: 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
|
||||
#+NAME:test_ee
|
||||
#+begin_src python :results output :exports none :noweb yes
|
||||
<<ee_dist>>
|
||||
for p in range(cord_num+1):
|
||||
# for i in range(elec_num):
|
||||
# for j in range(elec_num):
|
||||
for i in range(3):
|
||||
for j in range(3):
|
||||
print(f"assert( fabs(een_rescaled_e[0][{p}][{i}][{j}] - ({een_rescaled_e[p,i,j]})) < 1.e-10 );")
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :tangle (eval c_test)
|
||||
#+begin_src c :tangle (eval c_test) :noweb yes
|
||||
assert(qmckl_electron_provided(context));
|
||||
{
|
||||
|
||||
double een_rescaled_e[walk_num][(cord_num + 1)][elec_num][elec_num];
|
||||
rc = qmckl_get_jastrow_champ_een_distance_rescaled_e(context, &(een_rescaled_e[0][0][0][0]),elec_num*elec_num*(cord_num+1)*walk_num);
|
||||
double een_rescaled_e[walk_num][(cord_num + 1)][elec_num][elec_num];
|
||||
rc = qmckl_get_jastrow_champ_een_distance_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.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);
|
||||
<<test_ee()>>
|
||||
}
|
||||
|
||||
{
|
||||
@ -7236,13 +7225,13 @@ assert(fabs(een_rescaled_e[0][2][1][5]- 0.5257156022077619 ) < 1.e-12);
|
||||
double een_rescaled_e_doc[walk_num][cord_num+1][elec_num][elec_num];
|
||||
memset(&(een_rescaled_e_doc[0][0][0][0]), 0, sizeof(een_rescaled_e_doc));
|
||||
rc = qmckl_compute_een_rescaled_e(context, walk_num, elec_num, cord_num,
|
||||
0.6, &(ee_distance[0]), &(een_rescaled_e_doc[0][0][0][0]));
|
||||
rescale_factor_ee, &(ee_distance[0]), &(een_rescaled_e_doc[0][0][0][0]));
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
|
||||
double een_rescaled_e_hpc[walk_num][cord_num+1][elec_num][elec_num];
|
||||
memset(&(een_rescaled_e_hpc[0][0][0][0]), 0, sizeof(een_rescaled_e_hpc));
|
||||
rc = qmckl_compute_een_rescaled_e_hpc(context, walk_num, elec_num, cord_num,
|
||||
0.6, &(ee_distance[0]), &(een_rescaled_e_hpc[0][0][0][0]));
|
||||
rescale_factor_ee, &(ee_distance[0]), &(een_rescaled_e_hpc[0][0][0][0]));
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
|
||||
for (int64_t i = 0; i < walk_num; i++) {
|
||||
@ -7270,7 +7259,12 @@ assert(fabs(een_rescaled_e[0][2][1][5]- 0.5257156022077619 ) < 1.e-12);
|
||||
|
||||
\[ \frac{\partial}{\partial x} \left[ {g_\text{e}(r)}\right]^p =
|
||||
-\frac{x}{r} \kappa_\text{e}\, p\,\left[ {g_\text{e}(r)}\right]^p \]
|
||||
\[ \Delta \left[ {g_\text{e}(r)}\right]^p = \frac{2}{r} \kappa_\text{e}\, p\,\left[ {g_\text{e}(r)}\right]^p \right] + \left(\frac{\partial}{\partial x}\left[ {g_\text{e}(r)}\right]^p \right)^2 + \left(\frac{\partial}{\partial y}\left[ {g_\text{e}(r)}\right]^p \right)^2 + \left(\frac{\partial}{\partial z}\left[ {g_\text{e}(r)}\right]^p \right)^2 \]
|
||||
|
||||
\[ \Delta \left[ {g_\text{e}(r)}\right]^p = \kappa_\text{e}\,
|
||||
p\,\left[ - \frac{2}{r} + \kappa_\text{e}\, p \right] \left[
|
||||
{g_\text{e}(r)} \right]^p \]
|
||||
|
||||
Derivatives are set to zero at $r_{ii}$ to avoid NaNs.
|
||||
|
||||
**** Get
|
||||
|
||||
@ -7408,6 +7402,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_gl(qmckl_context context)
|
||||
:END:
|
||||
|
||||
#+NAME: qmckl_factor_een_rescaled_e_gl_args
|
||||
|---------------------+-------------------------------------------------------+--------+--------------------------------------|
|
||||
| Variable | Type | In/Out | Description |
|
||||
|---------------------+-------------------------------------------------------+--------+--------------------------------------|
|
||||
| ~context~ | ~qmckl_context~ | in | Global state |
|
||||
@ -7419,6 +7414,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_gl(qmckl_context context)
|
||||
| ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | in | Electron-electron distances |
|
||||
| ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | in | Electron-electron distances |
|
||||
| ~een_rescaled_e_gl~ | ~double[walk_num][0:cord_num][elec_num][4][elec_num]~ | out | Electron-electron rescaled distances |
|
||||
|---------------------+-------------------------------------------------------+--------+--------------------------------------|
|
||||
|
||||
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
||||
function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_doc( &
|
||||
@ -7433,7 +7429,7 @@ function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_doc( &
|
||||
integer(c_int64_t) , intent(in), value :: elec_num
|
||||
integer(c_int64_t) , intent(in), value :: cord_num
|
||||
real(c_double) , intent(in), value :: rescale_factor_ee
|
||||
real(c_double) , intent(in) :: coord_ee(elec_num,3,walk_num)
|
||||
real(c_double) , intent(in) :: coord_ee(elec_num,walk_num,3)
|
||||
real(c_double) , intent(in) :: ee_distance(elec_num,elec_num,walk_num)
|
||||
real(c_double) , intent(in) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
|
||||
real(c_double) , intent(out) :: een_rescaled_e_gl(elec_num,4,elec_num,0:cord_num,walk_num)
|
||||
@ -7445,7 +7441,10 @@ function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_doc( &
|
||||
|
||||
double precision :: rij_inv(elec_num)
|
||||
|
||||
|
||||
allocate(elec_dist_gl(elec_num, 4, elec_num))
|
||||
elec_dist_gl = 0.d0
|
||||
een_rescaled_e_gl = 0.d0
|
||||
|
||||
info = QMCKL_SUCCESS
|
||||
|
||||
@ -7482,7 +7481,7 @@ function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_doc( &
|
||||
enddo
|
||||
do i = 1, elec_num
|
||||
do ii = 1, 3
|
||||
elec_dist_gl(i, ii, j) = (coord_ee(i, ii, nw) - coord_ee(j, ii, nw)) * rij_inv(i)
|
||||
elec_dist_gl(i, ii, j) = (coord_ee(i, nw, ii) - coord_ee(j, nw, ii)) * rij_inv(i)
|
||||
end do
|
||||
elec_dist_gl(i, 4, j) = 2.0d0 * rij_inv(i)
|
||||
end do
|
||||
@ -7492,27 +7491,20 @@ function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_doc( &
|
||||
een_rescaled_e_gl(:,:,:,0,nw) = 0.d0
|
||||
|
||||
do l = 1, cord_num
|
||||
kappa_l = - dble(l) * rescale_factor_ee
|
||||
kappa_l = -dble(l) * rescale_factor_ee
|
||||
do j = 1, elec_num
|
||||
do i = 1, elec_num
|
||||
een_rescaled_e_gl(i, 1, j, l, nw) = kappa_l * elec_dist_gl(i, 1, j)
|
||||
een_rescaled_e_gl(i, 2, j, l, nw) = kappa_l * elec_dist_gl(i, 2, j)
|
||||
een_rescaled_e_gl(i, 3, j, l, nw) = kappa_l * elec_dist_gl(i, 3, j)
|
||||
een_rescaled_e_gl(i, 4, j, l, nw) = kappa_l * elec_dist_gl(i, 4, j)
|
||||
end do
|
||||
|
||||
do i = 1, elec_num
|
||||
een_rescaled_e_gl(i, 4, j, l, nw) = een_rescaled_e_gl(i, 4, j, l, nw) &
|
||||
+ een_rescaled_e_gl(i, 1, j, l, nw) * een_rescaled_e_gl(i, 1, j, l, nw) &
|
||||
+ een_rescaled_e_gl(i, 2, j, l, nw) * een_rescaled_e_gl(i, 2, j, l, nw) &
|
||||
+ een_rescaled_e_gl(i, 3, j, l, nw) * een_rescaled_e_gl(i, 3, j, l, nw)
|
||||
end do
|
||||
|
||||
do i = 1, elec_num
|
||||
een_rescaled_e_gl(i,1,j,l,nw) = een_rescaled_e_gl(i,1,j,l,nw) * een_rescaled_e(i,j,l,nw)
|
||||
een_rescaled_e_gl(i,2,j,l,nw) = een_rescaled_e_gl(i,2,j,l,nw) * een_rescaled_e(i,j,l,nw)
|
||||
een_rescaled_e_gl(i,3,j,l,nw) = een_rescaled_e_gl(i,3,j,l,nw) * een_rescaled_e(i,j,l,nw)
|
||||
een_rescaled_e_gl(i,4,j,l,nw) = een_rescaled_e_gl(i,4,j,l,nw) * een_rescaled_e(i,j,l,nw)
|
||||
if (i /= j) then
|
||||
een_rescaled_e_gl(i, 1, j, l, nw) = kappa_l * elec_dist_gl(i, 1, j) * een_rescaled_e(i,j,l,nw)
|
||||
een_rescaled_e_gl(i, 2, j, l, nw) = kappa_l * elec_dist_gl(i, 2, j) * een_rescaled_e(i,j,l,nw)
|
||||
een_rescaled_e_gl(i, 3, j, l, nw) = kappa_l * elec_dist_gl(i, 3, j) * een_rescaled_e(i,j,l,nw)
|
||||
een_rescaled_e_gl(i, 4, j, l, nw) = kappa_l * (elec_dist_gl(i, 4, j) + kappa_l) * een_rescaled_e(i,j,l,nw)
|
||||
else
|
||||
een_rescaled_e_gl(i, 1, j, l, nw) = 0.d0
|
||||
een_rescaled_e_gl(i, 2, j, l, nw) = 0.d0
|
||||
een_rescaled_e_gl(i, 3, j, l, nw) = 0.d0
|
||||
een_rescaled_e_gl(i, 4, j, l, nw) = 0.d0
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
@ -7599,10 +7591,10 @@ qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_hpc (const qmckl_context co
|
||||
#endif
|
||||
for (int64_t nw = 0; nw < walk_num; ++nw) {
|
||||
|
||||
double rij_inv[elec_num];
|
||||
|
||||
for (int64_t j=0; j<elec_num; ++j) {
|
||||
|
||||
double rij_inv[elec_num];
|
||||
|
||||
#ifdef HAVE_OPENMP
|
||||
#pragma omp simd
|
||||
#endif
|
||||
@ -7619,18 +7611,18 @@ qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_hpc (const qmckl_context co
|
||||
}
|
||||
rij_inv[j] = 0.0;
|
||||
|
||||
const double xj = coord_ee[j + nw * elec_num * 3];
|
||||
const double yj = coord_ee[j + elec_num + nw * elec_num * 3];
|
||||
const double zj = coord_ee[j + 2 * elec_num + nw * elec_num * 3];
|
||||
const double xj = coord_ee[j+elec_num*nw];
|
||||
const double yj = coord_ee[j+elec_num*(nw+walk_num) ];
|
||||
const double zj = coord_ee[j+elec_num*(nw+walk_num*2)];
|
||||
|
||||
#ifdef HAVE_OPENMP
|
||||
#pragma omp simd
|
||||
#endif
|
||||
for (int64_t i = 0; i < elec_num ; ++i) {
|
||||
|
||||
const double xi = coord_ee[i + nw * elec_num * 3];
|
||||
const double yi = coord_ee[i + elec_num + nw * elec_num * 3];
|
||||
const double zi = coord_ee[i + 2 * elec_num + nw * elec_num * 3];
|
||||
const double xi = coord_ee[i+elec_num*nw];
|
||||
const double yi = coord_ee[i+elec_num*(nw+walk_num) ];
|
||||
const double zi = coord_ee[i+elec_num*(nw+walk_num*2)];
|
||||
|
||||
elec_dist_gl0[i + j * elec_num] = rij_inv[i] * (xi-xj);
|
||||
elec_dist_gl1[i + j * elec_num] = rij_inv[i] * (yi-yj);
|
||||
@ -7639,47 +7631,18 @@ qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_hpc (const qmckl_context co
|
||||
}
|
||||
}
|
||||
|
||||
for (int64_t j = 0; j < elec_num; ++j) {
|
||||
|
||||
double* restrict eegl = &een_rescaled_e_gl[ elec_num * 4 * (j + elec_num * (cord_num + 1) * nw)];
|
||||
#ifdef HAVE_OPENMP
|
||||
#pragma omp simd
|
||||
#endif
|
||||
for (int64_t i = 0; i < 4*elec_num; ++i) {
|
||||
eegl[i] = 0.0;
|
||||
}
|
||||
|
||||
}
|
||||
double* restrict eegl = &een_rescaled_e_gl[elec_num*4*elec_num*(cord_num+1)*nw];
|
||||
memset(eegl, 0, 4*elec_num*elec_num*sizeof(double));
|
||||
|
||||
for (int64_t l=1; l<=cord_num; ++l) {
|
||||
|
||||
double kappa_l = - (double)l * rescale_factor_ee;
|
||||
const double kappa_l = -rescale_factor_ee * (double) l;
|
||||
|
||||
for (int64_t j=0; j<elec_num; ++j) {
|
||||
|
||||
double* restrict eegl =
|
||||
&een_rescaled_e_gl[elec_num*4*(j+elec_num*(l+(cord_num+1)*nw))];
|
||||
|
||||
#ifdef HAVE_OPENMP
|
||||
#pragma omp simd
|
||||
#endif
|
||||
for (int64_t i=0; i<elec_num; ++i) {
|
||||
eegl[i ] = kappa_l * elec_dist_gl0[i + j * elec_num];
|
||||
eegl[i + elec_num ] = kappa_l * elec_dist_gl1[i + j * elec_num];
|
||||
eegl[i + elec_num * 2] = kappa_l * elec_dist_gl2[i + j * elec_num];
|
||||
eegl[i + elec_num * 3] = kappa_l * elec_dist_gl3[i + j * elec_num];
|
||||
}
|
||||
|
||||
#ifdef HAVE_OPENMP
|
||||
#pragma omp simd
|
||||
#endif
|
||||
for (int64_t i=0; i<elec_num; ++i) {
|
||||
eegl[i + elec_num*3] = eegl[i + elec_num*3] +
|
||||
eegl[i] * eegl[i] +
|
||||
eegl[i + elec_num*1] * eegl[i + elec_num*1] +
|
||||
eegl[i + elec_num*2] * eegl[i + elec_num*2];
|
||||
}
|
||||
|
||||
const double* restrict ee =
|
||||
&een_rescaled_e[elec_num*(j+elec_num*(l+(cord_num+1)*nw))];
|
||||
|
||||
@ -7687,11 +7650,15 @@ qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_hpc (const qmckl_context co
|
||||
#pragma omp simd
|
||||
#endif
|
||||
for (int64_t i=0; i<elec_num; ++i) {
|
||||
eegl[i ] *= ee[i];
|
||||
eegl[i + elec_num * 1] *= ee[i];
|
||||
eegl[i + elec_num * 2] *= ee[i];
|
||||
eegl[i + elec_num * 3] *= ee[i];
|
||||
eegl[i ] = kappa_l * ee[i] * elec_dist_gl0[i+j*elec_num];
|
||||
eegl[i+elec_num ] = kappa_l * ee[i] * elec_dist_gl1[i+j*elec_num];
|
||||
eegl[i+elec_num*2] = kappa_l * ee[i] * elec_dist_gl2[i+j*elec_num];
|
||||
eegl[i+elec_num*3] = kappa_l * ee[i] * (elec_dist_gl3[i+j*elec_num] + kappa_l);
|
||||
}
|
||||
eegl[j ] = 0.0;
|
||||
eegl[j+elec_num ] = 0.0;
|
||||
eegl[j+elec_num*2] = 0.0;
|
||||
eegl[j+elec_num*3] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -7730,79 +7697,36 @@ qmckl_exit_code qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl (
|
||||
|
||||
**** Test :noexport:
|
||||
|
||||
#+name: een_e_gl
|
||||
#+name: ee_dist_gl
|
||||
#+begin_src python :results output :exports none :noweb yes
|
||||
import numpy as np
|
||||
<<ee_dist>>
|
||||
|
||||
<<jastrow_data>>
|
||||
|
||||
elec_coord = np.array(elec_coord)[0]
|
||||
elec_dist = np.zeros(shape=(elec_num, elec_num),dtype=float)
|
||||
een_rescaled_e_gl = np.zeros(shape=(4, cord_num+1, elec_num, elec_num), dtype=float)
|
||||
for i in range(elec_num):
|
||||
for j in range(elec_num):
|
||||
elec_dist[i, j] = np.linalg.norm(elec_coord[i] - elec_coord[j])
|
||||
|
||||
elec_dist_gl = 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_gl[ii, i, j] = (elec_coord[i][ii] - elec_coord[j][ii]) * rij_inv
|
||||
elec_dist_gl[3, i, j] = 2.0 * rij_inv
|
||||
elec_dist_gl[:, j, j] = 0.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
|
||||
|
||||
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
|
||||
|
||||
een_rescaled_e_gl = np.zeros(shape=(elec_num,4,elec_num,cord_num+1),dtype=float)
|
||||
for l in range(0,cord_num+1):
|
||||
kappa_l = -1.0 * kappa * l
|
||||
for j in range(0,elec_num):
|
||||
for i in range(0,elec_num):
|
||||
for ii in range(0,4):
|
||||
een_rescaled_e_gl[i,ii,j,l] = kappa_l * elec_dist_gl[ii,i,j]
|
||||
een_rescaled_e_gl[i,3,j,l] = een_rescaled_e_gl[i,3,j,l] + \
|
||||
een_rescaled_e_gl[i,0,j,l] * een_rescaled_e_gl[i,0,j,l] + \
|
||||
een_rescaled_e_gl[i,1,j,l] * een_rescaled_e_gl[i,1,j,l] + \
|
||||
een_rescaled_e_gl[i,2,j,l] * een_rescaled_e_gl[i,2,j,l]
|
||||
|
||||
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])
|
||||
if j != i:
|
||||
r = elec_dist[i,j]
|
||||
for p in range(cord_num+1):
|
||||
f = een_rescaled_e[p,i,j]
|
||||
for k in range(3):
|
||||
een_rescaled_e_gl[k,p,i,j] = -kappa_e * p * (elec_coord[j,k] - elec_coord[i,k])/r * f
|
||||
een_rescaled_e_gl[3,p,i,j] = kappa_e*p*(-2.0/r + kappa_e*p)*f
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :tangle (eval c_test)
|
||||
#+NAME:test_ee_gl
|
||||
#+begin_src python :results output :exports none :noweb yes
|
||||
<<ee_dist_gl>>
|
||||
for p in range(cord_num+1):
|
||||
# for i in range(elec_num):
|
||||
for i in range(3):
|
||||
for k in range(4):
|
||||
# for j in range(elec_num):
|
||||
for j in range(3):
|
||||
print(f"printf( \"%e %e\\n\", een_rescaled_e_gl[0][{p}][{i}][{k}][{j}], {een_rescaled_e_gl[k,p,i,j]});")
|
||||
print(f"assert( fabs(een_rescaled_e_gl[0][{p}][{i}][{k}][{j}] - ({een_rescaled_e_gl[k,p,i,j]})) < 1.e-10 );")
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :tangle (eval c_test) :noweb yes
|
||||
assert(qmckl_electron_provided(context));
|
||||
|
||||
{
|
||||
@ -7810,13 +7734,7 @@ assert(qmckl_electron_provided(context));
|
||||
size_max=walk_num*(cord_num + 1)*elec_num*4*elec_num;
|
||||
rc = qmckl_get_jastrow_champ_een_distance_rescaled_e_gl(context,
|
||||
&(een_rescaled_e_gl[0][0][0][0][0]),size_max);
|
||||
|
||||
assert(fabs(een_rescaled_e_gl[0][1][0][0][2] + 0.09831391870751387 ) < 1.e-8);
|
||||
assert(fabs(een_rescaled_e_gl[0][1][0][0][3] + 0.017204157459682526 ) < 1.e-8);
|
||||
assert(fabs(een_rescaled_e_gl[0][1][0][0][4] + 0.013345768421098641 ) < 1.e-8);
|
||||
assert(fabs(een_rescaled_e_gl[0][2][1][0][3] + 0.03733086358273962 ) < 1.e-8);
|
||||
assert(fabs(een_rescaled_e_gl[0][2][1][0][4] + 0.004922634822943517 ) < 1.e-8);
|
||||
assert(fabs(een_rescaled_e_gl[0][2][1][0][5] + 0.5416751547830984 ) < 1.e-8);
|
||||
<<test_ee_gl()>>
|
||||
}
|
||||
|
||||
|
||||
@ -7848,19 +7766,27 @@ assert(qmckl_electron_provided(context));
|
||||
een_rescaled_e_gl_hpc);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
|
||||
for (int64_t i = 0; i < walk_num*(cord_num + 1)*elec_num*4*elec_num; i++) {
|
||||
if (fabs(een_rescaled_e_gl_doc[i] - een_rescaled_e_gl_hpc[i]) > 1.e-12) {
|
||||
printf("i = %ld, doc = %f, hpc = %f\n", i, een_rescaled_e_gl_doc[i], een_rescaled_e_gl_hpc[i]);
|
||||
for (int nw=0 ; nw < walk_num ; nw++) {
|
||||
for (int c=0 ; c <= cord_num ; c++) {
|
||||
for (int i=0 ; i < elec_num ; i++) {
|
||||
for (int j=0 ; j < elec_num ; j++) {
|
||||
for (int k=0 ; k < 4 ; k++) {
|
||||
printf("nw=%d c=%d i=%d k=%d j=%d doc=%e hpc=%e\n", nw, c, i, k, j, een_rescaled_e_gl_doc[nw*(cord_num + 1)*elec_num*4*elec_num + c*elec_num*4*elec_num + i*4*elec_num + k*elec_num + j], een_rescaled_e_gl_hpc[nw*(cord_num + 1)*elec_num*4*elec_num + c*elec_num*4*elec_num + i*4*elec_num + k*elec_num + j]);
|
||||
if (fabs(een_rescaled_e_gl_doc[nw*(cord_num + 1)*elec_num*4*elec_num + c*elec_num*4*elec_num + i*4*elec_num + k*elec_num + j] - een_rescaled_e_gl_hpc[nw*(cord_num + 1)*elec_num*4*elec_num + c*elec_num*4*elec_num + i*4*elec_num + k*elec_num + j]) > 1.e-12) {
|
||||
printf("nw=%d c=%d i=%d k=%d j=%d doc=%e hpc=%e\n", nw, c, i, k, j, een_rescaled_e_gl_doc[nw*(cord_num + 1)*elec_num*4*elec_num + c*elec_num*4*elec_num + i*4*elec_num + k*elec_num + j], een_rescaled_e_gl_hpc[nw*(cord_num + 1)*elec_num*4*elec_num + c*elec_num*4*elec_num + i*4*elec_num + k*elec_num + j]);
|
||||
fflush(stdout);
|
||||
}
|
||||
assert(fabs(een_rescaled_e_gl_doc[i] - een_rescaled_e_gl_hpc[i]) < 1.e-8);
|
||||
assert(fabs(een_rescaled_e_gl_doc[nw*(cord_num + 1)*elec_num*4*elec_num + c*elec_num*4*elec_num + i*4*elec_num + k*elec_num + j] - een_rescaled_e_gl_hpc[nw*(cord_num + 1)*elec_num*4*elec_num + c*elec_num*4*elec_num + i*4*elec_num + k*elec_num + j]) < 1.e-8);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
{
|
||||
/* Finite difference test fails and I can't understand why... */
|
||||
/*
|
||||
|
||||
printf("een_distance_rescaled_e_gl\n");
|
||||
|
||||
@ -7876,11 +7802,11 @@ assert(qmckl_electron_provided(context));
|
||||
|
||||
qmckl_exit_code rc;
|
||||
|
||||
double elec_coord[walk_num][3][elec_num];
|
||||
rc = qmckl_get_electron_coord (context, 'T', &(elec_coord[0][0][0]), 3*walk_num*elec_num);
|
||||
double elec_coord[walk_num][elec_num][3];
|
||||
rc = qmckl_get_electron_coord (context, 'N', &(elec_coord[0][0][0]), 3*walk_num*elec_num);
|
||||
assert (rc == QMCKL_SUCCESS);
|
||||
|
||||
double temp_coord[walk_num][3][elec_num];
|
||||
double temp_coord[walk_num][elec_num][3];
|
||||
memcpy(&(temp_coord[0][0][0]), &(elec_coord[0][0][0]), sizeof(temp_coord));
|
||||
|
||||
double function_values[walk_num][cord_num+1][elec_num][elec_num];
|
||||
@ -7890,12 +7816,13 @@ assert(qmckl_electron_provided(context));
|
||||
for (int64_t k = 0; k < 3; k++) {
|
||||
for (int64_t m = -4; m <= 4; m++) { // Apply finite difference displacement
|
||||
|
||||
memcpy(&(temp_coord[0][0][0]), &(elec_coord[0][0][0]), sizeof(temp_coord));
|
||||
for (int64_t nw=0 ; nw<walk_num ; nw++) {
|
||||
temp_coord[nw][k][i] = elec_coord[nw][k][i] + (double) m * delta_x;
|
||||
temp_coord[nw][i][k] = elec_coord[nw][i][k] + (double) m * delta_x;
|
||||
}
|
||||
|
||||
// Update coordinates in the context
|
||||
rc = qmckl_set_electron_coord (context, 'T', walk_num,
|
||||
rc = qmckl_set_electron_coord (context, 'N', walk_num,
|
||||
&(temp_coord[0][0][0]),
|
||||
walk_num*3*elec_num);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
@ -7917,14 +7844,11 @@ assert(qmckl_electron_provided(context));
|
||||
}
|
||||
|
||||
}
|
||||
for (int64_t nw=0 ; nw<walk_num ; nw++) {
|
||||
temp_coord[nw][k][i] = elec_coord[nw][k][i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Reset coordinates in the context
|
||||
rc = qmckl_set_electron_coord (context, 'T', walk_num,
|
||||
rc = qmckl_set_electron_coord (context, 'N', walk_num,
|
||||
&(elec_coord[0][0][0]),
|
||||
walk_num*3*elec_num);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
@ -7953,6 +7877,7 @@ assert(qmckl_electron_provided(context));
|
||||
walk_num*(cord_num+1)*elec_num*4*elec_num)
|
||||
);
|
||||
|
||||
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
|
||||
for (int nw = 0; nw < walk_num; nw++){
|
||||
@ -7983,9 +7908,7 @@ assert(qmckl_electron_provided(context));
|
||||
}
|
||||
}
|
||||
printf("OK\n");
|
||||
*/
|
||||
}
|
||||
|
||||
#+end_src
|
||||
|
||||
*** Electron-nucleus rescaled distances in $J_\text{eeN}$
|
||||
@ -7994,7 +7917,7 @@ assert(qmckl_electron_provided(context));
|
||||
electrons and nuclei raised to the power \(p\) defined by ~cord_num~:
|
||||
|
||||
\[
|
||||
C_{i\alpha,p} = \left[ \exp\left(-\kappa_\alpha\, R_{i\alpha}\right) \right]^p
|
||||
[g_{\alpha}(R_{i\alpha})]^p = \exp\left(-p\, \kappa_\alpha\, R_{i\alpha}\right)
|
||||
\]
|
||||
|
||||
where \(R_{i\alpha}\) is the matrix of electron-nucleus distances.
|
||||
@ -8741,6 +8664,20 @@ assert(fabs( 0.005359281880312882 - een_rescaled_n_gl[0][2][1][0][5]) < 1.e-12
|
||||
calculation of the three-body jastrow ~factor_een~ and its derivative
|
||||
~factor_een_gl~.
|
||||
|
||||
The array ~tmp_c~ corresponds to the tensor $P$ defined at the
|
||||
beginning of this section:
|
||||
|
||||
\[ P_{i\alpha}^{km} =
|
||||
\sum_{j=1}^{N_{\text{elec}}}
|
||||
\left[ g_\text{e}({r}_{ij}) \right]^k
|
||||
\left[ g_\alpha({R}_{j\alpha}) \right]^{m}
|
||||
\]
|
||||
|
||||
\[ \nabla_i P_{i\alpha}^{km} =
|
||||
\sum_{j=1}^{N_{\text{elec}}}
|
||||
\nabla_i \left[ g_\text{e}({r}_{ij}) \right]^k
|
||||
\left[ g_\alpha({R}_{j\alpha}) \right]^{m}
|
||||
\]
|
||||
**** Compute dim_c_vector
|
||||
:PROPERTIES:
|
||||
:Name: qmckl_compute_dim_c_vector
|
||||
@ -9693,7 +9630,7 @@ qmckl_exit_code qmckl_compute_tmp_c_doc (
|
||||
double* const tmp_c );
|
||||
#+end_src
|
||||
|
||||
***** CPU :noexport:
|
||||
***** HPC :noexport:
|
||||
|
||||
#+begin_src c :comments org :tangle (eval c) :noweb yes
|
||||
qmckl_exit_code qmckl_compute_tmp_c_hpc (
|
||||
@ -9802,8 +9739,9 @@ qmckl_exit_code qmckl_compute_tmp_c_hpc (const qmckl_context context,
|
||||
:END:
|
||||
|
||||
#+NAME: qmckl_factor_dtmp_c_args
|
||||
|---------------------+---------------------------------------------------------------------+--------+-----------------------------------------------|
|
||||
| Variable | Type | In/Out | Description |
|
||||
|---------------------+------------------------------------------------------------------+--------+-----------------------------------------------|
|
||||
|---------------------+---------------------------------------------------------------------+--------+-----------------------------------------------|
|
||||
| ~context~ | ~qmckl_context~ | in | Global state |
|
||||
| ~cord_num~ | ~int64_t~ | in | Order of polynomials |
|
||||
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
|
||||
@ -9811,7 +9749,8 @@ qmckl_exit_code qmckl_compute_tmp_c_hpc (const qmckl_context context,
|
||||
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
|
||||
| ~een_rescaled_e_gl~ | ~double[walk_num][0:cord_num][elec_num][4][elec_num]~ | in | Electron-electron rescaled factor derivatives |
|
||||
| ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled factor |
|
||||
| ~dtmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | out | vector of non-zero coefficients |
|
||||
| ~dtmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | out | vector of non-zero coefficients |
|
||||
|---------------------+---------------------------------------------------------------------+--------+-----------------------------------------------|
|
||||
|
||||
|
||||
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
||||
@ -9991,85 +9930,58 @@ qmckl_exit_code qmckl_compute_dtmp_c_hpc (
|
||||
import numpy as np
|
||||
|
||||
<<jastrow_data>>
|
||||
|
||||
elec_coord = np.array(elec_coord)[0]
|
||||
nucl_coord = np.array(nucl_coord)
|
||||
elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float)
|
||||
for i in range(elec_num):
|
||||
for a in range(nucl_num):
|
||||
elnuc_dist[i, a] = np.linalg.norm(elec_coord[i] - nucl_coord[:,a])
|
||||
|
||||
kappa = 0.6
|
||||
|
||||
een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float)
|
||||
een_rescaled_n[:,:,0] = 1.0
|
||||
|
||||
for a in range(nucl_num):
|
||||
for i in range(elec_num):
|
||||
een_rescaled_n[a, i, 1] = np.exp(-kappa * elnuc_dist[i, a])
|
||||
|
||||
for l in range(2,cord_num+1):
|
||||
for a in range(nucl_num):
|
||||
for i in range(elec_num):
|
||||
een_rescaled_n[a, i, l] = een_rescaled_n[a, i, l - 1] * een_rescaled_n[a, i, 1]
|
||||
|
||||
elec_dist = np.zeros(shape=(elec_num, elec_num),dtype=float)
|
||||
for i in range(elec_num):
|
||||
for j in range(elec_num):
|
||||
elec_dist[i, j] = np.linalg.norm(elec_coord[i] - elec_coord[j])
|
||||
|
||||
kappa = 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
|
||||
|
||||
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
|
||||
<<en_dist>>
|
||||
<<ee_dist_gl>>
|
||||
#+end_src
|
||||
|
||||
#+RESULTS: helper_funcs
|
||||
|
||||
#+begin_src c :tangle (eval c_test)
|
||||
assert(qmckl_electron_provided(context));
|
||||
#+NAME: test_tmpc
|
||||
#+begin_src python :results output :exports none :noweb yes
|
||||
<<helper_funcs>>
|
||||
tmp_c = np.einsum("pij,qaj->pqai", een_rescaled_e, een_rescaled_n)
|
||||
for p in range(cord_num):
|
||||
for q in range(cord_num+1):
|
||||
for a in range(nucl_num):
|
||||
# for i in range(elec_num):
|
||||
for i in range(3):
|
||||
print(f"assert( fabs(tmp_c[0][{p}][{q}][{a}][{i}] - ({tmp_c[p,q,a,i]})) < 1.e-10 );")
|
||||
#+end_src
|
||||
|
||||
double tmp_c[walk_num][cord_num][cord_num+1][nucl_num][elec_num];
|
||||
rc = qmckl_get_jastrow_champ_tmp_c(context, &(tmp_c[0][0][0][0][0]), sizeof(tmp_c)/sizeof(double));
|
||||
#+NAME: test_dtmpc
|
||||
#+begin_src python :results output :exports none :noweb yes
|
||||
<<helper_funcs>>
|
||||
dtmp_c = np.einsum("lpji,qaj->pqali", een_rescaled_e_gl, een_rescaled_n)
|
||||
for p in range(cord_num):
|
||||
for q in range(cord_num+1):
|
||||
for a in range(nucl_num):
|
||||
for l in range(4):
|
||||
# for i in range(elec_num):
|
||||
for i in range(3):
|
||||
print(f"printf(\"%e %e\\n\", dtmp_c[0][{p}][{q}][{a}][{l}][{i}], {dtmp_c[p,q,a,l,i]}); fflush(stdout);")
|
||||
print(f"assert( fabs(dtmp_c[0][{p}][{q}][{a}][{l}][{i}] - ({dtmp_c[p,q,a,l,i]})) < 1.e-10 );")
|
||||
#+end_src
|
||||
|
||||
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]), sizeof(dtmp_c)/sizeof(double));
|
||||
#+RESULTS: test_dtmpc
|
||||
|
||||
printf("%e\n%e\n", tmp_c[0][0][1][0][0], 3.954384);
|
||||
fflush(stdout);
|
||||
assert(fabs(tmp_c[0][0][1][0][0] - 3.954384) < 1e-6);
|
||||
#+begin_src c :tangle (eval c_test) :noweb yes
|
||||
{
|
||||
assert(qmckl_electron_provided(context));
|
||||
|
||||
printf("%e\n%e\n", dtmp_c[0][1][0][0][0][0],3.278657e-01);
|
||||
fflush(stdout);
|
||||
assert(fabs(dtmp_c[0][1][0][0][0][0] - 3.278657e-01 ) < 1e-6);
|
||||
printf("tmp_c\n");
|
||||
double tmp_c[walk_num][cord_num][cord_num+1][nucl_num][elec_num];
|
||||
rc = qmckl_get_jastrow_champ_tmp_c(context, &(tmp_c[0][0][0][0][0]), sizeof(tmp_c)/sizeof(double));
|
||||
|
||||
<<test_tmpc()>>
|
||||
}
|
||||
|
||||
{
|
||||
printf("dtmp_c\n");
|
||||
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]), sizeof(dtmp_c)/sizeof(double));
|
||||
|
||||
<<test_dtmpc()>>
|
||||
}
|
||||
#+end_src
|
||||
|
||||
*** Electron-electron-nucleus Jastrow $f_{een}$
|
||||
@ -11386,8 +11298,6 @@ import numpy as np
|
||||
|
||||
<<jastrow_data>>
|
||||
|
||||
<<een_e_gl>>
|
||||
|
||||
<<helper_funcs>>
|
||||
|
||||
kappa = 0.6
|
||||
|
||||
Loading…
Reference in New Issue
Block a user