TREX/qmckl

Fork 0

mirror of https://github.com/TREX-CoE/qmckl.git synced 2025-04-30 20:35:10 +02:00

Emiel Slootman 6a33fbe901 2-body force done

2024-12-11 14:17:33 +01:00

58 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

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;
} 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]);
int64_t size_max;

/* 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 new_coords[3][nucl_num];

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, 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(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_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(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_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");

58 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

58 KiB

Raw Blame History