TREX/qmckl

Fork 0

mirror of https://github.com/TREX-CoE/qmckl.git synced 2025-04-30 12:24:49 +02:00

Anthony Scemama 77c4eb9702 Debug flags OK

2024-12-12 10:23:58 +01:00

143 KiB

Raw Blame History

Forces

Introduction
Context
- Data structure
- Test
Finite-difference function
Force of en jastrow value
- Get
- Compute
- Test
Force of en jastrow gradient
- Get
- Compute
- Test
Force of en jastrow laplacien
- Get
- Compute
- Test
Force of tmp_c matrix
- Get
- Compute
- Test
Force of dtmp_c matrix
- Get
- Compute
- Test
Force of een jastrow value
- Get
- Compute
- Test
Force of een_rescaled_n_gl
- Get
- Compute
- Test
Force of een jastrow gradient
- Get
- Compute
- Test
Force of een jastrow laplacian
- Get
- Compute
- Test

Introduction

Context

Data structure

typedef struct qmckl_forces_struct{
 double  * restrict forces_jastrow_en;
 uint64_t  forces_jastrow_en_date;
 double  * restrict forces_jastrow_en_g;
 uint64_t  forces_jastrow_en_g_date;
 double  * restrict forces_jastrow_en_l;
 uint64_t  forces_jastrow_en_l_date;
 double  * restrict forces_tmp_c;
 uint64_t  forces_tmp_c_date;
 double  * restrict forces_dtmp_c;
 uint64_t  forces_dtmp_c_date;
 double  * restrict forces_een_n;
 uint64_t  forces_een_n_date;
 double  * restrict forces_jastrow_een;
 uint64_t  forces_jastrow_een_date;
 double  * restrict forces_jastrow_een_g;
 uint64_t  forces_jastrow_een_g_date;
 double  * restrict forces_jastrow_een_l;
 uint64_t  forces_jastrow_een_l_date;
} qmckl_forces_struct;

Test

/* Reference input data */
int64_t walk_num      = n2_walk_num;
int64_t elec_num      = n2_elec_num;
int64_t elec_up_num   = n2_elec_up_num;
int64_t elec_dn_num   = n2_elec_dn_num;
int64_t nucl_num      = n2_nucl_num;
double  rescale_factor_ee   = 0.6;
double  rescale_factor_en[2] = { 0.6, 0.6 };
double* elec_coord    = &(n2_elec_coord[0][0][0]);

const double*   nucl_charge   = n2_charge;
double*  nucl_coord    = &(n2_nucl_coord[0][0]);

/* Provide Electron data */

qmckl_exit_code rc;

assert(!qmckl_electron_provided(context));

rc = qmckl_check(context,
              qmckl_set_electron_num (context, elec_up_num, elec_dn_num)
              );
assert(rc == QMCKL_SUCCESS);

assert(qmckl_electron_provided(context));

rc = qmckl_check(context,
              qmckl_set_electron_coord (context, 'N', walk_num, elec_coord, walk_num*3*elec_num)
              );
assert(rc == QMCKL_SUCCESS);

double elec_coord2[walk_num*3*elec_num];

rc = qmckl_check(context,
              qmckl_get_electron_coord (context, 'N', elec_coord2, walk_num*3*elec_num)
              );
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*elec_num ; ++i) {
assert( elec_coord[i] == elec_coord2[i] );
}


/* Provide Nucleus data */

assert(!qmckl_nucleus_provided(context));

rc = qmckl_check(context,
              qmckl_set_nucleus_num (context, nucl_num)
              );
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_nucleus_provided(context));

double nucl_coord2[3*nucl_num];

rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);

rc = qmckl_check(context,
              qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num)
              );
assert(rc == QMCKL_SUCCESS);

assert(!qmckl_nucleus_provided(context));

rc = qmckl_check(context,
              qmckl_get_nucleus_coord (context, 'N', nucl_coord2, nucl_num*3)
              );
assert(rc == QMCKL_SUCCESS);
for (int64_t k=0 ; k<3 ; ++k) {
for (int64_t i=0 ; i<nucl_num ; ++i) {
 assert( nucl_coord[nucl_num*k+i] == nucl_coord2[3*i+k] );
}
}

rc = qmckl_check(context,
              qmckl_get_nucleus_coord (context, 'T', nucl_coord2, nucl_num*3)
              );
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<3*nucl_num ; ++i) {
assert( nucl_coord[i] == nucl_coord2[i] );
}

double nucl_charge2[nucl_num];

rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);

rc = qmckl_check(context,
              qmckl_set_nucleus_charge(context, nucl_charge, nucl_num)
              );
assert(rc == QMCKL_SUCCESS);

rc = qmckl_check(context,
              qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num)
              );
assert(rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i<nucl_num ; ++i) {
assert( nucl_charge[i] == nucl_charge2[i] );
}
assert(qmckl_nucleus_provided(context));

assert(qmckl_electron_provided(context));

int64_t type_nucl_num = n2_type_nucl_num;
int64_t* type_nucl_vector = &(n2_type_nucl_vector[0]);
int64_t aord_num = n2_aord_num;
int64_t bord_num = n2_bord_num;
int64_t cord_num = n2_cord_num;
double* a_vector = &(n2_a_vector[0][0]);
double* b_vector = &(n2_b_vector[0]);
double* c_vector = &(n2_c_vector[0][0]);
int64_t dim_c_vector=0;

assert(!qmckl_jastrow_champ_provided(context));

/* Set the data */
rc = qmckl_check(context,
              qmckl_set_jastrow_champ_spin_independent(context, 0)
              );

rc = qmckl_check(context,
              qmckl_set_jastrow_champ_aord_num(context, aord_num)
              );
rc = qmckl_check(context,
              qmckl_set_jastrow_champ_bord_num(context, bord_num)
              );
rc = qmckl_check(context,
              qmckl_set_jastrow_champ_cord_num(context, cord_num)
              );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
              qmckl_set_jastrow_champ_type_nucl_num(context, type_nucl_num)
              );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
              qmckl_set_jastrow_champ_type_nucl_vector(context, type_nucl_vector, nucl_num)
              );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
              qmckl_set_jastrow_champ_a_vector(context, a_vector,(aord_num+1)*type_nucl_num)
              );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
              qmckl_set_jastrow_champ_b_vector(context, b_vector,(bord_num+1))
              );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
              qmckl_get_jastrow_champ_dim_c_vector(context, &dim_c_vector)
              );
assert(rc == QMCKL_SUCCESS);
rc = qmckl_check(context,
              qmckl_set_jastrow_champ_c_vector(context, c_vector, dim_c_vector*type_nucl_num)
              );
assert(rc == QMCKL_SUCCESS);

double k_ee = 0.;
double k_en[2] = { 0., 0. };
rc = qmckl_check(context,
              qmckl_set_jastrow_champ_rescale_factor_en(context, rescale_factor_en, type_nucl_num)
              );
assert(rc == QMCKL_SUCCESS);

rc = qmckl_check(context,
              qmckl_set_jastrow_champ_rescale_factor_ee(context, rescale_factor_ee)
              );
assert(rc == QMCKL_SUCCESS);

rc = qmckl_check(context,
              qmckl_get_jastrow_champ_rescale_factor_ee (context, &k_ee)
              );
assert(rc == QMCKL_SUCCESS);
assert(k_ee == rescale_factor_ee);

rc = qmckl_check(context,
              qmckl_get_jastrow_champ_rescale_factor_en (context, &(k_en[0]), type_nucl_num)
              );
assert(rc == QMCKL_SUCCESS);
for (int i=0 ; i<type_nucl_num ; ++i) {
assert(k_en[i] == rescale_factor_en[i]);
}

/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));

double delta_x = 0.00001;
double coef[9] = { 1.0/280.0, -4.0/105.0, 1.0/5.0, -4.0/5.0, 0.0, 4.0/5.0, -1.0/5.0, 4.0/105.0, -1.0/280.0 };

Finite-difference function

We introduce here a general function to compute the derivatives of any quantity with respect to nuclear coordinates. using finite-differences.

Force of en jastrow value

Get

qmckl_exit_code
qmckl_get_forces_jastrow_en(qmckl_context context,
                               double* const forces_jastrow_en,
                               const int64_t size_max);

interface
  integer(qmckl_exit_code) function qmckl_get_forces_jastrow_en (context, &
       forces_jastrow_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)               :: forces_jastrow_en(size_max)
  end function qmckl_get_forces_jastrow_en
end interface

Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nuclei
`type_nucl_num`	`int64_t`	in	Number of unique nuclei
`type_nucl_vector`	`int64_t[nucl_num]`	in	IDs of unique nuclei
`aord_num`	`int64_t`	in	Number of coefficients
`a_vector`	`double[type_nucl_num][aord_num+1]`	in	List of coefficients
`en_distance_rescaled`	`double[walk_num][nucl_num][elec_num]`	in	Electron-nucleus distances
`en_distance_rescaled_gl`	`double[walk_num][nucl_num][elec_num][4]`	in	Electron-nucleus distance derivatives
`forces_jastrow_en`	`double[walk_num][nucl_num][3]`	out	Electron-nucleus forces

function qmckl_compute_forces_jastrow_en_doc( &
context, walk_num, elec_num, nucl_num, type_nucl_num, &
type_nucl_vector, aord_num, a_vector, &
en_distance_rescaled, en_distance_rescaled_gl, forces_jastrow_en) &
bind(C) result(info)
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 :: 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(out)         :: forces_jastrow_en(3,nucl_num,walk_num)
integer(qmckl_exit_code)                   :: info

integer*8 :: i, a, k, nw, ii
double precision   :: x, x1, kf
double precision   :: denom, invdenom, invdenom2, f
double precision   :: dx(3)

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 (aord_num < 0) then
info = QMCKL_INVALID_ARG_7
return
endif

do nw =1, walk_num
forces_jastrow_en(:,:,nw) = 0.0d0
do a = 1, nucl_num
  do i = 1, elec_num

     x = en_distance_rescaled(i,a,nw)
     if(abs(x) < 1.d-12) continue

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

     dx(1) = -en_distance_rescaled_gl(1,i,a,nw)
     dx(2) = -en_distance_rescaled_gl(2,i,a,nw)
     dx(3) = -en_distance_rescaled_gl(3,i,a,nw)

     f = a_vector(1, type_nucl_vector(a)+1) * invdenom2

     forces_jastrow_en(1,a,nw) = forces_jastrow_en(1,a,nw) + f * dx(1)
     forces_jastrow_en(2,a,nw) = forces_jastrow_en(2,a,nw) + f * dx(2)
     forces_jastrow_en(3,a,nw) = forces_jastrow_en(3,a,nw) + f * dx(3)


      kf = 2.d0
      x1 = x
      x = 1.d0
      do k=2, aord_num
         f = a_vector(k+1,type_nucl_vector(a)+1) * kf * x
         forces_jastrow_en(1,a,nw) = forces_jastrow_en(1,a,nw) + f * x1 * dx(1)
         forces_jastrow_en(2,a,nw) = forces_jastrow_en(2,a,nw) + f * x1 * dx(2)
         forces_jastrow_en(3,a,nw) = forces_jastrow_en(3,a,nw) + f * x1 * dx(3)
         x = x*x1
         kf = kf + 1.d0
      end do

   end do
end do
end do


end function qmckl_compute_forces_jastrow_en_doc

Test

printf("Forces Jastrow en\n");

/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));

rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);

double forces_jastrow_en[walk_num][nucl_num][3];
rc = qmckl_get_forces_jastrow_en(context, &forces_jastrow_en[0][0][0], 3*nucl_num*walk_num);
assert(rc == QMCKL_SUCCESS);

double finite_difference_force_en[walk_num][nucl_num][3];
rc = qmckl_finite_difference_deriv_n(context, delta_x, &qmckl_get_jastrow_champ_factor_en, &(finite_difference_force_en[0][0][0]), 1);

for (int nw = 0; nw < walk_num; nw++){
for (int a = 0; a < nucl_num; a++) {
 for (int k = 0; k < 3; k++){
   //printf("%.10f\t", finite_difference_force_en[nw][a][k]);
   //printf("%.10f\n", forces_jastrow_en[nw][a][k]);
 }
}
}
for (int nw = 0; nw < walk_num; nw++){
for (int a = 0; a < nucl_num; a++) {
 for (int k = 0; k < 3; k++){
   assert(fabs(finite_difference_force_en[nw][a][k] - forces_jastrow_en[nw][a][k]) < 1.e-8);
 }
}
}
printf("OK\n");

Force of en jastrow gradient

Get

qmckl_exit_code
qmckl_get_forces_jastrow_en_g(qmckl_context context,
                               double* const forces_jastrow_en_g,
                               const int64_t size_max);

interface
  integer(qmckl_exit_code) function qmckl_get_forces_jastrow_en_g (context, &
       forces_jastrow_en_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)               :: forces_jastrow_en_g(size_max)
  end function qmckl_get_forces_jastrow_en_g
end interface

Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nuclei
`type_nucl_num`	`int64_t`	in	Number of unique nuclei
`type_nucl_vector`	`int64_t[nucl_num]`	in	IDs of unique nuclei
`aord_num`	`int64_t`	in	Number of coefficients
`a_vector`	`double[type_nucl_num][aord_num+1]`	in	List of coefficients
`rescale_factor_en`	`double[type_nucl_num]`	in	Rescale factor for electron-nucleus
`en_distance`	`double[elec_num][nucl_num]`	in	Electron-nucleus distances
`en_distance_rescaled`	`double[walk_num][nucl_num][elec_num]`	in	Electron-nucleus distances
`en_distance_rescaled_gl`	`double[walk_num][nucl_num][elec_num][4]`	in	Electron-nucleus distance derivatives
`forces_jastrow_en_g`	`double[walk_num][nucl_num][3][elec_num][3]`	out	Electron-nucleus forces

function qmckl_compute_forces_jastrow_en_g_doc( &
context, walk_num, elec_num, nucl_num, type_nucl_num, &
type_nucl_vector, aord_num, a_vector, rescale_factor_en, en_distance, &
en_distance_rescaled, en_distance_rescaled_gl, forces_jastrow_en_g) &
bind(C) result(info)
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 :: 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)          :: rescale_factor_en(type_nucl_num)
real    (c_double ) , intent(in)          :: en_distance(nucl_num, elec_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(out)         :: forces_jastrow_en_g(3,elec_num,3,nucl_num,walk_num)
integer(qmckl_exit_code)                   :: info

integer*8 :: i, a, k, nw, ii, m,l
double precision   :: x, x1, kf
double precision   :: denom, invdenom, invdenom2, f, f2, expk, invdist
double precision   :: dx(4)

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 (aord_num < 0) then
info = QMCKL_INVALID_ARG_7
return
endif

do nw =1, walk_num
forces_jastrow_en_g(:,:,:,:,nw) = 0.0d0
do a = 1, nucl_num
 do i = 1, elec_num
   expk = dexp(rescale_factor_en(type_nucl_vector(a)+1) * en_distance(a,i))
   invdist = 1.d0 / en_distance(a,i)
   x = en_distance_rescaled(i,a,nw)
   if(abs(x) < 1.d-12) continue
   denom = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x
   invdenom = 1.0d0 / denom
   invdenom2 = invdenom*invdenom
   f = a_vector(1, type_nucl_vector(a)+1) * invdenom2

   do m = 1, 3
     dx(m) = en_distance_rescaled_gl(m,i,a,nw)
   end do

   do m = 1, 3
     do l = 1,3
       if (m == l) then
         forces_jastrow_en_g(m,i,l,a,nw) = forces_jastrow_en_g(m,i,l,a,nw) - f  * invdist / expk
       end if

       forces_jastrow_en_g(m,i,l,a,nw) = forces_jastrow_en_g(m,i,l,a,nw) + &
          f * dx(m) * dx(l) * invdist * expk
       forces_jastrow_en_g(m,i,l,a,nw) = forces_jastrow_en_g(m,i,l,a,nw) + 2.d0 * f * invdenom * &
         a_vector(2, type_nucl_vector(a)+1) * dx(m) * dx(l)
     end do
   end do

   kf = 2.d0
   x1 = x
   x = 1.d0
   do k=2, aord_num
     f = a_vector(k+1,type_nucl_vector(a)+1) * kf * x
     f2 = a_vector(k+1,type_nucl_vector(a)+1) * kf * x * (kf-1.d0)

     do m = 1, 3
       do l = 1, 3
         if (m == l) then
           forces_jastrow_en_g(m,i,l,a,nw) = forces_jastrow_en_g(m,i,l,a,nw) - f * x1 * invdist / expk
         end if
         forces_jastrow_en_g(m,i,l,a,nw) = forces_jastrow_en_g(m,i,l,a,nw) - f2 * dx(m) * dx(l) &
             + f * x1 * dx(m) * dx(l) * rescale_factor_en(type_nucl_vector(a)+1) * expk &
             + f * x1 * dx(m) * dx(l) * invdist * expk
       end do
     end do
     x = x*x1
     kf = kf + 1.d0
   end do
 end do
end do
end do



end function qmckl_compute_forces_jastrow_en_g_doc

Test

printf("Forces Jastrow en G\n");

/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));

rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);

double forces_jastrow_en_g[walk_num][nucl_num][3][elec_num][3];
rc = qmckl_get_forces_jastrow_en_g(context, &forces_jastrow_en_g[0][0][0][0][0], 3*3*nucl_num*walk_num*elec_num);
assert(rc == QMCKL_SUCCESS);

double finite_difference_force_en_g[walk_num][nucl_num][3][4][elec_num];
rc = qmckl_finite_difference_deriv_n(context, delta_x, &qmckl_get_jastrow_champ_factor_en_gl, &finite_difference_force_en_g[0][0][0][0][0], 4*elec_num);

for (int nw = 0; nw < walk_num; nw++){
for (int a = 0; a < nucl_num; a++) {
 for (int k = 0; k < 3; k++){
   for (int i = 0; i < elec_num; i++){
     for (int l = 1; l < 3; l++){
       //printf("finite_difference_force_en_g[%i][%i][%i][%i][%i] %+3.10f \n", nw,a,k,l,i,finite_difference_force_en_g[nw][a][k][l][i]);
       //printf("forces_jastrow_en_g         [%i][%i][%i][%i][%i] %+3.10f\n", nw,a,k,i,l,forces_jastrow_en_g[nw][a][k][i][l]);
     }
   }
 }
}
}
for (int nw = 0; nw < walk_num; nw++){
for (int a = 0; a < nucl_num; a++) {
 for (int k = 0; k < 3; k++){
   for (int i = 0; i < elec_num; i++){
     for (int l = 1; l < 3; l++){
       assert(fabs(finite_difference_force_en_g[nw][a][k][l][i] - forces_jastrow_en_g[nw][a][k][i][l]) < 1.e-8);
     }
   }
 }
}
}
printf("OK\n");

Force of en jastrow laplacien

Get

qmckl_exit_code
qmckl_get_forces_jastrow_en_l(qmckl_context context,
                               double* const forces_jastrow_en_l,
                               const int64_t size_max);

interface
  integer(qmckl_exit_code) function qmckl_get_forces_jastrow_en_l (context, &
       forces_jastrow_en_l, 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)               :: forces_jastrow_en_l(size_max)
  end function qmckl_get_forces_jastrow_en_l
end interface

Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nuclei
`type_nucl_num`	`int64_t`	in	Number of unique nuclei
`type_nucl_vector`	`int64_t[nucl_num]`	in	IDs of unique nuclei
`aord_num`	`int64_t`	in	Number of coefficients
`a_vector`	`double[type_nucl_num][aord_num+1]`	in	List of coefficients
`rescale_factor_en`	`double[type_nucl_num]`	in	Rescale factor for electron-nucleus
`en_distance`	`double[elec_num][nucl_num]`	in	Electron-nucleus distances
`en_distance_rescaled`	`double[walk_num][nucl_num][elec_num]`	in	Electron-nucleus distances
`en_distance_rescaled_gl`	`double[walk_num][nucl_num][elec_num][4]`	in	Electron-nucleus distance derivatives
`forces_jastrow_en_l`	`double[walk_num][nucl_num][3]`	out	Electron-nucleus forces

function qmckl_compute_forces_jastrow_en_l_doc( &
context, walk_num, elec_num, nucl_num, type_nucl_num, &
type_nucl_vector, aord_num, a_vector, rescale_factor_en, en_distance, &
en_distance_rescaled, en_distance_rescaled_gl, forces_jastrow_en_l) &
bind(C) result(info)
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 :: 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)          :: rescale_factor_en(type_nucl_num)
real    (c_double ) , intent(in)          :: en_distance(nucl_num, elec_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(out)         :: forces_jastrow_en_l(3,nucl_num,walk_num)
integer(qmckl_exit_code)                   :: info

integer*8 :: i, a, k, nw, ii, m,l
double precision   :: x, x1, kf
double precision   :: denom, invdenom, invdenom2, f, f2, expk, invdist
double precision   :: dx(4)

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 (aord_num < 0) then
info = QMCKL_INVALID_ARG_7
return
endif

do nw =1, walk_num
forces_jastrow_en_l(:,:,nw) = 0.0d0
do a = 1, nucl_num
 do i = 1, elec_num
   expk = dexp(rescale_factor_en(type_nucl_vector(a)+1) * en_distance(a,i))
   invdist = 1.d0 / en_distance(a,i)
   x = en_distance_rescaled(i,a,nw)
   if(abs(x) < 1.d-12) continue
   denom = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x
   invdenom = 1.0d0 / denom
   invdenom2 = invdenom*invdenom
   f = a_vector(1, type_nucl_vector(a)+1) * invdenom2

   do m = 1, 4
     dx(m) = en_distance_rescaled_gl(m,i,a,nw)
   end do

   !do m = 1, 3
   !  do l = 1,3
   !    if (m == l) then
   !      forces_jastrow_en_g(l,a,nw) = forces_jastrow_en_g(l,a,nw) - f  * invdist / expk
   !    end if

   !    forces_jastrow_en_g(l,a,nw) = forces_jastrow_en_g(l,a,nw) + f * dx(m) * dx(l) * invdist * expk
   !    forces_jastrow_en_g(l,a,nw) = forces_jastrow_en_g(l,a,nw) + 2.d0 * f * invdenom * &
   !      a_vector(2, type_nucl_vector(a)+1) * dx(m) * dx(l)
   !  end do
   !end do


   kf = 2.d0
   x1 = x
   x = 1.d0
   do k=2, aord_num
     f = a_vector(k+1,type_nucl_vector(a)+1) * kf * x
     do m = 1, 3
       forces_jastrow_en_l(m,a,nw) = forces_jastrow_en_l(m,a,nw) &
         - f * dx(m) * dx(4) * (kf-1.d0) &
         - f / x1 * (kf-1.d0) * (kf-2.d0) * dx(m) /expk /expk &
         + f * x1 * rescale_factor_en(type_nucl_vector(a)+1) * dx(m) * dx(4) * expk &
         + f * x1 * 2 * dx(m) * invdist * invdist &
         + 2 * f * (kf-1.d0) * dx(m) * rescale_factor_en(type_nucl_vector(a)+1) / expk
     end do
     x = x*x1
     kf = kf + 1.d0
   end do

 end do
end do
end do



end function qmckl_compute_forces_jastrow_en_l_doc

Test

printf("Forces Jastrow en L\n");

/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));

rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);

double forces_jastrow_en_l[walk_num][nucl_num][3];
rc = qmckl_get_forces_jastrow_en_l(context, &forces_jastrow_en_l[0][0][0], 3*nucl_num*walk_num);
assert(rc == QMCKL_SUCCESS);

