1
0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-11-03 20:54:09 +01:00

Moved rescale factors in Jastrow, removed kappa from names

This commit is contained in:
Anthony Scemama 2022-09-09 10:36:38 +02:00
parent 9f4731ff94
commit 368604633d
7 changed files with 2134 additions and 2321 deletions

View File

@ -292,6 +292,9 @@ qmckl_context qmckl_context_create() {
rc = qmckl_init_determinant(context); rc = qmckl_init_determinant(context);
assert (rc == QMCKL_SUCCESS); assert (rc == QMCKL_SUCCESS);
rc = qmckl_init_jastrow(context);
assert (rc == QMCKL_SUCCESS);
} }
/* Allocate qmckl_memory_struct */ /* Allocate qmckl_memory_struct */

View File

@ -1111,21 +1111,21 @@ const double* nucl_charge = chbrclf_charge;
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]); const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
rc = qmckl_set_electron_num (context, chbrclf_elec_up_num, chbrclf_elec_dn_num); rc = qmckl_set_electron_num (context, chbrclf_elec_up_num, chbrclf_elec_dn_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
rc = qmckl_set_electron_coord (context, 'N', chbrclf_walk_num, elec_coord, chbrclf_walk_num*chbrclf_elec_num*3);
qmckl_check(context, rc);
assert(qmckl_electron_provided(context)); assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', chbrclf_walk_num, elec_coord, chbrclf_walk_num*chbrclf_elec_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num); rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), chbrclf_nucl_num*3); rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), chbrclf_nucl_num*3);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
rc = qmckl_set_nucleus_charge(context, nucl_charge, chbrclf_nucl_num); rc = qmckl_set_nucleus_charge(context, nucl_charge, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(qmckl_nucleus_provided(context)); assert(qmckl_nucleus_provided(context));
@ -1145,57 +1145,57 @@ const char typ = 'G';
assert(!qmckl_ao_basis_provided(context)); assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_type (context, typ); rc = qmckl_set_ao_basis_type (context, typ);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(!qmckl_ao_basis_provided(context)); assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_num (context, chbrclf_shell_num); rc = qmckl_set_ao_basis_shell_num (context, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(!qmckl_ao_basis_provided(context)); assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_prim_num (context, chbrclf_prim_num); rc = qmckl_set_ao_basis_prim_num (context, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(!qmckl_ao_basis_provided(context)); assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, chbrclf_nucl_num); rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(!qmckl_ao_basis_provided(context)); assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, chbrclf_nucl_num); rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(!qmckl_ao_basis_provided(context)); assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, chbrclf_shell_num); rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(!qmckl_ao_basis_provided(context)); assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_factor (context, shell_factor, chbrclf_shell_num); rc = qmckl_set_ao_basis_shell_factor (context, shell_factor, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(!qmckl_ao_basis_provided(context)); assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, chbrclf_shell_num); rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(!qmckl_ao_basis_provided(context)); assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, chbrclf_shell_num); rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(!qmckl_ao_basis_provided(context)); assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_exponent (context, exponent, chbrclf_prim_num); rc = qmckl_set_ao_basis_exponent (context, exponent, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(!qmckl_ao_basis_provided(context)); assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_coefficient (context, coefficient, chbrclf_prim_num); rc = qmckl_set_ao_basis_coefficient (context, coefficient, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(!qmckl_ao_basis_provided(context)); assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_prim_factor (context, prim_factor, chbrclf_prim_num); rc = qmckl_set_ao_basis_prim_factor (context, prim_factor, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
rc = qmckl_set_ao_basis_ao_num(context, chbrclf_ao_num); rc = qmckl_set_ao_basis_ao_num(context, chbrclf_ao_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
rc = qmckl_set_ao_basis_ao_factor (context, ao_factor, chbrclf_ao_num); rc = qmckl_set_ao_basis_ao_factor (context, ao_factor, chbrclf_ao_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(qmckl_ao_basis_provided(context)); assert(qmckl_ao_basis_provided(context));
@ -1203,22 +1203,22 @@ assert(qmckl_ao_basis_provided(context));
double ao_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_ao_num]; double ao_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_ao_num];
rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]), (int64_t) 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_ao_num); rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]), (int64_t) 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_ao_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
/* Set up MO data */ /* Set up MO data */
rc = qmckl_set_mo_basis_mo_num(context, chbrclf_mo_num); rc = qmckl_set_mo_basis_mo_num(context, chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
const double * mo_coefficient = &(chbrclf_mo_coef[0]); const double * mo_coefficient = &(chbrclf_mo_coef[0]);
rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient); rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(qmckl_mo_basis_provided(context)); assert(qmckl_mo_basis_provided(context));
double mo_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_mo_num]; double mo_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_mo_num];
rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_mo_num); rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
/* Set up determinant data */ /* Set up determinant data */
@ -1238,19 +1238,19 @@ for(k = 0; k < det_num_beta; ++k)
mo_index_beta[k][i][j] = j + 1; mo_index_beta[k][i][j] = j + 1;
rc = qmckl_set_determinant_type (context, typ); rc = qmckl_set_determinant_type (context, typ);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
rc = qmckl_set_determinant_det_num_alpha (context, det_num_alpha); rc = qmckl_set_determinant_det_num_alpha (context, det_num_alpha);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
rc = qmckl_set_determinant_det_num_beta (context, det_num_beta); rc = qmckl_set_determinant_det_num_beta (context, det_num_beta);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
rc = qmckl_set_determinant_mo_index_alpha (context, &(mo_index_alpha[0][0][0])); rc = qmckl_set_determinant_mo_index_alpha (context, &(mo_index_alpha[0][0][0]));
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
rc = qmckl_set_determinant_mo_index_beta (context, &(mo_index_beta[0][0][0])); rc = qmckl_set_determinant_mo_index_beta (context, &(mo_index_beta[0][0][0]));
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
// Get slater-determinant // Get slater-determinant
@ -1258,10 +1258,10 @@ double det_vgl_alpha[det_num_alpha][chbrclf_walk_num][5][chbrclf_elec_up_num][ch
double det_vgl_beta[det_num_beta][chbrclf_walk_num][5][chbrclf_elec_dn_num][chbrclf_elec_dn_num]; double det_vgl_beta[det_num_beta][chbrclf_walk_num][5][chbrclf_elec_dn_num][chbrclf_elec_dn_num];
rc = qmckl_get_det_vgl_alpha(context, &(det_vgl_alpha[0][0][0][0][0])); rc = qmckl_get_det_vgl_alpha(context, &(det_vgl_alpha[0][0][0][0][0]));
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
rc = qmckl_get_det_vgl_beta(context, &(det_vgl_beta[0][0][0][0][0])); rc = qmckl_get_det_vgl_beta(context, &(det_vgl_beta[0][0][0][0][0]));
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
#+end_src #+end_src
@ -2018,10 +2018,10 @@ double det_inv_matrix_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num
double det_inv_matrix_beta[det_num_beta][chbrclf_walk_num][chbrclf_elec_dn_num][chbrclf_elec_dn_num]; double det_inv_matrix_beta[det_num_beta][chbrclf_walk_num][chbrclf_elec_dn_num][chbrclf_elec_dn_num];
rc = qmckl_get_det_inv_matrix_alpha(context, &(det_inv_matrix_alpha[0][0][0][0])); rc = qmckl_get_det_inv_matrix_alpha(context, &(det_inv_matrix_alpha[0][0][0][0]));
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
rc = qmckl_get_det_inv_matrix_beta(context, &(det_inv_matrix_beta[0][0][0][0])); rc = qmckl_get_det_inv_matrix_beta(context, &(det_inv_matrix_beta[0][0][0][0]));
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
#+end_src #+end_src
@ -2034,7 +2034,7 @@ assert (rc == QMCKL_SUCCESS);
*** Test *** Test
#+begin_src c :tangle (eval c_test) #+begin_src c :tangle (eval c_test)
rc = qmckl_context_destroy(context); rc = qmckl_context_destroy(context);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
return 0; return 0;
} }

