1
0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-07-18 08:53:47 +02:00

Added electron and nucleus potential. #41

This commit is contained in:
v1j4y 2021-10-11 11:34:42 +02:00
parent 2959281a10
commit 88178c5efa
2 changed files with 443 additions and 39 deletions

View File

@ -63,29 +63,36 @@ int main() {
The following data stored in the context:
| ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data |
| ~num~ | ~int64_t~ | Total number of electrons |
| ~up_num~ | ~int64_t~ | Number of up-spin electrons |
| ~down_num~ | ~int64_t~ | Number of down-spin electrons |
| ~walk_num~ | ~int64_t~ | Number of walkers |
| ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor |
| ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor |
| ~provided~ | ~bool~ | If true, ~electron~ is valid |
| ~coord_new~ | ~double[walk_num][3][num]~ | New set of electron coordinates |
| ~coord_old~ | ~double[walk_num][3][num]~ | Old set of electron coordinates |
| ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates |
| ~ee_distance~ | ~double[walk_num][num][num]~ | Electron-electron distances |
| ~ee_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~en_distance~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances |
| ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~ee_distance_rescaled~ | ~double[walk_num][num][num]~ | Electron-electron rescaled distances |
| ~ee_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~ee_distance_rescaled_deriv_e~ | ~double[walk_num][4][num][num]~ | Electron-electron rescaled distances derivatives |
| ~ee_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives |
| ~en_distance_rescaled~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances |
| ~en_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~en_distance_rescaled_deriv_e~ | ~double[walk_num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives |
| ~en_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives |
|-------------------------------------+--------------------------------------+----------------------------------------------------------------------|
| ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data |
| ~num~ | ~int64_t~ | Total number of electrons |
| ~up_num~ | ~int64_t~ | Number of up-spin electrons |
| ~down_num~ | ~int64_t~ | Number of down-spin electrons |
| ~walk_num~ | ~int64_t~ | Number of walkers |
| ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor |
| ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor |
| ~provided~ | ~bool~ | If true, ~electron~ is valid |
| ~coord_new~ | ~double[walk_num][3][num]~ | New set of electron coordinates |
| ~coord_old~ | ~double[walk_num][3][num]~ | Old set of electron coordinates |
| ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates |
Computed data:
|-------------------------------------+--------------------------------------+----------------------------------------------------------------------|
| ~ee_distance~ | ~double[walk_num][num][num]~ | Electron-electron distances |
| ~ee_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~en_distance~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances |
| ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~ee_distance_rescaled~ | ~double[walk_num][num][num]~ | Electron-electron rescaled distances |
| ~ee_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~ee_distance_rescaled_deriv_e~ | ~double[walk_num][4][num][num]~ | Electron-electron rescaled distances derivatives |
| ~ee_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives |
| ~ee_pot~ | ~double[walk_num][4][num][num]~ | Electron-electron rescaled distances derivatives |
| ~ee_pot_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives |
| ~en_distance_rescaled~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances |
| ~en_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~en_distance_rescaled_deriv_e~ | ~double[walk_num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives |
| ~en_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives |
** Data structure
@ -100,6 +107,7 @@ typedef struct qmckl_electron_struct {
int64_t coord_new_date;
int64_t ee_distance_date;
int64_t en_distance_date;
int64_t ee_pot_date;
int64_t ee_distance_rescaled_date;
int64_t ee_distance_rescaled_deriv_e_date;
int64_t en_distance_rescaled_date;
@ -108,6 +116,7 @@ typedef struct qmckl_electron_struct {
double* coord_old;
double* ee_distance;
double* en_distance;
double* ee_pot;
double* ee_distance_rescaled;
double* ee_distance_rescaled_deriv_e;
double* en_distance_rescaled;
@ -1378,7 +1387,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance_rescaled_deriv_e(qmckl_context co
#+end_src
*** Provide :noexport:
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context);
#+end_src
@ -1565,6 +1574,194 @@ rc = qmckl_get_electron_ee_distance_rescaled_deriv_e(context, ee_distance_rescal
#+end_src
** Electron-electron potential
~ee_pot~ calculates the ~ee~ potential energy.
\[
\mathcal{V}_{ee} = \sum_{i=1}^{N_e}\sum_{j>i}^{N_e}\frac{1}{r_{ij}}
\]
where \(\mathcal{V}_{ee}\) is the ~ee~ potential and \[r_{ij}\] the ~ee~
distance.
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* const ee_pot);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* const ee_pot)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_ee_potential(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
size_t sze = ctx->electron.walk_num * sizeof(double);
memcpy(ee_pot, ctx->electron.ee_pot, sze);
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_ee_potential(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
if (!ctx->electron.provided) return QMCKL_NOT_PROVIDED;
/* Compute if necessary */
if (ctx->electron.coord_new_date > ctx->electron.ee_pot_date) {
/* Allocate array */
if (ctx->electron.ee_pot == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walk_num * sizeof(double);
double* ee_pot = (double*) qmckl_malloc(context, mem_info);
if (ee_pot == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_ee_potential",
NULL);
}
ctx->electron.ee_pot = ee_pot;
}
qmckl_exit_code rc =
qmckl_compute_ee_potential(context,
ctx->electron.num,
ctx->electron.walk_num,
ctx->electron.ee_distance,
ctx->electron.ee_pot);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->electron.ee_pot_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
*** Compute
:PROPERTIES:
:Name: qmckl_compute_ee_potential
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_ee_potential_args
| qmckl_context | context | in | Global state |
| int64_t | elec_num | in | Number of electrons |
| int64_t | walk_num | in | Number of walkers |
| double | ee_distance[walk_num][elec_num][elec_num] | in | Electron-electron rescaled distances |
| double | ee_pot[walk_num] | out | Electron-electron potential |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_ee_potential_f(context, elec_num, walk_num, &
ee_distance, ee_pot) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: elec_num
integer*8 , intent(in) :: walk_num
double precision , intent(in) :: ee_distance(elec_num,elec_num,walk_num)
double precision , intent(out) :: ee_pot(walk_num)
integer*8 :: nw, i, j
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif
if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif
ee_pot = 0.0d0
do nw=1,walk_num
do j=1,elec_num
do i=1,j-1
ee_pot(nw) = ee_pot(nw) + 1.0d0/(ee_distance(i,j,nw))
end do
end do
end do
end function qmckl_compute_ee_potential_f
#+end_src
#+CALL: generate_c_header(table=qmckl_ee_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_ee_potential (
const qmckl_context context,
const int64_t elec_num,
const int64_t walk_num,
const double* ee_distance,
double* const ee_pot );
#+end_src
#+CALL: generate_c_interface(table=qmckl_ee_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_ee_potential &
(context, elec_num, walk_num, ee_distance, ee_pot) &
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 :: walk_num
real (c_double ) , intent(in) :: ee_distance(elec_num,elec_num,walk_num)
real (c_double ) , intent(out) :: ee_pot(walk_num)
integer(c_int32_t), external :: qmckl_compute_ee_potential_f
info = qmckl_compute_ee_potential_f &
(context, elec_num, walk_num, ee_distance, ee_pot)
end function qmckl_compute_ee_potential
#+end_src
*** Test
** Electron-nucleus distances
*** Get

View File

@ -62,18 +62,26 @@ int main() {
The following data stored in the context:
| ~uninitialized~ | int32_t | Keeps bit set for uninitialized data |
| ~num~ | int64_t | Total number of nuclei |
| ~provided~ | bool | If true, ~nucleus~ is valid |
| ~charge~ | double[num] | Nuclear charges |
| ~coord~ | double[3][num] | Nuclear coordinates, in transposed format |
| ~nn_distance~ | double[num][num] | Nucleus-nucleus distances |
| ~nn_distance_date~ | int64_t | Date when Nucleus-nucleus distances were computed |
| ~nn_distance_rescaled~ | double[num][num] | Nucleus-nucleus rescaled distances |
| ~nn_distance_rescaled_date~ | int64_t | Date when Nucleus-nucleus rescaled distances were computed |
| ~repulsion~ | double | Nuclear repulsion energy |
| ~repulsion_date~ | int64_t | Date when the nuclear repulsion energy was computed |
| ~rescale_factor_kappa~ | double | The distance scaling factor |
|------------------------+----------------+-------------------------------------------|
| ~uninitialized~ | int32_t | Keeps bit set for uninitialized data |
| ~num~ | int64_t | Total number of nuclei |
| ~provided~ | bool | If true, ~nucleus~ is valid |
| ~charge~ | double[num] | Nuclear charges |
| ~coord~ | double[3][num] | Nuclear coordinates, in transposed format |
| ~coord_date~ | int64_t | Nuclear coordinates, date if modified |
| ~rescale_factor_kappa~ | double | The distance scaling factor |
Computed data:
|-----------------------------+------------------+--------------------------------------------------------------|
| ~nn_distance~ | double[num][num] | Nucleus-nucleus distances |
| ~nn_distance_date~ | int64_t | Date when Nucleus-nucleus distances were computed |
| ~nn_distance_rescaled~ | double[num][num] | Nucleus-nucleus rescaled distances |
| ~nn_distance_rescaled_date~ | int64_t | Date when Nucleus-nucleus rescaled distances were computed |
| ~repulsion~ | double | Nuclear repulsion energy |
| ~repulsion_date~ | int64_t | Date when the nuclear repulsion energy was computed |
| ~en_pot~ | double | Electron-nucleus potential energy |
| ~en_pot_date~ | int64_t | Date when the electron-nucleus potential energy was computed |
** Data structure
@ -83,6 +91,7 @@ typedef struct qmckl_nucleus_struct {
int64_t repulsion_date;
int64_t nn_distance_date;
int64_t nn_distance_rescaled_date;
int64_t coord_date;
double* coord;
double* charge;
double* nn_distance;
@ -91,6 +100,8 @@ typedef struct qmckl_nucleus_struct {
double rescale_factor_kappa;
int32_t uninitialized;
bool provided;
int64_t en_pot_date;
double* en_pot;
} qmckl_nucleus_struct;
#+end_src
@ -125,8 +136,6 @@ qmckl_exit_code qmckl_init_nucleus(qmckl_context context) {
}
#+end_src
** Access functions
#+begin_src c :comments org :tangle (eval h_func) :exports none
@ -785,7 +794,6 @@ assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
#+end_src
** Nucleus-nucleus rescaled distances
*** Get
@ -1132,6 +1140,205 @@ assert(rep - 318.2309879436158 < 1.e-10);
#+end_src
** Electron-nucleus potential
~en_potential~ stores the ~en~ potential energy
\[
\mathcal{V}_{en} = \sum_{i=1}^{N_e}\sum_{A=1}^{N_n}\frac{Z_A}{r_{iA}}
\]
where \(\mathcal{V}_{en}\) is the ~en~ potential, \[r_{iA}\] the ~en~
distance and \[Z_A\] is the nuclear charge.
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* const en_pot);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* const en_pot)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_en_potential(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
size_t sze = ctx->electron.walk_num * sizeof(double);
memcpy(en_pot, ctx->nucleus.en_pot, sze);
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_en_potential(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_en_potential(qmckl_context context)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
if (!ctx->electron.provided) return QMCKL_NOT_PROVIDED;
if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED;
/* Compute if necessary */
if (ctx->nucleus.coord_date > ctx->nucleus.en_pot_date) {
/* Allocate array */
if (ctx->nucleus.en_pot == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walk_num * sizeof(double);
double* en_pot = (double*) qmckl_malloc(context, mem_info);
if (en_pot == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_en_potential",
NULL);
}
ctx->nucleus.en_pot = en_pot;
}
qmckl_exit_code rc =
qmckl_compute_en_potential(context,
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.walk_num,
ctx->nucleus.charge,
ctx->electron.en_distance,
ctx->nucleus.en_pot);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->nucleus.en_pot_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
*** Compute
:PROPERTIES:
:Name: qmckl_compute_en_potential
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_en_potential_args
| qmckl_context | context | in | Global state |
| int64_t | elec_num | in | Number of electrons |
| int64_t | nucl_num | in | Number of nucleii |
| int64_t | walk_num | in | Number of walkers |
| double | charge[nucl_num] | in | charge of nucleus |
| double | en_distance[walk_num][nucl_num][elec_num] | in | Electron-electron rescaled distances |
| double | en_pot[walk_num] | out | Electron-electron potential |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_en_potential_f(context, elec_num, nucl_num, walk_num, &
charge, en_distance, en_pot) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: elec_num
integer*8 , intent(in) :: nucl_num
integer*8 , intent(in) :: walk_num
double precision , intent(in) :: charge(nucl_num)
double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num)
double precision , intent(out) :: en_pot(walk_num)
integer*8 :: nw, i, j
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif
if (walk_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif
en_pot = 0.0d0
do nw=1,walk_num
do j=1,nucl_num
do i=1,elec_num
en_pot(nw) = en_pot(nw) + charge(j)/(en_distance(i,j,nw))
end do
end do
end do
end function qmckl_compute_en_potential_f
#+end_src
#+CALL: generate_c_header(table=qmckl_en_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_en_potential (
const qmckl_context context,
const int64_t elec_num,
const int64_t nucl_num,
const int64_t walk_num,
const double* charge,
const double* en_distance,
double* const en_pot );
#+end_src
#+CALL: generate_c_interface(table=qmckl_en_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_en_potential &
(context, elec_num, nucl_num, walk_num, charge, en_distance, en_pot) &
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 :: nucl_num
integer (c_int64_t) , intent(in) , value :: walk_num
real (c_double ) , intent(in) :: charge(nucl_num)
real (c_double ) , intent(in) :: en_distance(elec_num,nucl_num,walk_num)
real (c_double ) , intent(out) :: en_pot(walk_num)
integer(c_int32_t), external :: qmckl_compute_en_potential_f
info = qmckl_compute_en_potential_f &
(context, elec_num, nucl_num, walk_num, charge, en_distance, en_pot)
end function qmckl_compute_en_potential
#+end_src
*** Test
* End of files :noexport:
#+begin_src c :tangle (eval h_private_type)