double finite_difference_force_en_l[walk_num][nucl_num][3][4][elec_num];
rc = qmckl_finite_difference_deriv_n(context, delta_x, &qmckl_get_jastrow_champ_factor_en_gl, &finite_difference_force_en_l[0][0][0][0][0], 4*elec_num);


double finite_difference_force_en_l_sum[walk_num][nucl_num][3];
for (int nw = 0; nw < walk_num; nw++){
for (int a = 0; a < nucl_num; a++) {
 for (int k = 0; k < 3; k++){
   finite_difference_force_en_l_sum[nw][a][k] = 0;
   for (int i = 0; i < elec_num; i++){
     finite_difference_force_en_l_sum[nw][a][k] = finite_difference_force_en_l_sum[nw][a][k] + finite_difference_force_en_l[nw][a][k][3][i];
   }
 }
}
}

for (int nw = 0; nw < walk_num; nw++){
for (int a = 0; a < nucl_num; a++) {
 for (int k = 0; k < 3; k++){
   //printf("finite_difference_force_en_l_sum[%i][%i][%i] %+3.10f \n", nw,a,k,finite_difference_force_en_l_sum[nw][a][k]);
   //printf("forces_jastrow_en_l             [%i][%i][%i] %+3.10f\n", nw,a,k,forces_jastrow_en_l[nw][a][k]);
   
 }
}
}
for (int nw = 0; nw < walk_num; nw++){
for (int a = 0; a < nucl_num; a++) {
 for (int k = 0; k < 3; k++){
       assert(fabs(finite_difference_force_en_l_sum[nw][a][k] - forces_jastrow_en_l[nw][a][k]) < 1.e-8);
 }
}
}
printf("OK\n");

Force of tmp_c matrix

Get

qmckl_exit_code
qmckl_get_forces_tmp_c(qmckl_context context,
                                double* const forces_tmp_c,
                                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
`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_e`	`double[walk_num][0:cord_num][elec_num][elec_num]`	in	Electron-nucleus rescaled
`een_rescaled_n_gl`	`double[walk_num][0:cord_num][nucl_num][4][elec_num]`	in	Electron-nucleus rescaled factor
`forces_tmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]`	out	vector of non-zero coefficients

integer(qmckl_exit_code) function qmckl_compute_forces_tmp_c( &
context, walk_num, elec_num, nucl_num, cord_num,&
een_rescaled_e, een_rescaled_n_gl, forces_tmp_c) &
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   :: walk_num, elec_num, cord_num, nucl_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_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num)
real    (c_double )   , intent(out)         :: forces_tmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1, walk_num)

integer*8 :: nw, i

integer*8 :: l, m, k, a,j

info = QMCKL_SUCCESS

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

do nw=1, walk_num
do i=0, cord_num-1
 info = qmckl_dgemm(context,'N','N',elec_num*1_8,&
       nucl_num*(cord_num+1)*4_8, elec_num*1_8, -1.0d0,     &
       een_rescaled_e(1,1,i,nw),1_8*elec_num,              &
       een_rescaled_n_gl(1,1,1,0,nw),elec_num*1_8,         &
       0.0d0,                                              &
       forces_tmp_c(1,1,1,0,i,nw),1_8*elec_num)
end do
end do



end function qmckl_compute_forces_tmp_c

Test

printf("Forces tmp_c\n");

/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));

rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);

double forces_tmp_c[walk_num][cord_num][cord_num+1][nucl_num][4][elec_num];
rc = qmckl_get_forces_tmp_c(context, &forces_tmp_c[0][0][0][0][0][0], 4*nucl_num*walk_num*elec_num*(cord_num+1)*cord_num);
assert(rc == QMCKL_SUCCESS);


double finite_difference_force_tmp_c[walk_num][cord_num][cord_num+1][nucl_num][3][elec_num];

double* nucleus_coord = (double*) malloc(3 * nucl_num * sizeof(double));
if (nucleus_coord == NULL) {
return QMCKL_ALLOCATION_FAILED;
}

rc = qmckl_get_nucleus_coord(context, 'N', nucleus_coord, 3*nucl_num);

double* temp_coord = (double*) malloc(3 * nucl_num * sizeof(double));
if (temp_coord == NULL) {
free(nucleus_coord);
return QMCKL_ALLOCATION_FAILED;
}
double output[walk_num][cord_num][cord_num+1][nucl_num][elec_num];

// Copy original coordinates
for (int i = 0; i < 3 * nucl_num; i++) {
temp_coord[i] = nucleus_coord[i];
}



for (int64_t a = 0; a < nucl_num; a++) {
for (int64_t k = 0; k < 3; k++) {
 for (int64_t m = -4; m <= 4; m++) {
   
   // Apply finite difference displacement
   temp_coord[k+a*3] = nucleus_coord[k+3*a] + (double) m * delta_x;

   // Update coordinates in the context
   rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num);
   assert(rc == QMCKL_SUCCESS);
   
   rc = qmckl_context_touch(context);
   assert(rc == QMCKL_SUCCESS);
   
   // Call the provided function
   rc = qmckl_get_jastrow_champ_tmp_c(context,&output[0][0][0][0][0]);
   assert(rc == QMCKL_SUCCESS);

   // Accumulate derivative using finite-difference coefficients
   for (int nw=0 ; nw<walk_num ; nw++) {
     for (int l = 0; l < cord_num; l++) {
       for (int mm = 0; mm <= cord_num; mm++) {
         for (int i = 0; i < elec_num; i++) {

           if (m == -4) {
             finite_difference_force_tmp_c[nw][l][mm][a][k][i] = 0.0;
           }
           finite_difference_force_tmp_c[nw][l][mm][a][k][i] += coef[m + 4] * output[nw][l][mm][a][i]/delta_x;
         }
       }
     }
   }
 }
 temp_coord[k+a*3] = nucleus_coord[k+3*a];
}
}

// Reset coordinates in the context
rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_context_touch(context);
assert(rc == QMCKL_SUCCESS);
 
free(nucleus_coord);
free(temp_coord);


for (int nw = 0; nw < walk_num; nw++){
for (int l = 0; l < cord_num; l++){
 for (int m = 0; m <= cord_num; m++){
   for (int a = 0; a < nucl_num; a++) {
     for (int k = 0; k < 3; k++){
       for (int i = 0; i < elec_num; i++){
         //printf("nw=%i l=%i m=%i a=%i k=%i i=%i\n",nw,l,m,a,k,i);
         //printf("%.10f\t", finite_difference_force_tmp_c[nw][l][m][a][k][i]);
         //printf("%.10f\n", forces_tmp_c[nw][l][m][a][k][i]);
       
       }
     }
   }
 }
}
}
for (int nw = 0; nw < walk_num; nw++){
for (int l = 0; l < cord_num; l++){
 for (int m = 0; m <= cord_num; m++){
   for (int a = 0; a < nucl_num; a++) {
     for (int k = 0; k < 3; k++){
       for (int i = 0; i < elec_num; i++){
         assert(fabs(finite_difference_force_tmp_c[nw][l][m][a][k][i] - forces_tmp_c[nw][l][m][a][k][i]) < 1.e-8);
       }
     }
   }
 }
}
}


printf("OK\n");

Force of dtmp_c matrix

Get

qmckl_exit_code
qmckl_get_forces_dtmp_c(qmckl_context context,
                                double* const forces_dtmp_c,
                                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
`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_e_gl`	`double[walk_num][0:cord_num][elec_num][4][elec_num]`	in	Electron-nucleus rescaled
`een_rescaled_n_gl`	`double[walk_num][0:cord_num][nucl_num][4][elec_num]`	in	Electron-nucleus rescaled factor
`forces_dtmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][4][nucl_num][4][elec_num]`	out	vector of non-zero coefficients

integer(qmckl_exit_code) function qmckl_compute_forces_dtmp_c( &
context, walk_num, elec_num, nucl_num, cord_num,&
een_rescaled_e_gl, een_rescaled_n_gl, forces_dtmp_c) &
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   :: walk_num, elec_num, cord_num, nucl_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_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num)
real    (c_double )   , intent(out)         :: forces_dtmp_c(elec_num, 4, nucl_num,4,0:cord_num, 0:cord_num-1, walk_num)

