TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2025-04-30 20:35:10 +02:00
Code Issues Projects Releases Wiki Activity
qmckl/org/qmckl_forces.org
Emiel Slootman 7b3a8fd558 small things
2025-02-03 16:00:01 +01:00

6551 lines
247 KiB
Org Mode
Raw Blame History

#+TITLE: Forces
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
* Introduction
* Headers :noexport:
#+begin_src elisp :noexport :results none
(org-babel-lob-ingest "../tools/lib.org")
#+end_src
#+begin_src c :tangle (eval h_private_func)
#ifndef QMCKL_FORCES_HPF
#define QMCKL_FORCES_HPF
#+end_src
#+begin_src c :tangle (eval h_private_type)
#ifndef QMCKL_FORCES_HPT
#define QMCKL_FORCES_HPT
#include <stdbool.h>
#+end_src
#+begin_src c :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include <assert.h>
#include <math.h>
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <stdbool.h>
#include <stdio.h>
#include "n2.h"
#include "chbrclf.h"
#include "qmckl_jastrow_champ_private_func.h"
#include "qmckl_jastrow_champ_single_private_func.h"
#include "qmckl_forces_private_func.h"
int main() {
qmckl_context context;
context = qmckl_context_create();
#+end_src
#+begin_src c :tangle (eval c)
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_STDINT_H
#include <stdint.h>
#elif HAVE_INTTYPES_H
#include <inttypes.h>
#endif
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include "qmckl.h"
#include "qmckl_context_private_type.h"
#include "qmckl_memory_private_type.h"
#include "qmckl_memory_private_func.h"
#include "qmckl_jastrow_champ_private_type.h"
#include "qmckl_jastrow_champ_private_func.h"
#include "qmckl_jastrow_champ_single_private_type.h"
#include "qmckl_jastrow_champ_single_private_func.h"
#include "qmckl_forces_private_type.h"
#include "qmckl_forces_private_func.h"
#+end_src
* Context
** Data structure
#+begin_src c :comments org :tangle (eval h_private_type)
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;
double * restrict forces_ao_value;
uint64_t forces_ao_value_date;
double * restrict forces_mo_value;
uint64_t forces_mo_value_date;
double * restrict forces_mo_g;
uint64_t forces_mo_g_date;
double * restrict forces_mo_l;
uint64_t forces_mo_l_date;
double * forces_jastrow_single_en;
uint64_t forces_jastrow_single_en_date;
double * forces_jastrow_single_een;
uint64_t forces_jastrow_single_een_date;
double * forces_delta_p;
uint64_t forces_delta_p_date;
} qmckl_forces_struct;
#+end_src
** Test
#+begin_src c :tangle (eval c_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 };
#+end_src
* Finite-difference function
We introduce here a general function to compute the derivatives of any quantity with respect to nuclear coordinates.
using finite-differences.
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
typedef qmckl_exit_code (*function_callback)(qmckl_context context, double* const output, const int64_t size);
qmckl_exit_code qmckl_finite_difference_deriv_n(
qmckl_context context,
const double delta_x, // Step size for finite difference
function_callback get_function, // Function to compute values
double* const derivative_output, // Output derivative array: nucl_num*3*size
int64_t const size); // Size of the object to differentiate
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
typedef qmckl_exit_code (*function_callback)(qmckl_context context, double* const output, const int64_t size);
qmckl_exit_code qmckl_finite_difference_deriv_n(
qmckl_context context,
const double delta_x,
function_callback get_function,
double* const derivative_output,
int64_t const size)
{
// Finite difference coefficients for a 9-point stencil
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 };
qmckl_exit_code rc;
int64_t walk_num;
rc = qmckl_get_electron_walk_num(context, &walk_num);
if (rc != QMCKL_SUCCESS) {
return rc;
}
int64_t nucl_num;
rc = qmckl_get_nucleus_num(context, &nucl_num);
if (rc != QMCKL_SUCCESS) {
return rc;
}
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* function_values = (double*) malloc(walk_num*size * sizeof(double));
if (function_values == NULL) {
free(nucleus_coord);
free(temp_coord);
return QMCKL_ALLOCATION_FAILED;
}
memset(derivative_output, 0, nucl_num*3*walk_num*size*sizeof(double));
// 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);
rc = qmckl_single_touch(context);
assert(rc == QMCKL_SUCCESS);
// Call the provided function
rc = get_function(context, function_values, size);
assert(rc == QMCKL_SUCCESS);
// Accumulate derivative using finite-difference coefficients
for (int64_t nw=0 ; nw<walk_num ; nw++) {
int64_t shift = nucl_num*3*size*nw + size*(k + 3*a);
for (int64_t i = 0; i < size; i++) {
derivative_output[i+shift] += coef[m + 4] * function_values[nw*size+i];
}
}
}
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);
// Normalize by the step size
for (int64_t i = 0; i < size*3*nucl_num*walk_num ; i++) {
derivative_output[i] /= delta_x;
}
free(nucleus_coord);
free(temp_coord);
free(function_values);
return QMCKL_SUCCESS;
}
#+end_src
* Force of en jastrow value
** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_forces_jastrow_en(qmckl_context context,
double* const forces_jastrow_en,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_forces_jastrow_en(qmckl_context context,
double* const forces_jastrow_en,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_forces_jastrow_en(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = 3 * ctx->nucleus.num * ctx->electron.walker.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_forces_jastrow_en",
"Array too small. Expected 3*nucl_num*walk_num");
}
memcpy(forces_jastrow_en, ctx->forces.forces_jastrow_en, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org
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
#+end_src
** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_jastrow_en(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_jastrow_en(qmckl_context context)
{
qmckl_exit_code rc;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_provide_forces_jastrow_en",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->jastrow_champ.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_provide_forces_jastrow_en",
NULL);
}
/* Check if en rescaled distance is provided */
rc = qmckl_provide_en_distance_rescaled(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_en_distance_rescaled_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Compute if necessary */
if (ctx->date > ctx->forces.forces_jastrow_en_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
if (ctx->forces.forces_jastrow_en != NULL) {
rc = qmckl_free(context, ctx->forces.forces_jastrow_en);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_provide_forces_jastrow_en",
"Unable to free ctx->forces.forces_jastrow_en");
}
ctx->forces.forces_jastrow_en = NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_jastrow_en == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walker.num * 3 * ctx->nucleus.num * sizeof(double);
double* forces_jastrow_en = (double*) qmckl_malloc(context, mem_info);
if (forces_jastrow_en == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_provide_forces_jastrow_en",
NULL);
}
ctx->forces.forces_jastrow_en = forces_jastrow_en;
}
rc = qmckl_compute_forces_jastrow_en(context,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow_champ.type_nucl_num,
ctx->jastrow_champ.type_nucl_vector,
ctx->jastrow_champ.aord_num,
ctx->jastrow_champ.a_vector,
ctx->jastrow_champ.en_distance_rescaled,
ctx->jastrow_champ.en_distance_rescaled_gl,
ctx->forces.forces_jastrow_en);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->forces.forces_jastrow_en_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
** Compute
:PROPERTIES:
:Name: qmckl_compute_forces_jastrow_en
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_forces_jastrow_en_args
| 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 |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
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
#+end_src
# #+CALL: generate_c_header(table=qmckl_factor_en_gl_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_forces_jastrow_en_doc (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
const int64_t* type_nucl_vector,
const int64_t aord_num,
const double* a_vector,
const double* en_distance_rescaled,
const double* en_distance_rescaled_gl,
double* const forces_jastrow_en );
qmckl_exit_code qmckl_compute_forces_jastrow_en (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
const int64_t* type_nucl_vector,
const int64_t aord_num,
const double* a_vector,
const double* en_distance_rescaled,
const double* en_distance_rescaled_gl,
double* const forces_jastrow_en );
qmckl_exit_code qmckl_compute_forces_jastrow_en (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
const int64_t* type_nucl_vector,
const int64_t aord_num,
const double* a_vector,
const double* en_distance_rescaled,
const double* en_distance_rescaled_gl,
double* const forces_jastrow_en );
#+end_src
#+begin_src c :tangle (eval c) :comments org :exports none
qmckl_exit_code
qmckl_compute_forces_jastrow_en (const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
const int64_t* type_nucl_vector,
const int64_t aord_num,
const double* a_vector,
const double* en_distance_rescaled,
const double* en_distance_rescaled_gl,
double* const forces_jastrow_en )
{
#ifdef HAVE_HPC
return qmckl_compute_forces_jastrow_en_doc
#else
return qmckl_compute_forces_jastrow_en_doc
#endif
(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);
}
#+end_src
** Test
#+begin_src c :tangle (eval c_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");
#+end_src
* Force of en jastrow gradient
** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_forces_jastrow_en_g(qmckl_context context,
double* const forces_jastrow_en_g,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_forces_jastrow_en_g(qmckl_context context,
double* const forces_jastrow_en_g,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_forces_jastrow_en_g(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = 3*3 * ctx->nucleus.num * ctx->electron.walker.num * ctx->electron.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_forces_jastrow_en_g",
"Array too small. Expected 3*3*nucl_num*walk_num_elec_num");
}
memcpy(forces_jastrow_en_g, ctx->forces.forces_jastrow_en_g, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org
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
#+end_src
** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_jastrow_en_g(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_jastrow_en_g(qmckl_context context)
{
qmckl_exit_code rc;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_provide_forces_jastrow_en_g",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->jastrow_champ.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_provide_forces_jastrow_en_g",
NULL);
}
/* Check if en rescaled distance is provided */
rc = qmckl_provide_en_distance_rescaled(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_en_distance_rescaled_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_en_distance(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Compute if necessary */
if (ctx->date > ctx->forces.forces_jastrow_en_g_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
if (ctx->forces.forces_jastrow_en_g != NULL) {
rc = qmckl_free(context, ctx->forces.forces_jastrow_en_g);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_provide_forces_jastrow_en",
"Unable to free ctx->forces.forces_jastrow_en_g");
}
ctx->forces.forces_jastrow_en_g = NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_jastrow_en_g == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walker.num * 3 * 3 * ctx->nucleus.num * ctx->electron.num * sizeof(double);
double* forces_jastrow_en_g = (double*) qmckl_malloc(context, mem_info);
if (forces_jastrow_en_g == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_provide_forces_jastrow_en_g",
NULL);
}
ctx->forces.forces_jastrow_en_g = forces_jastrow_en_g;
}
rc = qmckl_compute_forces_jastrow_en_g(context,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow_champ.type_nucl_num,
ctx->jastrow_champ.type_nucl_vector,
ctx->jastrow_champ.aord_num,
ctx->jastrow_champ.a_vector,
ctx->jastrow_champ.rescale_factor_en,
ctx->electron.en_distance,
ctx->jastrow_champ.en_distance_rescaled,
ctx->jastrow_champ.en_distance_rescaled_gl,
ctx->forces.forces_jastrow_en_g);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->forces.forces_jastrow_en_g_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
** Compute
:PROPERTIES:
:Name: qmckl_compute_forces_jastrow_en_g
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_forces_jastrow_en_g_args
| 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 |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
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
#+end_src
# #+CALL: generate_c_header(table=qmckl_factor_en_gl_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_forces_jastrow_en_g_doc (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
const int64_t* type_nucl_vector,
const int64_t aord_num,
const double* a_vector,
const double* rescale_factor_en,
const double* en_distance,
const double* en_distance_rescaled,
const double* en_distance_rescaled_gl,
double* const forces_jastrow_en_g );
qmckl_exit_code qmckl_compute_forces_jastrow_en_g (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
const int64_t* type_nucl_vector,
const int64_t aord_num,
const double* a_vector,
const double* rescale_factor_en,
const double* en_distance,
const double* en_distance_rescaled,
const double* en_distance_rescaled_gl,
double* const forces_jastrow_en_g );
qmckl_exit_code qmckl_compute_forces_jastrow_en_g (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
const int64_t* type_nucl_vector,
const int64_t aord_num,
const double* a_vector,
const double* rescale_factor_en,
const double* en_distance,
const double* en_distance_rescaled,
const double* en_distance_rescaled_gl,
double* const forces_jastrow_en_g );
#+end_src
#+begin_src c :tangle (eval c) :comments org :exports none
qmckl_exit_code
qmckl_compute_forces_jastrow_en_g (const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
const int64_t* type_nucl_vector,
const int64_t aord_num,
const double* a_vector,
const double* rescale_factor_en,
const double* en_distance,
const double* en_distance_rescaled,
const double* en_distance_rescaled_gl,
double* const forces_jastrow_en_g )
{
#ifdef HAVE_HPC
return qmckl_compute_forces_jastrow_en_g_doc
#else
return qmckl_compute_forces_jastrow_en_g_doc
#endif
(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);
}
#+end_src
** Test
#+begin_src c :tangle (eval c_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");
#+end_src
* Force of en jastrow laplacien
** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_forces_jastrow_en_l(qmckl_context context,
double* const forces_jastrow_en_l,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_forces_jastrow_en_l(qmckl_context context,
double* const forces_jastrow_en_l,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_forces_jastrow_en_l(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = 3 * ctx->nucleus.num * ctx->electron.walker.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_forces_jastrow_en_l",
"Array too small. Expected 3*nucl_num*walk_num");
}
memcpy(forces_jastrow_en_l, ctx->forces.forces_jastrow_en_l, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org
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
#+end_src
** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_jastrow_en_l(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_jastrow_en_l(qmckl_context context)
{
qmckl_exit_code rc;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_provide_forces_jastrow_en_l",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->jastrow_champ.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_provide_forces_jastrow_en_l",
NULL);
}
/* Check if en rescaled distance is provided */
rc = qmckl_provide_en_distance_rescaled(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_en_distance_rescaled_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_en_distance(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Compute if necessary */
if (ctx->date > ctx->forces.forces_jastrow_en_l_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
if (ctx->forces.forces_jastrow_en_l != NULL) {
rc = qmckl_free(context, ctx->forces.forces_jastrow_en_l);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_provide_forces_jastrow_en_l",
"Unable to free ctx->forces.forces_jastrow_en_l");
}
ctx->forces.forces_jastrow_en_l = NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_jastrow_en_l == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walker.num * 3 * ctx->nucleus.num * sizeof(double);
double* forces_jastrow_en_l = (double*) qmckl_malloc(context, mem_info);
if (forces_jastrow_en_l == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_provide_forces_jastrow_en_l",
NULL);
}
ctx->forces.forces_jastrow_en_l = forces_jastrow_en_l;
}
rc = qmckl_compute_forces_jastrow_en_l(context,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow_champ.type_nucl_num,
ctx->jastrow_champ.type_nucl_vector,
ctx->jastrow_champ.aord_num,
ctx->jastrow_champ.a_vector,
ctx->jastrow_champ.rescale_factor_en,
ctx->electron.en_distance,
ctx->jastrow_champ.en_distance_rescaled,
ctx->jastrow_champ.en_distance_rescaled_gl,
ctx->forces.forces_jastrow_en_l);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->forces.forces_jastrow_en_l_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
** Compute
:PROPERTIES:
:Name: qmckl_compute_forces_jastrow_en_l
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_forces_jastrow_en_l_args
| 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 |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
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
#+end_src
# #+CALL: generate_c_header(table=qmckl_factor_en_gl_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_forces_jastrow_en_l_doc (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
const int64_t* type_nucl_vector,
const int64_t aord_num,
const double* a_vector,
const double* rescale_factor_en,
const double* en_distance,
const double* en_distance_rescaled,
const double* en_distance_rescaled_gl,
double* const forces_jastrow_en_l );
qmckl_exit_code qmckl_compute_forces_jastrow_en_l (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
const int64_t* type_nucl_vector,
const int64_t aord_num,
const double* a_vector,
const double* rescale_factor_en,
const double* en_distance,
const double* en_distance_rescaled,
const double* en_distance_rescaled_gl,
double* const forces_jastrow_en_l );
qmckl_exit_code qmckl_compute_forces_jastrow_en_l (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
const int64_t* type_nucl_vector,
const int64_t aord_num,
const double* a_vector,
const double* rescale_factor_en,
const double* en_distance,
const double* en_distance_rescaled,
const double* en_distance_rescaled_gl,
double* const forces_jastrow_en_l );
#+end_src
#+begin_src c :tangle (eval c) :comments org :exports none
qmckl_exit_code
qmckl_compute_forces_jastrow_en_l (const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
const int64_t* type_nucl_vector,
const int64_t aord_num,
const double* a_vector,
const double* rescale_factor_en,
const double* en_distance,
const double* en_distance_rescaled,
const double* en_distance_rescaled_gl,
double* const forces_jastrow_en_l )
{
#ifdef HAVE_HPC
return qmckl_compute_forces_jastrow_en_l_doc
#else
return qmckl_compute_forces_jastrow_en_l_doc
#endif
(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);
}
#+end_src
** Test
#+begin_src c :tangle (eval c_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");
#+end_src
* Force of tmp_c matrix
** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_forces_tmp_c(qmckl_context context,
double* const forces_tmp_c,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_forces_tmp_c(qmckl_context context,
double* const forces_tmp_c,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_forces_tmp_c(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = 4 * ctx->electron.walker.num * ctx->electron.num * ctx->nucleus.num *
ctx->jastrow_champ.cord_num * (ctx->jastrow_champ.cord_num+1);
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_forces_tmp_c",
"Array too small. Expected 4*walk_num*elec_num*nucl_num*cord_num*(cord_num+1)");
}
memcpy(forces_tmp_c, ctx->forces.forces_tmp_c, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_tmp_c(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_tmp_c(qmckl_context context)
{
qmckl_exit_code rc;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (ctx->jastrow_champ.cord_num > 0) {
/* Check if en rescaled distance is provided */
rc = qmckl_provide_een_rescaled_e(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance gl derivatives is provided */
rc = qmckl_provide_een_rescaled_n_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
}
/* Compute if necessary */
if (ctx->date > ctx->forces.forces_tmp_c_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
if (ctx->forces.forces_tmp_c != NULL) {
rc = qmckl_free(context, ctx->forces.forces_tmp_c);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_provide_forces_tmp_c",
"Unable to free ctx->forces.forces_tmp_c");
}
ctx->forces.forces_tmp_c = NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_tmp_c == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 4 * ctx->electron.num * ctx->jastrow_champ.cord_num *
(ctx->jastrow_champ.cord_num+1) * ctx->nucleus.num * ctx->electron.walker.num * sizeof(double);
double* forces_tmp_c = (double*) qmckl_malloc(context, mem_info);
if (forces_tmp_c == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_provide_forces_tmp_c",
NULL);
}
ctx->forces.forces_tmp_c = forces_tmp_c;
}
rc = qmckl_compute_forces_tmp_c(context,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow_champ.cord_num,
ctx->jastrow_champ.een_rescaled_e,
ctx->jastrow_champ.een_rescaled_n_gl,
ctx->forces.forces_tmp_c);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->forces.forces_tmp_c_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
** Compute
:PROPERTIES:
:Name: qmckl_compute_forces_tmp_c
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_forces_tmp_c_args
| 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 |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
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
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_forces_tmp_c (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t cord_num,
const double* een_rescaled_e,
const double* een_rescaled_n_gl,
double* const forces_tmp_c );
#+end_src
** Test
#+begin_src c :tangle (eval 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],
4*nucl_num*walk_num*elec_num*(cord_num+1)*cord_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 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");
#+end_src
* Force of dtmp_c matrix
** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_forces_dtmp_c(qmckl_context context,
double* const forces_dtmp_c,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_forces_dtmp_c(qmckl_context context,
double* const forces_dtmp_c,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_forces_dtmp_c(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = 4 * 4 * ctx->electron.walker.num * ctx->electron.num * ctx->nucleus.num *
ctx->jastrow_champ.cord_num * (ctx->jastrow_champ.cord_num+1);
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_forces_dtmp_c",
"Array too small. Expected 4*4*walk_num*elec_num*nucl_num*cord_num*(cord_num+1)");
}
memcpy(forces_dtmp_c, ctx->forces.forces_dtmp_c, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_dtmp_c(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_dtmp_c(qmckl_context context)
{
qmckl_exit_code rc;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (ctx->jastrow_champ.cord_num > 0) {
/* Check if en rescaled distance is provided */
rc = qmckl_provide_een_rescaled_e_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance gl derivatives is provided */
rc = qmckl_provide_een_rescaled_n_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
}
/* Compute if necessary */
if (ctx->date > ctx->forces.forces_dtmp_c_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
if (ctx->forces.forces_dtmp_c != NULL) {
rc = qmckl_free(context, ctx->forces.forces_dtmp_c);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_provide_forces_dtmp_c",
"Unable to free ctx->forces.forces_dtmp_c");
}
ctx->forces.forces_dtmp_c = NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_dtmp_c == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 4* 4 * ctx->electron.num * ctx->jastrow_champ.cord_num *
(ctx->jastrow_champ.cord_num+1) * ctx->nucleus.num * ctx->electron.walker.num * sizeof(double);
double* forces_dtmp_c = (double*) qmckl_malloc(context, mem_info);
if (forces_dtmp_c == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_provide_forces_dtmp_c",
NULL);
}
ctx->forces.forces_dtmp_c = forces_dtmp_c;
}
rc = qmckl_compute_forces_dtmp_c(context,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow_champ.cord_num,
ctx->jastrow_champ.een_rescaled_e_gl,
ctx->jastrow_champ.een_rescaled_n_gl,
ctx->forces.forces_dtmp_c);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->forces.forces_dtmp_c_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
** Compute
:PROPERTIES:
:Name: qmckl_compute_forces_dtmp_c
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_forces_dtmp_c_args
| 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 |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
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
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_forces_dtmp_c (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t cord_num,
const double* een_rescaled_e_gl,
const double* een_rescaled_n_gl,
double* const forces_dtmp_c );
#+end_src
** Test
#+begin_src c :tangle (eval 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],
4*4*nucl_num*walk_num*elec_num*(cord_num+1)*cord_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 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");
#+end_src
* Force of een jastrow value
** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_forces_jastrow_een(qmckl_context context,
double* const forces_jastrow_een,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_forces_jastrow_een(qmckl_context context,
double* const forces_jastrow_een,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_forces_jastrow_een(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = 3 * ctx->nucleus.num * ctx->electron.walker.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_forces_jastrow_een",
"Array too small. Expected 3*nucl_num*walk_num");
}
memcpy(forces_jastrow_een, ctx->forces.forces_jastrow_een, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org
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
#+end_src
** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_jastrow_een(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_jastrow_een(qmckl_context context)
{
qmckl_exit_code rc;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (ctx->jastrow_champ.cord_num > 0) {
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_een_rescaled_n(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance gl derivatives is provided */
rc = qmckl_provide_een_rescaled_n_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if tmp_c is provided */
rc = qmckl_provide_tmp_c(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if forces_tmp_c is provided */
rc = qmckl_provide_forces_tmp_c(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_jastrow_champ_c_vector_full(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_lkpm_combined_index(context);
if(rc != QMCKL_SUCCESS) return rc;
}
/* Compute if necessary */
if (ctx->date > ctx->forces.forces_jastrow_een_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
if (ctx->forces.forces_jastrow_een != NULL) {
rc = qmckl_free(context, ctx->forces.forces_jastrow_een);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_provide_forces_jastrow_een",
"Unable to free ctx->forces.forces_jastrow_een");
}
ctx->forces.forces_jastrow_een = NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_jastrow_een == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 3 * ctx->nucleus.num * ctx->electron.walker.num * sizeof(double);
double* forces_jastrow_een = (double*) qmckl_malloc(context, mem_info);
if (forces_jastrow_een == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_provide_forces_jastrow_een",
NULL);
}
ctx->forces.forces_jastrow_een = forces_jastrow_een;
}
rc = qmckl_compute_forces_jastrow_een(context,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow_champ.cord_num,
ctx->jastrow_champ.dim_c_vector,
ctx->jastrow_champ.c_vector_full,
ctx->jastrow_champ.lkpm_combined_index,
ctx->jastrow_champ.een_rescaled_n,
ctx->jastrow_champ.een_rescaled_n_gl,
ctx->jastrow_champ.tmp_c,
ctx->forces.forces_tmp_c,
ctx->forces.forces_jastrow_een);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->forces.forces_jastrow_een_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
** Compute
:PROPERTIES:
:Name: qmckl_compute_forces_jastrow_een
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_forces_jastrow_een_args
| 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_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 |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
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_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_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
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_forces_jastrow_een (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t cord_num,
const int64_t dim_c_vector,
const double* c_vector_full,
const int64_t* lkpm_combined_index,
const double* een_rescaled_n,
const double* een_rescaled_n_gl,
const double* tmp_c,
const double* forces_tmp_c,
double* const forces_jastrow_een );
#+end_src
** Test
#+begin_src c :tangle (eval c_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");
#+end_src
* Force of een_rescaled_n_gl
** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_forces_een_rescaled_n_gl(qmckl_context context,
double* const forces_een_n,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_forces_een_rescaled_n_gl(qmckl_context context,
double* const forces_een_n,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_forces_een_rescaled_n_gl(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = 3* ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1);
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_forces_een_rescaled_n_gl",
"Array too small. Expected ctx->electron.num * 3 * 4 * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1)");
}
memcpy(forces_een_n, ctx->forces.forces_een_n, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_een_rescaled_n_gl(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_een_rescaled_n_gl(qmckl_context context)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
/* Check if ee distance is provided */
qmckl_exit_code rc = qmckl_provide_en_distance(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if ee distance is provided */
rc = qmckl_provide_een_rescaled_n(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if ee gl distance is provided */
rc = qmckl_provide_een_rescaled_n_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Compute if necessary */
if (ctx->date > ctx->forces.forces_een_n_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
if (ctx->forces.forces_een_n != NULL) {
rc = qmckl_free(context, ctx->forces.forces_een_n);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_provide_forces_een_rescaled_n_gl",
"Unable to free ctx->forces.forces_een_n");
}
ctx->forces.forces_een_n = NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_een_n == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.num * 3 * 4 * ctx->nucleus.num *
ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1) * sizeof(double);
double* forces_een_n = (double*) qmckl_malloc(context, mem_info);
if (forces_een_n == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_provide_forces_een_rescaled_n_gl",
NULL);
}
ctx->forces.forces_een_n = forces_een_n;
}
rc = qmckl_compute_forces_een_rescaled_n_gl(context,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow_champ.type_nucl_num,
ctx->jastrow_champ.type_nucl_vector,
ctx->jastrow_champ.cord_num,
ctx->jastrow_champ.rescale_factor_en,
ctx->electron.en_distance,
ctx->jastrow_champ.een_rescaled_n,
ctx->jastrow_champ.een_rescaled_n_gl,
ctx->forces.forces_een_n);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->forces.forces_een_n_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
** Compute
:PROPERTIES:
:Name: qmckl_compute_forces_een_rescaled_n_gl
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_compute_forces_een_rescaled_n_gl_args
| 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 |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
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
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_forces_een_rescaled_n_gl (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
int64_t* const type_nucl_vector,
const int64_t cord_num,
const double* rescale_factor_en,
const double* en_distance,
const double* een_rescaled_n,
const double* een_rescaled_n_gl,
double* const forces_een_n );
#+end_src
** Test
#+begin_src c :tangle (eval c_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");
#+end_src
* Force of een jastrow gradient
** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_forces_jastrow_een_g(qmckl_context context,
double* const forces_jastrow_een_g,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_forces_jastrow_een_g(qmckl_context context,
double* const forces_jastrow_een_g,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_forces_jastrow_een_g(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = ctx->electron.walker.num * 3 * ctx->electron.num * 3 * ctx->nucleus.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_forces_jastrow_een_g",
"Array too small. Expected 3*3*walk_num*elec_num*nucl_num");
}
memcpy(forces_jastrow_een_g, ctx->forces.forces_jastrow_een_g, sze*sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org
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
#+end_src
#
** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_jastrow_een_g(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_jastrow_een_g(qmckl_context context)
{
qmckl_exit_code rc;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (ctx->jastrow_champ.cord_num > 0) {
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_een_rescaled_n(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_een_rescaled_n_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_jastrow_champ_c_vector_full(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_lkpm_combined_index(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if tmp_c is provided */
rc = qmckl_provide_tmp_c(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if dtmp_c is provided */
rc = qmckl_provide_dtmp_c(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if forces_tmp_c is provided */
rc = qmckl_provide_forces_tmp_c(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if forces_dtmp_c is provided */
rc = qmckl_provide_forces_dtmp_c(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if ne distance is provided */
rc = qmckl_provide_en_distance(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_forces_een_rescaled_n_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
}
/* Compute if necessary */
if (ctx->date > ctx->forces.forces_jastrow_een_g_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
if (ctx->forces.forces_jastrow_een_g != NULL) {
rc = qmckl_free(context, ctx->forces.forces_jastrow_een_g);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_provide_forces_jastrow_een_g",
"Unable to free ctx->forces.forces_jastrow_een_g");
}
ctx->forces.forces_jastrow_een_g = NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_jastrow_een_g == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 3 * 3 * ctx->electron.num * ctx->electron.walker.num * ctx->nucleus.num * sizeof(double);
double* forces_jastrow_een_g = (double*) qmckl_malloc(context, mem_info);
if (forces_jastrow_een_g == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_provide_forces_jastrow_een_g",
NULL);
}
ctx->forces.forces_jastrow_een_g = forces_jastrow_een_g;
}
rc = qmckl_compute_forces_jastrow_een_g(context,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow_champ.cord_num,
ctx->jastrow_champ.dim_c_vector,
ctx->jastrow_champ.c_vector_full,
ctx->jastrow_champ.lkpm_combined_index,
ctx->electron.en_distance,
ctx->jastrow_champ.tmp_c,
ctx->jastrow_champ.dtmp_c,
ctx->forces.forces_tmp_c,
ctx->forces.forces_dtmp_c,
ctx->forces.forces_een_n,
ctx->jastrow_champ.een_rescaled_n,
ctx->jastrow_champ.een_rescaled_n_gl,
ctx->forces.forces_jastrow_een_g);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->forces.forces_jastrow_een_g_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
** Compute
:PROPERTIES:
:Name: qmckl_compute_forces_jastrow_een_g
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_forces_jastrow_een_g_args
| 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_n_gl~ | ~double[walk_num][0:cord_num][nucl_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 |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
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_n_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_n_gl(elec_num, 4, nucl_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
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_forces_jastrow_een_g(
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t cord_num,
const int64_t dim_c_vector,
const double* c_vector_full,
const int64_t* lkpm_combined_index,
const double* en_distance,
const double* tmp_c,
const double* dtmp_c,
const double* forces_tmp_c,
const double* forces_dtmp_c,
const double* forces_een_n,
const double* een_rescaled_n,
const double* een_rescaled_n_gl,
double* const forces_jastrow_een_g );
#+end_src
** Test
#+begin_src c :tangle (eval c_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");
#+end_src
* Force of een jastrow laplacian
** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_forces_jastrow_een_l(qmckl_context context,
double* const forces_jastrow_een_l,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_forces_jastrow_een_l(qmckl_context context,
double* const forces_jastrow_een_l,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_forces_jastrow_een_l(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = ctx->electron.walker.num * 3 * ctx->nucleus.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_forces_jastrow_een_l",
"Array too small. Expected 3*walk_num*nucl_num");
}
memcpy(forces_jastrow_een_l, ctx->forces.forces_jastrow_een_l, sze*sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org
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
#+end_src
#
** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_jastrow_een_l(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_jastrow_een_l(qmckl_context context)
{
qmckl_exit_code rc;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (ctx->jastrow_champ.cord_num > 0) {
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_een_rescaled_n(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_een_rescaled_n_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_jastrow_champ_c_vector_full(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_lkpm_combined_index(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if tmp_c is provided */
rc = qmckl_provide_tmp_c(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if dtmp_c is provided */
rc = qmckl_provide_dtmp_c(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if forces_tmp_c is provided */
rc = qmckl_provide_forces_tmp_c(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if forces_dtmp_c is provided */
rc = qmckl_provide_forces_dtmp_c(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if ne distance is provided */
rc = qmckl_provide_en_distance(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_forces_een_rescaled_n_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
}
/* Compute if necessary */
if (ctx->date > ctx->forces.forces_jastrow_een_l_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
if (ctx->forces.forces_jastrow_een_l != NULL) {
rc = qmckl_free(context, ctx->forces.forces_jastrow_een_l);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_provide_forces_jastrow_een_l",
"Unable to free ctx->forces.forces_jastrow_een_l");
}
ctx->forces.forces_jastrow_een_l = NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_jastrow_een_l == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 3 * ctx->electron.walker.num * ctx->nucleus.num * sizeof(double);
double* forces_jastrow_een_l = (double*) qmckl_malloc(context, mem_info);
if (forces_jastrow_een_l == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_provide_forces_jastrow_een_l",
NULL);
}
ctx->forces.forces_jastrow_een_l = forces_jastrow_een_l;
}
rc = qmckl_compute_forces_jastrow_een_l(context,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow_champ.cord_num,
ctx->jastrow_champ.dim_c_vector,
ctx->jastrow_champ.c_vector_full,
ctx->jastrow_champ.lkpm_combined_index,
ctx->electron.en_distance,
ctx->jastrow_champ.tmp_c,
ctx->jastrow_champ.dtmp_c,
ctx->forces.forces_tmp_c,
ctx->forces.forces_dtmp_c,
ctx->forces.forces_een_n,
ctx->jastrow_champ.een_rescaled_n,
ctx->jastrow_champ.een_rescaled_n_gl,
ctx->forces.forces_jastrow_een_l);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->forces.forces_jastrow_een_l_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
** Compute
:PROPERTIES:
:Name: qmckl_compute_forces_jastrow_een_l
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_forces_jastrow_een_l_args
| 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_n_gl~ | ~double[walk_num][0:cord_num][nucl_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 |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
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_n_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_n_gl(elec_num, 4, nucl_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
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_forces_jastrow_een_l(
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t cord_num,
const int64_t dim_c_vector,
const double* c_vector_full,
const int64_t* lkpm_combined_index,
const double* en_distance,
const double* tmp_c,
const double* dtmp_c,
const double* forces_tmp_c,
const double* forces_dtmp_c,
const double* forces_een_n,
const double* een_rescaled_n,
const double* een_rescaled_n_gl,
double* const forces_jastrow_een_l );
#+end_src
** Test
#+begin_src c :tangle (eval c_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");
#+end_src
* Force of single en Jastrow
** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_forces_jastrow_single_en(qmckl_context context,
double* const forces_jastrow_single_en,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_forces_jastrow_single_en(qmckl_context context,
double* const forces_jastrow_single_en,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_forces_jastrow_single_en(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->electron.walker.num * 3 * ctx->nucleus.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_forces_jastrow_single_en",
"input array too small");
}
memcpy(forces_jastrow_single_en, ctx->forces.forces_jastrow_single_en, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(qmckl_exit_code) function qmckl_get_forces_jastrow_single_en (context, &
forces_jastrow_single_en, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: size_max
real(c_double), intent(out) :: forces_jastrow_single_en(size_max)
end function qmckl_get_forces_jastrow_single_en
end interface
#+end_src
** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_jastrow_single_en(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :export none
qmckl_exit_code qmckl_provide_forces_jastrow_single_en(qmckl_context context)
{
qmckl_exit_code rc = QMCKL_SUCCESS;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_provide_forces_jastrow_single_en",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->jastrow_champ.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_provide_forces_jastrow_single_en",
NULL);
}
/* Check if en rescaled distance is provided */
rc = qmckl_provide_en_distance_rescaled(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_en_distance_rescaled_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance is provided */
rc = qmckl_provide_en_rescaled_single(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_en_rescaled_single_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Compute if necessary */
if (ctx->single_point.date > ctx->forces.forces_jastrow_single_en_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
if (ctx->forces.forces_jastrow_single_en != NULL) {
rc = qmckl_free(context, ctx->forces.forces_jastrow_single_en);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_provide_forces_jastrow_single_en",
"Unable to free ctx->forces.forces_jastrow_single_en");
}
ctx->forces.forces_jastrow_single_en = NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_jastrow_single_en == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walker.num * 3 * ctx->nucleus.num * sizeof(double);
double* forces_jastrow_single_en = (double*) qmckl_malloc(context, mem_info);
if (forces_jastrow_single_en == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_provide_forces_jastrow_single_en",
NULL);
}
ctx->forces.forces_jastrow_single_en = forces_jastrow_single_en;
}
rc = qmckl_compute_forces_jastrow_single_en(context,
ctx->single_point.num,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow_champ.type_nucl_num,
ctx->jastrow_champ.type_nucl_vector,
ctx->jastrow_champ.aord_num,
ctx->jastrow_champ.a_vector,
ctx->jastrow_champ.en_distance_rescaled,
ctx->jastrow_champ.en_distance_rescaled_gl,
ctx->single_point.en_rescaled_single,
ctx->single_point.en_rescaled_single_gl,
ctx->forces.forces_jastrow_single_en);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->forces.forces_jastrow_single_en_date = ctx->single_point.date;
}
return QMCKL_SUCCESS;
}
#+end_src
** Compute
:PROPERTIES:
:Name: qmckl_compute_forces_jastrow_single_en
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_forces_single_en_args
| Variable | Type | In/Out | Description |
|---------------------------+-------------------------------------------+--------+---------------------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~num~ | ~int64_t~ | in | Index of single electron |
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~type_nucl_num~ | ~int64_t~ | in | Number of unique nuclei |
| ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | IDs of unique nuclei |
| ~aord_num~ | ~int64_t~ | in | Number of coefficients |
| ~a_vector~ | ~double[type_nucl_num][aord_num+1]~ | in | List of coefficients |
| ~en_distance_rescaled~ | ~double[walk_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled distances |
| ~en_distance_rescaled_gl~ | ~double[walk_num][nucl_num][elec_num][4]~ | in | Electron-nucleus rescaled distance derivatives |
| ~en_rescaled_single~ | ~double[walk_num][nucl_num]~ | in | Electron-nucleus rescaled single distances |
| ~en_rescaled_single_gl~ | ~double[walk_num][nucl_num][4]~ | in | Electron-nucleus rescaled single distance derivatives |
| ~forces_jastrow_single_en~ | ~double[walk_num][nucl_num][3]~ | out | Single electron-nucleus forces |
|---------------------------+-------------------------------------------+--------+---------------------------------------------------------|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
function qmckl_compute_forces_jastrow_single_en_doc( &
context, num_in, walk_num, elec_num, nucl_num, type_nucl_num, &
type_nucl_vector, aord_num, a_vector, &
en_distance_rescaled, en_distance_rescaled_gl, en_rescaled_single, en_rescaled_single_gl, forces_jastrow_single_en) &
bind(C) result(info)
use qmckl
implicit none
integer (qmckl_context), intent(in), value :: context
integer (c_int64_t) , intent(in) , value :: num_in
integer (c_int64_t) , intent(in) , value :: walk_num
integer (c_int64_t) , intent(in) , value :: elec_num
integer (c_int64_t) , intent(in) , value :: nucl_num
integer (c_int64_t) , intent(in) , value :: type_nucl_num
integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num)
integer (c_int64_t) , intent(in) , value :: aord_num
real (c_double ) , intent(in) :: a_vector(aord_num+1,type_nucl_num)
real (c_double ) , intent(in) :: en_distance_rescaled(elec_num,nucl_num,walk_num)
real (c_double ) , intent(in) :: en_distance_rescaled_gl(4, elec_num,nucl_num,walk_num)
real (c_double ) , intent(in) :: en_rescaled_single(nucl_num,walk_num)
real (c_double ) , intent(in) :: en_rescaled_single_gl(4, nucl_num,walk_num)
real (c_double ) , intent(out) :: forces_jastrow_single_en(3,nucl_num,walk_num)
integer(qmckl_exit_code) :: info
integer*8 :: i, a, k, nw, ii, num
double precision :: x, x1, kf, x_old, x1_old
double precision :: denom, invdenom, invdenom2, f
double precision :: denom_old, invdenom_old, invdenom2_old, f_old
double precision :: dx(3), dx_old(3)
num = num_in + 1
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif
if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif
if (nucl_num <= 0) then
info = QMCKL_INVALID_ARG_5
return
endif
if (aord_num < 0) then
info = QMCKL_INVALID_ARG_8
return
endif
do nw =1, walk_num
forces_jastrow_single_en(:,:,nw) = 0.0d0
do a = 1, nucl_num
x_old = en_distance_rescaled(num,a,nw)
x = en_rescaled_single(a,nw)
denom = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x
invdenom = 1.0d0 / denom
invdenom2 = invdenom*invdenom
denom_old = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x_old
invdenom_old = 1.0d0 / denom_old
invdenom2_old = invdenom_old*invdenom_old
dx(1) = -en_rescaled_single_gl(1,a,nw)
dx(2) = -en_rescaled_single_gl(2,a,nw)
dx(3) = -en_rescaled_single_gl(3,a,nw)
dx_old(1) = -en_distance_rescaled_gl(1,num,a,nw)
dx_old(2) = -en_distance_rescaled_gl(2,num,a,nw)
dx_old(3) = -en_distance_rescaled_gl(3,num,a,nw)
f = a_vector(1, type_nucl_vector(a)+1) * invdenom2
f_old = a_vector(1, type_nucl_vector(a)+1) * invdenom2_old
forces_jastrow_single_en(1,a,nw) = forces_jastrow_single_en(1,a,nw) + f * dx(1) - f_old * dx_old(1)
forces_jastrow_single_en(2,a,nw) = forces_jastrow_single_en(2,a,nw) + f * dx(2) - f_old * dx_old(2)
forces_jastrow_single_en(3,a,nw) = forces_jastrow_single_en(3,a,nw) + f * dx(3) - f_old * dx_old(3)
kf = 2.d0
x1 = x
x = 1.d0
x1_old = x_old
x_old = 1.d0
do k=2, aord_num
f = a_vector(k+1,type_nucl_vector(a)+1) * kf * x
f_old = a_vector(k+1,type_nucl_vector(a)+1) * kf * x_old
forces_jastrow_single_en(1,a,nw) = forces_jastrow_single_en(1,a,nw) + f * x1 * dx(1) - f_old * x1_old * dx_old(1)
forces_jastrow_single_en(2,a,nw) = forces_jastrow_single_en(2,a,nw) + f * x1 * dx(2) - f_old * x1_old * dx_old(2)
forces_jastrow_single_en(3,a,nw) = forces_jastrow_single_en(3,a,nw) + f * x1 * dx(3) - f_old * x1_old * dx_old(3)
x = x*x1
x_old = x_old*x1_old
kf = kf + 1.d0
end do
end do
end do
end function qmckl_compute_forces_jastrow_single_en_doc
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_forces_jastrow_single_en_doc (
const qmckl_context context,
const int64_t num,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
const int64_t* type_nucl_vector,
const int64_t aord_num,
const double* a_vector,
const double* en_distance_rescaled,
const double* en_distance_rescaled_gl,
const double* en_rescaled_single,
const double* en_rescaled_single_gl,
double* const forces_jastrow_single_en );
qmckl_exit_code qmckl_compute_forces_jastrow_single_en (
const qmckl_context context,
const int64_t num,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
const int64_t* type_nucl_vector,
const int64_t aord_num,
const double* a_vector,
const double* en_distance_rescaled,
const double* en_distance_rescaled_gl,
const double* en_rescaled_single,
const double* en_rescaled_single_gl,
double* const forces_jastrow_single_en );
#+end_src
#+begin_src c :tangle (eval c) :comments org :exports none
qmckl_exit_code
qmckl_compute_forces_jastrow_single_en (const qmckl_context context,
const int64_t num,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t type_nucl_num,
const int64_t* type_nucl_vector,
const int64_t aord_num,
const double* a_vector,
const double* en_distance_rescaled,
const double* en_distance_rescaled_gl,
const double* en_rescaled_single,
const double* en_rescaled_single_gl,
double* const forces_jastrow_single_en )
{
#ifdef HAVE_HPC
return qmckl_compute_forces_jastrow_single_en_doc
#else
return qmckl_compute_forces_jastrow_single_en_doc
#endif
(context, num, walk_num, elec_num, nucl_num, type_nucl_num, type_nucl_vector, aord_num,
a_vector, en_distance_rescaled, en_distance_rescaled_gl, en_rescaled_single, en_rescaled_single_gl, forces_jastrow_single_en );
}
#+end_src
** Test
#+begin_src c :tangle (eval c_test)
printf("Forces Jastrow single 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 new_coords[6] = {1.0,2.0,3.0,4.0,5.0,6.0};
rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3*walk_num);
assert (rc == QMCKL_SUCCESS);
double forces_jastrow_single_en[walk_num][nucl_num][3];
rc = qmckl_get_forces_jastrow_single_en(context, &forces_jastrow_single_en[0][0][0], 3*nucl_num*walk_num);
assert(rc == QMCKL_SUCCESS);
double finite_difference_force_single_en[walk_num][nucl_num][3];
rc = qmckl_finite_difference_deriv_n(context, delta_x, &qmckl_get_jastrow_champ_single_en, &(finite_difference_force_single_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_single_en[nw][a][k]);
//printf("%.10f\n", forces_jastrow_single_en[nw][a][k]);
assert(fabs(finite_difference_force_single_en[nw][a][k] - forces_jastrow_single_en[nw][a][k]) < 1.e-8);
}
}
}
printf("OK\n");
#+end_src
* Force of delta_p matrix
** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_forces_jastrow_delta_p(qmckl_context context,
double* const forces_delta_p,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_forces_jastrow_delta_p(qmckl_context context,
double* const forces_delta_p,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_forces_jastrow_delta_p(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = 3 * ctx->electron.walker.num * ctx->jastrow_champ.cord_num *
(ctx->jastrow_champ.cord_num + 1) * ctx->nucleus.num * ctx->electron.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_forces_jastrow_delta_p",
"Array too small.");
}
memcpy(forces_delta_p, ctx->forces.forces_delta_p, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
** provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_jastrow_delta_p(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_jastrow_delta_p(qmckl_context context)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
qmckl_exit_code rc = qmckl_provide_een_rescaled_single_e(context);
if(rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_een_rescaled_n(context);
if(rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_een_rescaled_e(context);
if(rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_een_rescaled_n_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_een_rescaled_single_n(context);
if(rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_een_rescaled_single_e(context);
if(rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_een_rescaled_single_n_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Compute if necessary */
if (ctx->single_point.date > ctx->forces.forces_delta_p_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
if (ctx->forces.forces_delta_p != NULL) {
rc = qmckl_free(context, ctx->forces.forces_delta_p);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_provide_forces_jastrow_delta_p",
"Unable to free ctx->forces.forces_delta_p");
}
ctx->forces.forces_delta_p = NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_delta_p == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 3 * ctx->electron.walker.num * ctx->jastrow_champ.cord_num *
(ctx->jastrow_champ.cord_num + 1) * ctx->nucleus.num * ctx->electron.num * sizeof(double);
double* forces_delta_p = (double*) qmckl_malloc(context, mem_info);
if (forces_delta_p == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_provide_forces_jastrow_delta_p",
NULL);
}
ctx->forces.forces_delta_p = forces_delta_p;
}
rc = qmckl_compute_forces_jastrow_delta_p(context,
ctx->single_point.num,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow_champ.cord_num,
ctx->jastrow_champ.een_rescaled_n,
ctx->jastrow_champ.een_rescaled_e,
ctx->single_point.een_rescaled_single_n,
ctx->single_point.een_rescaled_single_e,
ctx->jastrow_champ.een_rescaled_n_gl,
ctx->single_point.een_rescaled_single_n_gl,
ctx->forces.forces_delta_p);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->forces.forces_delta_p_date = ctx->single_point.date;
}
return QMCKL_SUCCESS;
}
#+end_src
** Compute
:PROPERTIES:
:Name: qmckl_compute_forces_jastrow_delta_p_doc
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_forces_delta_p_args
| Variable | Type | In/Out | Description |
|----------------------------+---------------------------------------------------------------------+--------+---------------------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~num~ | ~int64_t~ | in | Index of single electron |
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~cord_num~ | ~int64_t~ | in | order of polynomials |
| ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled distances |
| ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | in | Electron-electron rescaled distances |
| ~een_rescaled_single_n~ | ~double[walk_num][0:cord_num][nucl_num]~ | in | Electron-nucleus single rescaled distances |
| ~een_rescaled_single_e~ | ~double[walk_num][0:cord_num][elec_num]~ | in | Electron-electron single rescaled distances |
| ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Electron-nucleus rescaled distances derivatives |
| ~een_rescaled_single_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4]~ | in | Electron-nucleus single rescaled distances derivatives |
| ~forces_delta_p~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][3][elec_num]~ | out | Delta P matrix gradient and Laplacian |
|----------------------------+---------------------------------------------------------------------+--------+---------------------------------------------------------|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_delta_p_doc( &
context, num_in, walk_num, elec_num, nucl_num, cord_num, &
een_rescaled_n, een_rescaled_e, een_rescaled_single_n, een_rescaled_single_e, &
een_rescaled_n_gl, een_rescaled_single_n_gl, forces_delta_p) &
result(info) bind(C)
use, intrinsic :: iso_c_binding
use qmckl
implicit none
integer(qmckl_context), intent(in) :: context
integer(c_int64_t) , intent(in), value :: num_in, walk_num, elec_num, cord_num, nucl_num
real(c_double) , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
real(c_double) , intent(in) :: een_rescaled_e(elec_num, elec_num, 0:cord_num, walk_num)
real(c_double) , intent(in) :: een_rescaled_single_n(nucl_num, 0:cord_num, walk_num)
real(c_double) , intent(in) :: een_rescaled_single_e(elec_num, 0:cord_num, walk_num)
real(c_double) , intent(in) :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num)
real(c_double) , intent(in) :: een_rescaled_single_n_gl(4, nucl_num, 0:cord_num, walk_num)
real(c_double) , intent(out) :: forces_delta_p(elec_num,3,nucl_num,0:cord_num, 0:cord_num-1, walk_num)
double precision :: delta_n_gl
double precision :: delta_e(elec_num)
integer*8 :: i, a, j, l, k, p, m, n, nw, num
double precision :: tmp
integer*8 :: LDA, LDB, LDC
num = num_in + 1
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (walk_num <= 0) info = QMCKL_INVALID_ARG_3
if (elec_num <= 0) info = QMCKL_INVALID_ARG_4
if (nucl_num <= 0) info = QMCKL_INVALID_ARG_5
if (cord_num < 0) info = QMCKL_INVALID_ARG_6
if (info /= QMCKL_SUCCESS) return
if (cord_num == 0) then
forces_delta_p = 0.d0
return
endif
do nw=1, walk_num
do m=0, cord_num-1
do j = 1, elec_num
delta_e(j) = een_rescaled_single_e(j,m,nw) - een_rescaled_e(j,num,m,nw)
end do
do l=0, cord_num
do k = 1, 3
do a = 1, nucl_num
delta_n_gl = een_rescaled_single_n_gl(k,a,l,nw) - een_rescaled_n_gl(num, k, a, l, nw)
do j = 1, elec_num
forces_delta_p(j,k,a,l,m,nw) = delta_e(j) * (-een_rescaled_single_n_gl(k,a,l,nw)) - &
een_rescaled_e(j,num,m,nw) * delta_n_gl
end do
do j = 1, elec_num
forces_delta_p(num,k,a,l,m,nw) = forces_delta_p(num,k,a,l,m,nw) - delta_e(j) * een_rescaled_n_gl(j,k,a,l,nw)
end do
end do
end do
end do
end do
end do
end function qmckl_compute_forces_jastrow_delta_p_doc
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code
qmckl_compute_forces_jastrow_delta_p_doc (const qmckl_context context,
const int64_t num,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t cord_num,
const double* een_rescaled_n,
const double* een_rescaled_e,
const double* een_rescaled_single_n,
const double* een_rescaled_single_e,
const double* een_rescaled_n_gl,
const double* een_rescaled_single_n_gl,
double* const forces_delta_p );
qmckl_exit_code
qmckl_compute_forces_jastrow_delta_p (const qmckl_context context,
const int64_t num,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t cord_num,
const double* een_rescaled_n,
const double* een_rescaled_e,
const double* een_rescaled_single_n,
const double* een_rescaled_single_e,
const double* een_rescaled_n_gl,
const double* een_rescaled_single_n_gl,
double* const forces_delta_p );
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_compute_forces_jastrow_delta_p (const qmckl_context context,
const int64_t num,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t cord_num,
const double* een_rescaled_n,
const double* een_rescaled_e,
const double* een_rescaled_single_n,
const double* een_rescaled_single_e,
const double* een_rescaled_n_gl,
const double* een_rescaled_single_n_gl,
double* const forces_delta_p )
{
#ifdef HAVE_HPC
return qmckl_compute_forces_jastrow_delta_p_doc
#else
return qmckl_compute_forces_jastrow_delta_p_doc
#endif
(context, num, walk_num, elec_num, nucl_num, cord_num,
een_rescaled_n, een_rescaled_e, een_rescaled_single_n, een_rescaled_single_e,
een_rescaled_n_gl, een_rescaled_single_n_gl, forces_delta_p);
}
#+end_src
** Test
#+begin_src c :tangle (eval c_test)
printf("Forces delta p\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);
rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3*walk_num);
assert (rc == QMCKL_SUCCESS);
double forces_delta_p[walk_num][cord_num][cord_num+1][nucl_num][3][elec_num];
rc = qmckl_get_forces_jastrow_delta_p(context, &forces_delta_p[0][0][0][0][0][0], 3*nucl_num*walk_num*elec_num*(cord_num+1)*cord_num);
assert(rc == QMCKL_SUCCESS);
double finite_difference_force_delta_p[walk_num][cord_num][cord_num+1][nucl_num][3][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 output_delta_p[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);
rc = qmckl_single_touch(context);
assert(rc == QMCKL_SUCCESS);
// Call the provided function
rc = qmckl_get_jastrow_champ_delta_p(context,
&output_delta_p[0][0][0][0][0],
nucl_num*walk_num*elec_num*(cord_num+1)*cord_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 mm = 0; mm <= cord_num; mm++) {
for (int i = 0; i < elec_num; i++) {
if (m == -4) {
finite_difference_force_delta_p[nw][l][mm][a][k][i] = 0.0;
}
finite_difference_force_delta_p[nw][l][mm][a][k][i] += coef[m + 4] * output_delta_p[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_delta_p[nw][l][m][a][k][i]);
//printf("%.10f\n", forces_delta_p[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_delta_p[nw][l][m][a][k][i] - forces_delta_p[nw][l][m][a][k][i]) < 1.e-8);
}
}
}
}
}
}
printf("OK\n");
#+end_src
* Force of single een Jastrow
** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_forces_jastrow_single_een(qmckl_context context,
double* const forces_jastrow_single_een,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_forces_jastrow_single_een(qmckl_context context,
double* const forces_jastrow_single_een,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_forces_jastrow_single_een(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->electron.walker.num * 3 * ctx->nucleus.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_forces_jastrow_single_een",
"input array too small");
}
memcpy(forces_jastrow_single_een, ctx->forces.forces_jastrow_single_een, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(qmckl_exit_code) function qmckl_get_forces_jastrow_single_een (context, &
forces_jastrow_single_een, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: size_max
real(c_double), intent(out) :: forces_jastrow_single_een(size_max)
end function qmckl_get_forces_jastrow_single_een
end interface
#+end_src
** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_jastrow_single_een(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :export none
qmckl_exit_code qmckl_provide_forces_jastrow_single_een(qmckl_context context)
{
qmckl_exit_code rc = QMCKL_SUCCESS;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_provide_forces_jastrow_single_een",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->jastrow_champ.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_provide_forces_jastrow_single_een",
NULL);
}
/* Check if en rescaled distance is provided */
rc = qmckl_provide_een_rescaled_n(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_een_rescaled_n_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance is provided */
rc = qmckl_provide_een_rescaled_single_n(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Check if en rescaled distance derivatives is provided */
rc = qmckl_provide_een_rescaled_single_n_gl(context);
if(rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_tmp_c(context);
if(rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_forces_tmp_c(context);
if(rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_jastrow_champ_delta_p(context);
if(rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_forces_jastrow_delta_p(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Compute if necessary */
if (ctx->single_point.date > ctx->forces.forces_jastrow_single_een_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
if (ctx->forces.forces_jastrow_single_een != NULL) {
rc = qmckl_free(context, ctx->forces.forces_jastrow_single_een);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_provide_forces_jastrow_single_een",
"Unable to free ctx->forces.forces_jastrow_single_een");
}
ctx->forces.forces_jastrow_single_een = NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_jastrow_single_een == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walker.num * 3 * ctx->nucleus.num * sizeof(double);
double* forces_jastrow_single_een = (double*) qmckl_malloc(context, mem_info);
if (forces_jastrow_single_een == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_provide_forces_jastrow_single_een",
NULL);
}
ctx->forces.forces_jastrow_single_een = forces_jastrow_single_een;
}
rc = qmckl_compute_forces_jastrow_single_een(context,
ctx->single_point.num,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->jastrow_champ.cord_num,
ctx->jastrow_champ.dim_c_vector,
ctx->jastrow_champ.c_vector_full,
ctx->jastrow_champ.lkpm_combined_index,
ctx->single_point.delta_p,
ctx->forces.forces_delta_p,
ctx->jastrow_champ.tmp_c,
ctx->forces.forces_tmp_c,
ctx->jastrow_champ.een_rescaled_n,
ctx->single_point.een_rescaled_single_n,
ctx->jastrow_champ.een_rescaled_n_gl,
ctx->single_point.een_rescaled_single_n_gl,
ctx->forces.forces_jastrow_single_een);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->forces.forces_jastrow_single_een_date = ctx->single_point.date;
}
return QMCKL_SUCCESS;
}
#+end_src
** Compute
:PROPERTIES:
:Name: qmckl_compute_forces_jastrow_single_een
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_forces_single_een_args
| Variable | Type | In/Out | Description |
|---------------------------+-------------------------------------------+--------+---------------------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~num~ | ~int64_t~ | in | Index of single electron |
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~cord_num~ | ~int64_t~ | in | order of polynomials |
| ~dim_c_vector~ | ~int64_t~ | in | dimension of full coefficient vector |
| ~c_vector_full~ | ~double[dim_c_vector][nucl_num]~ | in | full coefficient vector |
| ~lkpm_combined_index~ | ~int64_t[4][dim_c_vector]~ | in | combined indices |
| ~delta_p~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | in | Single electron P matrix |
| ~forces_delta_p~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][3][elec_num]~ | in | Single electron P matrix |
| ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | in | Single electron P matrix |
| ~forces_tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | in | Single electron P matrix |
| ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled distances |
| ~een_rescaled_single_n~ | ~double[walk_num][0:cord_num][nucl_num]~ | in | Electron-nucleus single rescaled distances |
| ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Electron-nucleus rescaled distances derivatives |
| ~een_rescaled_single_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4]~ | in | Electron-nucleus single rescaled distances derivatives |
| ~forces_jastrow_single_een~ | ~double[walk_num][nucl_num][3]~ | out | Single electron-nucleus forces |
|---------------------------+-------------------------------------------+--------+---------------------------------------------------------|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
function qmckl_compute_forces_jastrow_single_een_doc( &
context, num_in, walk_num, elec_num, nucl_num, cord_num, &
dim_c_vector, c_vector_full, lkpm_combined_index, &
delta_p, forces_delta_p, tmp_c, forces_tmp_c, &
een_rescaled_n, een_rescaled_single_n, een_rescaled_n_gl, een_rescaled_single_n_gl, forces_jastrow_single_een) &
bind(C) result(info)
use qmckl
implicit none
integer (qmckl_context), intent(in), value :: context
integer (c_int64_t) , intent(in) , value :: num_in
integer (c_int64_t) , intent(in) , value :: walk_num
integer (c_int64_t) , intent(in) , value :: elec_num
integer (c_int64_t) , intent(in) , value :: nucl_num
integer (c_int64_t) , intent(in) , value :: dim_c_vector
integer (c_int64_t) , intent(in) , value :: cord_num
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) :: delta_p(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num)
real (c_double ) , intent(in) :: forces_delta_p(elec_num, 3, nucl_num,0:cord_num, 0:cord_num-1, 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(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
real (c_double ) , intent(in) :: een_rescaled_single_n(nucl_num, 0:cord_num, walk_num)
real (c_double ) , intent(in) :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num)
real (c_double ) , intent(in) :: een_rescaled_single_n_gl(4, nucl_num, 0:cord_num, walk_num)
real (c_double ) , intent(out) :: forces_jastrow_single_een(3,nucl_num,walk_num)
integer(qmckl_exit_code) :: info
double precision :: een_rescaled_delta_n(nucl_num, 0:cord_num), een_rescaled_delta_n_gl(nucl_num, 0:cord_num)
integer*8 :: i, a, j, l, k, p, m, n, nw, num, kk
double precision :: accu, accu2, cn
integer*8 :: LDA, LDB, LDC
num = num_in + 1
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
if (walk_num <= 0) info = QMCKL_INVALID_ARG_3
if (elec_num <= 0) info = QMCKL_INVALID_ARG_4
if (nucl_num <= 0) info = QMCKL_INVALID_ARG_5
if (cord_num < 0) info = QMCKL_INVALID_ARG_6
if (info /= QMCKL_SUCCESS) return
forces_jastrow_single_een = 0.0d0
if (cord_num == 0) return
do nw =1, walk_num
do kk = 1, 3
een_rescaled_delta_n(:,:) = een_rescaled_single_n(:,:,nw) - een_rescaled_n(num,:,:,nw)
een_rescaled_delta_n_gl(:,:) = een_rescaled_single_n_gl(kk,:,:,nw) - een_rescaled_n_gl(num,kk,:,:,nw)
do n = 1, dim_c_vector
l = lkpm_combined_index(n, 1)
k = lkpm_combined_index(n, 2)
p = lkpm_combined_index(n, 3)
m = lkpm_combined_index(n, 4)
do a = 1, nucl_num
cn = c_vector_full(a, n)
if(cn == 0.d0) cycle
accu = 0.0d0
do j = 1, elec_num
accu = accu - een_rescaled_n_gl(j,kk,a,m,nw) * delta_p(j,a,m+l,k,nw) + &
een_rescaled_n(j,a,m,nw) * forces_delta_p(j,kk,a,m+l,k,nw)
end do
accu = accu - een_rescaled_delta_n_gl(a,m) * (tmp_c(num,a,m+l,k,nw) + delta_p(num,a,m+l,k,nw)) + &
een_rescaled_delta_n(a,m) * (forces_tmp_c(num,kk,a,m+l,k,nw) + forces_delta_p(num,kk,a,m+l,k,nw))
forces_jastrow_single_een(kk,a,nw) = forces_jastrow_single_een(kk,a,nw) + accu * cn
end do
end do
end do
end do
end function qmckl_compute_forces_jastrow_single_een_doc
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_forces_jastrow_single_een_doc (
const qmckl_context context,
const int64_t num,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t cord_num,
const int64_t dim_c_vector,
const double* c_vector_full,
const int64_t* lkpm_combined_index,
const double* delta_p,
const double* forces_delta_p,
const double* tmp_c,
const double* forces_tmp_c,
const double* een_rescaled_n,
const double* een_rescaled_single_n,
const double* een_rescaled_n_gl,
const double* een_rescaled_single_n_gl,
double* const forces_jastrow_single_een );
qmckl_exit_code qmckl_compute_forces_jastrow_single_een (
const qmckl_context context,
const int64_t num,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t cord_num,
const int64_t dim_c_vector,
const double* c_vector_full,
const int64_t* lkpm_combined_index,
const double* delta_p,
const double* forces_delta_p,
const double* tmp_c,
const double* forces_tmp_c,
const double* een_rescaled_n,
const double* een_rescaled_single_n,
const double* een_rescaled_n_gl,
const double* een_rescaled_single_n_gl,
double* const forces_jastrow_single_een );
#+end_src
#+begin_src c :tangle (eval c) :comments org :exports none
qmckl_exit_code
qmckl_compute_forces_jastrow_single_een (const qmckl_context context,
const int64_t num,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t cord_num,
const int64_t dim_c_vector,
const double* c_vector_full,
const int64_t* lkpm_combined_index,
const double* delta_p,
const double* forces_delta_p,
const double* tmp_c,
const double* forces_tmp_c,
const double* een_rescaled_n,
const double* een_rescaled_single_n,
const double* een_rescaled_n_gl,
const double* een_rescaled_single_n_gl,
double* const forces_jastrow_single_een )
{
#ifdef HAVE_HPC
return qmckl_compute_forces_jastrow_single_een_doc
#else
return qmckl_compute_forces_jastrow_single_een_doc
#endif
(context, num, walk_num, elec_num, nucl_num, cord_num, dim_c_vector, c_vector_full, lkpm_combined_index,
delta_p, forces_delta_p, tmp_c, forces_tmp_c, een_rescaled_n, een_rescaled_single_n,
een_rescaled_n_gl, een_rescaled_single_n_gl, forces_jastrow_single_een );
}
#+end_src
** Test
#+begin_src c :tangle (eval c_test)
printf("Forces Jastrow single 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);
rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3*walk_num);
assert (rc == QMCKL_SUCCESS);
double forces_jastrow_single_een[walk_num][nucl_num][3];
rc = qmckl_get_forces_jastrow_single_een(context, &forces_jastrow_single_een[0][0][0], 3*nucl_num*walk_num);
assert(rc == QMCKL_SUCCESS);
double finite_difference_force_single_een[walk_num][nucl_num][3];
rc = qmckl_finite_difference_deriv_n(context, delta_x, &qmckl_get_jastrow_champ_single_een, &(finite_difference_force_single_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_single_een[nw][a][k]);
//printf("%.10f\n", forces_jastrow_single_een[nw][a][k]);
assert(fabs(finite_difference_force_single_een[nw][a][k] - forces_jastrow_single_een[nw][a][k]) < 1.e-8);
}
}
}
printf("OK\n");
#+end_src
* Reset test for orbitals
#+begin_src c :tangle (eval c_test)
rc = qmckl_context_destroy(context);
assert (rc == QMCKL_SUCCESS);
context = qmckl_context_create();
nucl_num = chbrclf_nucl_num;
nucl_charge = chbrclf_charge;
nucl_coord = &(chbrclf_nucl_coord[0][0]);
rc = qmckl_set_nucleus_num (context, nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));
const int64_t shell_num = chbrclf_shell_num;
const int64_t prim_num = chbrclf_prim_num;
const int64_t ao_num = chbrclf_ao_num;
const int64_t * nucleus_index = &(chbrclf_basis_nucleus_index[0]);
const int64_t * nucleus_shell_num = &(chbrclf_basis_nucleus_shell_num[0]);
const int32_t * shell_ang_mom = &(chbrclf_basis_shell_ang_mom[0]);
const int64_t * shell_prim_num = &(chbrclf_basis_shell_prim_num[0]);
const int64_t * shell_prim_index = &(chbrclf_basis_shell_prim_index[0]);
const double * shell_factor = &(chbrclf_basis_shell_factor[0]);
const double * exponent = &(chbrclf_basis_exponent[0]);
const double * coefficient = &(chbrclf_basis_coefficient[0]);
const double * prim_factor = &(chbrclf_basis_prim_factor[0]);
const double * ao_factor = &(chbrclf_basis_ao_factor[0]);
const char typ = 'G';
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_type (context, typ);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_num (context, shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_prim_num (context, prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, nucl_num);
assert(rc == QMCKL_ALREADY_SET);
rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_factor (context, shell_factor, shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_exponent (context, exponent, prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_coefficient (context, coefficient, prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_prim_factor (context, prim_factor, prim_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_ao_basis_ao_num(context, ao_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_ao_basis_ao_factor (context, ao_factor, ao_num);
assert(rc == QMCKL_SUCCESS);
assert(qmckl_ao_basis_provided(context));
int64_t shell_num_test ;
int64_t prim_num_test ;
int64_t ao_num_test ;
int64_t * nucleus_index_test ;
int64_t * nucleus_shell_num_test;
int32_t * shell_ang_mom_test ;
int64_t * shell_prim_num_test ;
int64_t * shell_prim_index_test ;
double * shell_factor_test ;
double * exponent_test ;
double * coefficient_test ;
double * prim_factor_test ;
double * ao_factor_test ;
char typ_test ;
rc = qmckl_get_ao_basis_type (context, &typ_test);
assert (rc == QMCKL_SUCCESS);
assert(typ == typ_test);
rc = qmckl_get_ao_basis_shell_num (context, &shell_num_test);
assert (rc == QMCKL_SUCCESS);
assert(shell_num == shell_num_test);
rc = qmckl_get_ao_basis_prim_num (context, &prim_num_test);
assert (rc == QMCKL_SUCCESS);
assert(prim_num == prim_num_test);
nucleus_index_test = (int64_t*) malloc (nucl_num * sizeof(int64_t));
rc = qmckl_get_ao_basis_nucleus_index (context, nucleus_index_test, nucl_num);
assert (rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i < nucl_num ; ++i) {
assert(nucleus_index_test[i] == nucleus_index[i]);
}
free(nucleus_index_test);
nucleus_shell_num_test = (int64_t*) malloc ( nucl_num * sizeof(int64_t));
rc = qmckl_get_ao_basis_nucleus_shell_num (context, nucleus_shell_num_test, nucl_num);
assert (rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i < nucl_num ; ++i) {
assert(nucleus_shell_num_test[i] == nucleus_shell_num[i]);
}
free(nucleus_shell_num_test);
shell_ang_mom_test = (int32_t*) malloc ( shell_num * sizeof(int32_t));
rc = qmckl_get_ao_basis_shell_ang_mom (context, shell_ang_mom_test, shell_num);
assert (rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i < shell_num ; ++i) {
assert(shell_ang_mom_test[i] == shell_ang_mom[i]);
}
free(shell_ang_mom_test);
shell_factor_test = (double*) malloc ( shell_num * sizeof(double));
rc = qmckl_get_ao_basis_shell_factor (context, shell_factor_test, shell_num);
assert (rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i < shell_num ; ++i) {
assert(shell_factor_test[i] == shell_factor[i]);
}
free(shell_factor_test);
shell_prim_num_test = (int64_t*) malloc ( shell_num * sizeof(int64_t));
rc = qmckl_get_ao_basis_shell_prim_num (context, shell_prim_num_test, shell_num);
assert (rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i < shell_num ; ++i) {
assert(shell_prim_num_test[i] == shell_prim_num[i]);
}
free(shell_prim_num_test);
shell_prim_index_test = (int64_t*) malloc ( shell_num * sizeof(int64_t));
rc = qmckl_get_ao_basis_shell_prim_index (context, shell_prim_index_test, shell_num);
assert (rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i < shell_num ; ++i) {
assert(shell_prim_index_test[i] == shell_prim_index[i]);
}
free(shell_prim_index_test);
exponent_test = (double*) malloc ( prim_num * sizeof(double));
rc = qmckl_get_ao_basis_exponent(context, exponent_test, prim_num);
assert (rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i < prim_num ; ++i) {
assert(exponent_test[i] == exponent[i]);
}
free(exponent_test);
coefficient_test = (double*) malloc ( prim_num * sizeof(double));
rc = qmckl_get_ao_basis_coefficient(context, coefficient_test, prim_num);
assert (rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i < prim_num ; ++i) {
assert(coefficient_test[i] == coefficient[i]);
}
free(coefficient_test);
prim_factor_test = (double*) malloc ( prim_num * sizeof(double));
rc = qmckl_get_ao_basis_prim_factor (context, prim_factor_test, prim_num);
assert (rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i < prim_num ; ++i) {
assert(prim_factor_test[i] == prim_factor[i]);
}
free(prim_factor_test);
rc = qmckl_get_ao_basis_ao_num(context, &ao_num_test);
assert(ao_num == ao_num_test);
ao_factor_test = (double*) malloc ( ao_num * sizeof(double));
rc = qmckl_get_ao_basis_ao_factor (context, ao_factor_test, ao_num);
assert (rc == QMCKL_SUCCESS);
for (int64_t i=0 ; i < ao_num ; ++i) {
assert(ao_factor_test[i] == ao_factor[i]);
}
free(ao_factor_test);
#define walk_num 1 // chbrclf_walk_num
#define elec_num chbrclf_elec_num
#define prim_num chbrclf_prim_num
elec_up_num = chbrclf_elec_up_num;
elec_dn_num = chbrclf_elec_dn_num;
elec_coord = &(chbrclf_elec_coord[0][0][0]);
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
int64_t point_num = elec_num;
rc = qmckl_set_point(context, 'N', point_num, elec_coord, point_num*3);
assert(rc == QMCKL_SUCCESS);
int64_t mo_num = chbrclf_mo_num;
rc = qmckl_set_mo_basis_mo_num(context, mo_num);
assert (rc == QMCKL_SUCCESS);
const double * mo_coefficient = &(chbrclf_mo_coef[0]);
rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient, chbrclf_mo_num*chbrclf_ao_num);
assert (rc == QMCKL_SUCCESS);
#+end_src
* Force of AO value
Here we compute the forces of the AO value.
** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_forces_ao_value(qmckl_context context,
double* const forces_ao_value,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_forces_ao_value(qmckl_context context,
double* const forces_ao_value,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_forces_ao_value(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t sze = ctx->ao_basis.ao_num * ctx->nucleus.num * 3 * ctx->point.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_forces_ao_value",
"Array too small. Expected walk_num*nucl_num*point_num*3");
}
memcpy(forces_ao_value, ctx->forces.forces_ao_value, sze*sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_ao_value(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_ao_value(qmckl_context context)
{
qmckl_exit_code rc = QMCKL_SUCCESS;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_provide_forces_ao_value",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_provide_forces_ao_value",
NULL);
}
/* Compute if necessary */
if (ctx->point.date > ctx->forces.forces_ao_value_date) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.ao_num * 3 * ctx->nucleus.num * ctx->point.num * sizeof(double);
if (ctx->forces.forces_ao_value != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->forces.forces_ao_value, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->forces.forces_ao_value);
assert (rc == QMCKL_SUCCESS);
ctx->forces.forces_ao_value = NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_ao_value == NULL) {
double* forces_ao_value = (double*) qmckl_malloc(context, mem_info);
if (forces_ao_value == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_forces_ao_value",
NULL);
}
ctx->forces.forces_ao_value = forces_ao_value;
}
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_ao_vgl", NULL);
}
rc = qmckl_compute_forces_ao_value_doc(context,
ctx->ao_basis.ao_num,
ctx->ao_basis.shell_num,
ctx->point.num,
ctx->nucleus.num,
ctx->ao_basis.nucleus_index,
ctx->ao_basis.nucleus_shell_num,
ctx->ao_basis.shell_ang_mom,
ctx->ao_basis.ao_factor,
ctx->ao_basis.ao_vgl,
ctx->forces.forces_ao_value);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->forces.forces_ao_value_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
** Compute
:PROPERTIES:
:Name: qmckl_compute_forces_ao_value
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_forces_ao_value_args_doc
| Variable | Type | In/Out | Description |
|-----------------------+-----------------------------------+--------+----------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~ao_num~ | ~int64_t~ | in | Number of AOs |
| ~shell_num~ | ~int64_t~ | in | Number of shells |
| ~point_num~ | ~int64_t~ | in | Number of points |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus |
| ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus |
| ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell |
| ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs |
| ~ao_vgl~ | ~double[point_num][5][shell_num]~ | in | Value, gradients and Laplacian of the shells |
| ~forces_ao_value~ | ~double[nucl_num][3][point_num][ao_num]~ | out | Forces of the AOs |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
function qmckl_compute_forces_ao_value_doc(context, &
ao_num, shell_num, point_num, nucl_num, &
nucleus_index, nucleus_shell_num, &
shell_ang_mom, ao_factor, ao_vgl, forces_ao_value) &
bind(C) result(info)
use qmckl_constants
use qmckl, only : qmckl_ao_polynomial_vgl, qmckl_get_numprec_precision
implicit none
integer (qmckl_context), intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: ao_num
integer (c_int64_t) , intent(in) , value :: shell_num
integer (c_int64_t) , intent(in) , value :: point_num
integer (c_int64_t) , intent(in) , value :: nucl_num
integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num)
integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num)
integer (c_int32_t) , intent(in) :: shell_ang_mom(shell_num)
real (c_double ) , intent(in) :: ao_factor(ao_num)
real (c_double ) , intent(in) :: ao_vgl(ao_num,5,point_num)
real (c_double ) , intent(out) :: forces_ao_value(ao_num,point_num,3,nucl_num)
integer(qmckl_exit_code) :: info
integer :: l, il, k
integer*8 :: ipoint, inucl, ishell
integer*8 :: ishell_start, ishell_end
integer :: lstart(0:20)
double precision :: x, y, z, r2
double precision :: cutoff
integer , allocatable :: ao_index(:)
allocate(ao_index(ao_num))
forces_ao_value = 0.d0
! Pre-computed data
do l=0,20
lstart(l) = l*(l+1)*(l+2)/6 +1
end do
k=1
do inucl=1,nucl_num
ishell_start = nucleus_index(inucl) + 1
ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl)
do ishell = ishell_start, ishell_end
l = shell_ang_mom(ishell)
ao_index(ishell) = k
k = k + lstart(l+1) - lstart(l)
end do
end do
info = QMCKL_SUCCESS
do ipoint = 1, point_num
do inucl=1,nucl_num
! Loop over shells
ishell_start = nucleus_index(inucl) + 1
ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl)
do ishell = ishell_start, ishell_end
k = ao_index(ishell)
l = shell_ang_mom(ishell)
do il = lstart(l), lstart(l+1)-1
forces_ao_value(k,ipoint,1,inucl) = -ao_vgl(k,2,ipoint)
forces_ao_value(k,ipoint,2,inucl) = -ao_vgl(k,3,ipoint)
forces_ao_value(k,ipoint,3,inucl) = -ao_vgl(k,4,ipoint)
k = k+1
end do
end do
end do
end do
end function qmckl_compute_forces_ao_value_doc
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_compute_forces_ao_value_doc(
const qmckl_context context,
const int64_t ao_num,
const int64_t shell_num,
const int64_t point_num,
const int64_t nucl_num,
const int64_t* nucleus_index,
const int64_t* nucleus_shell_num,
const int32_t* shell_ang_mom,
const double* ao_factor,
const double* ao_vgl,
double* const forces_ao_value );
#+end_src
** Test
#+begin_src c :tangle (eval c_test)
printf("Forces AO value\n");
rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);
double * forces_ao_value = (double*) malloc(nucl_num * 3 *point_num * ao_num *sizeof(double));
rc = qmckl_get_forces_ao_value(context, &forces_ao_value[0], 3*nucl_num*ao_num*point_num);
assert(rc == QMCKL_SUCCESS);
double * finite_difference_force_ao_value = (double*) malloc(3 *nucl_num * point_num * ao_num *sizeof(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);
temp_coord = (double*) malloc(3 * nucl_num * sizeof(double));
if (temp_coord == NULL) {
free(nucleus_coord);
return QMCKL_ALLOCATION_FAILED;
}
double * ao_output = (double*) malloc(point_num*ao_num*sizeof(double));
// 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_ao_basis_ao_value(context,&ao_output[0], point_num*ao_num);
assert(rc == QMCKL_SUCCESS);
// Accumulate derivative using finite-difference coefficients
for (int i = 0; i < point_num; i++) {
for (int j = 0; j < ao_num; j++) {
if (m == -4) {
finite_difference_force_ao_value[k*ao_num*point_num*nucl_num + a*ao_num*point_num + i*ao_num + j] = 0.0;
}
finite_difference_force_ao_value[k*ao_num*point_num*nucl_num + a*ao_num*point_num + i*ao_num + j] += coef[m + 4] * ao_output[i*ao_num + j]/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);
free(ao_output);
for (int j = 0; j < ao_num; j++){
for (int i = 0; i < point_num; i++){
for (int a = 0; a < nucl_num; a++) {
for (int k = 0; k < 3; k++){
//printf("k=%i a=%i i=%i j=%i\n", k, a, i, j);
//printf("%.10f\t", finite_difference_force_ao_value[k*ao_num*point_num*nucl_num + a*ao_num*point_num + i*ao_num + j]);
//printf("%.10f\n", forces_ao_value[a*3*ao_num*point_num + k*ao_num*point_num + i*ao_num + j]);
assert(fabs(finite_difference_force_ao_value[k*ao_num*point_num*nucl_num + a*ao_num*point_num + i*ao_num + j] - forces_ao_value[a*3*ao_num*point_num + k*ao_num*point_num + i*ao_num + j]) < 1.e-9);
}
}
}
}
free(forces_ao_value);
free(finite_difference_force_ao_value);
printf("OK\n");
#+end_src
* Force of MO value
** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_forces_mo_value(qmckl_context context,
double* const forces_mo_value,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_forces_mo_value(qmckl_context context,
double* const forces_mo_value,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_forces_mo_value(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->point.num * 3 * ctx->mo_basis.mo_num * ctx->nucleus.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_forces_mo_value",
"input array too small");
}
memcpy(forces_mo_value, ctx->forces.forces_mo_value, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(qmckl_exit_code) function qmckl_get_forces_mo_value (context, &
forces_mo_value, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real(c_double), intent(out) :: forces_mo_value(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_forces_mo_value
end interface
#+end_src
** Provide
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_mo_value(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :export none
qmckl_exit_code qmckl_provide_forces_mo_value(qmckl_context context)
{
qmckl_exit_code rc = QMCKL_SUCCESS;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_provide_forces_mo_value",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->mo_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_provide_mo_basis_mo_vgl",
NULL);
}
/* Compute if necessary */
if (ctx->point.date > ctx->forces.forces_mo_value_date) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 3 * ctx->mo_basis.mo_num * ctx->point.num * ctx->nucleus.num * sizeof(double);
if (ctx->forces.forces_mo_value != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->forces.forces_mo_value, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->forces.forces_mo_value);
assert (rc == QMCKL_SUCCESS);
ctx->forces.forces_mo_value = NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_mo_value == NULL) {
double* forces_mo_value = (double*) qmckl_malloc(context, mem_info);
if (forces_mo_value == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_forces_mo_value",
NULL);
}
ctx->forces.forces_mo_value = forces_mo_value;
}
rc = qmckl_provide_forces_ao_value(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_forces_ao_value",
NULL);
}
rc = qmckl_compute_forces_mo_value_doc(context,
ctx->nucleus.num,
ctx->ao_basis.ao_num,
ctx->mo_basis.mo_num,
ctx->point.num,
ctx->mo_basis.coefficient_t,
ctx->forces.forces_ao_value,
ctx->forces.forces_mo_value);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->forces.forces_mo_value_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
** Compute
:PROPERTIES:
:Name: qmckl_compute_forces_mo_value
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_forces_mo_value_args
| Variable | Type | In/Out | Description |
|---------------------+--------------------------------+--------+-------------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~nucl_num~ | ~int64_t~ | in | Number of AOs |
| ~ao_num~ | ~int64_t~ | in | Number of AOs |
| ~mo_num~ | ~int64_t~ | in | Number of MOs |
| ~point_num~ | ~int64_t~ | in | Number of points |
| ~coefficient_t~ | ~double[mo_num][ao_num]~ | in | Transpose of the AO to MO transformation matrix |
| ~forces_ao_value~ | ~double[nucl_num][3][point_num][ao_num]~ | in | Value, gradients and Laplacian of the AOs |
| ~forces_mo_value~ | ~double[nucl_num][3][point_num][mo_num]~ | out | Value, gradients and Laplacian of the MOs |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer(qmckl_exit_code) function qmckl_compute_forces_mo_value_doc(context, &
nucl_num,ao_num, mo_num, point_num, &
coefficient_t, forces_ao_value, forces_mo_value) &
result(info) bind(C)
use, intrinsic :: iso_c_binding
use qmckl
implicit none
integer(qmckl_context), intent(in) :: context
integer(c_int64_t) , intent(in), value ::nucl_num, ao_num, mo_num, point_num
real(c_double) , intent(in) :: forces_ao_value(ao_num,point_num,3,nucl_num)
real(c_double) , intent(in) :: coefficient_t(mo_num,ao_num)
real(c_double) , intent(out) :: forces_mo_value(mo_num,point_num,3,nucl_num)
integer*8 :: i,j,k,a
double precision :: c1, c2, c3
info = QMCKL_SUCCESS
forces_mo_value = 0.0d0
do j=1,point_num
do k=1,ao_num
do a=1,nucl_num
c1 = forces_ao_value(k,j,1,a)
c2 = forces_ao_value(k,j,2,a)
c3 = forces_ao_value(k,j,3,a)
do i=1,mo_num
forces_mo_value(i,j,1,a) = forces_mo_value(i,j,1,a) + coefficient_t(i,k) * c1
forces_mo_value(i,j,2,a) = forces_mo_value(i,j,2,a) + coefficient_t(i,k) * c2
forces_mo_value(i,j,3,a) = forces_mo_value(i,j,3,a) + coefficient_t(i,k) * c3
end do
end do
end do
end do
end function qmckl_compute_forces_mo_value_doc
#+end_src
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_compute_forces_mo_value_doc (
const qmckl_context context,
const int64_t nucl_num,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const double* coefficient_t,
const double* forces_ao_value,
double* const forces_mo_value );
#+end_src
** Test
#+begin_src c :tangle (eval c_test)
printf("Forces MO value\n");
rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);
double * forces_mo_value = (double*) malloc(nucl_num * point_num* 3 * mo_num *sizeof(double));
//double forces_mo_value[nucl_num][3][point_num][mo_num];
rc = qmckl_get_forces_mo_value(context, &forces_mo_value[0], 3*nucl_num*mo_num*point_num);
assert(rc == QMCKL_SUCCESS);
//double finite_difference_force_mo_value[3][nucl_num][point_num][mo_num];
double * finite_difference_force_mo_value = (double*) malloc(nucl_num* 3 * point_num * mo_num * sizeof(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);
temp_coord = (double*) malloc(3 * nucl_num * sizeof(double));
if (temp_coord == NULL) {
free(nucleus_coord);
return QMCKL_ALLOCATION_FAILED;
}
double* mo_output = (double*) malloc(point_num * mo_num * sizeof(double));
if (mo_output == NULL) {
free(temp_coord);
free(nucleus_coord);
return QMCKL_ALLOCATION_FAILED;
}
// 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_mo_basis_mo_value(context,&mo_output[0], point_num*mo_num);
assert(rc == QMCKL_SUCCESS);
// Accumulate derivative using finite-difference coefficients
for (int i = 0; i < point_num; i++) {
for (int j = 0; j < mo_num; j++) {
if (m == -4) {
finite_difference_force_mo_value[k*mo_num*point_num*nucl_num + a*mo_num*point_num + i*mo_num + j]= 0.0;
}
finite_difference_force_mo_value[k*mo_num*point_num*nucl_num + a*mo_num*point_num + i*mo_num + j] += coef[m + 4] * mo_output[i*mo_num+ j]/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);
free(mo_output);
for (int j = 0; j < mo_num; j++){
for (int i = 0; i < point_num; i++){
for (int a = 0; a < nucl_num; a++) {
for (int k = 0; k < 3; k++){
//printf("k=%i a=%i i=%i j=%i\n", k, a, i, j);
//printf("%.10f\t", finite_difference_force_mo_value[k*mo_num*point_num*nucl_num + a*mo_num*point_num + i*mo_num + j]);
//printf("%.10f\n", forces_mo_value[a*3*mo_num*point_num + k*mo_num*point_num + i*mo_num + j]);
assert(fabs(finite_difference_force_mo_value[k*mo_num*point_num*nucl_num + a*mo_num*point_num + i*mo_num + j] - forces_mo_value[a*3*mo_num*point_num + k*mo_num*point_num + i*mo_num + j]) < 1.e-9);
}
}
}
}
free(forces_mo_value);
free(finite_difference_force_mo_value);
printf("OK\n");
#+end_src
* Force of MO gradient
** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_forces_mo_g(qmckl_context context,
double* const forces_mo_g,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_forces_mo_g(qmckl_context context,
double* const forces_mo_g,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_forces_mo_g(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->point.num * 3 * 3 * ctx->mo_basis.mo_num * ctx->nucleus.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_forces_mo_g",
"input array too small");
}
memcpy(forces_mo_g, ctx->forces.forces_mo_g, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(qmckl_exit_code) function qmckl_get_forces_mo_g (context, &
forces_mo_g, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real(c_double), intent(out) :: forces_mo_g(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_forces_mo_g
end interface
#+end_src
** Provide
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_mo_g(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :export none
qmckl_exit_code qmckl_provide_forces_mo_g(qmckl_context context)
{
qmckl_exit_code rc = QMCKL_SUCCESS;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_provide_forces_mo_g",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->mo_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_provide_mo_basis_mo_vgl",
NULL);
}
/* Compute if necessary */
if (ctx->point.date > ctx->forces.forces_mo_g_date) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 3 * 3 * ctx->mo_basis.mo_num * ctx->point.num * ctx->nucleus.num * sizeof(double);
if (ctx->forces.forces_mo_g != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->forces.forces_mo_g, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->forces.forces_mo_g);
assert (rc == QMCKL_SUCCESS);
ctx->forces.forces_mo_g = NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_mo_g == NULL) {
double* forces_mo_g = (double*) qmckl_malloc(context, mem_info);
if (forces_mo_g == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_forces_mo_g",
NULL);
}
ctx->forces.forces_mo_g = forces_mo_g;
}
rc = qmckl_provide_ao_basis_ao_hessian(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_basis_ao_hessian",
NULL);
}
rc = qmckl_compute_forces_mo_g_doc(context,
ctx->ao_basis.ao_num,
ctx->mo_basis.mo_num,
ctx->point.num,
ctx->nucleus.num,
ctx->mo_basis.coefficient_t,
ctx->ao_basis.ao_hessian,
ctx->forces.forces_mo_g);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->forces.forces_mo_g_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
** Compute
:PROPERTIES:
:Name: qmckl_compute_forces_mo_g
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_compute_forces_mo_g_args
| Variable | Type | In/Out | Description |
|-----------------+--------------------------------+--------+-------------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~ao_num~ | ~int64_t~ | in | Number of AOs |
| ~mo_num~ | ~int64_t~ | in | Number of MOs |
| ~point_num~ | ~int64_t~ | in | Number of points |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~coefficient_t~ | ~double[mo_num][ao_num]~ | in | Transpose of the AO to MO transformation matrix |
| ~ao_hessian~ | ~double[nucl_num][3][point_num][4][ao_num]~ | in | Value, gradients and Laplacian of the AOs |
| ~forces_mo_g~ | ~double[nucl_num][3][point_num][3][mo_num]~ | out | Value, gradients and Laplacian of the MOs |
|-----------------+--------------------------------+--------+-------------------------------------------------|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_forces_mo_g_doc(context, &
ao_num, mo_num, point_num, nucl_num, &
coefficient_t, ao_hessian, forces_mo_g) &
bind(C) result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in), value :: context
integer (c_int64_t) , intent(in), value :: nucl_num, ao_num, mo_num, point_num
real (c_double ) , intent(in) :: coefficient_t(mo_num,ao_num)
real (c_double ) , intent(in) :: ao_hessian(ao_num,4,point_num,3,nucl_num)
real (c_double ) , intent(out) :: forces_mo_g(mo_num,3,point_num,3,nucl_num)
integer*8 :: i,j,k, m, n,a
info = QMCKL_SUCCESS
do j=1,point_num
forces_mo_g(:,:,j,:,:) = 0.d0
do k=1,ao_num
do i=1,mo_num
do a=1, nucl_num
do m = 1, 3
do n = 1, 3
forces_mo_g(i, m, j, n, a) = forces_mo_g(i, m, j, n, a) - coefficient_t(i,k) * ao_hessian(k, m, j, n, a)
end do
end do
end do
end do
end do
end do
end function qmckl_compute_forces_mo_g_doc
#+end_src
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_forces_mo_g_doc (
const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const int64_t nucl_num,
const double* coefficient_t,
const double* ao_hessian,
double* const forces_mo_g );
#+end_src
** Test
#+begin_src c :tangle (eval c_test)
printf("Forces MO gradient\n");
rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);
double * forces_mo_g = (double*) malloc(nucl_num* 3 * point_num* 3 * mo_num *sizeof(double));
rc = qmckl_get_forces_mo_g(context, &forces_mo_g[0], 3*3*nucl_num*mo_num*point_num);
assert(rc == QMCKL_SUCCESS);
double * finite_difference_force_mo_g = (double*) malloc(nucl_num* 3 * point_num* 3 * mo_num * sizeof(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);
temp_coord = (double*) malloc(3 * nucl_num * sizeof(double));
if (temp_coord == NULL) {
free(nucleus_coord);
return QMCKL_ALLOCATION_FAILED;
}
mo_output = (double*) malloc(5 * point_num * mo_num * sizeof(double));
if (mo_output == NULL) {
free(temp_coord);
free(nucleus_coord);
return QMCKL_ALLOCATION_FAILED;
}
// 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_mo_basis_mo_vgl(context,&mo_output[0], 5*point_num*mo_num);
assert(rc == QMCKL_SUCCESS);
// Accumulate derivative using finite-difference coefficients
for (int i = 0; i < point_num; i++) {
for (int n = 0; n < 3; n++){
for (int j = 0; j < mo_num; j++) {
if (m == -4) {
finite_difference_force_mo_g[k*3*mo_num*point_num*nucl_num + a*3*mo_num*point_num + i*3*mo_num + n*mo_num + j] = 0.0;
}
finite_difference_force_mo_g[k*3*mo_num*point_num*nucl_num + a*3*mo_num*point_num + i*3*mo_num + n*mo_num + j] += coef[m + 4] * mo_output[i*mo_num*5 + (n+1)*mo_num + j]/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);
free(mo_output);
for (int j = 0; j < mo_num; j++){
for (int n = 0; n < 3; n++){
for (int i = 0; i < point_num; i++){
for (int a = 0; a < nucl_num; a++) {
for (int k = 0; k < 3; k++){
//printf("k=%i a=%i i=%i n=%i j=%i\n", k, a, i, n, j);
//printf("%.10f\t", finite_difference_force_mo_g[k*3*mo_num*point_num*nucl_num + a*3*mo_num*point_num + i*3*mo_num + n*mo_num + j]);
//printf("%.10f\n", forces_mo_g[a*9*mo_num*point_num + k*3*mo_num*point_num + i*3*mo_num + n*mo_num + j]);
assert(fabs(finite_difference_force_mo_g[k*3*mo_num*point_num*nucl_num + a*3*mo_num*point_num + i*3*mo_num + n*mo_num + j] - forces_mo_g[a*9*mo_num*point_num + k*3*mo_num*point_num + i*3*mo_num + n*mo_num + j]) < 1.e-9);
}
}
}
}
}
free(forces_mo_g);
free(finite_difference_force_mo_g);
printf("OK\n");
#+end_src
* Force of MO laplacian
** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_forces_mo_l(qmckl_context context,
double* const forces_mo_l,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_forces_mo_l(qmckl_context context,
double* const forces_mo_l,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_forces_mo_l(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->point.num * 3 * ctx->mo_basis.mo_num * ctx->nucleus.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_forces_mo_l",
"input array too small");
}
memcpy(forces_mo_l, ctx->forces.forces_mo_l, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(qmckl_exit_code) function qmckl_get_forces_mo_l (context, &
forces_mo_l, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real(c_double), intent(out) :: forces_mo_l(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_forces_mo_l
end interface
#+end_src
** Provide
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_forces_mo_l(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :export none
qmckl_exit_code qmckl_provide_forces_mo_l(qmckl_context context)
{
qmckl_exit_code rc = QMCKL_SUCCESS;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT,
"qmckl_provide_forces_mo_l",
NULL);
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->mo_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_provide_mo_basis_mo_vgl",
NULL);
}
/* Compute if necessary */
if (ctx->point.date > ctx->forces.forces_mo_l_date) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 3 * ctx->mo_basis.mo_num * ctx->point.num * ctx->nucleus.num * sizeof(double);
if (ctx->forces.forces_mo_l != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->forces.forces_mo_l, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->forces.forces_mo_l);
assert (rc == QMCKL_SUCCESS);
ctx->forces.forces_mo_l= NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_mo_l == NULL) {
double* forces_mo_l = (double*) qmckl_malloc(context, mem_info);
if (forces_mo_l == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_forces_mo_l",
NULL);
}
ctx->forces.forces_mo_l = forces_mo_l;
}
rc = qmckl_provide_ao_basis_ao_hessian(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_basis_ao_hessian",
NULL);
}
rc = qmckl_compute_forces_mo_l_doc(context,
ctx->ao_basis.ao_num,
ctx->mo_basis.mo_num,
ctx->point.num,
ctx->nucleus.num,
ctx->mo_basis.coefficient_t,
ctx->ao_basis.ao_hessian,
ctx->forces.forces_mo_l);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->forces.forces_mo_l_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
** Compute
:PROPERTIES:
:Name: qmckl_compute_forces_mo_l
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_compute_forces_mo_l_args
| Variable | Type | In/Out | Description |
|-----------------+--------------------------------+--------+-------------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~ao_num~ | ~int64_t~ | in | Number of AOs |
| ~mo_num~ | ~int64_t~ | in | Number of MOs |
| ~point_num~ | ~int64_t~ | in | Number of points |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~coefficient_t~ | ~double[mo_num][ao_num]~ | in | Transpose of the AO to MO transformation matrix |
| ~ao_hessian~ | ~double[nucl_num][3][point_num][4][ao_num]~ | in | Value, gradients and Laplacian of the AOs |
| ~forces_mo_l~ | ~double[nucl_num][3][point_num][mo_num]~ | out | Value, gradients and Laplacian of the MOs |
|-----------------+--------------------------------+--------+-------------------------------------------------|
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_forces_mo_l_doc(context, &
ao_num, mo_num, point_num, nucl_num, &
coefficient_t, ao_hessian, forces_mo_l) &
bind(C) result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in), value :: context
integer (c_int64_t) , intent(in), value :: nucl_num, ao_num, mo_num, point_num
real (c_double ) , intent(in) :: coefficient_t(mo_num,ao_num)
real (c_double ) , intent(in) :: ao_hessian(ao_num,4,point_num,3,nucl_num)
real (c_double ) , intent(out) :: forces_mo_l(mo_num,point_num,3,nucl_num)
integer*8 :: i,j,k, m, n,a
info = QMCKL_SUCCESS
do j=1,point_num
forces_mo_l(:,j,:,:) = 0.d0
do k=1,ao_num
do i=1,mo_num
do a=1, nucl_num
do n = 1, 3
forces_mo_l(i, j, n, a) = forces_mo_l(i, j, n, a) - coefficient_t(i,k) * ao_hessian(k, 4, j, n, a)
end do
end do
end do
end do
end do
end function qmckl_compute_forces_mo_l_doc
#+end_src
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_forces_mo_l_doc (
const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
const int64_t nucl_num,
const double* coefficient_t,
const double* ao_hessian,
double* const forces_mo_l );
#+end_src
** Test
#+begin_src c :tangle (eval c_test)
printf("Forces MO laplacian\n");
rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num);
assert(rc == QMCKL_SUCCESS);
double * forces_mo_l = (double*) malloc(nucl_num * point_num* 3 * mo_num *sizeof(double));
rc = qmckl_get_forces_mo_l(context, &forces_mo_l[0], 3*nucl_num*mo_num*point_num);
assert(rc == QMCKL_SUCCESS);
double * finite_difference_force_mo_l= (double*) malloc(nucl_num* 3 * point_num * mo_num * sizeof(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);
temp_coord = (double*) malloc(3 * nucl_num * sizeof(double));
if (temp_coord == NULL) {
free(nucleus_coord);
return QMCKL_ALLOCATION_FAILED;
}
mo_output = (double*) malloc(5 * point_num * mo_num * sizeof(double));
if (mo_output == NULL) {
free(temp_coord);
free(nucleus_coord);
return QMCKL_ALLOCATION_FAILED;
}
// 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_mo_basis_mo_vgl(context,&mo_output[0], 5*point_num*mo_num);
assert(rc == QMCKL_SUCCESS);
// Accumulate derivative using finite-difference coefficients
for (int i = 0; i < point_num; i++) {
for (int j = 0; j < mo_num; j++) {
if (m == -4) {
finite_difference_force_mo_l[k*mo_num*point_num*nucl_num + a*mo_num*point_num + i*mo_num + j] = 0.0;
}
finite_difference_force_mo_l[k*mo_num*point_num*nucl_num + a*mo_num*point_num + i*mo_num + j] += coef[m + 4] * mo_output[i*mo_num*5 + 4*mo_num + j]/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);
free(mo_output);
for (int j = 0; j < mo_num; j++){
for (int i = 0; i < point_num; i++){
for (int a = 0; a < nucl_num; a++) {
for (int k = 0; k < 3; k++){
//printf("k=%i a=%i i=%i j=%i\n", k, a, i, j);
//printf("%.10f\t", finite_difference_force_mo_l[k*mo_num*point_num*nucl_num + a*mo_num*point_num + i*mo_num + j]);
//printf("%.10f\n", forces_mo_l[a*3*mo_num*point_num + k*mo_num*point_num + i*mo_num + j]);
assert(fabs(finite_difference_force_mo_l[k*mo_num*point_num*nucl_num + a*mo_num*point_num + i*mo_num + j] - forces_mo_l[a*3*mo_num*point_num + k*mo_num*point_num + i*mo_num + j]) < 1.e-8);
}
}
}
}
free(forces_mo_l);
free(finite_difference_force_mo_l);
printf("OK\n");
#+end_src
* End of files :noexport:
#+begin_src c :tangle (eval h_private_type)
#endif
#+end_src
#+begin_src c :tangle (eval h_private_func)
#endif
#+end_src
** Test
#+begin_src c :tangle (eval c_test)
rc = qmckl_context_destroy(context);
assert (rc == QMCKL_SUCCESS);
return 0;
}
#+end_src
Reference in New Issue View Git Blame Copy Permalink