TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-11-19 20:42:50 +01:00
Code Issues Projects Releases Wiki Activity
15d07b89a4
qmckl/org/qmckl_trexio.org

9.2 KiB
Raw Blame History

TREXIO I/O library

The TREXIO library enables easy and efficient input/output of wave function parameters. In this section we provide high-level functions to prepare the context by reading the required data from a TREXIO file.

Local functions

Functions defined in this section are all local: they should not be exposed in the API. To identify them, we append _X to their name.

Open file

We first define a helper function to open a file by first trying to use the TEXT back end, and then the HDF5 back end. If both strategies fail, a NULL pointer is returned. This will allow to open only once the file and call multiple small functions to read groups of data by passing the trexio_t handle.

#ifdef HAVE_TREXIO
trexio_t* qmckl_trexio_open_X(const char* file_name, qmckl_exit_code* rc)
{
*rc = QMCKL_SUCCESS;
trexio_t* file = NULL;

file = trexio_open(file_name, 'r', TREXIO_TEXT, rc);
if (file != NULL) return file;

file = trexio_open(file_name, 'r', TREXIO_HDF5, rc);
if (file != NULL) return file;
 
*rc = QMCKL_FAILURE;
return NULL;
}
#endif

Electron

In this section we read all the data into the electron data structure.

#ifdef HAVE_TREXIO
qmckl_exit_code
qmckl_trexio_read_electron_X(qmckl_context context, trexio_t* const file) 
{
// Should not be possible by construction.
assert (context != (qmckl_context) 0);
assert (file != NULL);

int rcio = 0;

int64_t up_num = 0;
int64_t dn_num = 0;

rcio = trexio_read_electron_up_num_64(file, &up_num);
if (rcio != TREXIO_SUCCESS) {
 return qmckl_failwith( context,
                        QMCKL_FAILURE,
                        "trexio_read_electron_up_num",
                        trexio_string_of_error(rcio));
}

assert (up_num >= 0);

rcio = trexio_read_electron_dn_num_64(file, &dn_num);
if (rcio != TREXIO_SUCCESS) {
 return qmckl_failwith( context,
                        QMCKL_FAILURE,
                        "trexio_read_electron_dn_num",
                        trexio_string_of_error(rcio));
}

assert (dn_num >= 0);


qmckl_exit_code rc;
rc = qmckl_set_electron_num(context, up_num, dn_num);
return rc;
}
#endif

Nucleus

In this section we read the molecular geometry and nuclear charges.

#ifdef HAVE_TREXIO
qmckl_exit_code
qmckl_trexio_read_nucleus_X(qmckl_context context, trexio_t* const file) 
{
// Should not be possible by construction.
assert (context != (qmckl_context) 0);
assert (file != NULL);

qmckl_exit_code rc;
int rcio = 0;

// Nucleus num
int64_t num = 0;

rcio = trexio_read_nucleus_num_64(file, &num);
if (rcio != TREXIO_SUCCESS) {
 return qmckl_failwith( context,
                        QMCKL_FAILURE,
                        "trexio_read_nucleus_num",
                        trexio_string_of_error(rcio));
}

assert (num > 0);
rc = qmckl_set_nucleus_num(context, num);


// Nuclear charges
{
 qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
 mem_info.size = num * sizeof(double);

 double* nucl_charge = (double*) qmckl_malloc(context, mem_info);

 if (nucl_charge == NULL) {
   return qmckl_failwith( context,
                          QMCKL_ALLOCATION_FAILED,
                          "qmckl_trexio_read_nucleus_X",
                          NULL);
 }

 assert (nucl_charge != NULL);

 rcio = trexio_read_nucleus_charge_64(file, nucl_charge);
 if (rcio != TREXIO_SUCCESS) {
   return qmckl_failwith( context,
                          QMCKL_FAILURE,
                          "trexio_read_nucleus_charge",
                          trexio_string_of_error(rcio));
 }

 rc = qmckl_set_nucleus_charge(context, nucl_charge);

 qmckl_free(context, nucl_charge);
 nucl_charge = NULL;

 if (rc != QMCKL_SUCCESS) {
   return rc;
 }
}

// Nuclear coordinates
{
 qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
 mem_info.size = num * 3 * sizeof(double);

 double* nucl_coord = (double*) qmckl_malloc(context, mem_info);

 if (nucl_coord == NULL) {
   return qmckl_failwith( context,
                          QMCKL_ALLOCATION_FAILED,
                          "qmckl_trexio_read_nucleus_X",
                          NULL);
 }

 assert (nucl_coord != NULL);

 rcio = trexio_read_nucleus_coord_64(file, nucl_coord);
 if (rcio != TREXIO_SUCCESS) {
   return qmckl_failwith( context,
                          QMCKL_FAILURE,
                          "trexio_read_nucleus_charge",
                          trexio_string_of_error(rcio));
 }

 rc = qmckl_set_nucleus_coord(context, 'N', nucl_coord);

 qmckl_free(context, nucl_coord);
 nucl_coord = NULL;

 if (rc != QMCKL_SUCCESS) {
   return rc;
 }
}

return QMCKL_SUCCESS;
}
#endif

