qmckl/org/qmckl_jastrow_champ_single.org

CHAMP Jastrow Factor Single

• Introduction
• Context
  • Data structure
• Single point
  • Set
  • Touch
• Electron-electron and electron-nucleus distances for single point
  • Electron-electron distances
    • Get
    • Compute
  • Electron-nucleus distances
    • Get
    • Compute
• Electron-electron-nucleus Jastrow
  • Electron-electron rescaled distances
    • Get
    • Compute
  • Electron-nucleus rescaled distances
    • Get
    • Compute
  • $\delta P$ matrix
    • Get
    • Compute
  • Electron-electron-nucleus Jastrow value
    • Get
    • Compute
  • Electron-nucleus rescaled distance derivative
    • Get
    • Compute
  • Electron-electron rescaled distances derivative
    • Get
    • Compute
  • $\delta P$ matrix gradients and Laplacian
    • Get
    • Compute
  • Electron-electron-nucleus Jastrow gradients and Laplacian
    • Get
    • Compute
  • $\delta P$ matrix gradient only
    • Get
    • Compute
  • Electron-electron-nucleus Jastrow gradients
    • Get
    • Compute
• Electron-electron Jastrow
  • Electron-electron rescaled distance
    • Get
    • Compute
  • Electron-electron Jastrow value
    • Get
    • Compute
  • Electron-electron rescaled distances derivatives
    • Get
    • Compute
  • Electron-electron Jastrow gradients and Laplacian
    • Get
    • Compute
• Electron-nucleus Jastrow
  • Electron-nucleus rescaled distance
    • Get
    • Compute
  • Single electron-nucleus Jastrow value
    • Get
    • Compute
  • Electron-nucleus rescaled distances derivatives
    • Get
    • Compute
  • Electron-nucleus Jastrow gradients and Laplacian
    • Get
    • Compute
• Accept single electron move
  • Code

Introduction

Single-electron move version of the Jastrow factor functions. The single-electron move calculates the difference between the old Jastrow values, gradients and derivatives and the new ones (if the single electron would have been moved). That is to say, it calculates \[ \delta J = J(\mathbf{r}^\prime,\mathbf{R}) - J(\mathbf{r},\mathbf{R}) \] for all the neccessery quantities.

Context

Data structure

typedef struct qmckl_jastrow_champ_single_struct{
int64_t      num;
uint64_t     date;
qmckl_matrix coord;
double  *    een_rescaled_single_e;
uint64_t     een_rescaled_single_e_date;
uint64_t     een_rescaled_single_e_maxsize;
double  *    een_rescaled_single_n;
uint64_t     een_rescaled_single_n_date;
uint64_t     een_rescaled_single_n_maxsize;
double*      single_ee_distance;
uint64_t     single_ee_distance_date;
uint64_t     single_ee_distance_maxsize;
double*      single_en_distance;
uint64_t     single_en_distance_date;
uint64_t     single_en_distance_maxsize;
double*      delta_een;
uint64_t     delta_een_date;
uint64_t     delta_een_maxsize;
double*      delta_p;
uint64_t     delta_p_date;
uint64_t     delta_p_maxsize;
double*      ee_rescaled_single;
uint64_t     ee_rescaled_single_date;
uint64_t     ee_rescaled_single_maxsize;
double*      en_rescaled_single;
uint64_t     en_rescaled_single_date;
uint64_t     en_rescaled_single_maxsize;
double*      delta_en;
uint64_t     delta_en_date;
uint64_t     delta_en_maxsize;
double*      delta_ee;
uint64_t     delta_ee_date;
uint64_t     delta_ee_maxsize;
double  *    een_rescaled_single_e_gl;
uint64_t     een_rescaled_single_e_gl_date;
uint64_t     een_rescaled_single_e_gl_maxsize;
double  *    een_rescaled_single_n_gl;
uint64_t     een_rescaled_single_n_gl_date;
uint64_t     een_rescaled_single_n_gl_maxsize;
double*      delta_p_gl;
uint64_t     delta_p_gl_date;
uint64_t     delta_p_gl_maxsize;
double*      delta_p_g;
uint64_t     delta_p_g_date;
uint64_t     delta_p_g_maxsize;
double*      delta_een_gl;
uint64_t     delta_een_gl_date;
uint64_t     delta_een_gl_maxsize;
double*      delta_een_g;
uint64_t     delta_een_g_date;
uint64_t     delta_een_g_maxsize;
double*      ee_rescaled_single_gl;
uint64_t     ee_rescaled_single_gl_date;
uint64_t     ee_rescaled_single_gl_maxsize;
double*      en_rescaled_single_gl;
uint64_t     en_rescaled_single_gl_date;
uint64_t     en_rescaled_single_gl_maxsize;
double*      delta_en_gl;
uint64_t     delta_en_gl_date;
uint64_t     delta_en_gl_maxsize;
double*      delta_ee_gl;
uint64_t     delta_ee_gl_date;
uint64_t     delta_ee_gl_maxsize;

} qmckl_jastrow_champ_single_struct;

Single point

Set

We set the coordinates of the num-th electron for all walkers, where num is the electron which has to be moved. The dimension of coord is

• [walk_num][3] if transp is 'N'
• [3][walk_num] if transp is 'T'

Internally, the coordinates are stored in 'N' format as opposed to elec_coord. This function has to be called before any other functions from this module.

qmckl_exit_code qmckl_set_single_point (qmckl_context context,
                              const char transp,
                              const int64_t num,
                              const double* coord,
                              const int64_t size_max);

The Fortran function shifts the num by 1 because of 1-based indexing.

qmckl_exit_code qmckl_set_single_point_f (qmckl_context context,
                              const char transp,
                              const int64_t num,
                              const double* coord,
                              const int64_t size_max);

qmckl_exit_code
qmckl_set_single_point (qmckl_context context,
              const char transp,
              const int64_t num,
              const double* coord,
              const int64_t size_max)
{

if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
 return QMCKL_NULL_CONTEXT;
}

if (num < 0) {
   return qmckl_failwith( context,
                          QMCKL_INVALID_ARG_3,
                          "qmckl_set_single_point",
                          "Incorrect point number");
}

if (transp != 'N' && transp != 'T') {
 return qmckl_failwith( context,
                        QMCKL_INVALID_ARG_2,
                        "qmckl_set_single_point",
                        "transp should be 'N' or 'T'");
}

if (coord == NULL) {
 return qmckl_failwith( context,
                        QMCKL_INVALID_ARG_3,
                        "qmckl_set_single_point",
                        "coord is a NULL pointer");
}

qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);

int64_t walk_num = ctx->electron.walker.num;

if (size_max < 3*walk_num) {
   return qmckl_failwith( context,
                          QMCKL_INVALID_ARG_4,
                          "qmckl_set_single_point",
                          "Array too small");
}

//qmckl_exit_code rc;

//if (ctx->single_point.coord.data != NULL) {
//  rc = qmckl_matrix_free(context, &(ctx->single_point.coord));
//  assert (rc == QMCKL_SUCCESS);
//}

if (ctx->single_point.coord.data == NULL) {
 ctx->single_point.coord = qmckl_matrix_alloc(context, walk_num, 3);
 if (ctx->single_point.coord.data == NULL) {
   return qmckl_failwith( context,
                         QMCKL_ALLOCATION_FAILED,
                         "qmckl_set_single_point",
                         NULL);
 }
}

ctx->single_point.num = num;

if (transp == 'N') {
 double *a = ctx->single_point.coord.data;
 for (int64_t i=0 ; i<3*walk_num ; ++i) {
   a[i] = coord[i];
 }
} else {
 for (int64_t i=0 ; i<walk_num ; ++i) {
   qmckl_mat(ctx->single_point.coord, i, 0) = coord[i*walk_num + 0];
   qmckl_mat(ctx->single_point.coord, i, 1) = coord[i*walk_num + 1];
   qmckl_mat(ctx->single_point.coord, i, 2) = coord[i*walk_num + 2];
 }
}

/* Increment the date of the single point */
ctx->single_point.date += 1UL;

return QMCKL_SUCCESS;

}

qmckl_exit_code
qmckl_set_single_point_f (qmckl_context context,
              const char transp,
              const int64_t num,
              const double* coord,
              const int64_t size_max)
{
return qmckl_set_single_point(context, transp, num-1, coord, size_max);
}

interface
integer(qmckl_exit_code) function qmckl_set_single_point(context, &
    transp, num, coord, size_max) bind(C, name="qmckl_set_single_point_f")
 use, intrinsic :: iso_c_binding
 import
 implicit none

 integer (c_int64_t) , intent(in)  , value :: context
 character(c_char)   , intent(in)  , value :: transp
 integer (c_int64_t) , intent(in)  , value :: num
 real    (c_double ) , intent(in)          :: coord(*)
 integer (c_int64_t) , intent(in)  , value :: size_max
end function
end interface

Touch

qmckl_exit_code
qmckl_single_touch (const qmckl_context context);

Electron-electron and electron-nucleus distances for single point

In order to calculate the $\delta J$, the updated electron-electron and electron-nucleus distances are required.

Electron-electron distances

Electron-electron distance between the single electron and all electrons for all walkers. Dimension is [walk_num][elec_num].

Get

qmckl_exit_code qmckl_get_single_electron_ee_distance(qmckl_context context,
                                                  double* const distance,
                                                  const int64_t size_max);

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
num	int64_t	in	Index of single electron
elec_num	int64_t	in	Number of electrons
walk_num	int64_t	in	Number of walkers
coord	double[3][walk_num][elec_num]	in	Electron coordinates
single_coord	double[walk_num][3]	in	Single electron coordinates
single_ee_distance	double[walk_num][elec_num]	out	Electron-electron distances for single electron

integer(qmckl_exit_code) function qmckl_compute_single_ee_distance(context, &
 num_in, elec_num, walk_num, coord, single_coord, single_ee_distance) &
 result(info) bind(C)
use, intrinsic :: iso_c_binding
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer (c_int64_t) , intent(in)  , value :: elec_num, num_in
integer (c_int64_t) , intent(in)  , value :: walk_num
real    (c_double ) , intent(in)          :: coord(elec_num,walk_num,3)
real    (c_double ) , intent(in)          :: single_coord(3,walk_num)
real    (c_double ) , intent(out)         :: single_ee_distance(elec_num,walk_num)

integer*8 :: k, i, j, num
double precision :: x, y, z

info = QMCKL_SUCCESS

num = num_in + 1

if (context == QMCKL_NULL_CONTEXT) then
 info = QMCKL_INVALID_CONTEXT
 return
endif

if (elec_num <= 0) then
 info = QMCKL_INVALID_ARG_2
 return
endif

if (walk_num <= 0) then
 info = QMCKL_INVALID_ARG_3
 return
endif

do k=1,walk_num
 info = qmckl_distance(context, 'T', 'N', elec_num, 1_8, &
      coord(1,k,1), elec_num*walk_num, &
      single_coord(1,k), 3_8, &
      single_ee_distance(1,k), elec_num)
 if (info /= QMCKL_SUCCESS) then
    exit
 endif
 single_ee_distance(num,k) = 0.0d0
end do



end function qmckl_compute_single_ee_distance

Electron-nucleus distances

Electron-nucleus distance between the single electron and all nuclei for all walkers. Dimension is [walk_num][nucl_num].

Get

qmckl_exit_code
qmckl_get_single_electron_en_distance(qmckl_context context,
                                  double* distance,
                                  const int64_t size_max);

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
nucl_num	int64_t	in	Number of nuclei
walk_num	int64_t	in	Number of walkers
elec_coord	double[3][walk_num]	in	Electron coordinates
nucl_coord	double[3][nucl_num]	in	Nuclear coordinates
single_en_distance	double[walk_num][nucl_num]	out	Electron-nucleus distances for single-electron

integer function qmckl_compute_single_en_distance(context, nucl_num, walk_num,  &
 elec_coord, nucl_coord, single_en_distance) result(info) bind(C)
use, intrinsic :: iso_c_binding
use qmckl
implicit none
integer(qmckl_context), intent(in)        :: context
integer (c_int64_t) , intent(in)  , value :: nucl_num, walk_num
real    (c_double ) , intent(in)          :: elec_coord(3,walk_num)
real    (c_double ) , intent(in)          :: nucl_coord(nucl_num,3)
real    (c_double ) , intent(out)         :: single_en_distance(nucl_num, walk_num)

integer*8 :: k

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
 info = QMCKL_INVALID_CONTEXT
 return
endif

if (nucl_num <= 0) then
 info = QMCKL_INVALID_ARG_2
 return
endif

info = qmckl_distance(context, 'T', 'N', nucl_num, walk_num, &
      nucl_coord, nucl_num, &
      elec_coord, 3_8, &
      single_en_distance, nucl_num)

end function qmckl_compute_single_en_distance

Electron-electron-nucleus Jastrow

Calcules the single electron contributions to the electron-electron-nucleus Jastrow term. First, the rescaled distances for the single electron electron-nucleus and electron-electron distances are computed. These rescaled distances can be calculates by \[ \widetilde{R}_{i\alpha} = e^{-\kappa l R_{i\alpha}} \quad \text{and} \quad \widetilde{r}_{ij} = e^{-\kappa l r_{ij}}. \] Here, $\kappa$ is the rescaling factor and $l$ the power. The neccessery powers are stored in the same array. From these, we can calculate the $\delta \widetilde{R}$ and $\delta \widetilde{r}$ using

\begin{eqnarray*} \delta \widetilde{R}_{i\alpha} &=& \widetilde{R}_{i\alpha}^\text{new} - \widetilde{R}_{i\alpha}^\text{old} \quad \text{and}\\ \delta \widetilde{r}_{ij} &=& \widetilde{r}_{ij}^\text{new} - \widetilde{r}_{ij}^\text{old}. \end{eqnarray*}

With these, we can now calculate the single electron contribution to the $\delta P$ matrix, using the equation

\begin{eqnarray*} \delta P_{i,\alpha,k,l} &=& P^{\text{new}}_{i,\alpha,k,l} - P^{\text{old}}_{i,\alpha,k,l}\\ &=& \sum_{j=1}^{N_\text{elec}} \left(\widetilde{r}_{i,j,k} \delta \widetilde{R}_{j,\alpha,l}\delta_{j,\text{num}} + \delta \widetilde{r}_{i,j,k} \widetilde{R}_{j,\alpha,l} (\delta_{j,\text{num}} + \delta_{i,\text{num}}) + \delta \widetilde{r}_{i,j,k} \delta \widetilde{R}_{j,\alpha,l} \delta_{j,\text{num}}\right)\\ &=& \sum_{j=1}^{N_\text{elec}} \left( \delta \widetilde{r}_{i,j,k} \widetilde{R}_{j,\alpha,l} \delta_{i,\text{num}} \right) + \widetilde{r}_{i,\text{num},k} \delta \widetilde{R}_{\text{num},\alpha,l} + \delta \widetilde{r}_{i,\text{num},k} \left( \widetilde{R}_{\text{num},\alpha,l} + \delta \widetilde{R}_{\text{num},\alpha,l} \right). \end{eqnarray*}

Then, the electron-electron-nucleus Jastrow value can be calculated by

\begin{eqnarray*} \delta J_{een} &=& J_{een}^{\text{new}} - J_{een}^{\text{old}}\\ &=&\sum_{p=2}^{N_\text{nord}} \sum_{k=0}^{p-1} \sum_{l=0}^{p-k-2\delta_{k,0}}\sum_{\alpha=1}^{N_\text{nucl}} c_{l,k,p,\alpha} \sum_{i=1}^{N_\text{elec}} \left( \delta \widetilde{R}_{i,\alpha,(p-k-l)/2} P_{i,\alpha,k,(p-k+l)/2} \delta_{i,\text{num}} + \widetilde{R}_{i,\alpha,(p-k-l)/2} \delta P_{i,\alpha,k,(p-k+l)/2} + \delta \widetilde{R}_{i,\alpha,(p-k-l)/2} \delta P_{i,\alpha,k,(p-k+l)/2} \delta_{i,\text{num}} \right)\\ &=& \sum_{p=2}^{N_\text{nord}} \sum_{k=0}^{p-1} \sum_{l=0}^{p-k-2\delta_{k,0}}\sum_{\alpha=1}^{N_\text{nucl}} c_{l,k,p,\alpha} \left( \sum_{i=1}^{N_\text{elec}} \left( \widetilde{R}_{i,\alpha,(p-k-l)/2} \delta P_{i,\alpha,k,(p-k+l)/2} \right) + \delta \widetilde{R}_{\text{num},\alpha,(p-k-l)/2} \left(P_{\text{num},\alpha,k,(p-k+l)/2} + \delta P_{\text{num},\alpha,k,(p-k+l)/2} \right)\right) \end{eqnarray*}

To calculate the gradients and Laplacian of the electron-electron-nucleus Jastrow, we first have to calculate the gradients and Laplacian of the rescaled distances, \[ \partial_{i,m} \widetilde{R}_{i\alpha} = -\kappa l e^{-\kappa l R_{i\alpha}} \frac{r_{i,m} - R_{\alpha,m}}{R_{i\alpha}} \quad \text{and} \quad \partial_{i,4} \widetilde{R}_{i\alpha} = -\kappa l \left(\frac{2}{R_{i\alpha}}- \kappa l \right) e^{-\kappa l R_{i\alpha}}, \] where $i$ is the electron of which we are taking the derivative and $m=1:3$ are the gradients and $m=4$ is the Laplacian. The derivatives of the single electron rescaled electron-nucleus distances are only nonzero when $i=\text{num}$. Similarly for $r$ we get \[ \partial_{i,m} \widetilde{r}_{ij} = -\kappa l e^{-\kappa l r_{ij}} \frac{r_{i,m} - r_{j,m}}{r_{ij}} \quad \text{and} \quad \partial_{i,4} \widetilde{r}_{ij} = -\kappa l \left(\frac{2}{r_{ij}}- \kappa l \right) e^{-\kappa l r_{ij}}. \]

With these, we can now calculate the gradient and Laplacian of the $\delta P$ matrix, using the equation

\begin{eqnarray*} \partial_{i,m} \delta P_{i,\alpha,k,l} &=& \partial_{i,m} P_{i,\alpha,k,l}^\text{new} - \partial_{i,m} P_{i,\alpha,k,l}^\text{old}\\ &=& \partial_{i,m}\delta \widetilde{r}_{\text{num},i,k} \left(\partial_{i,m}\delta \widetilde{R}_{\text{num},\alpha,l} + \partial_{i,m}\widetilde{R}_{\text{num},\alpha,l} \right) g_m + \partial_{i,m}\widetilde{r}_{\text{num},i,k} \delta \widetilde{R}_{\text{num},\alpha,l} + \delta_{i,\text{num}} \sum_{j=1}^{N_\text{elec}} \left( \partial_{j,m} \delta \widetilde{r}_{\text{num},j,k} \widetilde{R}_{j,\alpha,l} \right), \end{eqnarray*}

where $g_m = \{-1,-1,-1,1\}$.

Then, the gradient and Laplacian of the electron-electron-nucleus Jastrow value can be calculated by

\begin{eqnarray*} \delta \partial_{i,m} J_{een} &=& \partial_{i,m} J_{een}^{\text{new}} - \partial_{i,m} J_{een}^{\text{old}}\\ &=&\sum_{p=2}^{N_\text{nord}} \sum_{k=0}^{p-1} \sum_{l=0}^{p-k-2\delta_{k,0}}\sum_{\alpha=1}^{N_\text{nucl}} c_{l,k,p,\alpha} \left( \widetilde{R}_{i,\alpha,(p-k-l)/2} \partial_{i,m} \delta P_{i,\alpha,k,(p-k+l)/2} + \widetilde{R}_{i,\alpha,(p-k+l)/2} \partial_{i,m} \delta P_{i,\alpha,k,(p-k-l)/2} \right. \\ & &\ + \partial_{i,m}\widetilde{R}_{i,\alpha,(p-k-l)/2} \delta P_{i,\alpha,k,(p-k+l)/2} + \partial_{i,m}\widetilde{R}_{i,\alpha,(p-k+l)/2} \delta P_{i,\alpha,k,(p-k-l)/2} \\ & &\ + \delta_{i,\text{num}} \left( \delta \widetilde{R}_{i,\alpha,(p-k-l)/2} \left ( \partial_{i,m} P_{i,\alpha,k,(p-k+l)/2} + \partial_{i,m} \delta P_{i,\alpha,k,(p-k+l)/2} \right) + \delta \widetilde{R}_{i,\alpha,(p-k+l)/2} \left ( \partial_{i,m} P_{i,\alpha,k,(p-k-l)/2} + \partial_{i,m} \delta P_{i,\alpha,k,(p-k-l)/2} \right) \right. \\ & &\ \left. + \partial_{i,m} \delta \widetilde{R}_{i,\alpha,(p-k-l)/2} \left ( P_{i,\alpha,k,(p-k+l)/2} + \delta P_{i,\alpha,k,(p-k+l)/2} \right) + \partial_{i,m} \delta \widetilde{R}_{i,\alpha,(p-k+l)/2} \left ( P_{i,\alpha,k,(p-k-l)/2} + \delta P_{i,\alpha,k,(p-k-l)/2} \right) \right)\\ & &\ + \delta_{m,4} \sum_{d=1}^3 \left( \partial_{i,d} \widetilde{R}_{i,\alpha,(p-k-l)/2} \partial_{i,d} \delta P_{i,\alpha,k,(p-k+l)/2} + \partial_{i,d} \widetilde{R}_{i,\alpha,(p-k+l)/2} \partial_{i,d} \delta P_{i,\alpha,k,(p-k-l)/2} \right)\\ & &\ \left. + \delta_{m,4}\delta_{i,\text{num}} \sum_{d=1}^3\left( \partial_{i,d}\delta \widetilde{R}_{i,\alpha,(p-k-l)/2} \left( \partial_{i,d} P_{i,\alpha,k,(p-k+l)/2} + \partial_{i,d}\delta P_{i,\alpha,k,(p-k+l)/2} \right) + \partial_{i,d}\delta \widetilde{R}_{i,\alpha,(p-k+l)/2} \left( \partial_{i,d} P_{i,\alpha,k,(p-k-l)/2} + \partial_{i,d} \delta P_{i,\alpha,k,(p-k-l)/2} \right) \right) \right) \end{eqnarray*}