integer*8 :: nw, i, k

integer*8 :: l, m, a,j, kk

info = QMCKL_SUCCESS

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

do nw=1, walk_num
do i=0, cord_num-1
do k =1, 4
 info = qmckl_dgemm(context,'N','N',elec_num*1_8,&
       nucl_num*(cord_num+1)*4_8, elec_num*1_8, -1.0d0,     &
       een_rescaled_e_gl(1,k,1,i,nw),4_8*elec_num,              &
       een_rescaled_n_gl(1,1,1,0,nw),elec_num*1_8,         &
       0.0d0,                                              &
       forces_dtmp_c(1,k,1,1,0,i,nw),4_8*elec_num)
 end do
end do
end do

if (.true.) then
do nw = 1, walk_num
do l = 0, cord_num-1
 do m = 0, cord_num
   do k = 1, 4
     do a = 1, nucl_num
       do kk = 1, 4
         do i = 1, elec_num
           forces_dtmp_c(i,kk,a,k,m,l,nw) = 0.0
           do j = 1, elec_num
             forces_dtmp_c(i,kk,a,k,m,l,nw) = forces_dtmp_c(i,kk,a,k,m,l,nw) - &
                   een_rescaled_e_gl(i,kk,j,l,nw) * een_rescaled_n_gl(j,k,a,m,nw)
           end do
         end do
       end do
     end do
   end do
 end do
end do
end do
end if

end function qmckl_compute_forces_dtmp_c

Test

printf("Forces dtmp_c\n");

/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));

rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);

double forces_dtmp_c[walk_num][cord_num][cord_num+1][4][nucl_num][4][elec_num];
rc = qmckl_get_forces_dtmp_c(context, &forces_dtmp_c[0][0][0][0][0][0][0], 4*4*nucl_num*walk_num*elec_num*(cord_num+1)*cord_num);
assert(rc == QMCKL_SUCCESS);


double finite_difference_force_dtmp_c[walk_num][cord_num][cord_num+1][3][nucl_num][4][elec_num];

nucleus_coord = (double*) malloc(3 * nucl_num * sizeof(double));
if (nucleus_coord == NULL) {
return QMCKL_ALLOCATION_FAILED;
}

rc = qmckl_get_nucleus_coord(context, 'N', nucleus_coord, 3*nucl_num);

temp_coord = (double*) malloc(3 * nucl_num * sizeof(double));
if (temp_coord == NULL) {
free(nucleus_coord);
return QMCKL_ALLOCATION_FAILED;
}
double doutput[walk_num][cord_num][cord_num+1][nucl_num][4][elec_num];

// Copy original coordinates
for (int i = 0; i < 3 * nucl_num; i++) {
temp_coord[i] = nucleus_coord[i];
}



for (int64_t a = 0; a < nucl_num; a++) {
for (int64_t k = 0; k < 3; k++) {
 for (int64_t m = -4; m <= 4; m++) {
   
   // Apply finite difference displacement
   temp_coord[k+a*3] = nucleus_coord[k+3*a] + (double) m * delta_x;

   // Update coordinates in the context
   rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num);
   assert(rc == QMCKL_SUCCESS);
   
   rc = qmckl_context_touch(context);
   assert(rc == QMCKL_SUCCESS);
   
   // Call the provided function
   rc = qmckl_get_jastrow_champ_dtmp_c(context,&doutput[0][0][0][0][0][0]);
   assert(rc == QMCKL_SUCCESS);

   // Accumulate derivative using finite-difference coefficients
   for (int nw=0 ; nw<walk_num ; nw++) {
     for (int l = 0; l < cord_num; l++) {
       for (int mm = 0; mm <= cord_num; mm++) {
         for (int i = 0; i < elec_num; i++) {
           for (int j = 0; j < 4; j++) {
             if (m == -4) {
               finite_difference_force_dtmp_c[nw][l][mm][k][a][j][i] = 0.0;
             }
             finite_difference_force_dtmp_c[nw][l][mm][k][a][j][i] += coef[m + 4] * doutput[nw][l][mm][a][j][i]/delta_x;
           }

         }
       }
     }
   }
 }
 temp_coord[k+a*3] = nucleus_coord[k+3*a];
}
}

// Reset coordinates in the context
rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_context_touch(context);
assert(rc == QMCKL_SUCCESS);
 
free(nucleus_coord);
free(temp_coord);


for (int nw = 0; nw < walk_num; nw++){
for (int l = 0; l < cord_num; l++){
 for (int m = 0; m <= cord_num; m++){
   for (int a = 0; a < nucl_num; a++) {
     for (int k = 0; k < 3; k++){
       for (int i = 0; i < elec_num; i++){
         for (int kk = 0; kk<4; kk++){
         //printf("nw=%i l=%i m=%i k=%i a=%i kk=%i i=%i\n",nw,l,m,k,a,kk,i);
         //printf("%.10f\t", finite_difference_force_dtmp_c[nw][l][m][k][a][kk][i]);
         //printf("%.10f\n", forces_dtmp_c[nw][l][m][k][a][kk][i]);
         //assert(fabs(finite_difference_force_dtmp_c[nw][l][m][k][a][kk][i] - forces_dtmp_c[nw][l][m][k][a][kk][i]) < 1.e-8);
         }
       
       }
     }
   }
 }
}
}
for (int nw = 0; nw < walk_num; nw++){
for (int l = 0; l < cord_num; l++){
 for (int m = 0; m <= cord_num; m++){
   for (int a = 0; a < nucl_num; a++) {
     for (int k = 0; k < 3; k++){
       for (int i = 0; i < elec_num; i++){
         for (int kk = 0; kk<4; kk++){
           assert(fabs(finite_difference_force_dtmp_c[nw][l][m][k][a][kk][i] - forces_dtmp_c[nw][l][m][k][a][kk][i]) < 1.e-8);
         }
       }
     }
   }
 }
}
}


printf("OK\n");

Force of een jastrow value

Get

qmckl_exit_code
qmckl_get_forces_jastrow_een(qmckl_context context,
                               double* const forces_jastrow_een,
                               const int64_t size_max);

interface
  integer(qmckl_exit_code) function qmckl_get_forces_jastrow_een (context, &
       forces_jastrow_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)               :: forces_jastrow_een(size_max)
  end function qmckl_get_forces_jastrow_een