Read everything 

qmckl_exit_code qmckl_trexio_read(const qmckl_context context, const char* file_name);
qmckl_exit_code
qmckl_trexio_read(const qmckl_context context, const char* file_name)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
  return false;
}

qmckl_exit_code rc;

#ifdef HAVE_TREXIO
trexio_t* file = qmckl_trexio_open_X(file_name, &rc);
if (file == NULL) {
  trexio_close(file);
  return qmckl_failwith( context,
                         QMCKL_INVALID_ARG_2,
                         "qmckl_trexio_read",
                         trexio_string_of_error(rc));
}

assert (file != NULL);

rc = qmckl_trexio_read_electron_X(context, file);
if (rc != QMCKL_SUCCESS) {
  trexio_close(file);
  return qmckl_failwith( context,
                         QMCKL_FAILURE,
                         "qmckl_trexio_read",
                         "Error reading electron");
}

rc = qmckl_trexio_read_nucleus_X(context, file);
if (rc != QMCKL_SUCCESS) {
  trexio_close(file);
  return qmckl_failwith( context,
                         QMCKL_FAILURE,
                         "qmckl_trexio_read",
                         "Error reading nucleus");
}

trexio_close(file);
file = NULL;
#else

rc = qmckl_failwith( context,
                     QMCKL_FAILURE,
                     "qmckl_trexio_read",
                     "QMCkl was not compiled without TREXIO");
#endif
return rc;
}

Test 

#ifdef HAVE_TREXIO

qmckl_exit_code rc;
char fname[256];
char message[256];
rc = qmckl_trexio_read(context, "share/qmckl/test_data/chbrclf");
if (rc != QMCKL_SUCCESS) {
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 );

printf("Electrons");
int64_t num;
rc = qmckl_get_electron_up_num(context, &num);
assert (rc == QMCKL_SUCCESS);
assert (num == chbrclf_elec_up_num);

rc = qmckl_get_electron_down_num(context, &num);
assert (rc == QMCKL_SUCCESS);
assert (num == chbrclf_elec_dn_num);

rc = qmckl_get_nucleus_num(context, &num);
assert (rc == QMCKL_SUCCESS);
assert (num == chbrclf_nucl_num);

printf("Nuclear charges");
double * charge = (double*) malloc (num * sizeof(double));
rc = qmckl_get_nucleus_charge(context, charge);
assert (rc == QMCKL_SUCCESS);
for (int i=0 ; i<num ; i++) {
assert (charge[i] == chbrclf_charge[i]);
}
free(charge);
charge = NULL;

printf("Nuclear coordinates");
double * coord  = (double*) malloc (num * 3 * sizeof(double));
rc = qmckl_get_nucleus_coord(context, 'T', coord);
assert (rc == QMCKL_SUCCESS);
int k=0;
for (int j=0 ; j<3 ; ++j) {
for (int i=0 ; i<num ; ++i) {
//    printf("%f  %f\n", coord[k], chbrclf_nucl_coord[j][i]);
  assert (coord[k] == chbrclf_nucl_coord[j][i]);
  ++k;
}
}
free(coord);
coord = NULL;


#endif