Electron-electron rescaled distances

Computes the new components of the $\widetilde{r}_{ij}$ matrix. Because only one electron is moved, the new matrix only has one dimension of size elec_num.

Get

qmckl_exit_code
qmckl_get_een_rescaled_single_e(qmckl_context context,
                            double* const distance_rescaled,
                            const int64_t size_max);

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
num	int64_t	in	Number of single electron
walk_num	int64_t	in	Number of walkers
elec_num	int64_t	in	Number of electrons
cord_num	int64_t	in	Order of polynomials
rescale_factor_ee	double	in	Factor to rescale ee distances
single_ee_distance	double[walk_num][elec_num]	in	Single electron-electron distances for each walker
een_rescaled_e	double[walk_num][0:cord_num][elec_num][elec_num]	in	Rescaled electron-electron distances for each walker
een_rescaled_single_e	double[walk_num][0:cord_num][elec_num]	out	Single electron-electron rescaled distances for each walker

integer function qmckl_compute_een_rescaled_single_e_doc( &
context, num_in, walk_num, elec_num, cord_num, rescale_factor_ee,  &
single_ee_distance, een_rescaled_e, een_rescaled_single_e) &
result(info) bind(C)
use, intrinsic :: iso_c_binding
use qmckl
implicit none
integer(qmckl_context), intent(in)         :: context
integer(c_int64_t)    , intent(in), value  :: num_in
integer(c_int64_t)    , intent(in), value  :: walk_num
integer(c_int64_t)    , intent(in), value  :: elec_num
integer(c_int64_t)    , intent(in), value  :: cord_num
real(c_double)        , intent(in), value  :: rescale_factor_ee
real(c_double)        , intent(in)         :: single_ee_distance(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_single_e(elec_num,0:cord_num,walk_num)

double precision,allocatable        :: een_rescaled_single_e_ij(:,:)
double precision                    :: x
integer*8                           :: i, j, k, l, nw, num

num = num_in + 1
info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif

if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif

if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif

if (cord_num < 0) then
info = QMCKL_INVALID_ARG_4
return
endif

allocate(een_rescaled_single_e_ij(elec_num, cord_num + 1))

! Prepare table of exponentiated distances raised to appropriate power
do nw = 1, walk_num
een_rescaled_single_e_ij(:, 1) = 1.0d0


do j = 1, elec_num
   een_rescaled_single_e_ij(j, 2) = dexp(-rescale_factor_ee * single_ee_distance(j, nw))
end do


do l = 2, cord_num
   do k = 1, elec_num
      een_rescaled_single_e_ij(k, l + 1) = een_rescaled_single_e_ij(k, l) * een_rescaled_single_e_ij(k, 2)
   end do
end do

! prepare the actual een table
een_rescaled_single_e(:,0,nw) = 1.0d0

do l = 1, cord_num
   do j = 1, elec_num
       x = een_rescaled_single_e_ij(j, l + 1)
       een_rescaled_single_e(j, l, nw) = x
   end do
end do

een_rescaled_single_e(num, :, :) = 0.0d0

end do

end function qmckl_compute_een_rescaled_single_e_doc

Electron-nucleus rescaled distances

Computes the new components of the $\widetilde{R}_{i\alpha}$ matrix.

Get

qmckl_exit_code
qmckl_get_een_rescaled_single_n(qmckl_context context,
                            double* const distance_rescaled,
                            const int64_t size_max);

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
num	int64_t	in	Number of single electron
walk_num	int64_t	in	Number of walkers
elec_num	int64_t	in	Number of atoms
nucl_num	int64_t	in	Number of atoms
type_nucl_num	int64_t	in	Number of atom types
type_nucl_vector	int64_t[nucl_num]	in	Types of atoms
cord_num	int64_t	in	Order of polynomials
rescale_factor_en	double[nucl_num]	in	Factor to rescale ee distances
single_en_distance	double[walk_num][nucl_num]	in	Electron-nucleus distances
een_rescaled_n	double[walk_num][0:cord_num][nucl_num][elec_num]	in	Electron-nucleus rescaled distances
een_rescaled_single_n	double[walk_num][0:cord_num][nucl_num]	out	Single electron-nucleus rescaled distances

integer function qmckl_compute_een_rescaled_single_n( &
context, num_in, walk_num, elec_num, nucl_num, &
type_nucl_num, type_nucl_vector, cord_num, rescale_factor_en,  &
single_en_distance, een_rescaled_n, een_rescaled_single_n) &
result(info) bind(C)
use, intrinsic :: iso_c_binding
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer(c_int64_t)    , intent(in), value  :: num_in
integer(c_int64_t)    , intent(in), value  :: walk_num
integer(c_int64_t)    , intent(in), value  :: elec_num
integer(c_int64_t)    , intent(in), value  :: nucl_num
integer(c_int64_t)    , intent(in), value  :: type_nucl_num
integer(c_int64_t)    , intent(in)         :: type_nucl_vector(nucl_num)
integer(c_int64_t)    , intent(in), value  :: cord_num
real(c_double)        , intent(in)         :: rescale_factor_en(type_nucl_num)
real(c_double)        , intent(in)         :: single_en_distance(nucl_num,walk_num)
real(c_double)        , intent(in)         :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
real(c_double)        , intent(out)        :: een_rescaled_single_n(nucl_num,0:cord_num,walk_num)

double precision                    :: x
integer*8                           :: i, a, k, l, nw, num

num = num_in + 1

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif

if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif

if (nucl_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif

if (cord_num < 0) then
info = QMCKL_INVALID_ARG_4
return
endif

do nw = 1, walk_num

! prepare the actual een table
een_rescaled_single_n(:, 0, nw) = 1.0d0

do a = 1, nucl_num
   een_rescaled_single_n(a, 1, nw) = dexp(-rescale_factor_en(type_nucl_vector(a)+1) * single_en_distance(a, nw))
end do

do l = 2, cord_num
   do a = 1, nucl_num
      een_rescaled_single_n(a, l, nw) = een_rescaled_single_n(a, l - 1, nw) * een_rescaled_single_n(a, 1, nw)
   end do
end do

end do

end function qmckl_compute_een_rescaled_single_n

$\delta P$ matrix

Computes the $\delta P$ matrix.

\begin{eqnarray*} \delta P_{i,\alpha,k,l} = \sum_{j=1}^{N_\text{elec}} \left( \delta \widetilde{r}_{i,j,k} \widetilde{R}_{j,\alpha,l} \delta_{i,\text{num}} \right) + \widetilde{r}_{i,\text{num},k} \delta \widetilde{R}_{\text{num},\alpha,l} + \delta \widetilde{r}_{i,\text{num},k} \left( \widetilde{R}_{\text{num},\alpha,l} + \delta \widetilde{R}_{\text{num},\alpha,l} \right). \end{eqnarray*}

Get

qmckl_exit_code
qmckl_get_jastrow_champ_delta_p(qmckl_context context,
                            double* const delta_p,
                            const int64_t size_max);

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
num	int64_t	in	Single point index
walk_num	int64_t	in	Number of walkers
elec_num	int64_t	in	Number of electrons
nucl_num	int64_t	in	Number of nuclei
cord_num	int64_t	in	order of polynomials
een_rescaled_n	double[walk_num][0:cord_num][nucl_num][elec_num]	in	Electron-nucleus rescaled distances
een_rescaled_e	double[walk_num][0:cord_num][elec_num][elec_num]	in	Electron-electron rescaled distances
een_rescaled_single_n	double[walk_num][0:cord_num][nucl_num]	in	Electron-nucleus single rescaled distances
een_rescaled_single_e	double[walk_num][0:cord_num][elec_num]	in	Electron-electron single rescaled distances
delta_p	double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]	out	Delta P matrix

integer function qmckl_compute_jastrow_champ_delta_p_doc( &
context, num_in, walk_num, elec_num, nucl_num, cord_num,   &
een_rescaled_n, een_rescaled_e, een_rescaled_single_n, een_rescaled_single_e, delta_p) &
result(info) bind(C)
use, intrinsic :: iso_c_binding
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer(c_int64_t), intent(in), value  :: num_in, walk_num, elec_num, cord_num, nucl_num
real(c_double) , intent(in)            :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
real(c_double) , intent(in)            :: een_rescaled_e(elec_num, elec_num, 0:cord_num, walk_num)
real(c_double) , intent(in)            :: een_rescaled_single_n(nucl_num, 0:cord_num, walk_num)
real(c_double) , intent(in)            :: een_rescaled_single_e(elec_num, 0:cord_num, walk_num)
real(c_double) , intent(out)           :: delta_p(elec_num, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)

double precision        :: een_rescaled_delta_e(elec_num)

integer*8 :: i, a, c, j, l, k, p, m, n, nw, num
double precision :: dn, dn2
integer*8                           :: LDA, LDB, LDC

num = num_in + 1

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (walk_num <= 0)                 info = QMCKL_INVALID_ARG_3
if (elec_num <= 0)                 info = QMCKL_INVALID_ARG_4
if (nucl_num <= 0)                 info = QMCKL_INVALID_ARG_5
if (cord_num <  0)                 info = QMCKL_INVALID_ARG_6
if (info /= QMCKL_SUCCESS)         return

if (cord_num == 0) return

do nw=1, walk_num

do i=0, cord_num-1

   een_rescaled_delta_e(:) = een_rescaled_single_e(:,i,nw) - een_rescaled_e(:,num,i,nw)

   do c=0,cord_num
      do a=1,nucl_num
         dn = een_rescaled_single_n(a,c,nw) - een_rescaled_n(num,a,c,nw)
         dn2 = een_rescaled_single_n(a,c,nw)
         do j=1,elec_num
            delta_p(j,a,c,i,nw) = een_rescaled_e(j,num,i,nw)*dn + een_rescaled_delta_e(j) * dn2
         enddo
      end do
   end do

   info = qmckl_dgemm(context, 'T', 'N', 1_8, nucl_num * (cord_num+1_8), elec_num, 1.0d0,     &
        een_rescaled_delta_e,elec_num, &
        een_rescaled_n(1,1,0,nw),elec_num,     &
        1.0d0,                                 &
        delta_p(num,1,0,i,nw),elec_num)

enddo

end do

end function qmckl_compute_jastrow_champ_delta_p_doc

Electron-electron-nucleus Jastrow value

Computes $\delta J_{een}$ by combining $\delta \widetilde{R}_{i\alpha}$ and $\delta P$.

\begin{eqnarray*} J_{een} = \sum_{p=2}^{N_\text{nord}} \sum_{k=0}^{p-1} \sum_{l=0}^{p-k-2\delta_{k,0}}\sum_{\alpha=1}^{N_\text{nucl}} c_{l,k,p,\alpha} \left( \sum_{i=1}^{N_\text{elec}} \left( \widetilde{R}_{i,\alpha,(p-k-l)/2} \delta P_{i,\alpha,k,(p-k+l)/2} \right) + \delta \widetilde{R}_{\text{num},\alpha,(p-k-l)/2} \left(P_{\text{num},\alpha,k,(p-k+l)/2} + \delta P_{\text{num},\alpha,k,(p-k+l)/2} \right)\right) \end{eqnarray*}

Get

qmckl_exit_code
qmckl_get_jastrow_champ_single_een(qmckl_context context,
                            double* const delta_een,
                            const int64_t size_max);

interface
  integer(qmckl_exit_code) function qmckl_get_jastrow_champ_single_een (context, &
       delta_een, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (qmckl_context) , intent(in), value :: context
    integer(c_int64_t), intent(in), value       :: size_max
    real(c_double),   intent(out)               :: delta_een(size_max)
  end function
end interface

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
num	int64_t	in	Single point number
walk_num	int64_t	in	Number of walkers
elec_num	int64_t	in	Number of electrons
nucl_num	int64_t	in	Number of nuclei
cord_num	int64_t	in	order of polynomials
dim_c_vector	int64_t	in	dimension of full coefficient vector
c_vector_full	double[dim_c_vector][nucl_num]	in	full coefficient vector
lkpm_combined_index	int64_t[4][dim_c_vector]	in	combined indices
tmp_c	double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]	in	P matrix
delta_p	double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]	in	Delta P matrix
een_rescaled_n	double[walk_num][0:cord_num][nucl_num][elec_num]	in	Electron-nucleus rescaled distances
een_rescaled_e	double[walk_num][0:cord_num][elec_num][elec_num]	in	Electron-electron rescaled distances
een_rescaled_single_n	double[walk_num][0:cord_num][nucl_num]	in	Electron-nucleus single rescaled distances
een_rescaled_single_e	double[walk_num][0:cord_num][elec_num]	in	Electron-electron single rescaled distances
delta_een	double[walk_num]	out	Single electron electron-electron-nucleus Jastrow