end interface

Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nuclei
`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
`een_rescaled_e`	`double[walk_num][0:cord_num][elec_num][elec_num]`	in	Electron-nucleus rescaled
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	in	Electron-nucleus rescaled factor
`een_rescaled_n_gl`	`double[walk_num][0:cord_num][nucl_num][4][elec_num]`	in	Electron-nucleus rescaled factor
`tmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]`	in	vector of non-zero coefficients
`forces_tmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]`	in	vector of non-zero coefficients
`forces_jastrow_een`	`double[walk_num][nucl_num][3]`	out	Electron-nucleus jastrow

integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_een( &
context, walk_num, elec_num, nucl_num, cord_num,&
dim_c_vector, c_vector_full, lkpm_combined_index, &
een_rescaled_e, een_rescaled_n, een_rescaled_n_gl, tmp_c, forces_tmp_c,forces_jastrow_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  :: 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)  :: een_rescaled_e(elec_num, elec_num, 0:cord_num, 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_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num)
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)  :: forces_tmp_c(elec_num,4, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
real    (c_double )   , intent(out) :: forces_jastrow_een(3, nucl_num, walk_num)

integer*8 :: i, a, j, l, k, p, m, n, nw, ii
double precision :: accu, accu2, cn

info = QMCKL_SUCCESS

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

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 a = 1, nucl_num
   cn = c_vector_full(a, n)
   if(cn == 0.d0) cycle
     do j = 1, elec_num
       do ii = 1, 3
         accu = een_rescaled_n(j,a,m,nw) * forces_tmp_c(j,ii,a,m+l,k,nw) 
         accu = accu - een_rescaled_n_gl(j,ii,a,m,nw) * tmp_c(j,a,m+l,k,nw)

         forces_jastrow_een(ii, a, nw) = forces_jastrow_een(ii, a, nw) + accu * cn
       end do
     end do
 end do
end do
end do

end function qmckl_compute_forces_jastrow_een

Test

printf("Forces Jastrow een\n");

/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));

rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);

double forces_jastrow_een[walk_num][nucl_num][3];
rc = qmckl_get_forces_jastrow_een(context, &forces_jastrow_een[0][0][0], 3*nucl_num*walk_num);
assert(rc == QMCKL_SUCCESS);

double finite_difference_force_een[walk_num][nucl_num][3];
rc = qmckl_finite_difference_deriv_n(context, delta_x, &qmckl_get_jastrow_champ_factor_een, &(finite_difference_force_een[0][0][0]), 1);

for (int nw = 0; nw < walk_num; nw++){
for (int a = 0; a < nucl_num; a++) {
 for (int k = 0; k < 3; k++){
   //printf("%.10f\t", finite_difference_force_een[nw][a][k]);
   //printf("%.10f\n", forces_jastrow_een[nw][a][k]);
 }
}
}
for (int nw = 0; nw < walk_num; nw++){
for (int a = 0; a < nucl_num; a++) {
 for (int k = 0; k < 3; k++){
   assert(fabs(finite_difference_force_een[nw][a][k] - forces_jastrow_een[nw][a][k]) < 1.e-8);
 }
}
}
printf("OK\n");

Force of een_rescaled_n_gl

Get

qmckl_exit_code
qmckl_get_forces_een_rescaled_n_gl(qmckl_context context,
                                    double* const forces_een_n,
                                    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
`elec_num`	`int64_t`	in	Number of electrons
`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
`en_distance`	`double[walk_num][elec_num][nucl_num]`	in	Electron-nucleus distances
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	in	Electron-nucleus distances
`een_rescaled_n_gl`	`double[walk_num][0:cord_num][nucl_num][4][elec_num]`	in	Electron-nucleus distances
`forces_een_n`	`double[walk_num][0:cord_num][3][nucl_num][4][elec_num]`	out	Electron-nucleus rescaled distances

integer(qmckl_exit_code) function qmckl_compute_forces_een_rescaled_n_gl( &
context, walk_num, elec_num, nucl_num, type_nucl_num, type_nucl_vector, &
cord_num, rescale_factor_en, &
en_distance, een_rescaled_n, een_rescaled_n_gl,forces_een_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  :: walk_num, elec_num, nucl_num, type_nucl_num, cord_num
integer (c_int64_t)   , intent(in)  :: type_nucl_vector(nucl_num)
real    (c_double )   , intent(in)  :: rescale_factor_en(type_nucl_num)
real    (c_double )   , intent(in)  :: en_distance(nucl_num,elec_num,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_n_gl(elec_num,4,nucl_num,0:cord_num,walk_num) 
real    (c_double )   , intent(out) :: forces_een_n(elec_num,4,nucl_num,3,0:cord_num,walk_num)

double precision                    :: x, ria_inv, kappa_l
integer*8                           :: i, a, k, l, nw, m, n


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 (cord_num < 0) then
info = QMCKL_INVALID_ARG_5
return
endif



forces_een_n = 0.0d0
do nw = 1, walk_num
do i = 1, elec_num
 do a = 1, nucl_num
   do l = 0, cord_num
     kappa_l = dble(l)*rescale_factor_en(type_nucl_vector(a)+1)
     do m = 1, 4
       do n = 1, 3
         if (l == 0) then
           forces_een_n(i,m,a,n,l,nw) = 0.0d0
         else
           if (m == n) then
             forces_een_n(i,m,a,n,l,nw) = forces_een_n(i,m,a,n,l,nw) + &
                  kappa_l*een_rescaled_n(i,a,l,nw)/en_distance(a,i,nw)
           end if
           if (m < 4) then
           forces_een_n(i,m,a,n,l,nw) = forces_een_n(i,m,a,n,l,nw) - &
             een_rescaled_n_gl(i,m,a,l,nw) * een_rescaled_n_gl(i,n,a,l,nw) / &
             (kappa_l*een_rescaled_n(i,a,l,nw)*en_distance(a,i,nw)) - &
              een_rescaled_n_gl(i,m,a,l,nw) * &
              een_rescaled_n_gl(i,n,a,l,nw)/een_rescaled_n(i,a,l,nw)
           else
           forces_een_n(i,m,a,n,l,nw) = forces_een_n(i,m,a,n,l,nw) + &
             2.0d0 * een_rescaled_n_gl(i,n,a,l,nw) / &
             (en_distance(a,i,nw)*en_distance(a,i,nw)) - &
             een_rescaled_n_gl(i,m,a,l,nw) * &
             een_rescaled_n_gl(i,n,a,l,nw)/een_rescaled_n(i,a,l,nw)
           end if

         end if
       end do
     end do
   end do
 end do
end do
end do


end function qmckl_compute_forces_een_rescaled_n_gl

Test

printf("Forces of een_rescaled_n_gl\n");

/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));

rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);

double forces_een_n[walk_num][cord_num+1][3][nucl_num][4][elec_num];
rc = qmckl_get_forces_een_rescaled_n_gl(context, &forces_een_n[0][0][0][0][0][0], 3*4*nucl_num*walk_num*elec_num*(cord_num+1));
assert(rc == QMCKL_SUCCESS);


double finite_difference_force_een_n[walk_num][cord_num+1][3][nucl_num][4][elec_num];

nucleus_coord = (double*) malloc(3 * nucl_num * sizeof(double));
if (nucleus_coord == NULL) {
return QMCKL_ALLOCATION_FAILED;
}

rc = qmckl_get_nucleus_coord(context, 'N', nucleus_coord, 3*nucl_num);

temp_coord = (double*) malloc(3 * nucl_num * sizeof(double));
if (temp_coord == NULL) {
free(nucleus_coord);
return QMCKL_ALLOCATION_FAILED;
}
double ddoutput[walk_num][cord_num+1][nucl_num][4][elec_num];

// Copy original coordinates
for (int i = 0; i < 3 * nucl_num; i++) {
temp_coord[i] = nucleus_coord[i];
}



for (int64_t a = 0; a < nucl_num; a++) {
for (int64_t k = 0; k < 3; k++) {
 for (int64_t m = -4; m <= 4; m++) {
   
   // Apply finite difference displacement
   temp_coord[k+a*3] = nucleus_coord[k+3*a] + (double) m * delta_x;

   // Update coordinates in the context
   rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num);
   assert(rc == QMCKL_SUCCESS);
   
   rc = qmckl_context_touch(context);
   assert(rc == QMCKL_SUCCESS);
   
   // Call the provided function
   rc = qmckl_get_jastrow_champ_een_rescaled_n_gl(context,&ddoutput[0][0][0][0][0], 4*elec_num*nucl_num*(cord_num+1)*walk_num);
   assert(rc == QMCKL_SUCCESS);

   // Accumulate derivative using finite-difference coefficients
   for (int nw=0 ; nw<walk_num ; nw++) {
     for (int l = 0; l <= cord_num; l++) {
         for (int i = 0; i < elec_num; i++) {
           for (int j = 0; j < 4; j++) {
             if (m == -4) {
               finite_difference_force_een_n[nw][l][k][a][j][i] = 0.0;
             }
             finite_difference_force_een_n[nw][l][k][a][j][i] += coef[m + 4] * ddoutput[nw][l][a][j][i]/delta_x;
           }


       }
     }
   }
 }
 temp_coord[k+a*3] = nucleus_coord[k+3*a];
}
}

// Reset coordinates in the context
rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_context_touch(context);
assert(rc == QMCKL_SUCCESS);
 
free(nucleus_coord);
free(temp_coord);