File diff suppressed because it is too large Load Diff

View File

@ -608,6 +608,62 @@ qmckl_last_error(qmckl_context context, char* buffer) {
end interface end interface
#+end_src #+end_src
* Helper functions for debugging
The following function prints to ~stderr~ an error message is the return code is
not ~QMCKL_SUCCESS~.
# Header
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_check(qmckl_context context, qmckl_exit_code rc);
#+end_src
# Source
#+begin_src c :tangle (eval c) :exports none
#include <stdio.h>
qmckl_exit_code
qmckl_check(qmckl_context context, qmckl_exit_code rc)
{
char fname[QMCKL_MAX_FUN_LEN];
char message[QMCKL_MAX_MSG_LEN];
if (rc != QMCKL_SUCCESS) {
fprintf(stderr, "===========\nQMCKL ERROR\n%s\n", qmckl_string_of_error(rc));
qmckl_get_error(context, &rc, fname, message);
fprintf(stderr, "Function: %s\nMessage: %s\n===========\n", fname, message);
}
return rc;
}
#+end_src
It should be used as:
#+begin_src c
rc = qmckl_check(context,
qmckl_...(context, ...)
);
assert (rc == QMCKL_SUCCESS);
#+end_src
** Fortran inteface
#+begin_src f90 :tangle (eval fh_func) :exports none :noweb yes
interface
function qmckl_check (context, rc) bind(C, name='qmckl_check')
use, intrinsic :: iso_c_binding
import
implicit none
integer(qmckl_exit_code) :: qmckl_check
integer (c_int64_t) , intent(in), value :: context
integer(qmckl_exit_code), intent(in) :: rc
end function qmckl_check
end interface
#+end_src
* End of files :noexport: * End of files :noexport:
#+begin_src c :comments link :tangle (eval h_private_type) #+begin_src c :comments link :tangle (eval h_private_type)
@ -623,12 +679,13 @@ qmckl_last_error(qmckl_context context, char* buffer) {
qmckl_exit_code exit_code; qmckl_exit_code exit_code;
exit_code = 1; exit_code = 1;
assert (qmckl_set_error(context, exit_code, "qmckl_transpose", "Success") == QMCKL_SUCCESS); assert (qmckl_set_error(context, exit_code, "my_function", "Message") == QMCKL_SUCCESS);
assert (qmckl_get_error(context, &exit_code, function_name, message) == QMCKL_SUCCESS); assert (qmckl_get_error(context, &exit_code, function_name, message) == QMCKL_SUCCESS);
assert (exit_code == 1); assert (exit_code == 1);
assert (strcmp(function_name,"qmckl_transpose") == 0);
assert (strcmp(message,"Success") == 0); assert (strcmp(function_name,"my_function") == 0);
assert (strcmp(message,"Message") == 0);
exit_code = qmckl_context_destroy(context); exit_code = qmckl_context_destroy(context);
assert(exit_code == QMCKL_SUCCESS); assert(exit_code == QMCKL_SUCCESS);

File diff suppressed because it is too large Load Diff

View File

@ -74,15 +74,12 @@ int main() {
| ~charge~ | qmckl_vector | Nuclear charges | | ~charge~ | qmckl_vector | Nuclear charges |
| ~coord~ | qmckl_matrix | Nuclear coordinates, in transposed format | | ~coord~ | qmckl_matrix | Nuclear coordinates, in transposed format |
| ~coord_date~ | int64_t | Nuclear coordinates, date if modified | | ~coord_date~ | int64_t | Nuclear coordinates, date if modified |
| ~rescale_factor_kappa~ | double | The distance scaling factor |
Computed data: Computed data:
|-----------------------------+--------------+------------------------------------------------------------| |-----------------------------+--------------+------------------------------------------------------------|
| ~nn_distance~ | qmckl_matrix | Nucleus-nucleus distances | | ~nn_distance~ | qmckl_matrix | Nucleus-nucleus distances |
| ~nn_distance_date~ | int64_t | Date when Nucleus-nucleus distances were computed | | ~nn_distance_date~ | int64_t | Date when Nucleus-nucleus distances were computed |
| ~nn_distance_rescaled~ | qmckl_matrix | Nucleus-nucleus rescaled distances |
| ~nn_distance_rescaled_date~ | int64_t | Date when Nucleus-nucleus rescaled distances were computed |
| ~repulsion~ | double | Nuclear repulsion energy | | ~repulsion~ | double | Nuclear repulsion energy |
| ~repulsion_date~ | int64_t | Date when the nuclear repulsion energy was computed | | ~repulsion_date~ | int64_t | Date when the nuclear repulsion energy was computed |
@ -93,14 +90,11 @@ typedef struct qmckl_nucleus_struct {
int64_t num; int64_t num;
int64_t repulsion_date; int64_t repulsion_date;
int64_t nn_distance_date; int64_t nn_distance_date;
int64_t nn_distance_rescaled_date;
int64_t coord_date; int64_t coord_date;
qmckl_vector charge; qmckl_vector charge;
qmckl_matrix coord; qmckl_matrix coord;
qmckl_matrix nn_distance; qmckl_matrix nn_distance;
qmckl_matrix nn_distance_rescaled;
double repulsion; double repulsion;
double rescale_factor_kappa;
int32_t uninitialized; int32_t uninitialized;
bool provided; bool provided;
} qmckl_nucleus_struct; } qmckl_nucleus_struct;
@ -131,7 +125,6 @@ qmckl_exit_code qmckl_init_nucleus(qmckl_context context) {
ctx->nucleus.uninitialized = (1 << 3) - 1; ctx->nucleus.uninitialized = (1 << 3) - 1;
/* Default values */ /* Default values */
ctx->nucleus.rescale_factor_kappa = 1.0;
return QMCKL_SUCCESS; return QMCKL_SUCCESS;
} }
@ -265,53 +258,6 @@ end interface
#+begin_src c :comments org :tangle (eval h_func) :exports none #+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code qmckl_exit_code
qmckl_get_nucleus_rescale_factor(const qmckl_context context,
double* const rescale_factor_kappa);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_nucleus_rescale_factor (const qmckl_context context,
double* const rescale_factor_kappa)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (rescale_factor_kappa == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_get_nucleus_rescale_factor",
"rescale_factor_kappa is a null pointer");
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
assert (ctx->nucleus.rescale_factor_kappa > 0.0);
(*rescale_factor_kappa) = ctx->nucleus.rescale_factor_kappa;
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_nucleus_rescale_factor(context, kappa) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double) , intent(out) :: kappa
end function qmckl_get_nucleus_rescale_factor
end interface
#+end_src
#+begin_src c :comments org :tangle (eval h_func) :exports none
qmckl_exit_code
qmckl_get_nucleus_coord(const qmckl_context context, qmckl_get_nucleus_coord(const qmckl_context context,
const char transp, const char transp,
double* const coord, double* const coord,
@ -631,55 +577,12 @@ interface
end interface end interface
#+end_src #+end_src
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code
qmckl_set_nucleus_rescale_factor(qmckl_context context,
const double kappa);
#+end_src
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_rescale_factor(qmckl_context context,
const double rescale_factor_kappa)
{
int32_t mask = 0; // Can be updated
<<pre2>>
if (rescale_factor_kappa <= 0.0) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_set_nucleus_rescale_factor",
"rescale_factor_kappa cannot be <= 0.");
}
ctx->nucleus.rescale_factor_kappa = rescale_factor_kappa;
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_set_nucleus_rescale_factor(context, kappa) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
real (c_double) , intent(in) , value :: kappa
end function qmckl_set_nucleus_rescale_factor
end interface
#+end_src
** Test ** Test
#+begin_src c :tangle (eval c_test) #+begin_src c :tangle (eval c_test)
const double* nucl_charge = chbrclf_charge; const double* nucl_charge = chbrclf_charge;
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]); const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
const double nucl_rescale_factor_kappa = 2.0;
/* --- */ /* --- */
@ -693,38 +596,25 @@ assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num); rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(!qmckl_nucleus_provided(context)); assert(!qmckl_nucleus_provided(context));
rc = qmckl_get_nucleus_num (context, &n); rc = qmckl_get_nucleus_num (context, &n);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(n == chbrclf_nucl_num); assert(n == chbrclf_nucl_num);
double k;
rc = qmckl_get_nucleus_rescale_factor (context, &k);
assert(rc == QMCKL_SUCCESS);
assert(k == 1.0);
rc = qmckl_set_nucleus_rescale_factor (context, nucl_rescale_factor_kappa);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_get_nucleus_rescale_factor (context, &k);
assert(rc == QMCKL_SUCCESS);
assert(k == nucl_rescale_factor_kappa);
double nucl_coord2[3*chbrclf_nucl_num]; double nucl_coord2[3*chbrclf_nucl_num];
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*chbrclf_nucl_num); rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*chbrclf_nucl_num);
assert(rc == QMCKL_NOT_PROVIDED); assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*chbrclf_nucl_num); rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert(!qmckl_nucleus_provided(context)); assert(!qmckl_nucleus_provided(context));
rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, 3*chbrclf_nucl_num); rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, 3*chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
for (size_t k=0 ; k<3 ; ++k) { for (size_t k=0 ; k<3 ; ++k) {
for (int64_t i=0 ; i<chbrclf_nucl_num ; ++i) { for (int64_t i=0 ; i<chbrclf_nucl_num ; ++i) {
assert( nucl_coord[chbrclf_nucl_num*k+i] == nucl_coord2[3*i+k] ); assert( nucl_coord[chbrclf_nucl_num*k+i] == nucl_coord2[3*i+k] );
@ -732,7 +622,7 @@ for (size_t k=0 ; k<3 ; ++k) {
} }
rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*chbrclf_nucl_num); rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
for (int64_t i=0 ; i<3*chbrclf_nucl_num ; ++i) { for (int64_t i=0 ; i<3*chbrclf_nucl_num ; ++i) {
assert( nucl_coord[i] == nucl_coord2[i] ); assert( nucl_coord[i] == nucl_coord2[i] );
} }
@ -743,10 +633,10 @@ rc = qmckl_get_nucleus_charge(context, nucl_charge2, chbrclf_nucl_num);
assert(rc == QMCKL_NOT_PROVIDED); assert(rc == QMCKL_NOT_PROVIDED);
rc = qmckl_set_nucleus_charge(context, nucl_charge, chbrclf_nucl_num); rc = qmckl_set_nucleus_charge(context, nucl_charge, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
rc = qmckl_get_nucleus_charge(context, nucl_charge2, chbrclf_nucl_num); rc = qmckl_get_nucleus_charge(context, nucl_charge2, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS); qmckl_check(context, rc);
for (int64_t i=0 ; i<chbrclf_nucl_num ; ++i) { for (int64_t i=0 ; i<chbrclf_nucl_num ; ++i) {
assert( nucl_charge[i] == nucl_charge2[i] ); assert( nucl_charge[i] == nucl_charge2[i] );
} }
@ -952,202 +842,6 @@ assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
#+end_src #+end_src
** Nucleus-nucleus rescaled distances
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code
qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context,
double* distance_rescaled,
const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code
qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context,
double* distance_rescaled,
const int64_t size_max)
{
/* Check input parameters */
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (char) 0;
}
qmckl_exit_code rc = qmckl_provide_nn_distance_rescaled(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->nucleus.num * ctx->nucleus.num;
if (sze > size_max) {
return qmckl_failwith(context,
QMCKL_INVALID_ARG_3,
"qmckl_get_nucleus_nn_distance_rescaled",
"Array too small");
}
rc = qmckl_double_of_matrix(context,
ctx->nucleus.nn_distance_rescaled,
distance_rescaled,
size_max);
return rc;
}
#+end_src
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_get_nucleus_nn_distance_rescaled(context, distance_rescaled, 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) :: distance_rescaled(*)
integer (c_int64_t) , intent(in) , value :: size_max
end function
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_nn_distance_rescaled(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context)
{
/* Check input parameters */
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return (char) 0;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED;
/* Allocate array */
if (ctx->nucleus.nn_distance_rescaled.data == NULL) {
ctx->nucleus.nn_distance_rescaled =
qmckl_matrix_alloc(context, ctx->nucleus.num, ctx->nucleus.num);
if (ctx->nucleus.nn_distance_rescaled.data == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_nn_distance_rescaled",
NULL);
}
}
qmckl_exit_code rc =
qmckl_compute_nn_distance_rescaled(context,
ctx->nucleus.num,
ctx->nucleus.rescale_factor_kappa,
ctx->nucleus.coord.data,
ctx->nucleus.nn_distance_rescaled.data);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->nucleus.nn_distance_rescaled_date = ctx->date;
return QMCKL_SUCCESS;
}
#+end_src
*** Compute
#+NAME: qmckl_nn_distance_rescaled_args
| qmckl_context | context | in | Global state |
| int64_t | nucl_num | in | Number of nuclei |
| double | coord[3][nucl_num] | in | Nuclear coordinates (au) |
| double | nn_distance_rescaled[nucl_num][nucl_num] | out | Nucleus-nucleus rescaled distances (au) |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_nn_distance_rescaled_f(context, nucl_num, rescale_factor_kappa, coord, nn_distance_rescaled) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: nucl_num
double precision , intent(in) :: rescale_factor_kappa
double precision , intent(in) :: coord(nucl_num,3)
double precision , intent(out) :: nn_distance_rescaled(nucl_num,nucl_num)
integer*8 :: k
info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
if (nucl_num <= 0) then
info = QMCKL_INVALID_ARG_2
return
endif
info = qmckl_distance_rescaled(context, 'T', 'T', nucl_num, nucl_num, &
coord, nucl_num, &
coord, nucl_num, &
nn_distance_rescaled, nucl_num, rescale_factor_kappa)
end function qmckl_compute_nn_distance_rescaled_f
#+end_src
#+begin_src c :tangle (eval h_private_func) :comments org :exports none
qmckl_exit_code qmckl_compute_nn_distance_rescaled (
const qmckl_context context,
const int64_t nucl_num,
const double rescale_factor_kappa,
const double* coord,
double* const nn_distance_rescaled );
#+end_src
#+CALL: generate_c_interface(table=qmckl_nn_distance_rescaled_args,rettyp="qmckl_exit_code",fname="qmckl_compute_nn_distance")
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_nn_distance_rescaled &
(context, nucl_num, rescale_factor_kappa, coord, nn_distance_rescaled) &
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 :: nucl_num
real (c_double ) , intent(in) , value :: rescale_factor_kappa
real (c_double ) , intent(in) :: coord(nucl_num,3)
real (c_double ) , intent(out) :: nn_distance_rescaled(nucl_num,nucl_num)
integer(c_int32_t), external :: qmckl_compute_nn_distance_rescaled_f
info = qmckl_compute_nn_distance_rescaled_f &
(context, nucl_num, rescale_factor_kappa, coord, nn_distance_rescaled)
end function qmckl_compute_nn_distance_rescaled
#+end_src
*** Test
#+begin_src c :tangle (eval c_test)
/* Reference input data */
/* TODO */
//assert(qmckl_nucleus_provided(context));
//
//double distance[nucl_num*nucl_num];
//rc = qmckl_get_nucleus_nn_distance(context, distance, nucl_num*nucl_num);
//assert(distance[0] == 0.);
//assert(distance[1] == distance[nucl_num]);
//assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
#+end_src
** Nuclear repulsion energy ** Nuclear repulsion energy
\[ \[

View File

@ -1097,37 +1097,25 @@ qmckl_trexio_read(const qmckl_context context, const char* file_name, const int6
rc = qmckl_trexio_read_electron_X(context, file); rc = qmckl_trexio_read_electron_X(context, file);
if (rc != QMCKL_SUCCESS) { if (rc != QMCKL_SUCCESS) {
trexio_close(file); trexio_close(file);
return qmckl_failwith( context, return rc;
QMCKL_FAILURE,
"qmckl_trexio_read",
"Error reading electron");
} }
rc = qmckl_trexio_read_nucleus_X(context, file); rc = qmckl_trexio_read_nucleus_X(context, file);
if (rc != QMCKL_SUCCESS) { if (rc != QMCKL_SUCCESS) {
trexio_close(file); trexio_close(file);
return qmckl_failwith( context, return rc;
QMCKL_FAILURE,
"qmckl_trexio_read",
"Error reading nucleus");
} }
rc = qmckl_trexio_read_ao_X(context, file); rc = qmckl_trexio_read_ao_X(context, file);
if (rc != QMCKL_SUCCESS) { if (rc != QMCKL_SUCCESS) {
trexio_close(file); trexio_close(file);
return qmckl_failwith( context, return rc;
QMCKL_FAILURE,
"qmckl_trexio_read",
"Error reading AOs");
} }
rc = qmckl_trexio_read_mo_X(context, file); rc = qmckl_trexio_read_mo_X(context, file);
if (rc != QMCKL_SUCCESS) { if (rc != QMCKL_SUCCESS) {
trexio_close(file); trexio_close(file);
return qmckl_failwith( context, return rc;
QMCKL_FAILURE,
"qmckl_trexio_read",
"Error reading MOs");
} }
trexio_close(file); trexio_close(file);
@ -1149,27 +1137,19 @@ qmckl_trexio_read(const qmckl_context context, const char* file_name, const int6
#ifdef HAVE_TREXIO #ifdef HAVE_TREXIO
qmckl_exit_code rc; qmckl_exit_code rc;
char fname[256]; char filename[256];
char message[256];
#ifndef QMCKL_TEST_DIR #ifndef QMCKL_TEST_DIR
#error "QMCKL_TEST_DIR is not defined" #error "QMCKL_TEST_DIR is not defined"
#endif #endif
strncpy(fname, QMCKL_TEST_DIR,255); strncpy(filename, QMCKL_TEST_DIR,255);
strncat(fname, "/chbrclf", 255); strncat(filename, "/chbrclf", 255);
printf("Test file: %s\n", fname); printf("Test file: %s\n", filename);
rc = qmckl_trexio_read(context, fname, 255); rc = qmckl_trexio_read(context, filename, 255);
if (rc != QMCKL_SUCCESS) { qmckl_check(context, rc);
printf("%s\n", qmckl_string_of_error(rc));
qmckl_get_error(context, &rc, fname, message);
printf("%s\n", fname);
printf("%s\n", message);
}
assert ( rc == QMCKL_SUCCESS );
#+end_src #+end_src
@ -1179,11 +1159,11 @@ assert ( rc == QMCKL_SUCCESS );
printf("Electrons\n"); printf("Electrons\n");
int64_t up_num, dn_num; int64_t up_num, dn_num;
rc = qmckl_get_electron_up_num(context, &up_num); rc = qmckl_get_electron_up_num(context, &up_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert (up_num == chbrclf_elec_up_num); assert (up_num == chbrclf_elec_up_num);
rc = qmckl_get_electron_down_num(context, &dn_num); rc = qmckl_get_electron_down_num(context, &dn_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert (dn_num == chbrclf_elec_dn_num); assert (dn_num == chbrclf_elec_dn_num);
#+end_src #+end_src
@ -1195,13 +1175,13 @@ printf("Nuclei\n");
int64_t nucl_num; int64_t nucl_num;
rc = qmckl_get_nucleus_num(context, &nucl_num); rc = qmckl_get_nucleus_num(context, &nucl_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert (nucl_num == chbrclf_nucl_num); assert (nucl_num == chbrclf_nucl_num);
printf("Nuclear charges\n"); printf("Nuclear charges\n");
double * charge = (double*) malloc (nucl_num * sizeof(double)); double * charge = (double*) malloc (nucl_num * sizeof(double));
rc = qmckl_get_nucleus_charge(context, charge, nucl_num); rc = qmckl_get_nucleus_charge(context, charge, nucl_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
for (int i=0 ; i<nucl_num ; i++) { for (int i=0 ; i<nucl_num ; i++) {
assert (charge[i] == chbrclf_charge[i]); assert (charge[i] == chbrclf_charge[i]);
} }
@ -1211,7 +1191,7 @@ charge = NULL;
printf("Nuclear coordinates\n"); printf("Nuclear coordinates\n");
double * coord = (double*) malloc (nucl_num * 3 * sizeof(double)); double * coord = (double*) malloc (nucl_num * 3 * sizeof(double));
rc = qmckl_get_nucleus_coord(context, 'T', coord, 3*nucl_num); rc = qmckl_get_nucleus_coord(context, 'T', coord, 3*nucl_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
int k=0; int k=0;
for (int j=0 ; j<3 ; ++j) { for (int j=0 ; j<3 ; ++j) {
for (int i=0 ; i<nucl_num ; ++i) { for (int i=0 ; i<nucl_num ; ++i) {
@ -1248,7 +1228,7 @@ assert (ao_num == chbrclf_ao_num);
int64_t* nucleus_index = (int64_t*) malloc (nucl_num * sizeof(int64_t)); int64_t* nucleus_index = (int64_t*) malloc (nucl_num * sizeof(int64_t));
rc = qmckl_get_ao_basis_nucleus_index(context, nucleus_index, nucl_num); rc = qmckl_get_ao_basis_nucleus_index(context, nucleus_index, nucl_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
for (int i=0 ; i<nucl_num ; i++) { for (int i=0 ; i<nucl_num ; i++) {
assert (nucleus_index[i] == chbrclf_basis_nucleus_index[i]); assert (nucleus_index[i] == chbrclf_basis_nucleus_index[i]);
} }
@ -1257,7 +1237,7 @@ nucleus_index = NULL;
int64_t* nucleus_shell_num = (int64_t*) malloc (nucl_num * sizeof(int64_t)); int64_t* nucleus_shell_num = (int64_t*) malloc (nucl_num * sizeof(int64_t));
rc = qmckl_get_ao_basis_nucleus_shell_num(context, nucleus_shell_num, nucl_num); rc = qmckl_get_ao_basis_nucleus_shell_num(context, nucleus_shell_num, nucl_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
for (int i=0 ; i<nucl_num ; i++) { for (int i=0 ; i<nucl_num ; i++) {
assert (nucleus_shell_num[i] == chbrclf_basis_nucleus_shell_num[i]); assert (nucleus_shell_num[i] == chbrclf_basis_nucleus_shell_num[i]);
} }
@ -1266,7 +1246,7 @@ nucleus_shell_num = NULL;
int32_t* shell_ang_mom = (int32_t*) malloc (shell_num * sizeof(int32_t)); int32_t* shell_ang_mom = (int32_t*) malloc (shell_num * sizeof(int32_t));
rc = qmckl_get_ao_basis_shell_ang_mom(context, shell_ang_mom, shell_num); rc = qmckl_get_ao_basis_shell_ang_mom(context, shell_ang_mom, shell_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
for (int i=0 ; i<shell_num ; i++) { for (int i=0 ; i<shell_num ; i++) {
assert (shell_ang_mom[i] == chbrclf_basis_shell_ang_mom[i]); assert (shell_ang_mom[i] == chbrclf_basis_shell_ang_mom[i]);
} }
@ -1275,7 +1255,7 @@ shell_ang_mom = NULL;
int64_t* shell_prim_num = (int64_t*) malloc (shell_num * sizeof(int64_t)); int64_t* shell_prim_num = (int64_t*) malloc (shell_num * sizeof(int64_t));
rc = qmckl_get_ao_basis_shell_prim_num(context, shell_prim_num, shell_num); rc = qmckl_get_ao_basis_shell_prim_num(context, shell_prim_num, shell_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
for (int i=0 ; i<shell_num ; i++) { for (int i=0 ; i<shell_num ; i++) {
assert (shell_prim_num[i] == chbrclf_basis_shell_prim_num[i]); assert (shell_prim_num[i] == chbrclf_basis_shell_prim_num[i]);
} }
@ -1284,7 +1264,7 @@ shell_prim_num = NULL;
int64_t* shell_prim_index = (int64_t*) malloc (shell_num * sizeof(int64_t)); int64_t* shell_prim_index = (int64_t*) malloc (shell_num * sizeof(int64_t));
rc = qmckl_get_ao_basis_shell_prim_index(context, shell_prim_index, shell_num); rc = qmckl_get_ao_basis_shell_prim_index(context, shell_prim_index, shell_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
for (int i=0 ; i<shell_num ; i++) { for (int i=0 ; i<shell_num ; i++) {
assert (shell_prim_index[i] == chbrclf_basis_shell_prim_index[i]); assert (shell_prim_index[i] == chbrclf_basis_shell_prim_index[i]);
} }
@ -1293,7 +1273,7 @@ shell_prim_index = NULL;
double* shell_factor = (double*) malloc (shell_num * sizeof(double)); double* shell_factor = (double*) malloc (shell_num * sizeof(double));
rc = qmckl_get_ao_basis_shell_factor(context, shell_factor, shell_num); rc = qmckl_get_ao_basis_shell_factor(context, shell_factor, shell_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
for (int i=0 ; i<shell_num ; i++) { for (int i=0 ; i<shell_num ; i++) {
assert (fabs(shell_factor[i] - chbrclf_basis_shell_factor[i]) < 1.e-6); assert (fabs(shell_factor[i] - chbrclf_basis_shell_factor[i]) < 1.e-6);
} }
@ -1302,7 +1282,7 @@ shell_factor = NULL;
double* exponent = (double*) malloc (prim_num * sizeof(double)); double* exponent = (double*) malloc (prim_num * sizeof(double));
rc = qmckl_get_ao_basis_exponent(context, exponent, prim_num); rc = qmckl_get_ao_basis_exponent(context, exponent, prim_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
for (int i=0 ; i<prim_num ; i++) { for (int i=0 ; i<prim_num ; i++) {
assert (fabs((exponent[i] - chbrclf_basis_exponent[i])/chbrclf_basis_exponent[i]) < 1.e-7); assert (fabs((exponent[i] - chbrclf_basis_exponent[i])/chbrclf_basis_exponent[i]) < 1.e-7);
} }
@ -1311,7 +1291,7 @@ exponent = NULL;
double* coefficient = (double*) malloc (prim_num * sizeof(double)); double* coefficient = (double*) malloc (prim_num * sizeof(double));
rc = qmckl_get_ao_basis_coefficient(context, coefficient, prim_num); rc = qmckl_get_ao_basis_coefficient(context, coefficient, prim_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
for (int i=0 ; i<prim_num ; i++) { for (int i=0 ; i<prim_num ; i++) {
assert (fabs((coefficient[i] - chbrclf_basis_coefficient[i])/chbrclf_basis_coefficient[i]) < 1.e-7); assert (fabs((coefficient[i] - chbrclf_basis_coefficient[i])/chbrclf_basis_coefficient[i]) < 1.e-7);
} }
@ -1320,7 +1300,7 @@ coefficient = NULL;
double* prim_factor = (double*) malloc (prim_num * sizeof(double)); double* prim_factor = (double*) malloc (prim_num * sizeof(double));
rc = qmckl_get_ao_basis_prim_factor(context, prim_factor, prim_num); rc = qmckl_get_ao_basis_prim_factor(context, prim_factor, prim_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
for (int i=0 ; i<prim_num ; i++) { for (int i=0 ; i<prim_num ; i++) {
assert (fabs((prim_factor[i] - chbrclf_basis_prim_factor[i])/chbrclf_basis_prim_factor[i]) < 1.e-7); assert (fabs((prim_factor[i] - chbrclf_basis_prim_factor[i])/chbrclf_basis_prim_factor[i]) < 1.e-7);
} }
@ -1337,13 +1317,14 @@ printf("MOs\n");
int64_t mo_num; int64_t mo_num;
rc = qmckl_get_mo_basis_mo_num(context, &mo_num); rc = qmckl_get_mo_basis_mo_num(context, &mo_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
assert (mo_num == chbrclf_mo_num); assert (mo_num == chbrclf_mo_num);
printf("MO coefs\n"); printf("MO coefs\n");
double * mo_coef = (double*) malloc (ao_num * mo_num * sizeof(double)); double * mo_coef = (double*) malloc (ao_num * mo_num * sizeof(double));
rc = qmckl_get_mo_basis_coefficient(context, mo_coef, mo_num*ao_num); rc = qmckl_get_mo_basis_coefficient(context, mo_coef, mo_num*ao_num);
assert (rc == QMCKL_SUCCESS); qmckl_check(context, rc);
for (int i=0 ; i<ao_num * mo_num ; i++) { for (int i=0 ; i<ao_num * mo_num ; i++) {
printf("%d %e %e %e\n", i, mo_coef[i], chbrclf_mo_coef[i], printf("%d %e %e %e\n", i, mo_coef[i], chbrclf_mo_coef[i],
( fabs(mo_coef[i] - chbrclf_mo_coef[i])/fabs(mo_coef[i])) ); ( fabs(mo_coef[i] - chbrclf_mo_coef[i])/fabs(mo_coef[i])) );