integer function qmckl_compute_jastrow_champ_factor_single_een_doc( &
context, num_in, walk_num, elec_num, nucl_num, cord_num,   &
dim_c_vector, c_vector_full, lkpm_combined_index, &
tmp_c, delta_p, een_rescaled_n, een_rescaled_e, een_rescaled_single_n, &
een_rescaled_single_e, delta_een) &
result(info) bind(C)
use, intrinsic :: iso_c_binding
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer(c_int64_t)    , intent(in), value  :: num_in, walk_num, elec_num, cord_num, nucl_num, dim_c_vector
integer(c_int64_t)    , intent(in)  :: lkpm_combined_index(dim_c_vector,4)
real(c_double)        , intent(in)  :: c_vector_full(nucl_num, dim_c_vector)
real(c_double)        , intent(in)  :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
real(c_double)        , intent(in)  :: delta_p(elec_num, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
real(c_double)        , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_e(elec_num, elec_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_n(nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_e(elec_num, 0:cord_num, walk_num)
real(c_double)        , intent(out) :: delta_een(walk_num)

double precision        :: delta_c(nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
double precision        :: delta_c2(elec_num, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)

double precision        :: een_rescaled_delta_n(nucl_num, 0:cord_num)

integer*8 :: i, a, j, l, k, p, m, n, nw, num
double precision :: accu, accu2, cn
integer*8                           :: LDA, LDB, LDC

num = num_in + 1

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (walk_num <= 0)                 info = QMCKL_INVALID_ARG_3
if (elec_num <= 0)                 info = QMCKL_INVALID_ARG_4
if (nucl_num <= 0)                 info = QMCKL_INVALID_ARG_5
if (cord_num <  0)                 info = QMCKL_INVALID_ARG_6
if (info /= QMCKL_SUCCESS)         return

delta_een = 0.0d0

if (cord_num == 0) return

do nw =1, walk_num
een_rescaled_delta_n(:,:) = een_rescaled_single_n(:,:,nw) - een_rescaled_n(num,:,:,nw)
do n = 1, dim_c_vector
   l = lkpm_combined_index(n, 1)
   k = lkpm_combined_index(n, 2)
   p = lkpm_combined_index(n, 3)
   m = lkpm_combined_index(n, 4)

   do a = 1, nucl_num
      cn = c_vector_full(a, n)
      if(cn == 0.d0) cycle

      accu = 0.0d0
      do j = 1, elec_num
         accu = accu + een_rescaled_n(j,a,m,nw) * delta_p(j,a,m+l,k,nw)
      end do
      accu = accu + een_rescaled_delta_n(a,m) * (tmp_c(num,a,m+l,k,nw) + delta_p(num,a,m+l,k,nw))
      delta_een(nw) = delta_een(nw) + accu * cn
   end do
end do
end do

end function qmckl_compute_jastrow_champ_factor_single_een_doc

Electron-nucleus rescaled distance derivative

Calculates $\partial_{i,m} \widetilde{R}_{i,\alpha}$, which is required to compute the derivative of $J_{een}$. \[ \partial_{i,m} \widetilde{R}_{i\alpha} = -\kappa l e^{-\kappa l R_{i\alpha}} \frac{r_{i,m} - R_{\alpha,m}}{R_{i\alpha}} \quad \text{and} \quad \partial_{i,4} \widetilde{R}_{i\alpha} = -\kappa l \left(\frac{2}{R_{i\alpha}}- \kappa l \right) e^{-\kappa l R_{i\alpha}}, \]

Get

qmckl_exit_code
qmckl_get_een_rescaled_single_n_gl(qmckl_context context,
                                    double* const distance_rescaled,
                                    const int64_t size_max);

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
walk_num	int64_t	in	Number of walkers
nucl_num	int64_t	in	Number of atoms
type_nucl_num	int64_t	in	Number of atom types
type_nucl_vector	int64_t[nucl_num]	in	Types of atoms
cord_num	int64_t	in	Order of polynomials
rescale_factor_en	double[nucl_num]	in	Factor to rescale ee distances
coord_ee	double[walk_num][3]	in	Electron coordinates
coord_n	double[3][nucl_num]	in	Nuclear coordinates
single_en_distance	double[walk_num][nucl_num]	in	Electron-nucleus single distances
een_rescaled_single_n	double[walk_num][0:cord_num][nucl_num]	in	Electron-nucleus rescaled single distances
een_rescaled_single_n_gl	double[walk_num][0:cord_num][nucl_num][4]	out	Electron-nucleus rescaled single distances derivative

integer function qmckl_compute_een_rescaled_single_n_gl( &
context, walk_num, nucl_num, type_nucl_num, type_nucl_vector, &
cord_num, rescale_factor_en, coord_ee, coord_n, single_en_distance, &
een_rescaled_single_n, een_rescaled_single_n_gl) &
result(info) bind(C)
use, intrinsic :: iso_c_binding
use qmckl
implicit none
integer(qmckl_context), intent(in), value  :: context
integer(c_int64_t)    , intent(in), value  :: walk_num
integer(c_int64_t)    , intent(in), value  :: nucl_num
integer(c_int64_t)    , intent(in), value  :: type_nucl_num
integer(c_int64_t)    , intent(in)         :: type_nucl_vector(nucl_num)
integer(c_int64_t)    , intent(in), value  :: cord_num
real(c_double)        , intent(in)         :: rescale_factor_en(type_nucl_num)
real(c_double)        , intent(in)         :: coord_ee(3,walk_num)
real(c_double)        , intent(in)         :: coord_n(nucl_num,3)
real(c_double)        , intent(in)         :: single_en_distance(nucl_num,walk_num)
real(c_double)        , intent(in)         :: een_rescaled_single_n(nucl_num,0:cord_num,walk_num)
real(c_double)        , intent(out)        :: een_rescaled_single_n_gl(4,nucl_num,0:cord_num,walk_num)

double precision,allocatable   :: elnuc_dist_gl(:,:)
double precision               :: x, ria_inv, kappa_l
integer*8                      :: i, a, k, l, nw, ii

allocate(elnuc_dist_gl(4, nucl_num))

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif

if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif


if (nucl_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif

if (cord_num < 0) then
info = QMCKL_INVALID_ARG_4
return
endif

! Prepare table of exponentiated distances raised to appropriate power
een_rescaled_single_n_gl             = 0.0d0
do nw = 1, walk_num
! prepare the actual een table
do a = 1, nucl_num
   ria_inv = 1.0d0 / single_en_distance(a, nw)
   do ii = 1, 3
      elnuc_dist_gl(ii, a) = (coord_ee(ii,nw) - coord_n(a, ii)) * ria_inv
   end do
   elnuc_dist_gl(4, a) = 2.0d0 * ria_inv
end do

do l = 0, cord_num
   do a = 1, nucl_num
      kappa_l = - dble(l) * rescale_factor_en(type_nucl_vector(a)+1)
      een_rescaled_single_n_gl(1, a, l, nw) = kappa_l * elnuc_dist_gl(1, a)
      een_rescaled_single_n_gl(2, a, l, nw) = kappa_l * elnuc_dist_gl(2, a)
      een_rescaled_single_n_gl(3, a, l, nw) = kappa_l * elnuc_dist_gl(3, a)
      een_rescaled_single_n_gl(4, a, l, nw) = kappa_l * (elnuc_dist_gl(4, a) + kappa_l)

      een_rescaled_single_n_gl(1, a, l, nw) = een_rescaled_single_n_gl(1, a, l, nw) * &
           een_rescaled_single_n(a, l, nw)
      een_rescaled_single_n_gl(2, a, l, nw) = een_rescaled_single_n_gl(2, a, l, nw) * &
           een_rescaled_single_n(a, l, nw)
      een_rescaled_single_n_gl(3, a, l, nw) = een_rescaled_single_n_gl(3, a, l, nw) * &
           een_rescaled_single_n(a, l, nw)
      een_rescaled_single_n_gl(4, a, l, nw) = een_rescaled_single_n_gl(4, a, l, nw) * &
           een_rescaled_single_n(a, l, nw)
   end do
end do
end do

end function qmckl_compute_een_rescaled_single_n_gl

Electron-electron rescaled distances derivative

Calculates the updated terms of $\partial_{i,m} \widetilde{r}_{ij}$, required for computing the derivative of $J_{een}$. \[ \partial_{i,m} \widetilde{r}_{ij} = -\kappa l e^{-\kappa l r_{ij}} \frac{r_{i,m} - r_{j,m}}{r_{ij}} \quad \text{and} \quad \partial_{i,4} \widetilde{r}_{ij} = -\kappa l \left(\frac{2}{r_{ij}}- \kappa l \right) e^{-\kappa l r_{ij}}. \]

Get

qmckl_exit_code
qmckl_get_een_rescaled_single_e_gl(qmckl_context context,
                                    double* const distance_rescaled,
                                    const int64_t size_max);

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
num	int64_t	in	Index of single electron
walk_num	int64_t	in	Number of walkers
elec_num	int64_t	in	Number of electrons
cord_num	int64_t	in	Order of polynomials
rescale_factor_ee	double	in	Factor to rescale ee distances
coord	double[walk_num][3]	in	Single electron coordinates
coord_ee	double[3][walk_num][elec_num]	in	Electron coordinates
single_ee_distance	double[walk_num][elec_num]	in	Electron-electron single distances
een_rescaled_single_e	double[walk_num][0:cord_num][elec_num]	in	Electron-electron rescaled single distances
een_rescaled_single_e_gl	double[walk_num][0:cord_num][elec_num][4]	out	Electron-electron rescaled single distances derivative

integer function qmckl_compute_een_rescaled_single_e_gl_doc( &
context, num_in, walk_num, elec_num, cord_num, rescale_factor_ee,  &
coord, coord_ee, single_ee_distance, een_rescaled_single_e, een_rescaled_single_e_gl) &
result(info) bind(C)
use, intrinsic :: iso_c_binding
use qmckl
implicit none
integer(qmckl_context), intent(in)        :: context
integer(c_int64_t)    , intent(in), value :: num_in
integer(c_int64_t)    , intent(in), value :: walk_num
integer(c_int64_t)    , intent(in), value :: elec_num
integer(c_int64_t)    , intent(in), value :: cord_num
real(c_double)        , intent(in), value :: rescale_factor_ee
real(c_double)        , intent(in)        :: coord(3,walk_num)
real(c_double)        , intent(in)        :: coord_ee(elec_num,walk_num,3)
real(c_double)        , intent(in)        :: single_ee_distance(elec_num,walk_num)
real(c_double)        , intent(in)        :: een_rescaled_single_e(elec_num,0:cord_num,walk_num)
real(c_double)        , intent(out)       :: een_rescaled_single_e_gl(4,elec_num,0:cord_num,walk_num)

double precision,allocatable   :: elec_dist_gl(:,:)
double precision               :: x, rij_inv, kappa_l
integer*8                      :: i, j, k, l, nw, ii, num

num = num_in + 1

allocate(elec_dist_gl(4, elec_num))

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif

if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif

if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif

if (cord_num < 0) then
info = QMCKL_INVALID_ARG_5
return
endif

!   Not necessary: should be set to zero by qmckl_malloc
!   een_rescaled_single_e_gl = 0.0d0

! Prepare table of exponentiated distances raised to appropriate power

do nw = 1, walk_num
do i = 1, elec_num
   if (i == num) cycle
   rij_inv = 1.0d0 / single_ee_distance(i, nw)
   do ii = 1, 3
      elec_dist_gl(ii, i) = (coord(ii, nw) - coord_ee(i, nw, ii)) * rij_inv
   end do
   elec_dist_gl(4, i) = 2.0d0 * rij_inv
end do

elec_dist_gl(:, num) = 0.0d0

do l = 1, cord_num
   kappa_l = - dble(l) * rescale_factor_ee
   do i = 1, elec_num
      een_rescaled_single_e_gl(1, i, l, nw) = kappa_l * elec_dist_gl(1, i)
      een_rescaled_single_e_gl(2, i, l, nw) = kappa_l * elec_dist_gl(2, i)
      een_rescaled_single_e_gl(3, i, l, nw) = kappa_l * elec_dist_gl(3, i)
      een_rescaled_single_e_gl(4, i, l, nw) = kappa_l * (elec_dist_gl(4, i) + kappa_l)

      een_rescaled_single_e_gl(1,i,l,nw) = een_rescaled_single_e_gl(1,i,l,nw) * een_rescaled_single_e(i,l,nw)
      een_rescaled_single_e_gl(2,i,l,nw) = een_rescaled_single_e_gl(2,i,l,nw) * een_rescaled_single_e(i,l,nw)
      een_rescaled_single_e_gl(3,i,l,nw) = een_rescaled_single_e_gl(3,i,l,nw) * een_rescaled_single_e(i,l,nw)
      een_rescaled_single_e_gl(4,i,l,nw) = een_rescaled_single_e_gl(4,i,l,nw) * een_rescaled_single_e(i,l,nw)

   end do
end do
end do



end function qmckl_compute_een_rescaled_single_e_gl_doc

$\delta P$ matrix gradients and Laplacian

Calculate $\partial_{i,m} \delta P$.

\begin{eqnarray*} \partial_{i,m} \delta P_{i,\alpha,k,l} = \partial_{i,m}\delta \widetilde{r}_{\text{num},i,k} \left(\partial_{i,m}\delta \widetilde{R}_{\text{num},\alpha,l} + \partial_{i,m}\widetilde{R}_{\text{num},\alpha,l} \right) g_m + \partial_{i,m}\widetilde{r}_{\text{num},i,k} \delta \widetilde{R}_{\text{num},\alpha,l} + \delta_{i,\text{num}} \sum_{j=1}^{N_\text{elec}} \left( \partial_{j,m} \delta \widetilde{r}_{\text{num},j,k} \widetilde{R}_{j,\alpha,l} \right), \end{eqnarray*}

where $g_m = \{-1,-1,-1,1\}$.

Get

qmckl_exit_code
qmckl_get_jastrow_champ_delta_p_gl(qmckl_context context,
                            double* const delta_p_gl,
                            const int64_t size_max);

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
num	int64_t	in	Index of single electron
walk_num	int64_t	in	Number of walkers
elec_num	int64_t	in	Number of electrons
nucl_num	int64_t	in	Number of nuclei
cord_num	int64_t	in	order of polynomials
een_rescaled_n	double[walk_num][0:cord_num][nucl_num][elec_num]	in	Electron-nucleus rescaled distances
een_rescaled_e	double[walk_num][0:cord_num][elec_num][elec_num]	in	Electron-electron rescaled distances
een_rescaled_single_n	double[walk_num][0:cord_num][nucl_num]	in	Electron-nucleus single rescaled distances
een_rescaled_single_e	double[walk_num][0:cord_num][elec_num]	in	Electron-electron single rescaled distances
een_rescaled_n_gl	double[walk_num][0:cord_num][nucl_num][4][elec_num]	in	Electron-nucleus rescaled distances derivatives
een_rescaled_e_gl	double[walk_num][0:cord_num][elec_num][4][elec_num]	in	Electron-electron rescaled distances derivatives
een_rescaled_single_n_gl	double[walk_num][0:cord_num][nucl_num][4]	in	Electron-nucleus single rescaled distances derivatives
een_rescaled_single_e_gl	double[walk_num][0:cord_num][elec_num][4]	in	Electron-electron single rescaled distances derivatives
delta_p_gl	double[walk_num][0:cord_num-1][0:cord_num][4][nucl_num][elec_num]	out	Delta P matrix gradient and Laplacian

integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_delta_p_gl_doc( &
context, num_in, walk_num, elec_num, nucl_num, cord_num,   &
een_rescaled_n, een_rescaled_e, een_rescaled_single_n, een_rescaled_single_e, &
een_rescaled_n_gl, een_rescaled_e_gl, een_rescaled_single_n_gl, een_rescaled_single_e_gl, delta_p_gl) &
result(info) bind(C)
use, intrinsic :: iso_c_binding
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer(c_int64_t)    , intent(in), value  :: num_in, walk_num, elec_num, cord_num, nucl_num
real(c_double)        , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_e(elec_num, elec_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_n(nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_e(elec_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_e_gl(elec_num, 4, elec_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_n_gl(4, nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_e_gl(4,elec_num, 0:cord_num, walk_num)
real(c_double)        , intent(out)  :: delta_p_gl(elec_num,nucl_num,4,0:cord_num, 0:cord_num-1,  walk_num)

double precision        :: delta_e_gl(elec_num,4)

double precision        :: een_rescaled_delta_n, een_re_n, een_re_single_n

integer*8 :: i, a, j, l, k, p, m, n, nw, num
double precision :: tmp, cummu
integer*8        :: LDA, LDB, LDC

num = num_in + 1

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (walk_num <= 0)                 info = QMCKL_INVALID_ARG_3
if (elec_num <= 0)                 info = QMCKL_INVALID_ARG_4
if (nucl_num <= 0)                 info = QMCKL_INVALID_ARG_5
if (cord_num <  0)                 info = QMCKL_INVALID_ARG_6
if (info /= QMCKL_SUCCESS)         return


if (cord_num == 0) then
delta_p_gl = 0.d0
return
endif

do nw=1, walk_num
do m=1, cord_num-1
   do k = 1, 4
     do j = 1, elec_num
       delta_e_gl(j,k) = een_rescaled_single_e_gl(k,j,m,nw) - een_rescaled_e_gl(num, k, j, m, nw)
     end do
   end do
   do k = 1, 4
     delta_e_gl(num, k) = 0.0d0
   end do

   do l=0, cord_num
     do k = 1, 3
       do a = 1, nucl_num

         een_re_n = een_rescaled_n(num, a, l, nw)
         een_re_single_n = een_rescaled_single_n(a,l,nw)

         cummu = 0.0d0
         do i = 1, elec_num

           delta_p_gl(i,a,k,l,m,nw) = -een_rescaled_e_gl(i,k,num,m,nw) * een_re_n &
             - een_rescaled_single_e_gl(k,i,m,nw) * een_re_single_n

           cummu = cummu + delta_e_gl(i,k) * een_rescaled_n(i,a,l,nw)
         end do
         delta_p_gl(num,a,k,l,m,nw) = delta_p_gl(num,a,k,l,m,nw) + cummu

       end do
     end do
     do a = 1, nucl_num
       een_rescaled_delta_n = een_rescaled_single_n(a,l,nw) - een_rescaled_n(num, a, l, nw)
       cummu = 0.0d0
       een_re_single_n = een_rescaled_single_n(a,l,nw)
       
       do i = 1, elec_num
         delta_p_gl(i,a,4,l,m,nw) = een_rescaled_e_gl(i,4,num,m,nw) * een_rescaled_delta_n &
           +delta_e_gl(i,4) * een_re_single_n

         cummu = cummu + delta_e_gl(i,4) * een_rescaled_n(i,a,l,nw)
       end do
       delta_p_gl(num,a,4,l,m,nw) = delta_p_gl(num,a,4,l,m,nw) + cummu

     end do
   end do
end do
end do

end function qmckl_compute_jastrow_champ_delta_p_gl_doc

Electron-electron-nucleus Jastrow gradients and Laplacian

Calculates the gradients and Laplacian of $\delta J_{een}$.

\begin{eqnarray*} \partial_{i,m} J_{een} &=& \sum_{p=2}^{N_\text{nord}} \sum_{k=0}^{p-1} \sum_{l=0}^{p-k-2\delta_{k,0}}\sum_{\alpha=1}^{N_\text{nucl}} c_{l,k,p,\alpha} \left( \widetilde{R}_{i,\alpha,(p-k-l)/2} \partial_{i,m} \delta P_{i,\alpha,k,(p-k+l)/2} + \widetilde{R}_{i,\alpha,(p-k+l)/2} \partial_{i,m} \delta P_{i,\alpha,k,(p-k-l)/2} \right. \\ & &\ + \partial_{i,m}\widetilde{R}_{i,\alpha,(p-k-l)/2} \delta P_{i,\alpha,k,(p-k+l)/2} + \partial_{i,m}\widetilde{R}_{i,\alpha,(p-k+l)/2} \delta P_{i,\alpha,k,(p-k-l)/2} \\ & &\ + \delta_{i,\text{num}} \left( \delta \widetilde{R}_{i,\alpha,(p-k-l)/2} \left ( \partial_{i,m} P_{i,\alpha,k,(p-k+l)/2} + \partial_{i,m} \delta P_{i,\alpha,k,(p-k+l)/2} \right) + \delta \widetilde{R}_{i,\alpha,(p-k+l)/2} \left ( \partial_{i,m} P_{i,\alpha,k,(p-k-l)/2} + \partial_{i,m} \delta P_{i,\alpha,k,(p-k-l)/2} \right) \right. \\ & &\ \left. + \partial_{i,m} \delta \widetilde{R}_{i,\alpha,(p-k-l)/2} \left ( P_{i,\alpha,k,(p-k+l)/2} + \delta P_{i,\alpha,k,(p-k+l)/2} \right) + \partial_{i,m} \delta \widetilde{R}_{i,\alpha,(p-k+l)/2} \left ( P_{i,\alpha,k,(p-k-l)/2} + \delta P_{i,\alpha,k,(p-k-l)/2} \right) \right)\\ & &\ + \delta_{m,4} \sum_{d=1}^3 \left( \partial_{i,d} \widetilde{R}_{i,\alpha,(p-k-l)/2} \partial_{i,d} \delta P_{i,\alpha,k,(p-k+l)/2} + \partial_{i,d} \widetilde{R}_{i,\alpha,(p-k+l)/2} \partial_{i,d} \delta P_{i,\alpha,k,(p-k-l)/2} \right)\\ & &\ \left. + \delta_{m,4}\delta_{i,\text{num}} \sum_{d=1}^3\left( \partial_{i,d}\delta \widetilde{R}_{i,\alpha,(p-k-l)/2} \left( \partial_{i,d} P_{i,\alpha,k,(p-k+l)/2} + \partial_{i,d}\delta P_{i,\alpha,k,(p-k+l)/2} \right) + \partial_{i,d}\delta \widetilde{R}_{i,\alpha,(p-k+l)/2} \left( \partial_{i,d} P_{i,\alpha,k,(p-k-l)/2} + \partial_{i,d} \delta P_{i,\alpha,k,(p-k-l)/2} \right) \right) \right) \end{eqnarray*}

Get

qmckl_exit_code
qmckl_get_jastrow_champ_single_een_gl(qmckl_context context,
                            double* const delta_een_gl,
                            const int64_t size_max);

interface
  integer(qmckl_exit_code) function qmckl_get_jastrow_champ_single_een_gl (context, &
       delta_een_gl, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (qmckl_context) , intent(in), value :: context
    integer(c_int64_t), intent(in), value       :: size_max
    real(c_double),   intent(out)               :: delta_een_gl(size_max)
  end function
end interface

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
num	int64_t	in	Index of single electron
walk_num	int64_t	in	Number of walkers
elec_num	int64_t	in	Number of electrons
nucl_num	int64_t	in	Number of nuclei
cord_num	int64_t	in	order of polynomials
dim_c_vector	int64_t	in	dimension of full coefficient vector
c_vector_full	double[dim_c_vector][nucl_num]	in	full coefficient vector
lkpm_combined_index	int64_t[4][dim_c_vector]	in	combined indices
tmp_c	double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]	in	P matrix
dtmp_c	double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]	in	P matrix derivative
delta_p	double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]	in	Delta P matrix
delta_p_gl	double[walk_num][0:cord_num-1][0:cord_num][4][nucl_num][elec_num]	in	Delta P matrix derivative
een_rescaled_n	double[walk_num][0:cord_num][nucl_num][elec_num]	in	Electron-nucleus rescaled distances
een_rescaled_single_n	double[walk_num][0:cord_num][nucl_num]	in	Electron-nucleus single rescaled distances
een_rescaled_n_gl	double[walk_num][0:cord_num][nucl_num][4][elec_num]	in	Electron-nucleus rescaled distances derivatives
een_rescaled_single_n_gl	double[walk_num][0:cord_num][nucl_num][4]	in	Electron-nucleus single rescaled distances derivatives
delta_een_gl	double[walk_num][4][elec_num]	out	Delta electron-electron-nucleus jastrow gradient and Laplacian

integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_factor_single_een_gl_doc( &
context, num_in, walk_num, elec_num, nucl_num, cord_num,   &
dim_c_vector, c_vector_full, lkpm_combined_index, &
tmp_c, dtmp_c, delta_p, delta_p_gl, een_rescaled_n, een_rescaled_single_n, &
een_rescaled_n_gl, een_rescaled_single_n_gl, delta_een_gl) &
result(info) bind(C)
use, intrinsic :: iso_c_binding
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer(c_int64_t)    , intent(in), value  :: num_in, walk_num, elec_num, cord_num, nucl_num, dim_c_vector
integer(c_int64_t)    , intent(in)  :: lkpm_combined_index(dim_c_vector,4)
real(c_double)        , intent(in)  :: c_vector_full(nucl_num, dim_c_vector)
real(c_double)        , intent(in)  :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
real(c_double)        , intent(in)  :: dtmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
real(c_double)        , intent(in)  :: delta_p(elec_num, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
real(c_double)        , intent(in)  :: delta_p_gl(elec_num, nucl_num, 4, 0:cord_num, 0:cord_num-1,  walk_num)
real(c_double)        , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_n(nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_n_gl(4, nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(out) :: delta_een_gl(elec_num, 4, walk_num)

integer*8 :: i, a, j, l, k, p, m, n, nw, kk, num
double precision :: accu, accu2, cn
integer*8                           :: LDA, LDB, LDC

double precision  :: een_rescaled_delta_n(nucl_num, 0:cord_num)
double precision  :: een_rescaled_delta_n_gl(4, nucl_num, 0:cord_num)
double precision :: dpg1_m, dpg1_ml, dp_m, dp_ml, een_r_m, een_r_ml, een_r_gl_m, een_r_gl_ml
num = num_in + 1

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (walk_num <= 0)                 info = QMCKL_INVALID_ARG_3
if (elec_num <= 0)                 info = QMCKL_INVALID_ARG_4
if (nucl_num <= 0)                 info = QMCKL_INVALID_ARG_5
if (cord_num <  0)                 info = QMCKL_INVALID_ARG_6
if (info /= QMCKL_SUCCESS)         return

delta_een_gl = 0.0d0

if (cord_num == 0) return


do nw =1, walk_num

een_rescaled_delta_n(:,:) = een_rescaled_single_n(:,:,nw) - een_rescaled_n(num, :, :, nw)
een_rescaled_delta_n_gl(:,:,:) = een_rescaled_single_n_gl(:,:,:,nw) - een_rescaled_n_gl(num, :,:,:,nw)

do n = 1, dim_c_vector
   l = lkpm_combined_index(n, 1)
   k = lkpm_combined_index(n, 2)
   p = lkpm_combined_index(n, 3)
   m = lkpm_combined_index(n, 4)

   do kk = 1, 4
      do a = 1, nucl_num
         cn = c_vector_full(a, n)
         if(cn == 0.d0) cycle
         do i = 1, elec_num
            delta_een_gl(i,kk,nw) = delta_een_gl(i,kk,nw) + ( &
                 delta_p_gl(i,a,kk,m  ,k,nw) * een_rescaled_n(i,a,m+l,nw) + &
                 delta_p_gl(i,a,kk,m+l,k,nw) * een_rescaled_n(i,a,m  ,nw) + &
                 delta_p(i,a,m  ,k,nw) * een_rescaled_n_gl(i,kk,a,m+l,nw) + &
                 delta_p(i,a,m+l,k,nw) * een_rescaled_n_gl(i,kk,a,m  ,nw) ) * cn
         end do

         delta_een_gl(num,kk,nw) = delta_een_gl(num,kk,nw) + ( &
              (dtmp_c(num,kk,a,m  ,k,nw) + delta_p_gl(num,a,kk,m  ,k,nw)) * een_rescaled_delta_n(a,m+l) + &
              (dtmp_c(num,kk,a,m+l,k,nw) + delta_p_gl(num,a,kk,m+l,k,nw)) * een_rescaled_delta_n(a,m  ) + &
              (tmp_c(num,a,m  ,k,nw) + delta_p(num,a,m  ,k,nw)) * een_rescaled_delta_n_gl(kk,a,m+l)  + &
              (tmp_c(num,a,m+l,k,nw) + delta_p(num,a,m+l,k,nw)) * een_rescaled_delta_n_gl(kk,a,m  ) )* cn
      end do
   end do
   do a = 1, nucl_num
      cn = c_vector_full(a, n)
      if(cn == 0.d0) cycle
      cn = cn + cn
      do i = 1, elec_num
         delta_een_gl(i,4,nw) = delta_een_gl(i,4,nw) + ( &
              delta_p_gl(i,a,1,m  ,k,nw) * een_rescaled_n_gl(i,1,a,m+l,nw) + &
              delta_p_gl(i,a,1,m+l,k,nw) * een_rescaled_n_gl(i,1,a,m  ,nw) + &
              delta_p_gl(i,a,2,m  ,k,nw) * een_rescaled_n_gl(i,2,a,m+l,nw) + &
              delta_p_gl(i,a,2,m+l,k,nw) * een_rescaled_n_gl(i,2,a,m  ,nw) + &
              delta_p_gl(i,a,3,m  ,k,nw) * een_rescaled_n_gl(i,3,a,m+l,nw) + &
              delta_p_gl(i,a,3,m+l,k,nw) * een_rescaled_n_gl(i,3,a,m  ,nw) ) * cn
      end do
     delta_een_gl(num,4,nw) = delta_een_gl(num,4,nw) + ( &
          (delta_p_gl(num,a,1,m  ,k,nw) + dtmp_c(num,1,a,m  ,k,nw)) * een_rescaled_delta_n_gl(1,a,m+l) + &
          (delta_p_gl(num,a,1,m+l,k,nw) + dtmp_c(num,1,a,m+l,k,nw)) * een_rescaled_delta_n_gl(1,a,m  ) + &
          (delta_p_gl(num,a,2,m  ,k,nw) + dtmp_c(num,2,a,m  ,k,nw)) * een_rescaled_delta_n_gl(2,a,m+l) + &
          (delta_p_gl(num,a,2,m+l,k,nw) + dtmp_c(num,2,a,m+l,k,nw)) * een_rescaled_delta_n_gl(2,a,m  ) + &
          (delta_p_gl(num,a,3,m  ,k,nw) + dtmp_c(num,3,a,m  ,k,nw)) * een_rescaled_delta_n_gl(3,a,m+l) + &
          (delta_p_gl(num,a,3,m+l,k,nw) + dtmp_c(num,3,a,m+l,k,nw)) * een_rescaled_delta_n_gl(3,a,m  ) ) * cn
   end do
end do
end do

end function qmckl_compute_jastrow_champ_factor_single_een_gl_doc

qmckl_exit_code
qmckl_compute_jastrow_champ_factor_single_een_gl_hpc (const qmckl_context context,
                                                 const int64_t num,
                                                 const int64_t walk_num,
                                                 const int64_t elec_num,
                                                 const int64_t nucl_num,
                                                 const int64_t cord_num,
                                                 const int64_t dim_c_vector,
                                                 const double* restrict c_vector_full,
                                                 const int64_t* restrict lkpm_combined_index,
                                                 const double* restrict tmp_c,
                                                 const double* restrict dtmp_c,
                                                 const double* restrict delta_p,
                                                 const double* restrict delta_p_gl,
                                                 const double* restrict een_rescaled_n,
                                                 const double* restrict een_rescaled_single_n,
                                                 const double* restrict een_rescaled_n_gl,
                                                 const double* restrict een_rescaled_single_n_gl,
                                                 double* restrict const delta_een_gl )
{


if (context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT;
if (walk_num <= 0)                 return QMCKL_INVALID_ARG_3;
if (elec_num <= 0)                 return QMCKL_INVALID_ARG_4;
if (nucl_num <= 0)                 return QMCKL_INVALID_ARG_5;
if (cord_num <  0)                 return QMCKL_INVALID_ARG_6;

if (cord_num == 0) {
#pragma omp parallel for
for (size_t i=0 ; i<walk_num*4*elec_num ; ++i) {
 delta_een_gl[i] = 0.;
}
return QMCKL_SUCCESS;
}

#ifdef HAVE_OPENMP
#pragma omp parallel for
#endif
for (int64_t nw=0 ; nw<walk_num ; nw++) {
for (int64_t i=0 ; i<4*elec_num ; ++i) {
 delta_een_gl[i+nw*4*elec_num] = 0.;
}
double een_rescaled_delta_n[nucl_num*(cord_num+1)];
double een_rescaled_delta_n_gl[4*nucl_num*(cord_num+1)];
for (int64_t k=0 ; k<nucl_num*(cord_num+1) ; ++k) {
 een_rescaled_delta_n[k] = een_rescaled_single_n[k+nucl_num*(cord_num+1)*nw] - een_rescaled_n[num+elec_num*(k+nucl_num*(cord_num+1)*nw)];
 een_rescaled_delta_n_gl[0+4*k] = een_rescaled_single_n_gl[0+4*(k+nucl_num*(cord_num+1)*nw)] - een_rescaled_n_gl[num+elec_num*(0+4*(k+nucl_num*(cord_num+1)*nw))];
 een_rescaled_delta_n_gl[1+4*k] = een_rescaled_single_n_gl[1+4*(k+nucl_num*(cord_num+1)*nw)] - een_rescaled_n_gl[num+elec_num*(1+4*(k+nucl_num*(cord_num+1)*nw))];
 een_rescaled_delta_n_gl[2+4*k] = een_rescaled_single_n_gl[2+4*(k+nucl_num*(cord_num+1)*nw)] - een_rescaled_n_gl[num+elec_num*(2+4*(k+nucl_num*(cord_num+1)*nw))];
 een_rescaled_delta_n_gl[3+4*k] = een_rescaled_single_n_gl[3+4*(k+nucl_num*(cord_num+1)*nw)] - een_rescaled_n_gl[num+elec_num*(3+4*(k+nucl_num*(cord_num+1)*nw))];
}

for (int64_t n=0 ; n<dim_c_vector ; ++n) {
 const int64_t l = lkpm_combined_index[n];
 const int64_t k = lkpm_combined_index[n+dim_c_vector];
 const int64_t m = lkpm_combined_index[n+3*dim_c_vector];

 for (int64_t kk=0 ; kk<4 ; ++kk) {
   double* restrict dgl = &delta_een_gl[elec_num*(kk+4*nw)];

   for (int64_t a=0 ; a<nucl_num ; ++a) {
     const double cn = c_vector_full[a+n*nucl_num];
     if (cn == 0.) continue;

     const double* restrict dpg1_m    = &delta_p_gl[elec_num*(a+nucl_num*(kk+4*(m+(cord_num+1)*(k+cord_num*nw))))];
     const double* restrict dpg1_ml   = &delta_p_gl[elec_num*(a+nucl_num*(kk+4*(m+l+(cord_num+1)*(k+cord_num*nw))))];
     const double* restrict dp_m      = &delta_p[elec_num*(a+nucl_num*(m+(cord_num+1)*(k+cord_num*nw)))];
     const double* restrict dp_ml     = &delta_p[elec_num*(a+nucl_num*(m+l+(cord_num+1)*(k+cord_num*nw)))];
     const double* restrict een_r_m   = &een_rescaled_n[elec_num*(a+nucl_num*(m+(cord_num+1)*nw))];
     const double* restrict een_r_ml  = &een_rescaled_n[elec_num*(a+nucl_num*(m+l+(cord_num+1)*nw))];
     const double* restrict een_r_gl_m  = &een_rescaled_n_gl[elec_num*(kk+4*(a+nucl_num*(m+(cord_num+1)*nw)))];
     const double* restrict een_r_gl_ml = &een_rescaled_n_gl[elec_num*(kk+4*(a+nucl_num*(m+l+(cord_num+1)*nw)))];

#ifdef HAVE_OPENMP
     #pragma omp simd
#endif
     for (int64_t i=0 ; i<elec_num ; ++i) {
       dgl[i] += cn * (dpg1_m[i] * een_r_ml[i] + dpg1_ml[i] * een_r_m[i] +
                       dp_m[i] * een_r_gl_ml[i] + dp_ml[i] * een_r_gl_m[i]);
     }
     dgl[num] += ( (dtmp_c[num+elec_num*(kk+4*(a+nucl_num*(m+(cord_num+1)*(k+cord_num*nw))))] +
                    dpg1_m[num]
                    ) * een_rescaled_delta_n[a+nucl_num*(m+l)] +
                   (dtmp_c[num+elec_num*(kk+4*(a+nucl_num*(m+l+(cord_num+1)*(k+cord_num*nw))))] +
                    dpg1_ml[num]
                    ) * een_rescaled_delta_n[a+nucl_num*m] +
                   (tmp_c[num+elec_num*(a+nucl_num*(m+(cord_num+1)*(k+cord_num*nw)))] +
                    dp_m[num]
                    ) * een_rescaled_delta_n_gl[kk+4*(a+nucl_num*(m+l))] +
                   (tmp_c[num+elec_num*(a+nucl_num*(m+l+(cord_num+1)*(k+cord_num*nw)))] +
                    dp_ml[num]
                    ) * een_rescaled_delta_n_gl[kk+4*(a+nucl_num*m)] 
                   ) * cn;
   }
 }

 for (int a=0 ; a<nucl_num ; ++a) {
   const double cn = 2. * c_vector_full[a+n*nucl_num];
   if (cn == 0.) continue;
   
   double* restrict dgl4 = &delta_een_gl[elec_num*(3+4*nw)];
   
   const double* restrict dpg1_m    = &delta_p_gl[elec_num*(a+nucl_num*(0+4*(m+(cord_num+1)*(k+cord_num*nw))))];
   const double* restrict dpg2_m    = &delta_p_gl[elec_num*(a+nucl_num*(1+4*(m+(cord_num+1)*(k+cord_num*nw))))];
   const double* restrict dpg3_m    = &delta_p_gl[elec_num*(a+nucl_num*(2+4*(m+(cord_num+1)*(k+cord_num*nw))))];
   const double* restrict dpg1_ml   = &delta_p_gl[elec_num*(a+nucl_num*(0+4*(m+l+(cord_num+1)*(k+cord_num*nw))))];
   const double* restrict dpg2_ml   = &delta_p_gl[elec_num*(a+nucl_num*(1+4*(m+l+(cord_num+1)*(k+cord_num*nw))))];
   const double* restrict dpg3_ml   = &delta_p_gl[elec_num*(a+nucl_num*(2+4*(m+l+(cord_num+1)*(k+cord_num*nw))))];
   const double* restrict een_r_gl1_m  = &een_rescaled_n_gl[elec_num*(0+4*(a+nucl_num*(m+(cord_num+1)*nw)))];
   const double* restrict een_r_gl2_m  = &een_rescaled_n_gl[elec_num*(1+4*(a+nucl_num*(m+(cord_num+1)*nw)))];
   const double* restrict een_r_gl3_m  = &een_rescaled_n_gl[elec_num*(2+4*(a+nucl_num*(m+(cord_num+1)*nw)))];
   const double* restrict een_r_gl1_ml = &een_rescaled_n_gl[elec_num*(0+4*(a+nucl_num*(m+l+(cord_num+1)*nw)))];
   const double* restrict een_r_gl2_ml = &een_rescaled_n_gl[elec_num*(1+4*(a+nucl_num*(m+l+(cord_num+1)*nw)))];
   const double* restrict een_r_gl3_ml = &een_rescaled_n_gl[elec_num*(2+4*(a+nucl_num*(m+l+(cord_num+1)*nw)))];

#ifdef HAVE_OPENMP
     #pragma omp simd
#endif
   for (int64_t i=0 ; i<elec_num ; ++i) {
     dgl4[i] += (dpg1_m[i] * een_r_gl1_ml[i] + dpg1_ml[i] * een_r_gl1_m[i] +
                 dpg2_m[i] * een_r_gl2_ml[i] + dpg2_ml[i] * een_r_gl2_m[i] +
                 dpg3_m[i] * een_r_gl3_ml[i] + dpg3_ml[i] * een_r_gl3_m[i] ) * cn;
   }
   dgl4[num] += ((dpg1_m[num] +
                  dtmp_c[num+elec_num*(0+4*(a+nucl_num*(m+(cord_num+1)*(k+cord_num*nw))))]
                  ) * een_rescaled_delta_n_gl[0+4*(a+nucl_num*(m+l))] +
                 (dpg1_ml[num] +
                  dtmp_c[num+elec_num*(0+4*(a+nucl_num*(m+l+(cord_num+1)*(k+cord_num*nw))))]
                  ) * een_rescaled_delta_n_gl[0+4*(a+nucl_num*m)] +
                 (dpg2_m[num] +
                  dtmp_c[num+elec_num*(1+4*(a+nucl_num*(m+(cord_num+1)*(k+cord_num*nw))))]
                  ) * een_rescaled_delta_n_gl[1+4*(a+nucl_num*(m+l))] +
                 (dpg2_ml[num] +
                  dtmp_c[num+elec_num*(1+4*(a+nucl_num*(m+l+(cord_num+1)*(k+cord_num*nw))))]
                  ) * een_rescaled_delta_n_gl[1+4*(a+nucl_num*m)] +
                 (dpg3_m[num] +
                  dtmp_c[num+elec_num*(2+4*(a+nucl_num*(m+(cord_num+1)*(k+cord_num*nw))))]
                  ) * een_rescaled_delta_n_gl[2+4*(a+nucl_num*(m+l))] +
                 (dpg3_ml[num] +
                  dtmp_c[num+elec_num*(2+4*(a+nucl_num*(m+l+(cord_num+1)*(k+cord_num*nw))))]
                  ) * een_rescaled_delta_n_gl[2+4*(a+nucl_num*m)] ) * cn;
   
 }
}
}

return QMCKL_SUCCESS;
}

$\delta P$ matrix gradient only

To allow for more efficient treatment of nonlocal pseudopotentials, we also include a routine for only calculating the gradient (and not the Laplacian) of $\delta J_{een}$. For this, we calculate the gradient of $\delta P$.

Get

qmckl_exit_code
qmckl_get_jastrow_champ_delta_p_g(qmckl_context context,
                            double* const delta_p_g,
                            const int64_t size_max);

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
num	int64_t	in	Index of single electron
walk_num	int64_t	in	Number of walkers
elec_num	int64_t	in	Number of electrons
nucl_num	int64_t	in	Number of nuclei
cord_num	int64_t	in	order of polynomials
een_rescaled_n	double[walk_num][0:cord_num][nucl_num][elec_num]	in	Electron-nucleus rescaled distances
een_rescaled_e	double[walk_num][0:cord_num][elec_num][elec_num]	in	Electron-electron rescaled distances
een_rescaled_single_n	double[walk_num][0:cord_num][nucl_num]	in	Electron-nucleus single rescaled distances
een_rescaled_single_e	double[walk_num][0:cord_num][elec_num]	in	Electron-electron single rescaled distances
een_rescaled_n_gl	double[walk_num][0:cord_num][nucl_num][4][elec_num]	in	Electron-nucleus rescaled distances derivatives
een_rescaled_e_gl	double[walk_num][0:cord_num][elec_num][4][elec_num]	in	Electron-electron rescaled distances derivatives
een_rescaled_single_n_gl	double[walk_num][0:cord_num][nucl_num][4]	in	Electron-nucleus single rescaled distances derivatives
een_rescaled_single_e_gl	double[walk_num][0:cord_num][elec_num][4]	in	Electron-electron single rescaled distances derivatives
delta_p_g	double[walk_num][0:cord_num-1][0:cord_num][4][nucl_num][elec_num]	out	Delta P matrix gradient

integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_delta_p_g_doc( &
context, num_in, walk_num, elec_num, nucl_num, cord_num,   &
een_rescaled_n, een_rescaled_e, een_rescaled_single_n, een_rescaled_single_e, &
een_rescaled_n_gl, een_rescaled_e_gl, een_rescaled_single_n_gl, een_rescaled_single_e_gl, delta_p_g) &
result(info) bind(C)
use, intrinsic :: iso_c_binding
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer(c_int64_t)    , intent(in), value  :: num_in, walk_num, elec_num, cord_num, nucl_num
real(c_double)        , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_e(elec_num, elec_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_n(nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_e(elec_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_e_gl(elec_num, 4, elec_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_n_gl(4, nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_e_gl(4,elec_num, 0:cord_num, walk_num)
real(c_double)        , intent(out)  :: delta_p_g(elec_num,nucl_num,4,0:cord_num, 0:cord_num-1,  walk_num)

double precision        :: delta_e_gl(elec_num,4)

double precision        :: een_rescaled_delta_n, een_re_n, een_re_single_n

integer*8 :: i, a, j, l, k, p, m, n, nw, num
double precision :: tmp, cummu
integer*8        :: LDA, LDB, LDC

num = num_in + 1

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (num <= 0)                      info = QMCKL_INVALID_ARG_2
if (walk_num <= 0)                 info = QMCKL_INVALID_ARG_3
if (elec_num <= 0)                 info = QMCKL_INVALID_ARG_4
if (nucl_num <= 0)                 info = QMCKL_INVALID_ARG_5
if (cord_num <  0)                 info = QMCKL_INVALID_ARG_6
if (info /= QMCKL_SUCCESS)         return


if (cord_num == 0) then
delta_p_g = 0.d0
return
endif

do nw=1, walk_num
do m=1, cord_num-1
   do j = 1, elec_num
     do k = 1, 4
       delta_e_gl(j,k) = een_rescaled_single_e_gl(k,j,m,nw) - een_rescaled_e_gl(num, k, j, m, nw)
     end do
   end do
   do k = 1, 4
     delta_e_gl(num, k) = 0.0d0
   end do

   do l=0, cord_num
     do k = 1, 3
       do a = 1, nucl_num

         cummu = 0.0d0
         do i = 1, elec_num
           cummu = cummu + delta_e_gl(i,k) * een_rescaled_n(i,a,l,nw)
         end do
         delta_p_g(num,a,k,l,m,nw) = cummu


       end do
     end do
   end do
end do
end do

end function qmckl_compute_jastrow_champ_delta_p_g_doc

integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_delta_p_g_hpc( &
context, num_in, walk_num, elec_num, nucl_num, cord_num,   &
een_rescaled_n, een_rescaled_e, een_rescaled_single_n, een_rescaled_single_e, &
een_rescaled_n_gl, een_rescaled_e_gl, een_rescaled_single_n_gl, een_rescaled_single_e_gl, delta_p_g) &
result(info) bind(C)
use, intrinsic :: iso_c_binding
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer(c_int64_t)    , intent(in), value  :: num_in, walk_num, elec_num, cord_num, nucl_num
real(c_double)        , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_e(elec_num, elec_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_n(nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_e(elec_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_e_gl(elec_num, 4, elec_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_n_gl(4, nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_e_gl(4,elec_num, 0:cord_num, walk_num)
real(c_double)        , intent(out)  :: delta_p_g(elec_num,nucl_num,4,0:cord_num, 0:cord_num-1,  walk_num)

double precision        :: delta_e_gl(3,elec_num)

double precision        :: een_rescaled_delta_n, een_re_n, een_re_single_n

integer*8 :: i, a, j, l, k, p, m, n, nw, num
double precision, allocatable :: tmp(:,:,:)
integer*8        :: LDA, LDB, LDC

num = num_in + 1

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (num <= 0)                      info = QMCKL_INVALID_ARG_2
if (walk_num <= 0)                 info = QMCKL_INVALID_ARG_3
if (elec_num <= 0)                 info = QMCKL_INVALID_ARG_4
if (nucl_num <= 0)                 info = QMCKL_INVALID_ARG_5
if (cord_num <  0)                 info = QMCKL_INVALID_ARG_6
if (info /= QMCKL_SUCCESS)         return


if (cord_num == 0) then
delta_p_g = 0.d0
return
endif

allocate( tmp(3,nucl_num,0:cord_num) )
do nw=1, walk_num
do m=1, cord_num-1
   delta_e_gl(1:3,1:elec_num) = een_rescaled_single_e_gl(1:3,1:elec_num,m,nw) - &
        een_rescaled_e_gl(num, 1:3, 1:elec_num, m, nw)
   delta_e_gl(1:3,num) = 0.0d0

   call dgemm('N','N', 3, nucl_num*(cord_num+1), elec_num, 1.d0, &
        delta_e_gl(1,1), 3, een_rescaled_n(1,1,0,nw), elec_num, 0.d0, &
        tmp, 3)

   delta_p_g(num,1:nucl_num,1,0:cord_num,m,nw) = tmp(1,1:nucl_num,0:cord_num)
   delta_p_g(num,1:nucl_num,2,0:cord_num,m,nw) = tmp(2,1:nucl_num,0:cord_num)
   delta_p_g(num,1:nucl_num,3,0:cord_num,m,nw) = tmp(3,1:nucl_num,0:cord_num)
end do
end do
deallocate(tmp)

end function qmckl_compute_jastrow_champ_delta_p_g_hpc

qmckl_exit_code
qmckl_compute_jastrow_champ_delta_p_g_hpc2 (const qmckl_context context,
                                      const int64_t num,
                                      const int64_t walk_num,
                                      const int64_t elec_num,
                                      const int64_t nucl_num,
                                      const int64_t cord_num,
                                      const double* een_rescaled_n,
                                      const double* een_rescaled_e,
                                      const double* een_rescaled_single_n,
                                      const double* een_rescaled_single_e,
                                      const double* een_rescaled_n_gl,
                                      const double* een_rescaled_e_gl,
                                      const double* een_rescaled_single_n_gl,
                                      const double* een_rescaled_single_e_gl,
                                      double* const delta_p_g )
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) return QMCKL_NULL_CONTEXT;
if (num < 0)       return QMCKL_INVALID_ARG_2;
if (walk_num <= 0) return QMCKL_INVALID_ARG_3;
if (elec_num <= 0) return QMCKL_INVALID_ARG_4;
if (nucl_num <= 0) return QMCKL_INVALID_ARG_5;
if (cord_num <  0) return QMCKL_INVALID_ARG_6;
if (een_rescaled_n == NULL) return QMCKL_INVALID_ARG_7;
if (een_rescaled_e == NULL) return QMCKL_INVALID_ARG_8;
if (een_rescaled_single_n == NULL) return QMCKL_INVALID_ARG_9;
if (een_rescaled_single_e == NULL) return QMCKL_INVALID_ARG_10;
if (een_rescaled_n_gl == NULL) return QMCKL_INVALID_ARG_11;
if (een_rescaled_e_gl == NULL) return QMCKL_INVALID_ARG_12;
if (een_rescaled_single_n_gl == NULL) return QMCKL_INVALID_ARG_13;
if (een_rescaled_single_e_gl == NULL) return QMCKL_INVALID_ARG_14;
if (delta_p_g == NULL) return QMCKL_INVALID_ARG_15;

if (cord_num == 0) {
#pragma omp parallel for
for (int64_t nw=0 ; nw<walk_num ; ++nw) {
 memset(&delta_p_g[elec_num*nucl_num*4*(cord_num+1)*cord_num*nw], 0,
        elec_num*nucl_num*4*(cord_num+1)*cord_num*sizeof(double));
}
return QMCKL_SUCCESS;
}

#pragma omp parallel for
for (int64_t nw=0 ; nw<walk_num ; ++nw) {
double delta_e_gl[elec_num*4];
for (int64_t m=1 ; m<=cord_num-1 ; ++m) {
 const double* een_rescaled_single_e_gl_ = &een_rescaled_single_e_gl[4*elec_num*(m+(cord_num+1)*nw)];
 const double* een_rescaled_e_gl_ = &een_rescaled_e_gl[num+elec_num*4*elec_num*(m+(cord_num+1)*nw)];

 #pragma omp simd
 for (int64_t kj=0 ; kj<4*elec_num ; ++kj) {
   delta_e_gl[kj] = een_rescaled_single_e_gl_[kj] - een_rescaled_e_gl_[kj+elec_num];
 }
 #pragma omp simd
 for (int64_t k=0 ; k<4; ++k) {
   delta_e_gl[4*num+k] = 0.;
 }

 for (int64_t l=0 ; l<=cord_num ; ++l) {
   for (int64_t k=0 ; k<3; ++k) {
     for (int64_t a=0 ; a<nucl_num; ++a) {
       double accu = 0.0;
       #pragma omp simd
       for (int64_t i=0 ; i<elec_num ; ++i) {
         accu += delta_e_gl[k+4*i] * een_rescaled_n[i+elec_num*(a+nucl_num*(l+(cord_num+1)*nw))];
       }
       delta_p_g[num+elec_num*(a+nucl_num*(k+4*(l+(cord_num+1)*(m+(cord_num-1+1)*nw))))] = accu;
     }
   }
 }
}
}

return QMCKL_SUCCESS;
}

Electron-electron-nucleus Jastrow gradients

Computes the gradient of the $\delta J_{een}$.

Get

qmckl_exit_code
qmckl_get_jastrow_champ_single_een_g(qmckl_context context,
                            double* const delta_een_g,
                            const int64_t size_max);

interface
  integer(qmckl_exit_code) function qmckl_get_jastrow_champ_single_een_g (context, &
       delta_een_g, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (qmckl_context) , intent(in), value :: context
    integer(c_int64_t), intent(in), value       :: size_max
    real(c_double),   intent(out)               :: delta_een_g(size_max)
  end function
end interface

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
num	int64_t	in	Index of single electron
walk_num	int64_t	in	Number of walkers
elec_num	int64_t	in	Number of electrons
nucl_num	int64_t	in	Number of nuclei
cord_num	int64_t	in	order of polynomials
dim_c_vector	int64_t	in	dimension of full coefficient vector
c_vector_full	double[dim_c_vector][nucl_num]	in	full coefficient vector
lkpm_combined_index	int64_t[4][dim_c_vector]	in	combined indices
tmp_c	double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]	in	P matrix
dtmp_c	double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]	in	P matrix derivative
delta_p	double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]	in	Delta P matrix
delta_p_gl	double[walk_num][0:cord_num-1][0:cord_num][4][nucl_num][elec_num]	in	Delta P matrix derivative
een_rescaled_n	double[walk_num][0:cord_num][nucl_num][elec_num]	in	Electron-nucleus rescaled distances
een_rescaled_single_n	double[walk_num][0:cord_num][nucl_num]	in	Electron-nucleus single rescaled distances
een_rescaled_n_gl	double[walk_num][0:cord_num][nucl_num][4][elec_num]	in	Electron-nucleus rescaled distances derivatives
een_rescaled_single_n_gl	double[walk_num][0:cord_num][nucl_num][4]	in	Electron-nucleus single rescaled distances derivatives
delta_een_g	double[walk_num][4][elec_num]	out	Delta electron-electron-nucleus jastrow gradient

integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_factor_single_een_g_doc( &
context, num_in, walk_num, elec_num, nucl_num, cord_num,   &
dim_c_vector, c_vector_full, lkpm_combined_index, &
tmp_c, dtmp_c, delta_p, delta_p_gl, een_rescaled_n, een_rescaled_single_n, &
een_rescaled_n_gl, een_rescaled_single_n_gl, delta_een_g) &
result(info) bind(C)
use, intrinsic :: iso_c_binding
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer(c_int64_t)    , intent(in), value  :: num_in, walk_num, elec_num, cord_num, nucl_num, dim_c_vector
integer(c_int64_t)    , intent(in)  :: lkpm_combined_index(dim_c_vector,4)
real(c_double)        , intent(in)  :: c_vector_full(nucl_num, dim_c_vector)
real(c_double)        , intent(in)  :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
real(c_double)        , intent(in)  :: dtmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
real(c_double)        , intent(in)  :: delta_p(elec_num, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
real(c_double)        , intent(in)  :: delta_p_gl(elec_num, nucl_num, 4, 0:cord_num, 0:cord_num-1,  walk_num)
real(c_double)        , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_n(nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(in)  :: een_rescaled_single_n_gl(4, nucl_num, 0:cord_num, walk_num)
real(c_double)        , intent(out) :: delta_een_g(elec_num, 4, walk_num)

integer*8 :: i, a, j, l, k, p, m, n, nw, kk, num
double precision :: accu, accu2, cn
integer*8                           :: LDA, LDB, LDC

double precision  :: een_rescaled_delta_n_gl(4, nucl_num, 0:cord_num, walk_num)
double precision  :: een_rescaled_delta_n(nucl_num, 0:cord_num, walk_num)

num = num_in + 1

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (walk_num <= 0)                 info = QMCKL_INVALID_ARG_3
if (elec_num <= 0)                 info = QMCKL_INVALID_ARG_4
if (nucl_num <= 0)                 info = QMCKL_INVALID_ARG_5
if (cord_num <  0)                 info = QMCKL_INVALID_ARG_6
if (info /= QMCKL_SUCCESS)         return

delta_een_g = 0.0d0

if (cord_num == 0) return

een_rescaled_delta_n(:,:,:) = een_rescaled_single_n(:,:,:) - een_rescaled_n(num, :, :, :)
een_rescaled_delta_n_gl(:,:,:,:) = een_rescaled_single_n_gl(:,:,:,:) - een_rescaled_n_gl(num, :,:,:,:)


do nw =1, walk_num
do n = 1, dim_c_vector
   l = lkpm_combined_index(n, 1)
   k = lkpm_combined_index(n, 2)
   p = lkpm_combined_index(n, 3)
   m = lkpm_combined_index(n, 4)

   do kk = 1, 3
      do a = 1, nucl_num
         cn = c_vector_full(a, n)
         if(cn == 0.d0) cycle

         !   delta_een_g(num,kk,nw) = delta_een_g(num,kk,nw) + ( &
         !        delta_p_gl(num,a,kk,m  ,k,nw) * een_rescaled_n(num,a,m+l,nw) + &
         !        delta_p_gl(num,a,kk,m+l,k,nw) * een_rescaled_n(num,a,m  ,nw) + &
         !        delta_p(num,a,m  ,k,nw) * een_rescaled_n_gl(num,kk,a,m+l,nw) + &
         !        delta_p(num,a,m+l,k,nw) * een_rescaled_n_gl(num,kk,a,m  ,nw) ) * cn

         !delta_een_g(num,kk,nw) = delta_een_g(num,kk,nw) + ( &
         !     (dtmp_c(num,kk,a,m  ,k,nw) + delta_p_gl(num,a,kk,m  ,k,nw)) * een_rescaled_delta_n(a,m+l,nw) + &
         !     (dtmp_c(num,kk,a,m+l,k,nw) + delta_p_gl(num,a,kk,m+l,k,nw)) * een_rescaled_delta_n(a,m  ,nw) + &
         !     (tmp_c(num,a,m  ,k,nw) + delta_p(num,a,m  ,k,nw)) * een_rescaled_delta_n_gl(kk,a,m+l,nw)  + &
         !     (tmp_c(num,a,m+l,k,nw) + delta_p(num,a,m+l,k,nw)) * een_rescaled_delta_n_gl(kk,a,m  ,nw) )* cn



         delta_een_g(num,kk,nw) = delta_een_g(num,kk,nw) + ( &
              dtmp_c(num,kk,a,m  ,k,nw) * een_rescaled_delta_n(a,m+l,nw) + &
              dtmp_c(num,kk,a,m+l,k,nw) * een_rescaled_delta_n(a,m  ,nw) + &
              tmp_c(num,a,m  ,k,nw) * een_rescaled_delta_n_gl(kk,a,m+l,nw)  + &
              tmp_c(num,a,m+l,k,nw) * een_rescaled_delta_n_gl(kk,a,m  ,nw) + &
              delta_p_gl(num,a,kk,m  ,k,nw) * een_rescaled_single_n(a,m+l,nw) + &
              delta_p_gl(num,a,kk,m+l,k,nw) * een_rescaled_single_n(a,m  ,nw) + &
              delta_p(num,a,m  ,k,nw) * een_rescaled_single_n_gl(kk,a,m+l,nw) + & 
              delta_p(num,a,m+l,k,nw) * een_rescaled_single_n_gl(kk,a,m  ,nw) )* cn
      end do
   end do
end do
end do

end function qmckl_compute_jastrow_champ_factor_single_een_g_doc

Electron-electron Jastrow

Calculates the single-electron move contribution to the $J_{ee}$ Jastrow factor. Similar to the 3-body case above, we need to compute the single-electron contributions to the rescaled distances. The rescaled electron-electron distance for the 2-body term is defined as \[ \widetilde{r}_{ij} = \frac{1-e^{-\kappa r_{ij}}}{\kappa}. \] For the electron-electro Jastrow, no intermediate terms have to be calculated, as we can immidiatly calculate \[ \delta J_{ee} = \sum_i \sum_{j\lt i} \delta_{j,\text{num}} \left(\frac{b_1 \widetilde{r}_{ij}^{\text{new}}}{1+b_2 \widetilde{r}_{ij}^{\text{new}}} - \frac{b_1 \widetilde{r}_{ij}^{\text{old}}}{1+b_2 \widetilde{r}_{ij}^{\text{old}}} + \sum_{k=2}^{N_{\text{bord}}} b_k ({\widetilde{r}_{ij}^{\text{new}}}^k - {\widetilde{r}_{ij}^{\text{old}}}^k )\right), \] where num is the electron that is being moved. The gradient and Laplacian of the rescaled distances are given by \[ \partial_{i,m} \widetilde{r}_{ij} = \frac{r_{i,m} - r_{j,m}}{r_{ij}} e^{-\kappa r_{ij}} \quad \text{and} \quad \partial_{i,4} \widetilde{r}_{ij} = \left(\frac{2}{r_{ij}} - \kappa \right)e^{-\kappa r_{ij}} \] Then,

\begin{eqnarray*} \partial_{i,m}\delta J_{ee} = & \sum_j \delta_{j,\text{num}} \left( \frac{b_1 \partial_{i,m}\widetilde{r}^{\text{new}}_{ij}}{(1+b_2 \widetilde{r}^{\text{new}}_{ij})^2} - \frac{b_1 \partial_{i,m}\widetilde{r}^{\text{old}}_{ij}}{(1+b_2 \widetilde{r}^{\text{old}}_{ij})^2} + \sum_{k=2}^{N_{\text{bord}}} b_k k \left({\widetilde{r}_{ij}^{\text{new}}}^{k-1} \partial_{i,m}\widetilde{r}^{\text{new}}_{ij} - {\widetilde{r}_{ij}^{\text{old}}}^{k-1} \partial_{i,m}\widetilde{r}^{\text{old}}_{ij} \right) \right)\\ \partial_{i,4}\delta J_{ee} = & \sum_j \delta_{j,\text{num}} \left( \frac{b_1}{(1+b_2 \widetilde{r}^{\text{new}}_{ij})^2} \left(\partial_{i,4}\widetilde{r}^{\text{new}}_{ij} - 2 b_2 \frac{\nabla\widetilde{r}^{\text{new}}_{ij} \cdot \nabla \widetilde{r}^{\text{new}}_{ij}}{1+b_2 \widetilde{r}^{\text{new}}_{ij}} \right) - \frac{b_1}{(1+b_2 \widetilde{r}^{\text{old}}_{ij})^2} \left(\partial_{i,4}\widetilde{r}^{\text{old}}_{ij} - 2 b_2 \frac{\nabla\widetilde{r}^{\text{old}}_{ij} \cdot \nabla \widetilde{r}^{\text{old}}_{ij}}{1+b_2 \widetilde{r}^{\text{old}}_{ij}} \right) \right.\\ & \left.+ \sum_{k=2}^{N_{\text{bord}}} b_k k \left({\widetilde{r}_{ij}^{\text{new}}}^{k-1}\partial_{i,4}\widetilde{r}_{ij}^{\text{new}} + (k-1)(\nabla\widetilde{r}_{ij}^{\text{new}} \cdot \nabla \widetilde{r}_{ij}^{\text{new}}) {\widetilde{r}_{ij}^{\text{new}}}^{k-2} - {\widetilde{r}_{ij}^{\text{old}}}^{k-1}\partial_{i,4}\widetilde{r}_{ij}^{\text{old}} - (k-1)(\nabla\widetilde{r}_{ij}^{\text{old}} \cdot \nabla \widetilde{r}_{ij}^{\text{old}}) {\widetilde{r}_{ij}^{\text{old}}}^{k-2} \right) \right) \end{eqnarray*}

In the special case that $i = \text{num}$, the gradient and Laplacian change to

\begin{eqnarray*} \partial_{\text{num},m}\delta J_{ee} = & -\sum_i \sum_j \delta_{j,\text{num}} \left(\frac{b_1 \partial_{i,m}\widetilde{r}^{\text{new}}_{ij}}{(1+b_2 \widetilde{r}^{\text{new}}_{ij})^2} - \frac{b_1 \partial_{i,m}\widetilde{r}^{\text{old}}_{ij}}{(1+b_2 \widetilde{r}^{\text{old}}_{ij})^2} + \sum_{k=2}^{N_{\text{bord}}} b_k k \left({\widetilde{r}_{ij}^{\text{new}}}^{k-1} \partial_{i,m}\widetilde{r}^{\text{new}}_{ij} - {\widetilde{r}_{ij}^{\text{old}}}^{k-1} \partial_{i,m}\widetilde{r}^{\text{old}}_{ij} \right) \right) \\ \partial_{\text{num},4}\delta J_{ee} = & -\sum_i \sum_j \delta_{j,\text{num}} \left( \frac{b_1}{(1+b_2 \widetilde{r}^{\text{new}}_{ij})^2} \left(\partial_{i,4}\widetilde{r}^{\text{new}}_{ij} - 2 b_2 \frac{\nabla\widetilde{r}^{\text{new}}_{ij} \cdot \nabla \widetilde{r}^{\text{new}}_{ij}}{1+b_2 \widetilde{r}^{\text{new}}_{ij}} \right) - \frac{b_1}{(1+b_2 \widetilde{r}^{\text{old}}_{ij})^2} \left(\partial_{i,4}\widetilde{r}^{\text{old}}_{ij} - 2 b_2 \frac{\nabla\widetilde{r}^{\text{old}}_{ij} \cdot \nabla \widetilde{r}^{\text{old}}_{ij}}{1+b_2 \widetilde{r}^{\text{old}}_{ij}} \right)\right. \\ & \left. + \sum_{k=2}^{N_{\text{bord}}} b_k k \left({\widetilde{r}_{ij}^{\text{new}}}^{k-1}\partial_{i,4}\widetilde{r}_{ij}^{\text{new}} + (k-1)(\nabla\widetilde{r}_{ij}^{\text{new}} \cdot \nabla \widetilde{r}_{ij}^{\text{new}}) {\widetilde{r}_{ij}^{\text{new}}}^{k-2} - {\widetilde{r}_{ij}^{\text{old}}}^{k-1}\partial_{i,4}\widetilde{r}_{ij}^{\text{old}} - (k-1)(\nabla\widetilde{r}_{ij}^{\text{old}} \cdot \nabla \widetilde{r}_{ij}^{\text{old}}) {\widetilde{r}_{ij}^{\text{old}}}^{k-2} \right) \right) \end{eqnarray*}

Electron-electron rescaled distance

Calculates the rescaled electron-electron distances \[ \widetilde{r}_{ij} = \frac{1-e^{-\kappa r_{ij}}}{\kappa}. \]

Get

qmckl_exit_code
qmckl_get_ee_rescaled_single(qmckl_context context,
                            double* const distance_rescaled,
                            const int64_t size_max);

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
elec_num	int64_t	in	Number of electrons
rescale_factor_ee	double	in	Factor to rescale ee distances
walk_num	int64_t	in	Number of walkers
single_ee_distance	double[walk_num][elec_num]	in	Single electron-electron distances
ee_rescaled_single	double[walk_num][elec_num]	out	Electron-electron rescaled distances

function qmckl_compute_ee_rescaled_single_doc(context, &
elec_num, rescale_factor_ee, walk_num, &
single_ee_distance, ee_rescaled_single) &
bind(C) result(info)
use qmckl
implicit none

integer(qmckl_context), intent(in), value :: context
integer (c_int64_t) , intent(in)  , value :: elec_num
real    (c_double ) , intent(in)  , value :: rescale_factor_ee
integer (c_int64_t) , intent(in)  , value :: walk_num
real    (c_double ) , intent(in)          :: single_ee_distance(elec_num,walk_num)
real    (c_double ) , intent(out)         :: ee_rescaled_single(elec_num,walk_num)
integer(qmckl_exit_code)                  :: info

integer*8 :: k, i
real (c_double) :: inverse_rescale_factor_ee

inverse_rescale_factor_ee = 1.0d0 / rescale_factor_ee

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif

if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif

if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif

do k=1,walk_num
 do i=1,elec_num
ee_rescaled_single(i,k) = (1.0d0 - dexp(-rescale_factor_ee * single_ee_distance(i,k))) * inverse_rescale_factor_ee
enddo
end do

end function qmckl_compute_ee_rescaled_single_doc

Electron-electron Jastrow value

Calculate the single-electron move contribution to the $J_{ee}$ Jastrow factor. \[ \delta J_{ee} = \sum_i \sum_{j\lt i} \delta_{j,\text{num}} \left(\frac{b_1 \widetilde{r}_{ij}^{\text{new}}}{1+b_2 \widetilde{r}_{ij}^{\text{new}}} - \frac{b_1 \widetilde{r}_{ij}^{\text{old}}}{1+b_2 \widetilde{r}_{ij}^{\text{old}}} + \sum_{k=2}^{N_{\text{bord}}} b_k ({\widetilde{r}_{ij}^{\text{new}}}^k - {\widetilde{r}_{ij}^{\text{old}}}^k )\right), \]

Get

qmckl_exit_code
qmckl_get_jastrow_champ_single_ee(qmckl_context context,
                       double* const delta_ee,
                       const int64_t size_max);

interface
  integer(qmckl_exit_code) function qmckl_get_jastrow_champ_single_ee (context, &
       delta_ee, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (qmckl_context) , intent(in), value :: context
    integer(c_int64_t), intent(in), value       :: size_max
    real(c_double),   intent(out)               :: delta_ee(size_max)
  end function qmckl_get_jastrow_champ_single_ee
end interface

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
num	int64_t	in	Index of single point
walk_num	int64_t	in	Number of walkers
elec_num	int64_t	in	Number of electrons
up_num	int64_t	in	Number of alpha electrons
bord_num	int64_t	in	Number of coefficients
b_vector	double[bord_num+1]	in	List of coefficients
ee_distance_rescaled	double[walk_num][elec_num][elec_num]	in	Electron-electron rescaled distances
ee_rescaled_single	double[walk_num][elec_num]	in	Electron-electron rescaled single distances
delta_ee	double[walk_num]	out	Single electron-electron Jastrow

function qmckl_compute_jastrow_champ_single_ee_doc(context, &
num_in, walk_num, elec_num, up_num, bord_num, b_vector, &
ee_distance_rescaled, ee_rescaled_single, spin_independent, delta_ee) &
bind(C) result(info)
use qmckl
implicit none

integer (qmckl_context), intent(in), value :: context
integer (c_int64_t)    , intent(in), value :: num_in
integer (c_int64_t)    , intent(in), value :: walk_num
integer (c_int64_t)    , intent(in), value :: elec_num
integer (c_int64_t)    , intent(in), value :: up_num
integer (c_int64_t)    , intent(in), value :: bord_num
real    (c_double )    , intent(in)        :: b_vector(bord_num+1)
real    (c_double )    , intent(in)        :: ee_distance_rescaled(elec_num,elec_num,walk_num)
real    (c_double )    , intent(in)        :: ee_rescaled_single(elec_num,walk_num)
integer (c_int32_t)    , intent(in), value :: spin_independent
real    (c_double )    , intent(out)       :: delta_ee(walk_num)
integer(qmckl_exit_code)                   :: info

integer*8 :: i, j, k, nw, num
double precision   :: x, xk, y, yk

num = num_in + 1

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif

if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif

if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif

if (bord_num < 0) then
info = QMCKL_INVALID_ARG_5
return
endif

do nw =1, walk_num

delta_ee(nw) = 0.0d0
do i=1,elec_num
    !print *,i, ee_rescaled_single(i,nw)
    !print *, i, ee_distance_rescaled(i,num,nw)
    !print *, '   '
   if (i.ne.num) then
      x = ee_distance_rescaled(i,num,nw)
      y = ee_rescaled_single(i,nw)
      if (spin_independent == 1) then
         delta_ee(nw) = delta_ee(nw) - (b_vector(1) * x / (1.d0 + b_vector(2) * x)) &
              + (b_vector(1) * y / (1.d0 + b_vector(2) * y))
      else
         if ((i <= up_num .and. num <= up_num ) .or. (i >  up_num .and. num >  up_num)) then
            delta_ee(nw) = delta_ee(nw) - (0.5d0 * b_vector(1) * x / (1.d0 + b_vector(2) * x)) &
                 + (0.5d0 * b_vector(1) * y / (1.d0 + b_vector(2) * y))
         else
            delta_ee(nw) = delta_ee(nw) - (b_vector(1) * x / (1.d0 + b_vector(2) * x)) &
                 + (b_vector(1) * y / (1.d0 + b_vector(2) * y))
         endif
      endif

      xk = x
      yk = y
      do k=2,bord_num
         xk = xk * x
         yk = yk * y
         delta_ee(nw) = delta_ee(nw) - (b_vector(k+1) * xk) + (b_vector(k+1) * yk)
      end do
   endif
end do

end do

end function qmckl_compute_jastrow_champ_single_ee_doc

Electron-electron rescaled distances derivatives

Calculate the derivative of the rescaled electron-electron distances. \[ \partial_{i,m} \widetilde{r}_{ij} = \frac{r_{i,m} - r_{j,m}}{r_{ij}} e^{-\kappa r_{ij}} \quad \text{and} \quad \partial_{i,4} \widetilde{r}_{ij} = \left(\frac{2}{r_{ij}} - \kappa \right)e^{-\kappa r_{ij}} \]

Get

qmckl_exit_code qmckl_get_ee_rescaled_single_gl(qmckl_context context,
                                           double* const distance_rescaled_gl,
                                           const int64_t size_max);

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
num	int64_t	in	Index of single electron
elec_num	int64_t	in	Number of electrons
rescale_factor_ee	double	in	Factor to rescale ee distances
walk_num	int64_t	in	Number of walkers
single_ee_distance	double[elec_num][walk_num]	in	Single electron-electron distances
elec_coord	double[3][walk_num][elec_num]	in	Electron coordinates
coord	double[walk_num][3]	in	Single electron coordinates
ee_rescaled_single_gl	double[walk_num][elec_num][4]	out	Electron-electron rescaled single distance derivatives

function qmckl_compute_ee_rescaled_single_gl_doc(context, num_in,  &
elec_num, rescale_factor_ee, walk_num, single_ee_distance, elec_coord, coord, ee_rescaled_single_gl) &
bind(C) result(info)
use qmckl
implicit none

integer(qmckl_context), intent(in), value :: context
integer (c_int64_t) , intent(in)  , value :: num_in
integer (c_int64_t) , intent(in)  , value :: elec_num
real    (c_double ) , intent(in)  , value :: rescale_factor_ee
integer (c_int64_t) , intent(in)  , value :: walk_num
real    (c_double ) , intent(in)          :: single_ee_distance(elec_num,walk_num)
real    (c_double ) , intent(in)          :: elec_coord(elec_num,walk_num,3)
real    (c_double ) , intent(in)          :: coord(3,walk_num)
real    (c_double ) , intent(out)         :: ee_rescaled_single_gl(4,elec_num,walk_num)
integer(qmckl_exit_code)                  :: info

integer*8 :: nw, i, ii, num
double precision :: rij_inv, elel_dist_gl(4, elec_num), kappa_l

num = num_in + 1

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif

if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif

if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif


ee_rescaled_single_gl = 0.0d0
do nw = 1, walk_num

! prepare the actual een table
do i = 1, elec_num
   rij_inv = 1.0d0 / single_ee_distance(i, nw)
   do ii = 1, 3
     elel_dist_gl(ii, i) = (elec_coord(i,nw, ii) - coord(ii,nw)) * rij_inv
   end do
   elel_dist_gl(4, i) = 2.0d0 * rij_inv
end do

 do i = 1, elec_num
   kappa_l = -1 * rescale_factor_ee
   ee_rescaled_single_gl(1, i, nw) = elel_dist_gl(1, i)
   ee_rescaled_single_gl(2, i, nw) = elel_dist_gl(2, i)
   ee_rescaled_single_gl(3, i, nw) = elel_dist_gl(3, i)
   ee_rescaled_single_gl(4, i, nw) = elel_dist_gl(4, i)

   ee_rescaled_single_gl(4, i, nw) = ee_rescaled_single_gl(4, i, nw) + kappa_l

   ee_rescaled_single_gl(1, i, nw) = ee_rescaled_single_gl(1, i, nw)  * dexp(kappa_l * single_ee_distance(i,nw))
   ee_rescaled_single_gl(2, i, nw) = ee_rescaled_single_gl(2, i, nw)  * dexp(kappa_l * single_ee_distance(i,nw))
   ee_rescaled_single_gl(3, i, nw) = ee_rescaled_single_gl(3, i, nw)  * dexp(kappa_l * single_ee_distance(i,nw))
   ee_rescaled_single_gl(4, i, nw) = ee_rescaled_single_gl(4, i, nw)  * dexp(kappa_l * single_ee_distance(i,nw))
end do

ee_rescaled_single_gl(1, num, nw) = 0.0d0
ee_rescaled_single_gl(2, num, nw) = 0.0d0
ee_rescaled_single_gl(3, num, nw) = 0.0d0
ee_rescaled_single_gl(4, num, nw) = 0.0d0
end do


end function qmckl_compute_ee_rescaled_single_gl_doc

Electron-electron Jastrow gradients and Laplacian

Calculates the gradient and Laplacian of the electron-electron Jastrow factor.

\begin{eqnarray*} \partial_{i,m}\delta J_{ee} = & \sum_j \delta_{j,\text{num}} \left( \frac{b_1 \partial_{i,m}\widetilde{r}^{\text{new}}_{ij}}{(1+b_2 \widetilde{r}^{\text{new}}_{ij})^2} - \frac{b_1 \partial_{i,m}\widetilde{r}^{\text{old}}_{ij}}{(1+b_2 \widetilde{r}^{\text{old}}_{ij})^2} + \sum_{k=2}^{N_{\text{bord}}} b_k k \left({\widetilde{r}_{ij}^{\text{new}}}^{k-1} \partial_{i,m}\widetilde{r}^{\text{new}}_{ij} - {\widetilde{r}_{ij}^{\text{old}}}^{k-1} \partial_{i,m}\widetilde{r}^{\text{old}}_{ij} \right) \right)\\ \partial_{i,4}\delta J_{ee} = & \sum_j \delta_{j,\text{num}} \left( \frac{b_1}{(1+b_2 \widetilde{r}^{\text{new}}_{ij})^2} \left(\partial_{i,4}\widetilde{r}^{\text{new}}_{ij} - 2 b_2 \frac{\nabla\widetilde{r}^{\text{new}}_{ij} \cdot \nabla \widetilde{r}^{\text{new}}_{ij}}{1+b_2 \widetilde{r}^{\text{new}}_{ij}} \right) - \frac{b_1}{(1+b_2 \widetilde{r}^{\text{old}}_{ij})^2} \left(\partial_{i,4}\widetilde{r}^{\text{old}}_{ij} - 2 b_2 \frac{\nabla\widetilde{r}^{\text{old}}_{ij} \cdot \nabla \widetilde{r}^{\text{old}}_{ij}}{1+b_2 \widetilde{r}^{\text{old}}_{ij}} \right) \right.\\ & \left.+ \sum_{k=2}^{N_{\text{bord}}} b_k k \left({\widetilde{r}_{ij}^{\text{new}}}^{k-1}\partial_{i,4}\widetilde{r}_{ij}^{\text{new}} + (k-1)(\nabla\widetilde{r}_{ij}^{\text{new}} \cdot \nabla \widetilde{r}_{ij}^{\text{new}}) {\widetilde{r}_{ij}^{\text{new}}}^{k-2} - {\widetilde{r}_{ij}^{\text{old}}}^{k-1}\partial_{i,4}\widetilde{r}_{ij}^{\text{old}} - (k-1)(\nabla\widetilde{r}_{ij}^{\text{old}} \cdot \nabla \widetilde{r}_{ij}^{\text{old}}) {\widetilde{r}_{ij}^{\text{old}}}^{k-2} \right) \right) \end{eqnarray*}

In the special case that $i = \text{num}$, the gradient and Laplacian change to

\begin{eqnarray*} \partial_{\text{num},m}\delta J_{ee} = & -\sum_i \sum_j \delta_{j,\text{num}} \left(\frac{b_1 \partial_{i,m}\widetilde{r}^{\text{new}}_{ij}}{(1+b_2 \widetilde{r}^{\text{new}}_{ij})^2} - \frac{b_1 \partial_{i,m}\widetilde{r}^{\text{old}}_{ij}}{(1+b_2 \widetilde{r}^{\text{old}}_{ij})^2} + \sum_{k=2}^{N_{\text{bord}}} b_k k \left({\widetilde{r}_{ij}^{\text{new}}}^{k-1} \partial_{i,m}\widetilde{r}^{\text{new}}_{ij} - {\widetilde{r}_{ij}^{\text{old}}}^{k-1} \partial_{i,m}\widetilde{r}^{\text{old}}_{ij} \right) \right) \\ \partial_{\text{num},4}\delta J_{ee} = & -\sum_i \sum_j \delta_{j,\text{num}} \left( \frac{b_1}{(1+b_2 \widetilde{r}^{\text{new}}_{ij})^2} \left(\partial_{i,4}\widetilde{r}^{\text{new}}_{ij} - 2 b_2 \frac{\nabla\widetilde{r}^{\text{new}}_{ij} \cdot \nabla \widetilde{r}^{\text{new}}_{ij}}{1+b_2 \widetilde{r}^{\text{new}}_{ij}} \right) - \frac{b_1}{(1+b_2 \widetilde{r}^{\text{old}}_{ij})^2} \left(\partial_{i,4}\widetilde{r}^{\text{old}}_{ij} - 2 b_2 \frac{\nabla\widetilde{r}^{\text{old}}_{ij} \cdot \nabla \widetilde{r}^{\text{old}}_{ij}}{1+b_2 \widetilde{r}^{\text{old}}_{ij}} \right)\right. \\ & \left. + \sum_{k=2}^{N_{\text{bord}}} b_k k \left({\widetilde{r}_{ij}^{\text{new}}}^{k-1}\partial_{i,4}\widetilde{r}_{ij}^{\text{new}} + (k-1)(\nabla\widetilde{r}_{ij}^{\text{new}} \cdot \nabla \widetilde{r}_{ij}^{\text{new}}) {\widetilde{r}_{ij}^{\text{new}}}^{k-2} - {\widetilde{r}_{ij}^{\text{old}}}^{k-1}\partial_{i,4}\widetilde{r}_{ij}^{\text{old}} - (k-1)(\nabla\widetilde{r}_{ij}^{\text{old}} \cdot \nabla \widetilde{r}_{ij}^{\text{old}}) {\widetilde{r}_{ij}^{\text{old}}}^{k-2} \right) \right) \end{eqnarray*}

Get

qmckl_exit_code
qmckl_get_jastrow_champ_single_ee_gl(qmckl_context context,
                               double* const delta_ee_gl,
                               const int64_t size_max);

interface
  integer(qmckl_exit_code) function qmckl_get_jastrow_champ_single_ee_gl (context, &
       delta_ee_gl, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (qmckl_context) , intent(in), value :: context
    integer(c_int64_t), intent(in), value       :: size_max
    real(c_double),   intent(out)               :: delta_ee_gl(size_max)
  end function qmckl_get_jastrow_champ_single_ee_gl
end interface

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
num	int64_t	in	Index of single electron
walk_num	int64_t	in	Number of walkers
elec_num	int64_t	in	Number of electrons
up_num	int64_t	in	Number of alpha electrons
bord_num	int64_t	in	Number of coefficients
b_vector	double[bord_num+1]	in	List of coefficients
ee_distance_rescaled	double[walk_num][elec_num][elec_num]	in	Electron-electron rescaled distances
ee_distance_rescaled_gl	double[walk_num][4][elec_num][elec_num]	in	Electron-electron rescaled distances derivatives
ee_rescaled_single	double[walk_num][elec_num]	in	Electron-electron rescaled single distances
ee_rescaled_single_gl	double[walk_num][4][elec_num]	in	Electron-electron rescaled single distances derivatives
spin_independent	int32_t	in	If 1, same parameters for parallel and antiparallel spins
delta_ee_gl	double[walk_num][elec_num][4]	out	Single electron-electron jastrow gradients and Laplacian

function qmckl_compute_jastrow_champ_single_ee_gl_doc( &
context, num_in, walk_num, elec_num, up_num, bord_num, &
b_vector, ee_distance_rescaled, ee_distance_rescaled_gl,  &
ee_rescaled_single, ee_rescaled_single_gl,  &
spin_independent, delta_ee_gl) &
bind(C) result(info)
use qmckl
implicit none

integer (qmckl_context), intent(in), value :: context
integer (c_int64_t) , intent(in)  , value :: num_in
integer (c_int64_t) , intent(in)  , value :: walk_num
integer (c_int64_t) , intent(in)  , value :: elec_num
integer (c_int64_t) , intent(in)  , value :: up_num
integer (c_int64_t) , intent(in)  , value :: bord_num
real    (c_double ) , intent(in)          :: b_vector(bord_num+1)
real    (c_double ) , intent(in)          :: ee_distance_rescaled(elec_num,elec_num,walk_num)
real    (c_double ) , intent(in)          :: ee_distance_rescaled_gl(4,elec_num,elec_num,walk_num)
real    (c_double ) , intent(in)          :: ee_rescaled_single(elec_num,walk_num)
real    (c_double ) , intent(in)          :: ee_rescaled_single_gl(4,elec_num,walk_num)
integer (c_int32_t) , intent(in)  , value :: spin_independent
real    (c_double ) , intent(out)         :: delta_ee_gl(4,elec_num,walk_num)
integer(qmckl_exit_code)                  :: info

integer*8 :: i, j, k, nw, ii, num
double precision   :: x, x1, kf, x_old, x1_old
double precision   :: denom, invdenom, invdenom2, f
double precision   :: denom_old, invdenom_old, invdenom2_old, f_old
double precision   :: grad_c2, grad_c2_old
double precision   :: dx(4), dx_old(4)

num = num_in + 1

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif

if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif

if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif

if (bord_num < 0) then
info = QMCKL_INVALID_ARG_5
return
endif

if ((spin_independent < 0).or.(spin_independent > 1)) then
info = QMCKL_INVALID_ARG_8
return
endif

do nw =1, walk_num
delta_ee_gl(:,:,nw) = 0.0d0
   do i = 1, elec_num
      if (i == num) cycle

      x = ee_rescaled_single(i,nw)
      x_old = ee_distance_rescaled(i,num,nw)

      denom         = 1.0d0 + b_vector(2) * x
      invdenom      = 1.0d0 / denom
      invdenom2     = invdenom * invdenom

      denom_old         = 1.0d0 + b_vector(2) * x_old
      invdenom_old      = 1.0d0 / denom_old
      invdenom2_old     = invdenom_old * invdenom_old

      dx(1) = ee_rescaled_single_gl(1, i, nw)
      dx(2) = ee_rescaled_single_gl(2, i, nw)
      dx(3) = ee_rescaled_single_gl(3, i, nw)
      dx(4) = ee_rescaled_single_gl(4, i, nw)

      dx_old(1) = ee_distance_rescaled_gl(1, i, num, nw)
      dx_old(2) = ee_distance_rescaled_gl(2, i, num, nw)
      dx_old(3) = ee_distance_rescaled_gl(3, i, num, nw)
      dx_old(4) = ee_distance_rescaled_gl(4, i, num, nw)

      grad_c2 = dx(1)*dx(1) + dx(2)*dx(2) + dx(3)*dx(3)
      grad_c2_old = dx_old(1)*dx_old(1) + dx_old(2)*dx_old(2) + dx_old(3)*dx_old(3)

      if (spin_independent == 1) then
         f = b_vector(1) * invdenom2
         f_old = b_vector(1) * invdenom2_old
      else
         if((i <= up_num .and. num <= up_num ) .or. (i >  up_num .and. num >  up_num)) then
            f = 0.5d0 * b_vector(1) * invdenom2
            f_old = 0.5d0 * b_vector(1) * invdenom2_old
         else
            f = b_vector(1) * invdenom2
            f_old = b_vector(1) * invdenom2_old
         end if
      end if

      delta_ee_gl(1,i,nw) = delta_ee_gl(1,i,nw) + f * dx(1) - f_old * dx_old(1)
      delta_ee_gl(2,i,nw) = delta_ee_gl(2,i,nw) + f * dx(2) - f_old * dx_old(2)
      delta_ee_gl(3,i,nw) = delta_ee_gl(3,i,nw) + f * dx(3) - f_old * dx_old(3)
      delta_ee_gl(4,i,nw) = delta_ee_gl(4,i,nw) &
           + f * (dx(4) - 2.d0 * b_vector(2) * grad_c2 * invdenom) &
           - f_old * (dx_old(4) - 2.d0 * b_vector(2) * grad_c2_old * invdenom_old)

      delta_ee_gl(1,num,nw) = delta_ee_gl(1,num,nw) - f * dx(1) + f_old * dx_old(1)
      delta_ee_gl(2,num,nw) = delta_ee_gl(2,num,nw) - f * dx(2) + f_old * dx_old(2)
      delta_ee_gl(3,num,nw) = delta_ee_gl(3,num,nw) - f * dx(3) + f_old * dx_old(3)
      delta_ee_gl(4,num,nw) = delta_ee_gl(4,num,nw) &
           + f * (dx(4) - 2.d0 * b_vector(2) * grad_c2 * invdenom) &
           - f_old * (dx_old(4) - 2.d0 * b_vector(2) * grad_c2_old * invdenom_old)


      kf = 2.d0
      x1 = x
      x1_old = x_old
      x = 1.d0
      x_old = 1.d0
      do k=2, bord_num
         f = b_vector(k+1) * kf * x
         f_old = b_vector(k+1) * kf * x_old
         delta_ee_gl(1,i,nw) = delta_ee_gl(1,i,nw) + f * x1 * dx(1) - f_old * x1_old * dx_old(1)
         delta_ee_gl(2,i,nw) = delta_ee_gl(2,i,nw) + f * x1 * dx(2) - f_old * x1_old * dx_old(2)
         delta_ee_gl(3,i,nw) = delta_ee_gl(3,i,nw) + f * x1 * dx(3) - f_old * x1_old * dx_old(3)
         delta_ee_gl(4,i,nw) = delta_ee_gl(4,i,nw) &
              + f * (x1 * dx(4) + (kf-1.d0) * grad_c2) &
              - f_old * (x1_old * dx_old(4) + (kf-1.d0) * grad_c2_old)

         delta_ee_gl(1,num,nw) = delta_ee_gl(1,num,nw) - f * x1 * dx(1) + f_old * x1_old * dx_old(1)
         delta_ee_gl(2,num,nw) = delta_ee_gl(2,num,nw) - f * x1 * dx(2) + f_old * x1_old * dx_old(2)
         delta_ee_gl(3,num,nw) = delta_ee_gl(3,num,nw) - f * x1 * dx(3) + f_old * x1_old * dx_old(3)
         delta_ee_gl(4,num,nw) = delta_ee_gl(4,num,nw) &
              + f * (x1 * dx(4) + (kf-1.d0) * grad_c2) &
              - f_old * (x1_old * dx_old(4) + (kf-1.d0) * grad_c2_old)
         x = x*x1
         x_old = x_old*x1_old
         kf = kf + 1.d0
      end do


   end do

end do

end function qmckl_compute_jastrow_champ_single_ee_gl_doc

Electron-nucleus Jastrow

Here we calculate the single-electron move contribution to the $J_{en}$ Jastrow factor. First, we need to compute the single-electron contributions to the rescaled distances. The rescaled electron-electron distance for the 2-body term is defined as \[ \widetilde{R}_{i\alpha} = \frac{1-e^{-\kappa R_{i\alpha}}}{\kappa}. \] For the electron-nucleus Jastrow, no intermediate terms have to be calculated, as we can immidiatly calculate \[ \delta J_{en} = \sum_i \delta_{i,\text{num}} \left( \frac{a_1 \widetilde{R}_{i\alpha}^{\text{new}}}{1+a_2 \widetilde{R}_{i\alpha}^{\text{new}}} - \frac{a_1 \widetilde{R}_{i\alpha}^{\text{old}}}{1+a_2 \widetilde{R}_{i\alpha}^{\text{old}}} + \sum_{k=2}^{N_\text{aord}} a_k\left({\widetilde{R}_{i\alpha}^{\text{new}}}^k - {\widetilde{R}_{i\alpha}^{\text{old}}}^k \right)\right) \] where num is the electron that is being moved. The gradient and Laplacian of the rescaled distances are given by \[ \partial_{i,m} \widetilde{R}_{i\alpha} = \frac{r_{i,m} - R_{\alpha,m}}{R_{i\alpha}} e^{-\kappa R_{i\alpha}} \quad \text{and} \quad \partial_{i,4} \widetilde{R}_{i\alpha} = \left(\frac{2}{R_{i\alpha}} - \kappa \right)e^{-\kappa R_{i\alpha}} \] Then,

\begin{eqnarray*} \partial_{i,m}\delta J_{en} = & \sum_\alpha \delta_{i,\text{num}} \left( \frac{a_1 \partial_{i,m}\widetilde{R}^{\text{new}}_{i\alpha}}{(1+a_2 \widetilde{R}^{\text{new}}_{i\alpha})^2} - \frac{a_1 \partial_{i,m}\widetilde{R}^{\text{old}}_{i\alpha}}{(1+a_2 \widetilde{R}^{\text{old}}_{i\alpha})^2} + \sum_{k=2}^{N_{\text{aord}}} a_k k \left({\widetilde{R}_{i\alpha}^{\text{new}}}^{k-1} \partial_{i,m}\widetilde{R}^{\text{new}}_{i\alpha} - {\widetilde{R}_{i\alpha}^{\text{old}}}^{k-1} \partial_{i,m}\widetilde{R}^{\text{old}}_{i\alpha} \right) \right)\\ \partial_{i,4}\delta J_{en} = & \sum_\alpha \delta_{i,\text{num}} \left( \frac{a_1}{(1+a_2 \widetilde{R}^{\text{new}}_{i\alpha})^2} \left(\partial_{i,4}\widetilde{R}^{\text{new}}_{i\alpha} - 2 a_2 \frac{\nabla\widetilde{R}^{\text{new}}_{i\alpha} \cdot \nabla \widetilde{R}^{\text{new}}_{i\alpha}}{1+a_2 \widetilde{R}^{\text{new}}_{i\alpha}} \right) - \frac{a_1}{(1+a_2 \widetilde{R}^{\text{old}}_{i\alpha})^2} \left(\partial_{i,4}\widetilde{R}^{\text{old}}_{i\alpha} - 2 a_2 \frac{\nabla\widetilde{R}^{\text{old}}_{i\alpha} \cdot \nabla \widetilde{R}^{\text{old}}_{i\alpha}}{1+b_2 \widetilde{R}^{\text{old}}_{i\alpha}} \right) \right.\\ & \left.+ \sum_{k=2}^{N_{\text{aord}}} a_k k \left({\widetilde{R}_{i\alpha}^{\text{new}}}^{k-1}\partial_{i,4}\widetilde{R}_{i\alpha}^{\text{new}} + (k-1)(\nabla\widetilde{R}_{i\alpha}^{\text{new}} \cdot \nabla \widetilde{R}_{i\alpha}^{\text{new}}) {\widetilde{R}_{i\alpha}^{\text{new}}}^{k-2} - {\widetilde{R}_{i\alpha}^{\text{old}}}^{k-1}\partial_{i,4}\widetilde{R}_{i\alpha}^{\text{old}} - (k-1)(\nabla\widetilde{R}_{i\alpha}^{\text{old}} \cdot \nabla \widetilde{R}_{i\alpha}^{\text{old}}) {\widetilde{R}_{i\alpha}^{\text{old}}}^{k-2} \right) \right) \end{eqnarray*}

Electron-nucleus rescaled distance

Calculate the updated rescaled electron-nucleus distances \[ \widetilde{R}_{i\alpha} = \frac{1-e^{-\kappa R_{i\alpha}}}{\kappa}. \]

Get

qmckl_exit_code
qmckl_get_en_rescaled_single(qmckl_context context,
                            double* const distance_rescaled,
                            const int64_t size_max);

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
elec_num	int64_t	in	Number of electrons
nucl_num	int64_t	in	Number of nuclei
type_nucl_num	int64_t	in	Number of types of nuclei
type_nucl_vector	int64_t[nucl_num]	in	Number of types of nuclei
rescale_factor_en	double[type_nucl_num]	in	The factor for rescaled distances
walk_num	int64_t	in	Number of walkers
single_en_distance	double[walk_num][nucl_num]	in	Single electron-nucleus distances
en_rescaled_single	double[walk_num][nucl_num]	out	Electron-nucleus rescaled distances

function qmckl_compute_en_rescaled_single_doc(context, &
nucl_num, type_nucl_num, type_nucl_vector, rescale_factor_en, &
walk_num, single_en_distance, en_rescaled_single) &
bind(C) result(info)
use qmckl
implicit none
integer (qmckl_context), intent(in), value :: context
integer (c_int64_t) , intent(in)  , value  :: nucl_num
integer (c_int64_t) , intent(in)  , value  :: type_nucl_num
integer (c_int64_t) , intent(in)           :: type_nucl_vector(nucl_num)
real    (c_double ) , intent(in)           :: rescale_factor_en(type_nucl_num)
integer (c_int64_t) , intent(in)  , value  :: walk_num
real    (c_double ) , intent(in)           :: single_en_distance(nucl_num,walk_num)
real    (c_double ) , intent(out)          :: en_rescaled_single(nucl_num,walk_num)
integer(qmckl_exit_code)                   :: info

integer*8 :: i, k
double precision      :: coord(3)

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif


if (nucl_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif

if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif

do i=1, nucl_num
do k=1,walk_num
   en_rescaled_single(i,k) = (1.0d0 - dexp(-rescale_factor_en(type_nucl_vector(i)+1) * &
        single_en_distance(i,k))) / rescale_factor_en(type_nucl_vector(i)+1)
end do
end do

end function qmckl_compute_en_rescaled_single_doc

Single electron-nucleus Jastrow value

Calculates $\delta J_{en}$ \[ \delta J_{en} = \sum_i \delta_{i,\text{num}} \left( \frac{a_1 \widetilde{R}_{i\alpha}^{\text{new}}}{1+a_2 \widetilde{R}_{i\alpha}^{\text{new}}} - \frac{a_1 \widetilde{R}_{i\alpha}^{\text{old}}}{1+a_2 \widetilde{R}_{i\alpha}^{\text{old}}} + \sum_{k=2}^{N_\text{aord}} a_k\left({\widetilde{R}_{i\alpha}^{\text{new}}}^k - {\widetilde{R}_{i\alpha}^{\text{old}}}^k \right)\right) \]

Get

qmckl_exit_code
qmckl_get_jastrow_champ_single_en(qmckl_context context,
                       double* const delta_en,
                       const int64_t size_max);

interface
  integer(qmckl_exit_code) function qmckl_get_jastrow_champ_single_en (context, &
       delta_en, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (qmckl_context) , intent(in), value :: context
    integer(c_int64_t), intent(in), value       :: size_max
    real(c_double),   intent(out)               :: delta_en(size_max)
  end function qmckl_get_jastrow_champ_single_en
end interface

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
num	int64_t	in	Index of single point
walk_num	int64_t	in	Number of walkers
elec_num	int64_t	in	Number of electrons
nucl_num	int64_t	in	Number of nuclei
type_nucl_num	int64_t	in	Number of unique nuclei
type_nucl_vector	int64_t[nucl_num]	in	IDs of unique nuclei
aord_num	int64_t	in	Number of coefficients
a_vector	double[type_nucl_num][aord_num+1]	in	List of coefficients
en_distance_rescaled	double[walk_num][nucl_num][elec_num]	in	Electron-nucleus rescaled distances
~en_rescaled_single ~	double[walk_num][nucl_num]	in	Electron-nucleus rescaled single distances
delta_en	double[walk_num]	out	Single electron-nucleus jastrow

function qmckl_compute_jastrow_champ_single_en_doc( &
context, num_in, walk_num, elec_num, nucl_num, type_nucl_num, &
type_nucl_vector, aord_num, a_vector, &
en_distance_rescaled, en_rescaled_single, delta_en) &
bind(C) result(info)
use qmckl
implicit none

integer (qmckl_context), intent(in), value :: context
integer (c_int64_t) , intent(in)  , value :: num_in
integer (c_int64_t) , intent(in)  , value :: walk_num
integer (c_int64_t) , intent(in)  , value :: elec_num
integer (c_int64_t) , intent(in)  , value :: nucl_num
integer (c_int64_t) , intent(in)  , value :: type_nucl_num
integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
integer (c_int64_t) , intent(in)  , value :: aord_num
real    (c_double ) , intent(in)          :: a_vector(aord_num+1,type_nucl_num)
real    (c_double ) , intent(in)          :: en_distance_rescaled(elec_num,nucl_num,walk_num)
real    (c_double ) , intent(in)          :: en_rescaled_single(nucl_num,walk_num)
real    (c_double ) , intent(out)         :: delta_en(walk_num)
integer(qmckl_exit_code)                  :: info

integer*8 :: i, a, p, nw, num
double precision   :: x, power_ser, y

num = num_in + 1

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif

if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif

if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif

if (nucl_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif

if (type_nucl_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif

if (aord_num < 0) then
info = QMCKL_INVALID_ARG_7
return
endif


do nw =1, walk_num
delta_en(nw) = 0.0d0
do a = 1, nucl_num
     x = en_distance_rescaled(num, a, nw)
     y = en_rescaled_single(a, nw)

     delta_en(nw) = delta_en(nw) - a_vector(1, type_nucl_vector(a)+1) * x / (1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x)
     delta_en(nw) = delta_en(nw) + a_vector(1, type_nucl_vector(a)+1) * y / (1.0d0 + a_vector(2, type_nucl_vector(a)+1) * y)

     do p = 2, aord_num
       x = x * en_distance_rescaled(num, a, nw)
       y = y * en_rescaled_single(a, nw)
       delta_en(nw) = delta_en(nw) - a_vector(p + 1, type_nucl_vector(a)+1) * x + a_vector(p + 1, type_nucl_vector(a)+1) * y
     end do

end do
end do

end function qmckl_compute_jastrow_champ_single_en_doc

Electron-nucleus rescaled distances derivatives

Calculate the gradients and Laplacians of the rescaled electron-nucleus distances \[ \partial_{i,m} \widetilde{R}_{i\alpha} = \frac{r_{i,m} - R_{\alpha,m}}{R_{i\alpha}} e^{-\kappa R_{i\alpha}} \quad \text{and} \quad \partial_{i,4} \widetilde{R}_{i\alpha} = \left(\frac{2}{R_{i\alpha}} - \kappa \right)e^{-\kappa R_{i\alpha}} \]

Get

qmckl_exit_code qmckl_get_en_rescaled_single_gl(qmckl_context context,
                                           double* distance_rescaled_gl,
                                           const size_t size_max);

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
nucl_num	int64_t	in	Number of nuclei
type_nucl_num	int64_t	in	Number of nucleus types
type_nucl_vector	int64_t[nucl_num]	in	Array of nucleus types
rescale_factor_en	double[nucl_num]	in	The factors for rescaled distances
walk_num	int64_t	in	Number of walkers
single_en_distance	double[walk_num][nucl_num]	in	Single electorn distances
coord	double[walk_num][3]	in	Single electron coordinates
nucl_coord	double[3][nucl_num]	in	Nucleus coordinates
en_rescaled_single_gl	double[walk_num][nucl_num][4]	out	Electron-nucleus rescaled single distance derivatives

integer function qmckl_compute_en_rescaled_single_gl_doc_f(context, nucl_num, &
type_nucl_num, type_nucl_vector, rescale_factor_en, walk_num, &
single_en_distance, coord, nucl_coord, en_rescaled_single_gl) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: nucl_num
integer*8             , intent(in)  :: type_nucl_num
integer*8             , intent(in)  :: type_nucl_vector(nucl_num)
double precision      , intent(in)  :: rescale_factor_en(nucl_num)
integer*8             , intent(in)  :: walk_num
double precision      , intent(in)  :: single_en_distance(nucl_num, walk_num)
double precision      , intent(in)  :: coord(3,walk_num)
double precision      , intent(in)  :: nucl_coord(nucl_num,3)
double precision      , intent(out) :: en_rescaled_single_gl(4,nucl_num,walk_num)

integer*8 :: nw, a, ii
double precision :: ria_inv, elnuc_dist_gl(4, nucl_num), kappa_l

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif


if (nucl_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif

if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif


en_rescaled_single_gl = 0.0d0
do nw = 1, walk_num

! prepare the actual een table
do a = 1, nucl_num
   ria_inv = 1.0d0 / single_en_distance(a, nw)
   do ii = 1, 3
     elnuc_dist_gl(ii, a) = (coord(ii,nw) - nucl_coord(a, ii)) * ria_inv
   end do
   elnuc_dist_gl(4, a) = 2.0d0 * ria_inv
end do

 do a = 1, nucl_num
   kappa_l = -1 * rescale_factor_en(type_nucl_vector(a)+1)
   en_rescaled_single_gl(1, a, nw) = elnuc_dist_gl(1, a)
   en_rescaled_single_gl(2, a, nw) = elnuc_dist_gl(2, a)
   en_rescaled_single_gl(3, a, nw) = elnuc_dist_gl(3, a)
   en_rescaled_single_gl(4, a, nw) = elnuc_dist_gl(4, a)

   en_rescaled_single_gl(4, a, nw) = en_rescaled_single_gl(4, a, nw) + kappa_l

   en_rescaled_single_gl(1, a, nw) = en_rescaled_single_gl(1, a, nw)  * dexp(kappa_l * single_en_distance(a,nw))
   en_rescaled_single_gl(2, a, nw) = en_rescaled_single_gl(2, a, nw)  * dexp(kappa_l * single_en_distance(a,nw))
   en_rescaled_single_gl(3, a, nw) = en_rescaled_single_gl(3, a, nw)  * dexp(kappa_l * single_en_distance(a,nw))
   en_rescaled_single_gl(4, a, nw) = en_rescaled_single_gl(4, a, nw)  * dexp(kappa_l * single_en_distance(a,nw))
end do

end do


end function qmckl_compute_en_rescaled_single_gl_doc_f

Electron-nucleus Jastrow gradients and Laplacian

Calculate the gradients and Laplacian of the $\delta J_{en}$.

\begin{eqnarray*} \partial_{i,m}\delta J_{en} = & \sum_\alpha \delta_{i,\text{num}} \left( \frac{a_1 \partial_{i,m}\widetilde{R}^{\text{new}}_{i\alpha}}{(1+a_2 \widetilde{R}^{\text{new}}_{i\alpha})^2} - \frac{a_1 \partial_{i,m}\widetilde{R}^{\text{old}}_{i\alpha}}{(1+a_2 \widetilde{R}^{\text{old}}_{i\alpha})^2} + \sum_{k=2}^{N_{\text{aord}}} a_k k \left({\widetilde{R}_{i\alpha}^{\text{new}}}^{k-1} \partial_{i,m}\widetilde{R}^{\text{new}}_{i\alpha} - {\widetilde{R}_{i\alpha}^{\text{old}}}^{k-1} \partial_{i,m}\widetilde{R}^{\text{old}}_{i\alpha} \right) \right)\\ \partial_{i,4}\delta J_{en} = & \sum_\alpha \delta_{i,\text{num}} \left( \frac{a_1}{(1+a_2 \widetilde{R}^{\text{new}}_{i\alpha})^2} \left(\partial_{i,4}\widetilde{R}^{\text{new}}_{i\alpha} - 2 a_2 \frac{\nabla\widetilde{R}^{\text{new}}_{i\alpha} \cdot \nabla \widetilde{R}^{\text{new}}_{i\alpha}}{1+a_2 \widetilde{R}^{\text{new}}_{i\alpha}} \right) - \frac{a_1}{(1+a_2 \widetilde{R}^{\text{old}}_{i\alpha})^2} \left(\partial_{i,4}\widetilde{R}^{\text{old}}_{i\alpha} - 2 a_2 \frac{\nabla\widetilde{R}^{\text{old}}_{i\alpha} \cdot \nabla \widetilde{R}^{\text{old}}_{i\alpha}}{1+b_2 \widetilde{R}^{\text{old}}_{i\alpha}} \right) \right.\\ & \left.+ \sum_{k=2}^{N_{\text{aord}}} a_k k \left({\widetilde{R}_{i\alpha}^{\text{new}}}^{k-1}\partial_{i,4}\widetilde{R}_{i\alpha}^{\text{new}} + (k-1)(\nabla\widetilde{R}_{i\alpha}^{\text{new}} \cdot \nabla \widetilde{R}_{i\alpha}^{\text{new}}) {\widetilde{R}_{i\alpha}^{\text{new}}}^{k-2} - {\widetilde{R}_{i\alpha}^{\text{old}}}^{k-1}\partial_{i,4}\widetilde{R}_{i\alpha}^{\text{old}} - (k-1)(\nabla\widetilde{R}_{i\alpha}^{\text{old}} \cdot \nabla \widetilde{R}_{i\alpha}^{\text{old}}) {\widetilde{R}_{i\alpha}^{\text{old}}}^{k-2} \right) \right) \end{eqnarray*}

Get

qmckl_exit_code
qmckl_get_jastrow_champ_single_en_gl(qmckl_context context,
                               double* const delta_en_gl,
                               const int64_t size_max);

interface
  integer(qmckl_exit_code) function qmckl_get_jastrow_champ_single_en_gl (context, &
       delta_en_gl, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (qmckl_context) , intent(in), value :: context
    integer(c_int64_t), intent(in), value       :: size_max
    real(c_double),   intent(out)               :: delta_en_gl(size_max)
  end function qmckl_get_jastrow_champ_single_en_gl
end interface

Compute

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
num	int64_t	in	Index of single electron
walk_num	int64_t	in	Number of walkers
elec_num	int64_t	in	Number of electrons
nucl_num	int64_t	in	Number of nuclei
type_nucl_num	int64_t	in	Number of unique nuclei
type_nucl_vector	int64_t[nucl_num]	in	IDs of unique nuclei
aord_num	int64_t	in	Number of coefficients
a_vector	double[type_nucl_num][aord_num+1]	in	List of coefficients
en_distance_rescaled	double[walk_num][nucl_num][elec_num]	in	Electron-nucleus rescaled distances
en_distance_rescaled_gl	double[walk_num][nucl_num][elec_num][4]	in	Electron-nucleus rescaled distance derivatives
en_rescaled_single	double[walk_num][nucl_num]	in	Electron-nucleus rescaled single distances
en_rescaled_single_gl	double[walk_num][nucl_num][4]	in	Electron-nucleus rescaled single distance derivatives
delta_en_gl	double[walk_num][elec_num][4]	out	Single electron-nucleus Jastrow gradients and Laplacian

function qmckl_compute_jastrow_champ_single_en_gl_doc( &
context, num_in, walk_num, elec_num, nucl_num, type_nucl_num, &
type_nucl_vector, aord_num, a_vector, &
en_distance_rescaled, en_distance_rescaled_gl, en_rescaled_single, en_rescaled_single_gl, delta_en_gl) &
bind(C) result(info)
use qmckl
implicit none

integer (qmckl_context), intent(in), value :: context
integer (c_int64_t) , intent(in)  , value :: num_in
integer (c_int64_t) , intent(in)  , value :: walk_num
integer (c_int64_t) , intent(in)  , value :: elec_num
integer (c_int64_t) , intent(in)  , value :: nucl_num
integer (c_int64_t) , intent(in)  , value :: type_nucl_num
integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
integer (c_int64_t) , intent(in)  , value :: aord_num
real    (c_double ) , intent(in)          :: a_vector(aord_num+1,type_nucl_num)
real    (c_double ) , intent(in)          :: en_distance_rescaled(elec_num,nucl_num,walk_num)
real    (c_double ) , intent(in)          :: en_distance_rescaled_gl(4, elec_num,nucl_num,walk_num)
real    (c_double ) , intent(in)          :: en_rescaled_single(nucl_num,walk_num)
real    (c_double ) , intent(in)          :: en_rescaled_single_gl(4, nucl_num,walk_num)
real    (c_double ) , intent(out)         :: delta_en_gl(4,elec_num,walk_num)
integer(qmckl_exit_code)                   :: info

integer*8 :: i, a, k, nw, ii, num
double precision   :: x, x1, kf, x_old, x1_old
double precision   :: denom, invdenom, invdenom2, f
double precision   :: denom_old, invdenom_old, invdenom2_old, f_old
double precision   :: grad_c2, grad_c2_old
double precision   :: dx(4), dx_old(4)

num = num_in + 1

info = QMCKL_SUCCESS

if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif

if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif

if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif

if (nucl_num <= 0) then
info = QMCKL_INVALID_ARG_5
return
endif

if (aord_num < 0) then
info = QMCKL_INVALID_ARG_8
return
endif


do nw =1, walk_num
delta_en_gl(:,:,nw) = 0.0d0

do a = 1, nucl_num

     x_old = en_distance_rescaled(num,a,nw)
     x = en_rescaled_single(a,nw)


     denom = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x
     invdenom = 1.0d0 / denom
     invdenom2 = invdenom*invdenom

     denom_old = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x_old
     invdenom_old = 1.0d0 / denom_old
     invdenom2_old = invdenom_old*invdenom_old

     dx(1) = en_rescaled_single_gl(1,a,nw)
     dx(2) = en_rescaled_single_gl(2,a,nw)
     dx(3) = en_rescaled_single_gl(3,a,nw)
     dx(4) = en_rescaled_single_gl(4,a,nw)

     dx_old(1) = en_distance_rescaled_gl(1,num,a,nw)
     dx_old(2) = en_distance_rescaled_gl(2,num,a,nw)
     dx_old(3) = en_distance_rescaled_gl(3,num,a,nw)
     dx_old(4) = en_distance_rescaled_gl(4,num,a,nw)

     f = a_vector(1, type_nucl_vector(a)+1) * invdenom2
     grad_c2 = dx(1)*dx(1) + dx(2)*dx(2) + dx(3)*dx(3)

     f_old = a_vector(1, type_nucl_vector(a)+1) * invdenom2_old
     grad_c2_old = dx_old(1)*dx_old(1) + dx_old(2)*dx_old(2) + dx_old(3)*dx_old(3)

     delta_en_gl(1,num,nw) = delta_en_gl(1,num,nw) + f * dx(1) - f_old * dx_old(1)
     delta_en_gl(2,num,nw) = delta_en_gl(2,num,nw) + f * dx(2) - f_old * dx_old(2)
     delta_en_gl(3,num,nw) = delta_en_gl(3,num,nw) + f * dx(3) - f_old * dx_old(3)
     delta_en_gl(4,num,nw) = delta_en_gl(4,num,nw) &
          + f * (dx(4) - 2.d0 * a_vector(2, type_nucl_vector(a)+1) * grad_c2 * invdenom) &
          - f_old * (dx_old(4) - 2.d0 * a_vector(2, type_nucl_vector(a)+1) * grad_c2_old * invdenom_old)


      kf = 2.d0
      x1 = x
      x = 1.d0
      x1_old = x_old
      x_old = 1.d0
      do k=2, aord_num
         f = a_vector(k+1,type_nucl_vector(a)+1) * kf * x
         f_old = a_vector(k+1,type_nucl_vector(a)+1) * kf * x_old
         delta_en_gl(1,num,nw) = delta_en_gl(1,num,nw) + f * x1 * dx(1) - f_old * x1_old * dx_old(1)
         delta_en_gl(2,num,nw) = delta_en_gl(2,num,nw) + f * x1 * dx(2) - f_old * x1_old * dx_old(2)
         delta_en_gl(3,num,nw) = delta_en_gl(3,num,nw) + f * x1 * dx(3) - f_old * x1_old * dx_old(3)
         delta_en_gl(4,num,nw) = delta_en_gl(4,num,nw) &
              + f * (x1 * dx(4) + (kf-1.d0) * grad_c2) &
              - f_old * (x1_old * dx_old(4) + (kf-1.d0) * grad_c2_old)
         x = x*x1
         x_old = x_old*x1_old
         kf = kf + 1.d0
      end do
end do
end do

end function qmckl_compute_jastrow_champ_single_en_gl_doc

Accept single electron move

This function 'accepts' a proposed single-electron move. It overrides (or adds to) the required arrays and matrices to bring the positions of the electron and the full Jastrow factor up to date. For instances,

\begin{eqnarray} J_{en} = & J_{en} + \delta J_{en}\\ J_{ee} = & J_{ee} + \delta J_{ee}\\ J_{een} = & J_{een} + \delta J_{een}\\ \partial_{i,m} J_{en} = & \partial_{i,m} J_{en} + \partial_{i,m} \delta J_{en}\\ \partial_{i,m} J_{ee} = & \partial_{i,m} J_{ee} + \partial_{i,m} \delta J_{ee}\\ \partial_{i,m} J_{een} = & \partial_{i,m} J_{een} + \partial_{i,m} \delta J_{een}\\ \end{eqnarray}

The function has similar functionaly for the intermediate products, such as $P$, $\widetilde{r}$, $\widetilde{R}$, and the derivative of these. Furthermore, it also updates the dates, such that immidiatly calling a Jastrow function after an accepted single-electron move will not trigger a re-calculation of the Jastrow factors.

Code

qmckl_exit_code
qmckl_get_jastrow_champ_single_accept(qmckl_context context);

qmckl_exit_code
qmckl_get_jastrow_champ_single_accept(qmckl_context context)
{

if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}

qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);

qmckl_exit_code rc;

rc = qmckl_provide_jastrow_champ_single_en_gl(context);
if (rc != QMCKL_SUCCESS) return rc;

rc = qmckl_provide_jastrow_champ_single_ee_gl(context);
if (rc != QMCKL_SUCCESS) return rc;

if (ctx->jastrow_champ.cord_num > 0) {
rc = qmckl_provide_jastrow_champ_single_een_gl(context);
if (rc != QMCKL_SUCCESS) return rc;
}

double metric[4] = {-1.0, -1.0, -1.0, 1.0};

int shift1, shift2, shift3, shift4, shift5, shift6, shift7;

int shift8, shift9, shift10, shift11, shift12, shift13, shift14;

shift2 = ctx->electron.num*ctx->electron.num;
shift5 = ctx->electron.num*4*ctx->electron.num;
shift6 = ctx->electron.num*4;
shift11 = ctx->nucleus.num*4*ctx->electron.num;
shift13 = ctx->nucleus.num*4;
shift14 = ctx->nucleus.num*ctx->electron.num;


if (ctx->jastrow_champ.cord_num > 0) {

shift1 = (ctx->jastrow_champ.cord_num+1)*ctx->electron.num*ctx->electron.num;

shift3 = (ctx->jastrow_champ.cord_num+1)*ctx->electron.num;

shift4 = (ctx->jastrow_champ.cord_num+1)*ctx->electron.num*4*ctx->electron.num;

shift7 = (ctx->jastrow_champ.cord_num+1)*ctx->electron.num*4;

shift8 = (ctx->jastrow_champ.cord_num+1)*ctx->nucleus.num*ctx->electron.num;
shift9 = (ctx->jastrow_champ.cord_num+1)*ctx->nucleus.num;

shift10 = (ctx->jastrow_champ.cord_num+1)*ctx->nucleus.num*4*ctx->electron.num;

shift12 = (ctx->jastrow_champ.cord_num+1)*ctx->nucleus.num*4;

for (int nw = 0; nw < ctx->electron.walker.num; nw++) {
 ctx->jastrow_champ.factor_een[nw] = ctx->jastrow_champ.factor_een[nw] + ctx->single_point.delta_een[nw];
 for (int l = 0; l <= ctx->jastrow_champ.cord_num; l++){
   for (int i = 0; i < ctx->electron.num; i++) {
     ctx->jastrow_champ.een_rescaled_e[nw*shift1
       + l*shift2 
       + i*ctx->electron.num 
       + ctx->single_point.num] = 
     ctx->single_point.een_rescaled_single_e[nw*shift3
       + l*ctx->electron.num
       + i];

     ctx->jastrow_champ.een_rescaled_e[nw*shift1
       + l*shift2
       + ctx->single_point.num*ctx->electron.num 
       + i] = 
     ctx->single_point.een_rescaled_single_e[nw*shift3
       + l*ctx->electron.num
       + i];
     for (int k = 0; k < 4; k++){
       ctx->jastrow_champ.een_rescaled_e_gl[nw*shift4
         + l*shift5
         + i*shift6
         + k*ctx->electron.num 
         + ctx->single_point.num] = 
       ctx->single_point.een_rescaled_single_e_gl[nw*shift7
         + l*shift6
         + i*4 
         + k]; 
       ctx->jastrow_champ.een_rescaled_e_gl[nw*shift4
         + l*shift5
         + ctx->single_point.num*shift6
         + k*ctx->electron.num 
         + i] = 
       metric[k] * ctx->single_point.een_rescaled_single_e_gl[nw*shift7
         + l*shift6
         + i*4
         + k]; 
     }
   }
   for (int a = 0; a < ctx->nucleus.num; a++){
     ctx->jastrow_champ.een_rescaled_n[nw*shift8
       + l*shift14
       + a*ctx->electron.num
       + ctx->single_point.num] = 
     ctx->single_point.een_rescaled_single_n[nw*shift9
       + l*ctx->nucleus.num
       + a];
     for (int k = 0; k < 4; k++){
       ctx->jastrow_champ.een_rescaled_n_gl[nw*shift10
         + l*shift11
         + a*shift6
         + k*ctx->electron.num
         + ctx->single_point.num] = 
       ctx->single_point.een_rescaled_single_n_gl[nw*shift12
         + l*shift13
         + a*4
         + k];
     }
   }
 }
}
for (int i = 0; i < ctx->electron.walker.num * 4 * ctx->electron.num; i++) {
 ctx->jastrow_champ.factor_een_gl[i] = ctx->jastrow_champ.factor_een_gl[i] + ctx->single_point.delta_een_gl[i];
}  
for (int i = 0; i < (ctx->electron.walker.num*(ctx->jastrow_champ.cord_num+1)*ctx->nucleus.num*ctx->electron.num*ctx->jastrow_champ.cord_num); i++) {
ctx->jastrow_champ.tmp_c[i] = ctx->jastrow_champ.tmp_c[i] + ctx->single_point.delta_p[i];
}

/*
for (int nw = 0; nw < ctx->electron.walker.num; nw++) {
 for (int m = 0; m < ctx->jastrow_champ.cord_num; m++){
   for (int l = 0; l <= ctx->jastrow_champ.cord_num; l++) {
     for (int a = 0; a < ctx->nucleus.num; a++) {
       for (int k = 0; k < 4; k++){
         for (int i = 0; i < ctx->electron.num; i++) {
           ctx->jastrow_champ.dtmp_c[nw*ctx->electron.num*4*ctx->nucleus.num*ctx->jastrow_champ.cord_num*(ctx->jastrow_champ.cord_num+1)
             + m*ctx->electron.num*4*ctx->nucleus.num*(ctx->jastrow_champ.cord_num+1)
             + l*ctx->electron.num*4*ctx->nucleus.num
             + a*4*ctx->electron.num
             + k*ctx->electron.num
             + i] =
           ctx->jastrow_champ.dtmp_c[nw*ctx->electron.num*4*ctx->nucleus.num*ctx->jastrow_champ.cord_num*(ctx->jastrow_champ.cord_num+1)
             + m*ctx->electron.num*4*ctx->nucleus.num*(ctx->jastrow_champ.cord_num+1)
             + l*ctx->electron.num*4*ctx->nucleus.num
             + a*4*ctx->electron.num
             + k*ctx->electron.num
             + i] +
           ctx->single_point.delta_p_gl[nw*ctx->electron.num*4*ctx->nucleus.num*ctx->jastrow_champ.cord_num*(ctx->jastrow_champ.cord_num+1)
             + m*ctx->electron.num*4*ctx->nucleus.num*(ctx->jastrow_champ.cord_num+1)
             + l*ctx->electron.num*4*ctx->nucleus.num
             + k*ctx->electron.num*ctx->nucleus.num
             + a*ctx->electron.num
             + i];
         }
       }
     }
   }
 }
}

*/

for (int nw = 0; nw < ctx->electron.walker.num*(ctx->jastrow_champ.cord_num+1)*ctx->jastrow_champ.cord_num; nw++) {
 for (int a = 0; a < ctx->nucleus.num; a++) {
   for (int k = 0; k < 4; k++){
     for (int i = 0; i < ctx->electron.num; i++) {
       ctx->jastrow_champ.dtmp_c[nw*shift11
         + a*shift6
         + k*ctx->electron.num
         + i] =
       ctx->jastrow_champ.dtmp_c[nw*shift11
         + a*shift6
         + k*ctx->electron.num
         + i] +
       ctx->single_point.delta_p_gl[nw*shift11
         + k*shift14
         + a*ctx->electron.num
         + i];
     }
   }
 }
}
}


for (int nw = 0; nw < ctx->electron.walker.num; nw++) {
ctx->jastrow_champ.factor_en[nw] = ctx->jastrow_champ.factor_en[nw] + ctx->single_point.delta_en[nw];
ctx->jastrow_champ.factor_ee[nw] = ctx->jastrow_champ.factor_ee[nw] + ctx->single_point.delta_ee[nw];
for (int a = 0; a < ctx->nucleus.num; a++) {
 ctx->electron.en_distance[nw*shift14
   + ctx->single_point.num*ctx->nucleus.num
   + a] =
 ctx->single_point.single_en_distance[nw*ctx->nucleus.num
   + a];

 ctx->jastrow_champ.en_distance_rescaled[nw*shift14
   + a*ctx->electron.num
   + ctx->single_point.num] =
 ctx->single_point.en_rescaled_single[nw*ctx->nucleus.num
   + a];
 for (int k = 0; k < 4; k++){
   ctx->jastrow_champ.en_distance_rescaled_gl[nw*shift11
     + a*shift6
     + ctx->single_point.num*4
     + k] = 
   ctx->single_point.en_rescaled_single_gl[nw*shift13
     + a*4
     + k];
 }
}
for (int i = 0; i < ctx->electron.num; i++) {
 ctx->jastrow_champ.ee_distance_rescaled[nw*shift2
   + i*ctx->electron.num
   + ctx->single_point.num] =
 ctx->single_point.ee_rescaled_single[nw*ctx->electron.num
   + i];

 ctx->jastrow_champ.ee_distance_rescaled[nw*shift2
   + ctx->single_point.num*ctx->electron.num
   + i] =
 ctx->single_point.ee_rescaled_single[nw*ctx->electron.num
   + i];

 ctx->electron.ee_distance[nw*shift2
   + i*ctx->electron.num
   + ctx->single_point.num] =
 ctx->single_point.single_ee_distance[nw*ctx->electron.num
   + i];

 ctx->electron.ee_distance[nw*shift2
   + ctx->single_point.num*ctx->electron.num
   + i] =
 ctx->single_point.single_ee_distance[nw*ctx->electron.num
   + i];

 for (int k = 0; k < 4; k++){
   ctx->jastrow_champ.ee_distance_rescaled_gl[nw*shift5
     + i*shift6
     + ctx->single_point.num*4 
     + k] = 
   metric[k] * ctx->single_point.ee_rescaled_single_gl[nw*shift6
     + i*4 
     + k];
   ctx->jastrow_champ.ee_distance_rescaled_gl[nw*shift5
     + ctx->single_point.num*shift6
     + i*4 
     + k] =  
   ctx->single_point.ee_rescaled_single_gl[nw*shift6
     + i*4 
     + k];
 }
}
for (int k = 0; k < 4; k++){
 for (int i = 0; i < ctx->electron.num; i++) {

   ctx->jastrow_champ.factor_ee_gl[nw*shift6
     + k*ctx->electron.num
     + i] = 
   ctx->jastrow_champ.factor_ee_gl[nw*shift6
     + k*ctx->electron.num
     + i] + 
   ctx->single_point.delta_ee_gl[nw*shift6
     + i*4
     + k];

   ctx->jastrow_champ.factor_en_gl[nw*shift6
     + k*ctx->electron.num
     + i] = 
   ctx->jastrow_champ.factor_en_gl[nw*shift6
     + k*ctx->electron.num
     + i] + 
   ctx->single_point.delta_en_gl[nw*shift6
     + i*4
     + k];
 }
}
}

for (int nw = 0; nw < ctx->electron.walker.num; nw++) {
for (int k = 0; k < 3; k++) {
 ctx->point.coord.data[k*ctx->electron.walker.num*ctx->electron.num + nw*ctx->electron.num + ctx->single_point.num] = ctx->single_point.coord.data[nw*3 + k];
}
}

rc = qmckl_context_touch(context);
if (rc != QMCKL_SUCCESS) return rc;


if (ctx->jastrow_champ.cord_num > 0){
ctx->jastrow_champ.dtmp_c_date = ctx->date;
ctx->jastrow_champ.een_rescaled_e_date = ctx->date;
ctx->jastrow_champ.een_rescaled_e_gl_date = ctx->date;
ctx->jastrow_champ.een_rescaled_n_date = ctx->date;
ctx->jastrow_champ.een_rescaled_n_gl_date = ctx->date;
ctx->jastrow_champ.factor_een_date = ctx->date;
ctx->jastrow_champ.factor_een_gl_date = ctx->date;
ctx->jastrow_champ.tmp_c_date = ctx->date;

}
ctx->jastrow_champ.ee_distance_rescaled_date = ctx->date;
ctx->jastrow_champ.ee_distance_rescaled_gl_date = ctx->date;
ctx->jastrow_champ.en_distance_rescaled_date = ctx->date;
ctx->jastrow_champ.en_distance_rescaled_gl_date = ctx->date;
ctx->jastrow_champ.factor_ee_date = ctx->date;
ctx->jastrow_champ.factor_ee_gl_date = ctx->date;
ctx->jastrow_champ.factor_en_date = ctx->date;
ctx->jastrow_champ.factor_en_gl_date = ctx->date;

ctx->electron.ee_distance_date = ctx->date;
ctx->electron.en_distance_date = ctx->date;

ctx->single_point.date = ctx->date;

if (ctx->jastrow_champ.cord_num > 0){
ctx->single_point.een_rescaled_single_e_date = ctx->single_point.date;
ctx->single_point.een_rescaled_single_n_date = ctx->single_point.date;
ctx->single_point.delta_een_date = ctx->single_point.date;
ctx->single_point.delta_p_date = ctx->single_point.date;
ctx->single_point.een_rescaled_single_e_gl_date = ctx->single_point.date;
ctx->single_point.een_rescaled_single_n_gl_date = ctx->single_point.date;
ctx->single_point.delta_p_gl_date = ctx->single_point.date;
ctx->single_point.delta_een_gl_date = ctx->single_point.date;
ctx->single_point.delta_p_g_date = ctx->single_point.date;
ctx->single_point.delta_een_g_date = ctx->single_point.date;
ctx->forces.forces_delta_p_date = ctx->single_point.date;
ctx->forces.forces_jastrow_single_een_date = ctx->single_point.date;
}
ctx->single_point.single_ee_distance_date = ctx->single_point.date;
ctx->single_point.single_en_distance_date = ctx->single_point.date;
ctx->single_point.ee_rescaled_single_date = ctx->single_point.date;
ctx->single_point.en_rescaled_single_date = ctx->single_point.date;
ctx->single_point.delta_en_date = ctx->single_point.date;
ctx->single_point.delta_ee_date = ctx->single_point.date;
ctx->single_point.ee_rescaled_single_gl_date = ctx->single_point.date;
ctx->single_point.en_rescaled_single_gl_date = ctx->single_point.date;
ctx->single_point.delta_en_gl_date = ctx->single_point.date;
ctx->single_point.delta_ee_gl_date = ctx->single_point.date;
ctx->forces.forces_jastrow_single_en_date = ctx->single_point.date;


return QMCKL_SUCCESS;
}

interface
  integer(qmckl_exit_code) function qmckl_get_jastrow_champ_single_accept (context) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (qmckl_context) , intent(in), value :: context
  end function qmckl_get_jastrow_champ_single_accept
end interface