From ef51934c7c1b69ad689c2d9d89f2cb647814e7fb Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Wed, 26 May 2021 11:59:03 +0530 Subject: [PATCH] Added kappa (distance rescale factor) to electron and nucleus structs along with tests. #15 --- org/qmckl_context.org | 6 +- org/qmckl_electron.org | 178 ++++++++++++++++++++++++++++++++++++----- org/qmckl_nucleus.org | 111 +++++++++++++++++++------ 3 files changed, 247 insertions(+), 48 deletions(-) diff --git a/org/qmckl_context.org b/org/qmckl_context.org index 70f7610..6472365 100644 --- a/org/qmckl_context.org +++ b/org/qmckl_context.org @@ -117,7 +117,6 @@ typedef struct qmckl_context_struct { qmckl_ao_basis_struct ao_basis; /* To be implemented: - qmckl_nucleus_struct nucleus; qmckl_mo_struct mo; qmckl_determinant_struct det; ,*/ @@ -174,6 +173,7 @@ qmckl_context qmckl_context_check(const qmckl_context context) { To create a new context, ~qmckl_context_create()~ should be used. - Upon success, it returns a pointer to a new context with the ~qmckl_context~ type - It returns ~QMCKL_NULL_CONTEXT~ upon failure to allocate the internal data structure + - A new context always has all its members initialized with a NULL value # Header #+begin_src c :comments org :tangle (eval h_func) :exports none @@ -222,8 +222,8 @@ qmckl_context qmckl_context_create() { ctx->numprec.range = QMCKL_DEFAULT_RANGE; ctx->ao_basis.uninitialized = (1 << 10) - 1; - ctx->nucleus.uninitialized = (1 << 3) - 1; - ctx->electron.uninitialized = (1 << 2) - 1; + ctx->nucleus.uninitialized = (1 << 4) - 1; + ctx->electron.uninitialized = (1 << 3) - 1; /* Allocate qmckl_memory_struct */ { diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index 67474b0..aa6b600 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -63,19 +63,25 @@ 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 | - | ~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 | + | ~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 | + | ~kappa_ee~ | double | The distance scaling 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 distances | + | ~ee_distance_rescaled_date~ | uint64_t | Last modification date of the electron-electron distances | + | ~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 | ** Data structure @@ -85,13 +91,19 @@ typedef struct qmckl_electron_struct { int64_t up_num; int64_t down_num; int64_t walk_num; + double kappa_ee; + double kappa_en; int64_t coord_new_date; int64_t ee_distance_date; int64_t en_distance_date; + int64_t ee_distance_rescaled_date; + int64_t en_distance_rescaled_date; double* coord_new; double* coord_old; double* ee_distance; double* en_distance; + double* ee_distance_rescaled; + double* en_distance_rescaled; int32_t uninitialized; bool provided; } qmckl_electron_struct; @@ -165,7 +177,7 @@ qmckl_get_electron_num (const qmckl_context context, int64_t* const num) { qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1; + int32_t mask = 1 << 0; if ( (ctx->electron.uninitialized & mask) != 0) { return QMCKL_NOT_PROVIDED; @@ -193,7 +205,7 @@ qmckl_get_electron_up_num (const qmckl_context context, int64_t* const up_num) { qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1; + int32_t mask = 1 << 0; if ( (ctx->electron.uninitialized & mask) != 0) { return QMCKL_NOT_PROVIDED; @@ -221,7 +233,7 @@ qmckl_get_electron_down_num (const qmckl_context context, int64_t* const down_nu qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1; + int32_t mask = 1 << 0; if ( (ctx->electron.uninitialized & mask) != 0) { return QMCKL_NOT_PROVIDED; @@ -260,7 +272,7 @@ qmckl_get_electron_walk_num (const qmckl_context context, int64_t* const walk_nu qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 2; + int32_t mask = 1 << 1; if ( (ctx->electron.uninitialized & mask) != 0) { return QMCKL_NOT_PROVIDED; @@ -272,6 +284,71 @@ qmckl_get_electron_walk_num (const qmckl_context context, int64_t* const walk_nu } #+end_src +*** Scaling factors Kappa + + #+begin_src c :comments org :tangle (eval h_func) :exports none +qmckl_exit_code qmckl_get_kappa_ee (const qmckl_context context, double* const kappa_ee); +qmckl_exit_code qmckl_get_kappa_en (const qmckl_context context, double* const kappa_en); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_get_kappa_ee (const qmckl_context context, double* const kappa_ee) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + if (kappa_ee == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_kappa_ee", + "kappa_ee is a null pointer"); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int32_t mask = 1 << 2; + + if ( (ctx->electron.uninitialized & mask) != 0) { + return QMCKL_NOT_PROVIDED; + } + + // TODO: assert (ctx->electron.kappa_ee > (double) 0); + ,*kappa_ee = ctx->electron.kappa_ee; + return QMCKL_SUCCESS; +} + + +qmckl_exit_code +qmckl_get_kappa_en (const qmckl_context context, double* const kappa_en) { + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + if (kappa_en == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_kappa_en", + "kappa_en is a null pointer"); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int32_t mask = 1 << 2; + + if ( (ctx->electron.uninitialized & mask) != 0) { + return QMCKL_NOT_PROVIDED; + } + + // TODO: assert (ctx->electron.kappa_en > (double) 0); + ,*kappa_en = ctx->electron.kappa_en; + return QMCKL_SUCCESS; +} + #+end_src + *** Electron coordinates Returns the current electron coordinates. The pointer is assumed @@ -359,6 +436,7 @@ qmckl_get_electron_coord (const qmckl_context context, const char transp, double #+begin_src c :comments org :tangle (eval h_func) qmckl_exit_code qmckl_set_electron_num (qmckl_context context, const int64_t up_num, const int64_t down_num); +qmckl_exit_code qmckl_set_kappa (qmckl_context context, const double kappa_ee, const double kappa_en); qmckl_exit_code qmckl_set_electron_walk_num (qmckl_context context, const int64_t walk_num); qmckl_exit_code qmckl_set_electron_coord (qmckl_context context, const char transp, const double* coord); #+end_src @@ -437,7 +515,7 @@ qmckl_set_electron_num(qmckl_context context, "down_num <= 0"); } - int32_t mask = 1; + int32_t mask = 1 << 0; ctx->electron.up_num = up_num; ctx->electron.down_num = down_num; @@ -447,6 +525,7 @@ qmckl_set_electron_num(qmckl_context context, } #+end_src + The following function sets the number of walkers. #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code @@ -461,13 +540,46 @@ qmckl_set_electron_walk_num(qmckl_context context, const int64_t walk_num) { "walk_num <= 0"); } - int32_t mask = 2; + int32_t mask = 1 << 1; ctx->electron.walk_num = walk_num; <> } #+end_src + Next we set the rescale parameter for the rescaled distance metric. + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_set_kappa(qmckl_context context, + const double kappa_ee, + const double kappa_en) { + <> + + // TODO: Check for 0 values + //if (kappa_ee != 0) { + // return qmckl_failwith( context, + // QMCKL_INVALID_ARG_2, + // "qmckl_set_kappa", + // "kappa_ee == 0"); + //} + + //if (kappa_en <= 0) { + // return qmckl_failwith( context, + // QMCKL_INVALID_ARG_3, + // "qmckl_set_kappa", + // "kappa_en == 0"); + //} + + int32_t mask = 1 << 2; + + ctx->electron.kappa_ee = kappa_ee; + ctx->electron.kappa_en = kappa_en; + + <> +} + #+end_src + The following function sets the electron coordinates of all the walkers. When this is done, the pointers to the old and new sets of coordinates are swapped, and the new coordinates are @@ -568,11 +680,14 @@ int64_t walk_num = chbrclf_walk_num; int64_t elec_num = chbrclf_elec_num; int64_t elec_up_num = chbrclf_elec_up_num; int64_t elec_dn_num = chbrclf_elec_dn_num; +double kappa_ee = 1.0; // TODO Get kappa_ee from chbrclf +double kappa_en = 1.0; // TODO Get kappa_en from chbrclf double* elec_coord = &(chbrclf_elec_coord[0][0][0]); int64_t nucl_num = chbrclf_nucl_num; double* charge = chbrclf_charge; double* nucl_coord = &(chbrclf_nucl_coord[0][0]); +double nucl_kappa = 1.0; // TODO Change get kappa from chbrclf example /* --- */ @@ -607,6 +722,26 @@ rc = qmckl_get_electron_num (context, &n); assert(rc == QMCKL_SUCCESS); assert(n == elec_num); +double k_ee; +double k_en; +rc = qmckl_get_kappa_ee (context, &k_ee); +assert(rc == QMCKL_NOT_PROVIDED); + +rc = qmckl_get_kappa_en (context, &k_en); +assert(rc == QMCKL_NOT_PROVIDED); + +rc = qmckl_set_kappa (context, kappa_ee, kappa_en); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_electron_provided(context)); + +rc = qmckl_get_kappa_ee (context, &k_ee); +assert(rc == QMCKL_SUCCESS); +assert(k_ee == kappa_ee); + +rc = qmckl_get_kappa_en (context, &k_en); +assert(rc == QMCKL_SUCCESS); +assert(k_en == kappa_en); + int64_t w; rc = qmckl_get_electron_walk_num (context, &w); @@ -1104,6 +1239,9 @@ assert(qmckl_electron_provided(context)); rc = qmckl_set_nucleus_num (context, nucl_num); assert(rc == QMCKL_SUCCESS); +rc = qmckl_set_nucleus_kappa (context, nucl_kappa); +assert(rc == QMCKL_SUCCESS); + rc = qmckl_set_nucleus_charge (context, charge); assert (rc == QMCKL_SUCCESS); @@ -1152,7 +1290,7 @@ assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12); } #+end_src -**✸ Compute file names +*** Compute file names #+begin_src emacs-lisp ; The following is required to compute the file names diff --git a/org/qmckl_nucleus.org b/org/qmckl_nucleus.org index b09cf6f..418c368 100644 --- a/org/qmckl_nucleus.org +++ b/org/qmckl_nucleus.org @@ -62,15 +62,18 @@ 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 | - | ~repulsion~ | double | Nuclear repulsion energy | - | ~repulsion_date~ | int64_t | Date when the nuclear repulsion energy was computed | + | ~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 | + | ~kappa~ | double | The distance scaling factor | ** Data structure @@ -85,21 +88,23 @@ typedef struct qmckl_nucleus_struct { double* nn_distance; double* nn_distance_rescaled; double repulsion; + double kappa; int32_t uninitialized; bool provided; } qmckl_nucleus_struct; #+end_src The ~uninitialized~ integer contains one bit set to one for each - initialization function which has not bee called. It becomes equal + initialization function which has not been called. It becomes equal to zero after all initialization functions have been called. The struct is then initialized and ~provided == true~. ** Access functions - + #+begin_src c :comments org :tangle (eval h_func) :exports none qmckl_exit_code qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num); qmckl_exit_code qmckl_get_nucleus_charge (const qmckl_context context, double* const charge); +qmckl_exit_code qmckl_get_nucleus_kappa (const qmckl_context context, double* const kappa); qmckl_exit_code qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double* const coord); #+end_src @@ -181,6 +186,31 @@ qmckl_get_nucleus_charge (const qmckl_context context, double* const charge) { } +qmckl_exit_code +qmckl_get_nucleus_kappa (const qmckl_context context, double* const kappa) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int32_t mask = 1 << 2; + + if ( (ctx->nucleus.uninitialized & mask) != 0) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_get_nucleus_kappa", + "nucleus data is not provided"); + } + + (*kappa) = ctx->nucleus.kappa; + + return QMCKL_SUCCESS; +} + + qmckl_exit_code qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double* const coord) { @@ -205,7 +235,7 @@ qmckl_get_nucleus_coord (const qmckl_context context, const char transp, double* qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - int32_t mask = 1 << 2; + int32_t mask = 1 << 3; if ( (ctx->nucleus.uninitialized & mask) != 0) { return qmckl_failwith( context, @@ -265,6 +295,7 @@ bool qmckl_nucleus_provided(const qmckl_context context) { #+begin_src c :comments org :tangle (eval h_func) qmckl_exit_code qmckl_set_nucleus_num (qmckl_context context, const int64_t num); qmckl_exit_code qmckl_set_nucleus_charge (qmckl_context context, const double* charge); +qmckl_exit_code qmckl_set_nucleus_kappa (qmckl_context context, const double kappa); qmckl_exit_code qmckl_set_nucleus_coord (qmckl_context context, const char transp, const double* coord); #+end_src @@ -349,6 +380,32 @@ qmckl_set_nucleus_charge(qmckl_context context, const double* charge) { ctx->nucleus.charge= memcpy(ctx->nucleus.charge, charge, num*sizeof(double)); assert (ctx->nucleus.charge != NULL); + <> +} + #+end_src + + The following function sets the rescale parameter for the nuclear distances. + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_set_nucleus_kappa(qmckl_context context, const double kappa) { + <> + + //TODO: Check for small values of kappa + //if (kappa == 0) { + // return qmckl_failwith( context, + // QMCKL_INVALID_ARG_2, + // "qmckl_set_nucleus_kappa", + // "kappa cannot be 0"); + //} + + int32_t mask = 1 << 2; + + int64_t num; + qmckl_exit_code rc; + + ctx->nucleus.kappa = kappa; + <> } #+end_src @@ -364,7 +421,7 @@ qmckl_set_nucleus_coord(qmckl_context context, const char transp, const double* int64_t nucl_num = (int64_t) 0; qmckl_exit_code rc; - int32_t mask = 1 << 2; + int32_t mask = 1 << 3; rc = qmckl_get_nucleus_num(context, &nucl_num); if (rc != QMCKL_SUCCESS) return rc; @@ -405,6 +462,7 @@ qmckl_set_nucleus_coord(qmckl_context context, const char transp, const double* const int64_t nucl_num = chbrclf_nucl_num; const double* nucl_charge = chbrclf_charge; const double* nucl_coord = &(chbrclf_nucl_coord[0][0]); +const double nucl_kappa = 1.0; // TODO Change get kappa from chbrclf example /* --- */ @@ -425,6 +483,19 @@ rc = qmckl_get_nucleus_num (context, &n); assert(rc == QMCKL_SUCCESS); assert(n == nucl_num); +double k; +rc = qmckl_get_nucleus_kappa (context, &k); +assert(rc == QMCKL_NOT_PROVIDED); + + +rc = qmckl_set_nucleus_kappa (context, nucl_kappa); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_nucleus_provided(context)); + +rc = qmckl_get_nucleus_kappa (context, &k); +assert(rc == QMCKL_SUCCESS); +assert(k == nucl_kappa); + double nucl_coord2[3*nucl_num]; rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2); @@ -476,7 +547,7 @@ assert(qmckl_nucleus_provided(context)); current date is stored. ** Nucleus-nucleus distances - + *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes @@ -816,7 +887,7 @@ qmckl_exit_code qmckl_compute_nn_distance_rescaled ( \] *** Get - + #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* energy); #+end_src @@ -960,16 +1031,6 @@ qmckl_exit_code qmckl_compute_nucleus_repulsion ( *** Test - #+begin_src c :tangle (eval c_test) -/* Reference input data */ - -assert(qmckl_nucleus_provided(context)); - -double rep; -rc = qmckl_get_nucleus_repulsion(context, &rep); -assert(rep - 318.2309879436158 < 1.e-10); - - #+end_src * End of files :noexport: