TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-06-13 16:55:35 +02:00
Code Issues Projects Releases Wiki Activity

Improved en_distance

Browse Source
This commit is contained in:
Anthony Scemama 2023-03-14 19:13:49 +01:00
parent 21336e0178
commit c0131d5da4
4 changed files with 317 additions and 122 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
org/qmckl_blas.org
View File

@ -248,7 +248,7 @@ qmckl_matrix_free( qmckl_context context,
| Variable | Type | Description |
|----------+-----------------------------------+-----------------------------|
| ~order~ | ~int64_t~ | Order of the tensor |
| ~order~ | ~int32_t~ | Order of the tensor |
| ~size~ | ~int64_t[QMCKL_TENSOR_ORDER_MAX]~ | Dimension of each component |
| ~data~ | ~double*~ | Elements |
@ -260,7 +260,8 @@ qmckl_matrix_free( qmckl_context context,
typedef struct qmckl_tensor {
double* restrict data;
int64_t order;
int32_t order;
int32_t __space__;
int64_t size[QMCKL_TENSOR_ORDER_MAX];
} qmckl_tensor;
#+end_src
@ -269,7 +270,7 @@ typedef struct qmckl_tensor {
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_tensor
qmckl_tensor_alloc( qmckl_context context,
const int64_t order,
const int32_t order,
const int64_t* size);
#+end_src
@ -279,7 +280,7 @@ qmckl_tensor_alloc( qmckl_context context,
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_tensor
qmckl_tensor_alloc( qmckl_context context,
const int64_t order,
const int32_t order,
const int64_t* size)
{
/* Should always be true by contruction */
@ -288,10 +289,11 @@ qmckl_tensor_alloc( qmckl_context context,
assert (size != NULL);
qmckl_tensor result;
memset(&result, 0, sizeof(qmckl_tensor));
result.order = order;
int64_t prod_size = (int64_t) 1;
for (int64_t i=0 ; i<order ; ++i) {
for (int32_t i=0 ; i<order ; ++i) {
assert (size[i] > (int64_t) 0);
result.size[i] = size[i];
prod_size *= size[i];
@ -383,7 +385,7 @@ qmckl_matrix_of_vector(const qmckl_vector vector,
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_tensor
qmckl_tensor_of_vector(const qmckl_vector vector,
const int64_t order,
const int32_t order,
const int64_t* size);
#+end_src
@ -392,13 +394,13 @@ qmckl_tensor_of_vector(const qmckl_vector vector,
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_tensor
qmckl_tensor_of_vector(const qmckl_vector vector,
const int64_t order,
const int32_t order,
const int64_t* size)
{
qmckl_tensor result;
int64_t prod_size = 1;
for (int64_t i=0 ; i<order ; ++i) {
for (int32_t i=0 ; i<order ; ++i) {
result.size[i] = size[i];
prod_size *= size[i];
}
@ -437,7 +439,7 @@ qmckl_vector_of_matrix(const qmckl_matrix matrix)
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_tensor
qmckl_tensor_of_matrix(const qmckl_matrix matrix,
const int64_t order,
const int32_t order,
const int64_t* size);
#+end_src
@ -446,13 +448,13 @@ qmckl_tensor_of_matrix(const qmckl_matrix matrix,
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_tensor
qmckl_tensor_of_matrix(const qmckl_matrix matrix,
const int64_t order,
const int32_t order,
const int64_t* size)
{
qmckl_tensor result;
int64_t prod_size = 1;
for (int64_t i=0 ; i<order ; ++i) {
for (int32_t i=0 ; i<order ; ++i) {
result.size[i] = size[i];
prod_size *= size[i];
}
@ -478,7 +480,7 @@ qmckl_vector
qmckl_vector_of_tensor(const qmckl_tensor tensor)
{
int64_t prod_size = (int64_t) tensor.size[0];
for (int64_t i=1 ; i<tensor.order ; i++) {
for (int32_t i=1 ; i<tensor.order ; i++) {
prod_size *= tensor.size[i];
}
@ -510,7 +512,7 @@ qmckl_matrix_of_tensor(const qmckl_tensor tensor,
{
/* Always true */
int64_t prod_size = (int64_t) 1;
for (int64_t i=0 ; i<tensor.order ; i++) {
for (int32_t i=0 ; i<tensor.order ; i++) {
prod_size *= tensor.size[i];
}
assert (prod_size == size1 * size2);
@ -602,7 +604,7 @@ qmckl_tensor
qmckl_tensor_set(qmckl_tensor tensor, double value)
{
qmckl_vector vector = qmckl_vector_of_tensor(tensor);
for (int64_t i=0 ; i<vector.size ; ++i) {
for (int32_t i=0 ; i<vector.size ; ++i) {
qmckl_vec(vector, i) = value;
}
return qmckl_tensor_of_vector(vector, tensor.order, tensor.size);

111
org/qmckl_electron.org
View File

@ -14,6 +14,8 @@ the components is =[ (x,y,z), (x,y,z), ... ]=
Using the ='T'= flage, it is =[ (x,x,x,...), (y,y,y,...), (z,z,z,...) ]=
# TODO: replace the qmckl_matrix by qmckl_point data structures.
# TODO: in provide funcions, replace the way to check that dimensions
# have changed
* Headers :noexport:
#+begin_src elisp :noexport :results none
@ -91,16 +93,16 @@ int main() {
Computed data:
| Variable | Type | Description |
|-------------------------------------+----------------------------------------+----------------------------------------------------------------------|
| ~ee_distance~ | ~double[walker.num][num][num]~ | Electron-electron distances |
| ~ee_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~en_distance~ | ~double[walker.num][num][nucl_num]~ | Electron-nucleus distances |
| ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~ee_potential~ | ~double[walker.num]~ | Electron-electron potential energy |
| ~ee_potential_date~ | ~uint64_t~ | Last modification date of the electron-electron potential |
| ~en_potential~ | ~double[walker.num]~ | Electron-nucleus potential energy |
| ~en_potential_date~ | ~int64_t~ | Date when the electron-nucleus potential energy was computed |
| Variable | Type | Description |
|---------------------+--------------------------------+--------------------------------------------------------------|
| ~ee_distance~ | ~double[walker.num][num][num]~ | Electron-electron distances |
| ~ee_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~en_distance~ | ~double[num][nucl_num]~ | Electron-nucleus distances |
| ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~ee_potential~ | ~double[walker.num]~ | Electron-electron potential energy |
| ~ee_potential_date~ | ~uint64_t~ | Last modification date of the electron-electron potential |
| ~en_potential~ | ~double[walker.num]~ | Electron-nucleus potential energy |
| ~en_potential_date~ | ~int64_t~ | Date when the electron-nucleus potential energy was computed |
** Data structure
@ -715,7 +717,7 @@ qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context)
/* Compute if necessary */
if (ctx->electron.walker.point.date > ctx->electron.ee_distance_date) {
if (ctx->point.date > ctx->electron.ee_distance_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->electron.ee_distance);
@ -965,7 +967,7 @@ qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context)
if (rc != QMCKL_SUCCESS) return rc;
/* Compute if necessary */
if (ctx->electron.walker.point.date > ctx->electron.ee_potential_date) {
if (ctx->point.date > ctx->electron.ee_potential_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->electron.ee_potential);
@ -1144,7 +1146,7 @@ qmckl_exit_code qmckl_get_electron_en_distance(qmckl_context context, double* di
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num;
size_t sze = ctx->point.num * ctx->nucleus.num;
memcpy(distance, ctx->electron.en_distance, sze * sizeof(double));
return QMCKL_SUCCESS;
@ -1176,19 +1178,28 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)
}
/* Compute if necessary */
if (ctx->electron.walker.point.date > ctx->electron.en_distance_date) {
if (ctx->point.date > ctx->electron.en_distance_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->electron.en_distance);
ctx->electron.en_distance = NULL;
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->point.num * ctx->nucleus.num * sizeof(double);
if (ctx->electron.en_distance != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
qmckl_exit_code rc = qmckl_get_malloc_info(context, ctx->electron.en_distance, &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->electron.en_distance);
assert (rc == QMCKL_SUCCESS);
ctx->electron.en_distance = NULL;
}
}
/* Allocate array */
if (ctx->electron.en_distance == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.num * ctx->nucleus.num *
ctx->electron.walker.num * sizeof(double);
double* en_distance = (double*) qmckl_malloc(context, mem_info);
if (en_distance == NULL) {
@ -1202,10 +1213,9 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)
qmckl_exit_code rc =
qmckl_compute_en_distance(context,
ctx->electron.num,
ctx->point.num,
ctx->nucleus.num,
ctx->electron.walker.num,
ctx->electron.walker.point.coord.data,
ctx->point.coord.data,
ctx->nucleus.coord.data,
ctx->electron.en_distance);
if (rc != QMCKL_SUCCESS) {
@ -1227,28 +1237,26 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)
:END:
#+NAME: qmckl_en_distance_args
| Variable | Type | In/Out | Description |
|---------------+----------------------------------------+--------+----------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~elec_num~ | ~int64_t~ | in | Number of electrons |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~walk_num~ | ~int64_t~ | in | Number of walkers |
| ~elec_coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates |
| ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates |
| ~en_distance~ | ~double[walk_num][elec_num][nucl_num]~ | out | Electron-nucleus distances |
| Variable | Type | In/Out | Description |
|---------------+-------------------------------+--------+----------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~point_num~ | ~int64_t~ | in | Number of points |
| ~nucl_num~ | ~int64_t~ | in | Number of nuclei |
| ~elec_coord~ | ~double[3][point_num]~ | in | Electron coordinates |
| ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates |
| ~en_distance~ | ~double[point_num][nucl_num]~ | out | Electron-nucleus distances |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_num, elec_coord, nucl_coord, en_distance) &
integer function qmckl_compute_en_distance_f(context, point_num, nucl_num, elec_coord, nucl_coord, en_distance) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: elec_num
integer*8 , intent(in) :: point_num
integer*8 , intent(in) :: nucl_num
integer*8 , intent(in) :: walk_num
double precision , intent(in) :: elec_coord(elec_num,walk_num,3)
double precision , intent(in) :: elec_coord(point_num,3)
double precision , intent(in) :: nucl_coord(nucl_num,3)
double precision , intent(out) :: en_distance(nucl_num,elec_num,walk_num)
double precision , intent(out) :: en_distance(nucl_num,point_num)
integer*8 :: k
@ -1259,7 +1267,7 @@ integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_n
return
endif
if (elec_num <= 0) then
if (point_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif
@ -1269,14 +1277,9 @@ integer function qmckl_compute_en_distance_f(context, elec_num, nucl_num, walk_n
return
endif
if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif
info = qmckl_distance(context, 'T', 'T', nucl_num, elec_num * walk_num, &
info = qmckl_distance(context, 'T', 'T', nucl_num, point_num, &
nucl_coord, nucl_num, &
elec_coord, elec_num * walk_num, &
elec_coord, point_num, &
en_distance, nucl_num)
end function qmckl_compute_en_distance_f
@ -1285,9 +1288,8 @@ end function qmckl_compute_en_distance_f
#+begin_src c :tangle (eval h_private_func) :comments org :exports none
qmckl_exit_code qmckl_compute_en_distance (
const qmckl_context context,
const int64_t elec_num,
const int64_t point_num,
const int64_t nucl_num,
const int64_t walk_num,
const double* elec_coord,
const double* nucl_coord,
double* const en_distance );
@ -1298,23 +1300,22 @@ qmckl_exit_code qmckl_compute_en_distance (
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_en_distance &
(context, elec_num, nucl_num, walk_num, elec_coord, nucl_coord, en_distance) &
(context, point_num, nucl_num, elec_coord, nucl_coord, en_distance) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: elec_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) , value :: walk_num
real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3)
real (c_double ) , intent(in) :: nucl_coord(elec_num,3)
real (c_double ) , intent(out) :: en_distance(nucl_num,elec_num,walk_num)
real (c_double ) , intent(in) :: elec_coord(point_num,3)
real (c_double ) , intent(in) :: nucl_coord(nucl_num,3)
real (c_double ) , intent(out) :: en_distance(nucl_num,point_num)
integer(c_int32_t), external :: qmckl_compute_en_distance_f
info = qmckl_compute_en_distance_f &
(context, elec_num, nucl_num, walk_num, elec_coord, nucl_coord, en_distance)
(context, point_num, nucl_num, elec_coord, nucl_coord, en_distance)
end function qmckl_compute_en_distance
#+end_src
@ -1451,7 +1452,7 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context)
if (rc != QMCKL_SUCCESS) return rc;
/* Compute if necessary */
if (ctx->electron.walker.point.date > ctx->electron.en_potential_date) {
if (ctx->point.date > ctx->electron.en_potential_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->electron.en_potential);

296
org/qmckl_mo.org
View File

@ -157,7 +157,7 @@ typedef struct qmckl_mo_basis_struct {
double * restrict mo_vgl;
double * restrict mo_value;
double * restrict r_cusp_param;
qmckl_tensor r_cusp_param;
uint64_t mo_vgl_date;
uint64_t mo_value_date;
@ -188,6 +188,8 @@ qmckl_exit_code qmckl_init_mo_basis(qmckl_context context) {
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
ctx->mo_basis.r_cusp = NULL;
ctx->mo_basis.uninitialized = (1 << 2) - 1;
return QMCKL_SUCCESS;
@ -252,7 +254,9 @@ qmckl_exit_code qmckl_set_mo_basis_mo_num(qmckl_context context, const int64_t m
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_set_mo_basis_coefficient(qmckl_context context,
const double* coefficient,
const int64_t size_max)
@ -292,55 +296,6 @@ qmckl_exit_code qmckl_set_mo_basis_coefficient(qmckl_context context,
<<post>>
}
qmckl_exit_code
qmckl_set_mo_basis_r_cusp(qmckl_context context,
const double* r_cusp,
const int64_t size_max)
{
int32_t mask = 0;
<<pre>>
if (r_cusp == NULL) {
return qmckl_failwith( context, QMCKL_INVALID_ARG_2,
"qmckl_set_mo_basis_r_cusp",
"r_cusp: Null pointer");
}
if (size_max < ctx->nucleus.num) {
return qmckl_failwith( context, QMCKL_INVALID_ARG_3,
"qmckl_set_mo_basis_r_cusp",
"Array too small");
}
if (ctx->mo_basis.r_cusp != NULL) {
qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.r_cusp);
if (rc != QMCKL_SUCCESS) {
return rc;
}
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->nucleus.num * sizeof(double);
double* new_array = (double*) qmckl_malloc(context, mem_info);
if (new_array == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_mo_basis_r_cusp",
NULL);
}
memcpy(new_array, r_cusp, mem_info.size);
ctx->mo_basis.r_cusp = new_array;
return QMCKL_SUCCESS;
}
#+end_src
When the basis set is completely entered, other data structures are
@ -406,10 +361,246 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
ctx->mo_basis.mo_value = NULL;
ctx->mo_basis.mo_value_date = 0;
}
return qmckl_context_touch(context);
}
#+end_src
** Cusp adjsutment functions
To activate the cusp adjustment, the user must enter the radius of
the fitting function for each atom.
This function requires the computation of the value and gradients
of the $s$ AOs at the distance equal to the radius, and the values
of the non-$s$ AOs at the center.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_set_mo_basis_r_cusp(qmckl_context context,
const double* r_cusp,
const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
if (r_cusp == NULL) {
return qmckl_failwith( context, QMCKL_INVALID_ARG_2,
"qmckl_set_mo_basis_r_cusp",
"r_cusp: Null pointer");
}
if (size_max < ctx->nucleus.num) {
return qmckl_failwith( context, QMCKL_INVALID_ARG_3,
"qmckl_set_mo_basis_r_cusp",
"Array too small");
}
// Nullify r_cusp
if (ctx->mo_basis.r_cusp != NULL) {
rc = qmckl_free(context, ctx->mo_basis.r_cusp);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->mo_basis.r_cusp = NULL;
}
// Save old points
int64_t old_point_num = ctx->point.num;
double* old_coord = NULL;
if (old_point_num > 0) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = old_point_num * 3 * sizeof(double);
old_coord = (double*) qmckl_malloc(context, mem_info);
rc = qmckl_get_point(context, 'T', old_coord, (old_point_num * 3));
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_FAILURE,
"qmckl_set_mo_basis_r_cusp",
"Unable to get old coordinates");
}
}
double* coord = (double*) malloc(ctx->nucleus.num * 3 * sizeof(double));
// Evaluate MO vgl at r_cusp
qmckl_tensor mo_vgl_at_r_cusp;
{
int64_t sze[3] = { ctx->mo_basis.mo_num, 5, ctx->nucleus.num };
mo_vgl_at_r_cusp = qmckl_tensor_alloc(context, 3, &(sze[0]));
rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, ctx->nucleus.num * 3);
if (rc != QMCKL_SUCCESS) return rc;
int64_t i=0;
for (int64_t inucl=0 ; inucl<ctx->nucleus.num ; ++inucl) {
coord[i] += r_cusp[inucl];
i+=3;
}
rc = qmckl_set_point(context, 'T', ctx->nucleus.num, coord, (ctx->nucleus.num * 3));
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_FAILURE,
"qmckl_set_mo_basis_r_cusp",
"Unable to set coordinates at r_cusp");
}
rc = qmckl_get_mo_basis_mo_vgl(context,
&(qmckl_ten3(mo_vgl_at_r_cusp,0,0,0)),
ctx->mo_basis.mo_num * 5 * ctx->nucleus.num);
if (rc != QMCKL_SUCCESS) return rc;
}
// Set r_cusp
{
assert (ctx->mo_basis.r_cusp == NULL);
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->nucleus.num * sizeof(double);
ctx->mo_basis.r_cusp = (double*) qmckl_malloc(context, mem_info);
if (ctx->mo_basis.r_cusp == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_set_mo_basis_r_cusp",
NULL);
}
memcpy(ctx->mo_basis.r_cusp, r_cusp, mem_info.size);
}
// Allocate cusp parameters and set them to zero
{
if (ctx->mo_basis.r_cusp_param.size[0] != 0) {
rc = qmckl_tensor_free(context, &(ctx->mo_basis.r_cusp_param));
if (rc != QMCKL_SUCCESS) return rc;
}
int64_t sze[3] = { 4, ctx->mo_basis.mo_num, ctx->nucleus.num };
ctx->mo_basis.r_cusp_param = qmckl_tensor_alloc(context, 3, &(sze[0]));
ctx->mo_basis.r_cusp_param = qmckl_tensor_set(ctx->mo_basis.r_cusp_param, 0.);
}
// Evaluate MO value at nucleus without s components
qmckl_matrix mo_value_at_nucl_no_s;
{
mo_value_at_nucl_no_s = qmckl_matrix_alloc(context, ctx->mo_basis.mo_num, ctx->nucleus.num);
rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, ctx->nucleus.num * 3);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_set_point(context, 'T', ctx->nucleus.num, coord, (ctx->nucleus.num * 3));
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_FAILURE,
"qmckl_set_mo_basis_r_cusp",
"Unable to set coordinates at the nuclei");
}
rc = qmckl_get_mo_basis_mo_value(context,
&(qmckl_mat(mo_value_at_nucl_no_s,0,0)),
ctx->mo_basis.mo_num * ctx->nucleus.num);
if (rc != QMCKL_SUCCESS) return rc;
}
// Evaluate MO vgl at r_cusp without s components
qmckl_tensor mo_vgl_at_r_cusp_no_s;
{
int64_t sze[3] = { ctx->mo_basis.mo_num, 5, ctx->nucleus.num };
mo_vgl_at_r_cusp_no_s = qmckl_tensor_alloc(context, 3, &(sze[0]));
rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, ctx->nucleus.num * 3);
if (rc != QMCKL_SUCCESS) return rc;
int64_t i=0;
for (int64_t inucl=0 ; inucl<ctx->nucleus.num ; ++inucl) {
coord[i] += r_cusp[inucl];
i+=3;
}
rc = qmckl_set_point(context, 'T', ctx->nucleus.num, coord, (ctx->nucleus.num * 3));
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_FAILURE,
"qmckl_set_mo_basis_r_cusp",
"Unable to set coordinates at r_cusp");
}
rc = qmckl_get_mo_basis_mo_vgl(context,
&(qmckl_ten3(mo_vgl_at_r_cusp_no_s,0,0,0)),
ctx->mo_basis.mo_num * 5 * ctx->nucleus.num);
if (rc != QMCKL_SUCCESS) return rc;
}
// Compute parameters
{
for (int64_t inucl=0 ; inucl < ctx->nucleus.num ; ++inucl) {
const double Z = qmckl_vec(ctx->nucleus.charge,inucl);
if (Z < 0.1) continue; // Avoid dummy atoms
const double R = r_cusp[inucl];
for (int64_t i=0 ; i<ctx->mo_basis.mo_num ; ++i) {
const double phi = qmckl_ten3(mo_vgl_at_r_cusp,i,0,inucl) - qmckl_ten3(mo_vgl_at_r_cusp_no_s,i,0,inucl);
const double grad_phi = qmckl_ten3(mo_vgl_at_r_cusp,i,1,inucl) - qmckl_ten3(mo_vgl_at_r_cusp_no_s,i,1,inucl);
const double lap_phi = qmckl_ten3(mo_vgl_at_r_cusp,i,4,inucl) - qmckl_ten3(mo_vgl_at_r_cusp_no_s,i,4,inucl);
const double eta = qmckl_mat(mo_value_at_nucl_no_s,i,inucl);
qmckl_ten3(ctx->mo_basis.r_cusp_param,0,i,inucl) =
-(R*(2.*eta*Z-6.*grad_phi)+lap_phi*R*R+6.*phi)/(2.*R*Z-6.);
qmckl_ten3(ctx->mo_basis.r_cusp_param,1,i,inucl) =
(lap_phi*R*R*Z-6.*grad_phi*R*Z+6.*phi*Z+6.*eta*Z)/(2.*R*Z-6.);
qmckl_ten3(ctx->mo_basis.r_cusp_param,2,i,inucl) =
-(R*(-5.*grad_phi*Z-1.5*lap_phi)+lap_phi*R*R*Z+3.*phi*Z+3.*eta*Z+6.*grad_phi)/(R*R*Z-3.*R);
qmckl_ten3(ctx->mo_basis.r_cusp_param,3,i,inucl) =
(R*(-2.*grad_phi*Z-lap_phi)+0.5*lap_phi*R*R*Z+phi*Z+eta*Z+3.*grad_phi)/(R*R*R*Z-3.*R*R);
printf("%ld %ld %f %f %f %f\n",i, inucl,
qmckl_ten3(ctx->mo_basis.r_cusp_param,0,i,inucl),
qmckl_ten3(ctx->mo_basis.r_cusp_param,1,i,inucl),
qmckl_ten3(ctx->mo_basis.r_cusp_param,2,i,inucl),
qmckl_ten3(ctx->mo_basis.r_cusp_param,3,i,inucl));
}
}
}
free(coord);
qmckl_matrix_free(context, &mo_value_at_nucl_no_s);
qmckl_tensor_free(context, &mo_vgl_at_r_cusp_no_s);
qmckl_tensor_free(context, &mo_vgl_at_r_cusp);
// Restore old points
if (old_point_num > 0) {
rc = qmckl_set_point(context, 'T', old_point_num, old_coord, (old_point_num * 3));
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_FAILURE,
"qmckl_set_mo_basis_r_cusp",
"Unable to set old coordinates");
}
rc = qmckl_free(context, old_coord);
if (rc != QMCKL_SUCCESS) return rc;
old_coord = NULL;
}
return QMCKL_SUCCESS;
}
#+end_src
** Access functions
#+begin_src c :comments org :tangle (eval h_func) :exports none
@ -908,10 +1099,11 @@ qmckl_exit_code qmckl_provide_mo_basis_mo_value(qmckl_context context)
} else {
rc = qmckl_provide_en_distance(context);
if (rc != QMCKL_SUCCESS) {
return rc;
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron_en_distance",
NULL);
"qmckl_provide_mo_basis_mo_value",
"en_distance");
}
rc = qmckl_compute_mo_basis_mo_value_cusp(context,

2
org/qmckl_trexio.org
View File

@ -267,7 +267,7 @@ qmckl_trexio_read_nucleus_X(qmckl_context context, trexio_t* const file)
#+begin_src c :tangle (eval c)
assert ( qmckl_nucleus_provided(context) );
return QMCKL_SUCCESS;
}
#endif