for (int nw = 0; nw < walk_num; nw++){
for (int l = 0; l <=cord_num; l++){
   for (int a = 0; a < nucl_num; a++) {
     for (int k = 0; k < 3; k++){
       for (int i = 0; i < elec_num; i++){
         for (int kk = 0; kk<4; kk++){
         //printf("nw=%i l=%i k=%i a=%i kk=%i i=%i\n",nw,l,k,a,kk,i);
         //printf("%.10f\t", finite_difference_force_een_n[nw][l][k][a][kk][i]);
         //printf("%.10f\n", forces_een_n[nw][l][k][a][kk][i]);
         //assert(fabs(finite_difference_force_een_n[nw][l][k][a][kk][i] - forces_een_n[nw][l][k][a][kk][i]) < 1.e-6);
         }
       
       }
     }
 }
}
}
for (int nw = 0; nw < walk_num; nw++){
for (int l = 0; l <= cord_num; l++){
   for (int a = 0; a < nucl_num; a++) {
     for (int k = 0; k < 3; k++){
       for (int i = 0; i < elec_num; i++){
         for (int kk = 0; kk<4; kk++){
           assert(fabs(finite_difference_force_een_n[nw][l][k][a][kk][i] - forces_een_n[nw][l][k][a][kk][i]) < 1.e-6);
         }
       }
   }
 }
}
}


printf("OK\n");

Force of een jastrow gradient

Get

qmckl_exit_code
qmckl_get_forces_jastrow_een_g(qmckl_context context,
                                double* const forces_jastrow_een_g,
                                const int64_t size_max);

interface
   integer(qmckl_exit_code) function qmckl_get_forces_jastrow_een_g (context, &
        forces_jastrow_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)               :: forces_jastrow_een_g(size_max)
   end function qmckl_get_forces_jastrow_een_g
end interface

Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nuclei
`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
`en_distance`	`double[elec_num][nucl_num]`	in	En distance
`tmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]`	in	Temporary intermediate tensor
`dtmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]`	in	vector of non-zero coefficients
`forces_tmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]`	in	Temporary intermediate tensor
`forces_dtmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][4][nucl_num][4][elec_num]`	in	vector of non-zero coefficients
`forces_een_n`	`double[walk_num][0:cord_num][3][nucl_num][4][elec_num]`	in	Derivative of Electron-nucleus rescaled factor
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	in	Electron-nucleus rescaled factor
`een_rescaled_e`	`double[walk_num][0:cord_num][elec_num][elec_num]`	in	Electron-nucleus rescaled factor
`een_rescaled_n_gl`	`double[walk_num][0:cord_num][nucl_num][4][elec_num]`	in	Derivative of Electron-nucleus rescaled factor
`een_rescaled_e_gl`	`double[walk_num][0:cord_num][elec_num][4][elec_num]`	in	Derivative of Electron-nucleus rescaled factor
`forces_jastrow_een_g`	`double[walk_num][3][nucl_num][3][elec_num]`	out	Derivative of Electron-nucleus jastrow

integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_een_g( &
context, walk_num, elec_num, nucl_num, &
cord_num, dim_c_vector, c_vector_full, lkpm_combined_index, &
en_distance, tmp_c, dtmp_c, forces_tmp_c, forces_dtmp_c, forces_een_n, een_rescaled_n, &
een_rescaled_e, een_rescaled_n_gl, een_rescaled_e_gl, forces_jastrow_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  :: 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)  :: en_distance(nucl_num, elec_num)
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)  :: forces_tmp_c(elec_num, 4,nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
real    (c_double )   , intent(in)  :: forces_dtmp_c(elec_num, 4, nucl_num,4,0:cord_num, 0:cord_num-1,  walk_num)
real    (c_double )   , intent(in)  :: forces_een_n(elec_num, 4, nucl_num, 3, 0:cord_num, 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_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(out) :: forces_jastrow_een_g(elec_num,3,nucl_num,3,walk_num)

integer*8 :: i, a, j, l, k, m, n, nw, ii, jj
double precision :: accu, accu2, cn


info = QMCKL_SUCCESS

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



if (cord_num == 0) return
forces_jastrow_een_g = 0.0d0
do nw =1, walk_num
do n = 1, dim_c_vector
   l = lkpm_combined_index(n, 1)
   k = lkpm_combined_index(n, 2)
   m = lkpm_combined_index(n, 4)

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

      do ii = 1, 3
         do i = 1, elec_num
           do jj = 1, 3
           forces_jastrow_een_g(i, ii, a, jj, nw) = forces_jastrow_een_g(i, ii, a, jj, nw) + (&
               tmp_c(i,a,              m,  k,nw) *  forces_een_n(i,ii,a,jj  ,m+l,nw)  &
             + forces_tmp_c(i,jj,a,    m,  k,nw) *  een_rescaled_n_gl(i,ii,a,m+l,nw)  &
             + tmp_c(i,a,              m+l,k,nw) *  forces_een_n(i,ii,a,jj,  m,  nw)  &
             + forces_tmp_c(i,jj,a,    m+l,k,nw) *  een_rescaled_n_gl(i,ii,a,m,  nw)  &
             - dtmp_c(i,ii,a,          m,  k,nw) *  een_rescaled_n_gl(i,jj,a,m+l,nw)  &
             + forces_dtmp_c(i,ii,a,jj,m,  k,nw) *  een_rescaled_n(i,a,      m+l,nw)  &
             - dtmp_c(i,ii,a,          m+l,k,nw) *  een_rescaled_n_gl(i,jj,a,m,  nw)  &
             + forces_dtmp_c(i,ii,a,jj,m+l,k,nw) *  een_rescaled_n(i,a,      m,  nw)  &
            ) * cn

           end do
         end do
      end do

   end do
end do
end do

end function qmckl_compute_forces_jastrow_een_g

Test

printf("Forces Jastrow een G\n");

/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));

rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);

double forces_jastrow_een_g[walk_num][3][nucl_num][3][elec_num];
rc = qmckl_get_forces_jastrow_een_g(context, &forces_jastrow_een_g[0][0][0][0][0], 3*3*nucl_num*walk_num*elec_num);
assert(rc == QMCKL_SUCCESS);

double finite_difference_force_een_g[walk_num][nucl_num][3][4][elec_num];
rc = qmckl_finite_difference_deriv_n(context, delta_x, &qmckl_get_jastrow_champ_factor_een_gl, &finite_difference_force_een_g[0][0][0][0][0], 4*elec_num);

for (int nw = 0; nw < walk_num; nw++){
for (int a = 0; a < nucl_num; a++) {
 for (int k = 0; k < 3; k++){
   for (int i = 0; i < elec_num; i++){
     for (int l = 0; l < 3; l++){
       //printf("finite_difference_force_een_g[%i][%i][%i][%i][%i] %+3.10f \n", nw,a,k,l,i,finite_difference_force_een_g[nw][a][k][l][i]);
       //printf("forces_jastrow_een_g         [%i][%i][%i][%i][%i] %+3.10f\n", nw,k,a,l,i,forces_jastrow_een_g[nw][k][a][l][i]);
       //assert(fabs(finite_difference_force_een_g[nw][a][k][l][i] - forces_jastrow_een_g[nw][k][a][l][i]) < 1.e-8);
     }
   }
 }
}
}
for (int nw = 0; nw < walk_num; nw++){
for (int a = 0; a < nucl_num; a++) {
 for (int k = 0; k < 3; k++){
   for (int i = 0; i < elec_num; i++){
     for (int l = 1; l < 3; l++){
       assert(fabs(finite_difference_force_een_g[nw][a][k][l][i] - forces_jastrow_een_g[nw][k][a][l][i]) < 1.e-8);
     }
   }
 }
}
}
printf("OK\n");

Force of een jastrow laplacian

Get

qmckl_exit_code
qmckl_get_forces_jastrow_een_l(qmckl_context context,
                                double* const forces_jastrow_een_l,
                                const int64_t size_max);

interface
   integer(qmckl_exit_code) function qmckl_get_forces_jastrow_een_l (context, &
        forces_jastrow_een_l, 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)               :: forces_jastrow_een_l(size_max)
   end function qmckl_get_forces_jastrow_een_l
end interface

Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`walk_num`	`int64_t`	in	Number of walkers
`elec_num`	`int64_t`	in	Number of electrons
`nucl_num`	`int64_t`	in	Number of nuclei
`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
`en_distance`	`double[elec_num][nucl_num]`	in	En distance
`tmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]`	in	Temporary intermediate tensor
`dtmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]`	in	vector of non-zero coefficients
`forces_tmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]`	in	Temporary intermediate tensor
`forces_dtmp_c`	`double[walk_num][0:cord_num-1][0:cord_num][4][nucl_num][4][elec_num]`	in	vector of non-zero coefficients
`forces_een_n`	`double[walk_num][0:cord_num][3][nucl_num][4][elec_num]`	in	Derivative of Electron-nucleus rescaled factor
`een_rescaled_n`	`double[walk_num][0:cord_num][nucl_num][elec_num]`	in	Electron-nucleus rescaled factor
`een_rescaled_e`	`double[walk_num][0:cord_num][elec_num][elec_num]`	in	Electron-nucleus rescaled factor
`een_rescaled_n_gl`	`double[walk_num][0:cord_num][nucl_num][4][elec_num]`	in	Derivative of Electron-nucleus rescaled factor
`een_rescaled_e_gl`	`double[walk_num][0:cord_num][elec_num][4][elec_num]`	in	Derivative of Electron-nucleus rescaled factor
`forces_jastrow_een_l`	`double[walk_num][3][nucl_num]`	out	Derivative of Electron-nucleus jastrow

integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_een_l( &
context, walk_num, elec_num, nucl_num, &
cord_num, dim_c_vector, c_vector_full, lkpm_combined_index, &
en_distance, tmp_c, dtmp_c, forces_tmp_c, forces_dtmp_c, forces_een_n, een_rescaled_n, &
een_rescaled_e, een_rescaled_n_gl, een_rescaled_e_gl, forces_jastrow_een_l)&
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  :: 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)  :: en_distance(nucl_num, elec_num)
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)  :: forces_tmp_c(elec_num, 4,nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
real    (c_double )   , intent(in)  :: forces_dtmp_c(elec_num, 4, nucl_num,4,0:cord_num, 0:cord_num-1,  walk_num)
real    (c_double )   , intent(in)  :: forces_een_n(elec_num, 4, nucl_num, 3, 0:cord_num, 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_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(out) :: forces_jastrow_een_l(nucl_num,3,walk_num)

integer*8 :: i, a, j, l, k, m, n, nw, ii, jj
double precision :: accu, accu2, cn


info = QMCKL_SUCCESS

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



if (cord_num == 0) return
forces_jastrow_een_l = 0.0d0
do nw =1, walk_num
do n = 1, dim_c_vector
   l = lkpm_combined_index(n, 1)
   k = lkpm_combined_index(n, 2)
   m = lkpm_combined_index(n, 4)

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

      do ii = 1, 3
         do i = 1, elec_num
           forces_jastrow_een_l(a, ii, nw) = forces_jastrow_een_l(a, ii, nw) + (&
                tmp_c(i,a,             m,  k,nw) *  forces_een_n(i,4,a,ii  ,m+l,nw) &
              + forces_tmp_c(i,ii,a,   m,  k,nw) *  een_rescaled_n_gl(i,4,a,m+l,nw) &
              + tmp_c(i,a,             m+l,k,nw) *  forces_een_n(i,4,a,ii,  m,  nw) &
              + forces_tmp_c(i,ii,a,   m+l,k,nw) *  een_rescaled_n_gl(i,4,a,m,  nw) &
              - dtmp_c(i,4,a,          m,  k,nw) *  een_rescaled_n_gl(i,ii,a,m+l,nw) &
              + forces_dtmp_c(i,4,a,ii,m,  k,nw) *  een_rescaled_n(i,a,      m+l,nw) &
              - dtmp_c(i,4,a,          m+l,k,nw) *  een_rescaled_n_gl(i,ii,a,m,  nw) &
              + forces_dtmp_c(i,4,a,ii,m+l,k,nw) *  een_rescaled_n(i,a,      m,  nw) &
            )   * cn
         end do
      end do

      cn = cn + cn
       do ii = 1, 3
         do i = 1, elec_num
           do jj = 1, 3
             forces_jastrow_een_l(a, ii, nw) = forces_jastrow_een_l(a, ii, nw) + (&
                 dtmp_c(i,jj,a,m, k,nw) * forces_een_n(i,jj,a,ii ,m+l,nw) + &
                 dtmp_c(i,jj,a,m+l, k,nw) * forces_een_n(i,jj,a,ii ,m,nw) + &
                 forces_dtmp_c(i,jj,a,ii,m, k,nw) * een_rescaled_n_gl(i,jj,a,m+l,nw) + &
                 forces_dtmp_c(i,jj,a,ii,m+l, k,nw) * een_rescaled_n_gl(i,jj,a,m,nw) &
               ) * cn
           end do
         end do
       end do
   end do
end do
end do

end function qmckl_compute_forces_jastrow_een_l

Test

printf("Forces Jastrow een l\n");

/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));

rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);

double forces_jastrow_een_l[walk_num][3][nucl_num];
rc = qmckl_get_forces_jastrow_een_l(context, &forces_jastrow_een_l[0][0][0], 3*nucl_num*walk_num);
assert(rc == QMCKL_SUCCESS);

double finite_difference_force_een_l[walk_num][nucl_num][3][4][elec_num];
rc = qmckl_finite_difference_deriv_n(context, delta_x, &qmckl_get_jastrow_champ_factor_een_gl, &finite_difference_force_een_l[0][0][0][0][0], 4*elec_num);

double finite_difference_force_een_l_sum[walk_num][nucl_num][3];
for (int nw = 0; nw < walk_num; nw++){
for (int a = 0; a < nucl_num; a++) {
 for (int k = 0; k < 3; k++){
   finite_difference_force_een_l_sum[nw][a][k] = 0;
   for (int i = 0; i < elec_num; i++){
     finite_difference_force_een_l_sum[nw][a][k] = finite_difference_force_een_l_sum[nw][a][k] + finite_difference_force_een_l[nw][a][k][3][i];
   }
 }
}
}

for (int nw = 0; nw < walk_num; nw++){
for (int a = 0; a < nucl_num; a++) {
 for (int k = 0; k < 3; k++){
   printf("finite_difference_force_een_l_sum[%i][%i][%i] %+3.10f \n", nw,a,k,finite_difference_force_een_l_sum[nw][a][k]);
   printf("forces_jastrow_een_l         [%i][%i][%i] %+3.10f\n", nw,k,a,forces_jastrow_een_l[nw][k][a]);
   assert(fabs(finite_difference_force_een_l_sum[nw][a][k] - forces_jastrow_een_l[nw][k][a]) < 1.e-8);
 }
}
}
for (int nw = 0; nw < walk_num; nw++){
for (int a = 0; a < nucl_num; a++) {
 for (int k = 0; k < 3; k++){
   assert(fabs(finite_difference_force_een_l_sum[nw][a][k] - forces_jastrow_een_l[nw][k][a]) < 1.e-8);
 }
}
}
printf("OK\n");

143 KiB Raw Blame History

Forces

Introduction

Context

Data structure

Test

Finite-difference function

Force of en jastrow value

Get

Compute

Test

Force of en jastrow gradient

Get

Compute

Test

Force of en jastrow laplacien

Get

Compute

Test

Force of tmp_c matrix

Get

Compute

Test

Force of dtmp_c matrix

Get

Compute

Test

Force of een jastrow value

Get

Compute

Test

Force of een_rescaled_n_gl

Get

Compute

Test

Force of een jastrow gradient

Get

Compute

Test

Force of een jastrow laplacian

Get

Compute

Test

143 KiB

Raw Blame History