hdf5 backend working

src/trio.c

@@ -129,11 +129,11 @@ trio_exit_code trio_read_nucleus_num(trio_t* file, uint64_t* num) {
   case TRIO_TEXT:
     return trio_text_read_nucleus_num(file, num);
     break;
-/*
+
   case TRIO_HDF5:
     return trio_hdf5_read_nucleus_num(file, num);
     break;
-
+/*
   case TRIO_JSON:
     return trio_json_read_nucleus_num(file, num);
     break;
@@ -151,11 +151,11 @@ trio_exit_code trio_write_nucleus_num(trio_t* file, uint64_t num) {
   case TRIO_TEXT:
     return trio_text_write_nucleus_num(file, num);
     break;
-/*
+
   case TRIO_HDF5:
     return trio_hdf5_write_nucleus_num(file, num);
     break;
-
+/*
   case TRIO_JSON:
     return trio_json_write_nucleus_num(file, num);
     break;
@@ -195,11 +195,11 @@ trio_exit_code trio_write_nucleus_coord(trio_t* file, double* coord) {
   case TRIO_TEXT:
     return trio_text_write_nucleus_coord(file, coord);
     break;
-/*
+
   case TRIO_HDF5:
     return trio_hdf5_write_nucleus_coord(file, coord);
     break;
-
+/*
   case TRIO_JSON:
     return trio_json_write_nucleus_coord(file, coord);
     break;

src/trio.org

@@ -251,7 +251,7 @@ trio_t* trio_create(const char* file_name, back_end_t back_end) {
   case TRIO_JSON:
     rc = trio_json_init(result);
     break;
-*/      
+,*/      
   default:
     assert (1 == 0); /* Impossible case */
   }
@@ -324,11 +324,11 @@ trio_exit_code trio_read_nucleus_num(trio_t* file, uint64_t* num) {
   case TRIO_TEXT:
     return trio_text_read_nucleus_num(file, num);
     break;
-/*
+
   case TRIO_HDF5:
     return trio_hdf5_read_nucleus_num(file, num);
     break;
-
+/*
   case TRIO_JSON:
     return trio_json_read_nucleus_num(file, num);
     break;
@@ -346,11 +346,11 @@ trio_exit_code trio_write_nucleus_num(trio_t* file, uint64_t num) {
   case TRIO_TEXT:
     return trio_text_write_nucleus_num(file, num);
     break;
-/*
+
   case TRIO_HDF5:
     return trio_hdf5_write_nucleus_num(file, num);
     break;
-
+/*
   case TRIO_JSON:
     return trio_json_write_nucleus_num(file, num);
     break;
@@ -398,11 +398,11 @@ trio_exit_code trio_write_nucleus_coord(trio_t* file, double* coord) {
   case TRIO_TEXT:
     return trio_text_write_nucleus_coord(file, coord);
     break;
-/*
+
   case TRIO_HDF5:
     return trio_hdf5_write_nucleus_coord(file, coord);
     break;
-
+/*
   case TRIO_JSON:
     return trio_json_write_nucleus_coord(file, coord);
     break;
@@ -839,6 +839,16 @@ trio_exit_code trio_text_write_nucleus_charge(const trio_t* file, const double*
     #+end_src
 ** HDF5 Back end
 
+*** HDF5 definitions
+
+#+begin_src c :tangle trio_hdf5.c
+  #define NUCLEUS_GROUP_NAME  "nucleus"
+  #define NUCLEUS_NUM_NAME    "nucleus_num"
+  #define NUCLEUS_CHARGE_NAME "nucleus_charge"
+  #define NUCLEUS_COORD_NAME  "nucleus_coord"
+  
+#+end_src
+
 *** HDF5 structures
 
   #+begin_src c :tangle trio_hdf5.h
@@ -892,6 +902,58 @@ trio_exit_code trio_hdf5_init(trio_t* file);
 
 
   #+begin_src c :tangle trio_hdf5.c
+
+/* 
+ * Currently H5LTread_dataset_ is used instead of this function 
+ * but keep it for later if we decide to get rid of the H5LT API 
+ */
+dset_t* trio_hdf5_read_dset_low(const trio_hdf5_t* file, const char *dset_name, void *buf) {
+
+  assert (file != NULL);
+  assert (dset_name != NULL);
+  assert (buf != NULL);
+  /*
+   * Low-level implementation. Involves dealing with all HDF5 handles and dimensions
+   */
+  dset_t* dset = (dset_t*) malloc(sizeof(dset_t));
+  assert (dset != NULL);
+
+  dset->dset_id = H5Dopen(file->nucleus_group, 
+                          dset_name, 
+                          H5P_DEFAULT);
+  assert (dset->dset_id > 0);
+  /*
+   * Get dataspace, datatype and dimensions
+   * dspace and dtype handles created below have to be closed when not used
+   */    
+  dset->dspace_id = H5Dget_space(dset->dset_id);
+  assert (dset->dspace_id > 0);
+  
+  dset->dtype_id = H5Dget_type(dset->dset_id);
+  assert (dset->dtype_id > 0);
+  
+  /* Check dimensions. Usefull, but then additional parameters 
+   * ranks and dims[] have to be passed to the function
+  int rrank;
+  const int rank = 1;
+  hsize_t dims[1] = {0};
+  rrank = H5Sget_simple_extent_dims(nucleus->h5_charge->dspace_id, 
+                                    dims, NULL);
+  assert (rrank == rank);
+  for (int i=0; i<rank; i++){
+    assert (dims[i] > 0);
+  }
+  */
+  herr_t status;
+  status = H5Dread(dset->dset_id, dset->dtype_id,
+               H5S_ALL, H5S_ALL, H5P_DEFAULT,
+	       buf);
+  assert (status >= 0);
+
+  return dset;
+
+}
+
 trio_exit_code trio_hdf5_init(trio_t* file) {
 
   trio_hdf5_t* f = (trio_hdf5_t*) file;
@@ -911,18 +973,13 @@ trio_exit_code trio_hdf5_init(trio_t* file) {
   }  
 
  /* Create groups in the hdf5 file */
-  const char* nucleus_group_name = "nucleus";
-  //const char* electron_group_name = "electron";
- 
   if (f_ishere == 0){
-    f->nucleus_group = H5Gcreate(f->file_id, nucleus_group_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+    f->nucleus_group = H5Gcreate(f->file_id, NUCLEUS_GROUP_NAME, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-    //f->electron_group = H5Gcreate(f->file_id, electron_group_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+    //f->electron_group = H5Gcreate(f->file_id, ELECTRON_GROUP_NAME, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
   } else { 
-    f->nucleus_group = H5Gopen(f->file_id, nucleus_group_name, H5P_DEFAULT);
+    f->nucleus_group = H5Gopen(f->file_id, NUCLEUS_GROUP_NAME, H5P_DEFAULT);
-    //f->electron_group = H5Gopen(f->file_id, electron_group_name, H5P_DEFAULT);   
+    //f->electron_group = H5Gopen(f->file_id, ELECTRON_GROUP_NAME, H5P_DEFAULT);   
   }
-  /* not sure if assert statement here makes sence 
-     H5Gcreate will raise its own H5 error if somethings is wrong*/
   assert (f->nucleus_group > 0L);
   //assert (f->electron_group > 0L);
 
@@ -948,6 +1005,10 @@ trio_exit_code trio_hdf5_finalize(trio_t* file) {
   H5Gclose(f->electron_group);
   f->electron_group = 0; 
 */
+
+  H5Fclose(f->file_id);
+  f->file_id = 0;
+
   return TRIO_SUCCESS;
 }
   #+end_src
@@ -969,14 +1030,16 @@ h5nucleus_t* trio_hdf5_read_nucleus(const trio_hdf5_t* file) {
   nucleus->h5_coord  = NULL;
   nucleus->h5_charge = NULL;
 
- /* Try to open the file. If HDF5 cannot open, return */ 
+  /* Check that the file was opened/created correctly, return */ 
   if (file->file_id < 0) return nucleus;
 
+ /* Quit if the dimensioning attribute is missing in the file */
+  if (H5Aexists(file->nucleus_group, NUCLEUS_NUM_NAME) == 0) return nucleus;
+
   herr_t status;  
  /* Read the nucleus_num attribute of nucleus group */
-  const char *num_name = "nucleus_num";
   hid_t num_id;
-  num_id = H5Aopen(file->nucleus_group, num_name, H5P_DEFAULT);
+  num_id = H5Aopen(file->nucleus_group, NUCLEUS_NUM_NAME, H5P_DEFAULT);
   assert (num_id > 0);
 
   status = H5Aread(num_id, H5T_NATIVE_ULLONG, &(nucleus->num));
@@ -986,85 +1049,152 @@ h5nucleus_t* trio_hdf5_read_nucleus(const trio_hdf5_t* file) {
   nucleus->charge = (double*) calloc(nucleus->num, sizeof(double));
   assert (nucleus->charge != NULL);
 
-  /*
+ /* High-level H5LT API. No need to deal with dataspaces and datatypes */
-   * High-level H5LT API. No need to deal with dataspaces and datatypes
-   */
-  /*  
   status = H5LTread_dataset_double(file->nucleus_group, 
-                                "nucleus_charge", 
+                                NUCLEUS_CHARGE_NAME, 
-				nucleus->charge)
-  */
-  /*
-   * Low-level implementation. Involves dealing with all HDF5 handles and dimensions
-   */
-  nucleus->h5_charge = (dset_t*) malloc(sizeof(dset_t));
-  assert (nucleus->h5_charge != NULL);
-
-  nucleus->h5_charge->dset_id = H5Dopen(file->nucleus_group, 
-                                        "nucleus_charge", 
-                                        H5P_DEFAULT);
-  assert (nucleus->h5_charge->dset_id > 0);
-  /*
-   * Get dataspace, datatype and dimensions
-   * dspace and dtype handles created below have to be closed when not used
-   */    
-  nucleus->h5_charge->dspace_id = H5Dget_space(nucleus->h5_charge->dset_id);
-  assert (nucleus->h5_charge->dspace_id > 0);
-  
-  nucleus->h5_charge->dtype_id = H5Dget_type(nucleus->h5_charge->dset_id);
-  assert (nucleus->h5_charge->dtype_id > 0);
-  
-  int rrank;
-  const int rank = 1;
-  hsize_t dims[1] = {0};
-  rrank = H5Sget_simple_extent_dims(nucleus->h5_charge->dspace_id, 
-                                    dims, NULL);
-  assert (rrank == rank);
-  for (int i=0; i<rank; i++){
-    assert (dims[i] > 0);
-  }
-
-  status = H5Dread(nucleus->h5_charge->dset_id, nucleus->h5_charge->dtype_id,
-               H5S_ALL, H5S_ALL, H5P_DEFAULT,
 				nucleus->charge);
-  assert (status >= 0);
 
  /* Allocate and read nucleus_coord array */
   nucleus->coord = (double*) calloc(3 * nucleus->num, sizeof(double));
   assert (nucleus->coord != NULL);
 
+ /* High-level H5LT API. No need to deal with dataspaces and datatypes */
   status = H5LTread_dataset_double(file->nucleus_group, 
-                                    "nucleus_coord", 
+                                    NUCLEUS_COORD_NAME, 
 				    nucleus->coord);
   assert (status >= 0);
- /* Print arrays */
-  /*
-  for (size_t i=0 ; i<nucleus->num ; i++) {
-    printf("%lf \n", nucleus->charge[i]);
-  }
 
-  for (size_t i=0 ; i<3*nucleus->num ; i++) {
+  /* Low-level read. Do not forget to close the associated IDs (dset,dtype,dspace) 
-      printf("%lf \n", nucleus->coord[i]);
+   *  when not used anymore, see below. Note how this function is similar to H5LTread_dataset_double
-  }
    */
+  /*
+  nucleus->h5_coord = trio_hdf5_read_dset_low(file, NUCLEUS_COORD_NAME, 
+  		  		nucleus->coord); 
+
+  H5Sclose(nucleus->h5_coord->dspace_id);
+  H5Tclose(nucleus->h5_coord->dtype_id);
+  H5Dclose(nucleus->h5_coord->dset_id);
+  */
+
   H5Aclose(num_id);
-  H5Sclose(nucleus->h5_charge->dspace_id);
-  H5Tclose(nucleus->h5_charge->dtype_id);
-  H5Dclose(nucleus->h5_charge->dset_id);
-  H5Fclose(file->file_id);
 
   return nucleus;
 }
 
 
 trio_exit_code trio_hdf5_write_nucleus(const trio_hdf5_t* file, h5nucleus_t* nucleus) {
-  assert (nucleus != NULL);
 
   assert (file != NULL);
+  assert (nucleus != NULL);
 
-  // TODO
+  herr_t status;
-  return TRIO_FAILURE;
+  hid_t dspace, dtype; 
+  hid_t attr_id;
+    
+  dtype = H5Tcopy(H5T_NATIVE_ULLONG);
+ /* Write the dimensioning variables */
+  if (H5Aexists(file->nucleus_group, NUCLEUS_NUM_NAME) == 0) {
+    dspace = H5Screate(H5S_SCALAR);
+    attr_id = H5Acreate(file->nucleus_group, NUCLEUS_NUM_NAME, dtype, dspace,
+                     H5P_DEFAULT, H5P_DEFAULT);
+    assert (attr_id > 0);
+
+    /* High-level routine does not work for some reason
+     * status = H5LTset_attribute_ulong (file->nucleus_group, "nucleus", NUCLEUS_NUM_NAME, 
+     *                                  &(nucleus->num), 1);
+     */
+  } else {
+    attr_id = H5Aopen(file->nucleus_group, NUCLEUS_NUM_NAME, H5P_DEFAULT);
+    assert (attr_id > 0);
+  }
+
+  status = H5Awrite(attr_id, dtype, &(nucleus->num));
+  assert (status >= 0); 
+
+  H5Aclose(attr_id);
+  
+ /* Write arrays */
+  hid_t dset_id;
+  int charge_rank = 1;
+  const hsize_t charge_dims[1] = {nucleus->num};
+
+  if ( H5LTfind_dataset(file->nucleus_group, NUCLEUS_CHARGE_NAME) != 1) { 
+    
+    status = H5LTmake_dataset_double (file->nucleus_group, NUCLEUS_CHARGE_NAME, 
+                                    charge_rank, charge_dims, nucleus->charge);
+    assert (status >= 0);
+
+  } else {   
+
+    dset_id = H5Dopen(file->nucleus_group, NUCLEUS_CHARGE_NAME, H5P_DEFAULT);
+    assert (dset_id > 0);
+    
+    dspace = H5Dget_space(dset_id);
+    assert (dspace > 0);
+
+    dtype = H5Dget_type(dset_id);
+    assert (dtype > 0);
+
+    
+    int rrank;
+    hsize_t dims[1] = {0};
+    rrank = H5Sget_simple_extent_dims(dspace,
+                                    dims, NULL);
+    assert (rrank == charge_rank);
+    // disabling asserts like this allows to overwrite _num variable
+    for (int i=0; i<rrank; i++){
+      assert (dims[i] == charge_dims[i]);
+    }
+
+    status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, nucleus->charge);
+    assert (status >= 0);
+
+    H5Sclose(dspace);
+    H5Tclose(dtype);
+    H5Dclose(dset_id);
+
+  }
+
+  int coord_rank = 2;
+  const hsize_t coord_dims[2] = {nucleus->num, 3};
+  if ( H5LTfind_dataset(file->nucleus_group, NUCLEUS_COORD_NAME) != 1) { 
+    status = H5LTmake_dataset_double (file->nucleus_group, NUCLEUS_COORD_NAME, 
+                                    coord_rank, coord_dims, nucleus->coord);
+    assert (status >= 0);
+
+  } else {   
+
+    dset_id = H5Dopen(file->nucleus_group, NUCLEUS_COORD_NAME, H5P_DEFAULT);
+    assert (dset_id > 0);
+    
+    dspace = H5Dget_space(dset_id);
+    assert (dspace > 0);
+
+    dtype = H5Dget_type(dset_id);
+    assert (dtype > 0);
+
+    
+    int rrank;
+    hsize_t dims[2] = {0, 0};
+    rrank = H5Sget_simple_extent_dims(dspace,
+                                    dims, NULL);
+    assert (rrank == coord_rank);
+    for (int i=0; i<rrank; i++){
+      assert (dims[i] == coord_dims[i]);
+    }
+
+    status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, nucleus->coord);
+    assert (status >= 0);
+
+    H5Sclose(dspace);
+    H5Tclose(dtype);
+    H5Dclose(dset_id);
+
+  }
+
+  return TRIO_SUCCESS;
 }
+
   #+end_src
 
 *** Free memory
@@ -1094,7 +1224,77 @@ trio_exit_code trio_hdf5_free_nucleus(h5nucleus_t* nucleus) {
 }
    #+end_src
 
-*** TODO Read/Write the num attribute
+*** Read/Write the num attribute
+
+   #+begin_src c :tangle trio_hdf5.h
+trio_exit_code trio_hdf5_read_nucleus_num(const trio_t* file, uint64_t* num);
+trio_exit_code trio_hdf5_write_nucleus_num(const trio_t* file, const uint64_t num);
+  #+end_src
+
+   #+begin_src c :tangle trio_hdf5.c
+trio_exit_code trio_hdf5_read_nucleus_num(const trio_t* file, uint64_t* num) {
+
+  assert (file != NULL);
+  assert (num  != NULL);
+
+  h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file);
+  
+  if (nucleus == NULL) {
+    return TRIO_FAILURE;
+  }
+
+ /**/ *num = nucleus->num;
+
+  trio_hdf5_free_nucleus(nucleus);
+  return TRIO_SUCCESS;
+}
+
+ 
+trio_exit_code trio_hdf5_write_nucleus_num(const trio_t* file, const uint64_t num) {
+
+  assert (file != NULL);
+  assert (num > 0L);
+  
+  h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file);
+
+  assert (nucleus != NULL);
+  
+  if (nucleus->num != num) {
+
+    if (nucleus->num != 0) {
+	printf("%ld -> %ld %s \n", num, nucleus->num, 
+	       "This variable alreasy exists. Overwriting it is not supported");
+        trio_hdf5_free_nucleus(nucleus);
+    	return TRIO_FAILURE;
+    }
+
+    nucleus->num = num;
+
+    if (nucleus->charge != NULL) free(nucleus->charge);
+    nucleus->charge = NULL;
+
+    nucleus->charge = (double*) calloc(num, sizeof(double));
+    assert (nucleus->charge != NULL);
+    
+    if (nucleus->coord  != NULL) free(nucleus->coord );
+    nucleus->coord = NULL;
+
+    nucleus->coord = (double*) calloc(3*num, sizeof(double));
+    assert (nucleus->coord != NULL);
+
+  } else {
+    nucleus->num = num;
+  }
+  
+  trio_exit_code rc = trio_hdf5_write_nucleus((trio_hdf5_t*) file, nucleus);
+  assert (rc == TRIO_SUCCESS);
+
+  trio_hdf5_free_nucleus(nucleus);
+  
+  return TRIO_SUCCESS;
+}
+    #+end_src
+
 
 *** Read/Write the coord attribute
 
@@ -1102,19 +1302,19 @@ trio_exit_code trio_hdf5_free_nucleus(h5nucleus_t* nucleus) {
     
    #+begin_src c :tangle trio_hdf5.h
 trio_exit_code trio_hdf5_read_nucleus_coord(const trio_t* file, double* coord);
-// TODO
+trio_exit_code trio_hdf5_write_nucleus_coord(const trio_t* file, const double* coord);
-//trio_exit_code trio_hdf5_write_nucleus_coord(const trio_t* file, const double* coord);
   #+end_src
 
    #+begin_src c :tangle trio_hdf5.c
 trio_exit_code trio_hdf5_read_nucleus_coord(const trio_t* file, double* coord) {
 
   assert (file != NULL);
+  assert (coord != NULL);
+
   h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file);
 
   if (nucleus == NULL) return TRIO_FAILURE;
-
+  assert (nucleus->coord != NULL);
-  assert (coord != NULL);
   
   for (size_t i=0 ; i<3*nucleus->num ; i++) {
     coord[i] = nucleus->coord[i];
@@ -1126,11 +1326,27 @@ trio_exit_code trio_hdf5_read_nucleus_coord(const trio_t* file, double* coord) {
 
  
 trio_exit_code trio_hdf5_write_nucleus_coord(const trio_t* file, const double* coord) {
+
   assert (file != NULL);
   assert (coord != NULL);
-  // TODO
+  
-  return TRIO_FAILURE;
+  h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file);
+
+  if (nucleus == NULL) return TRIO_FAILURE;
+  assert (nucleus->coord != NULL); 
+ 
+  for (size_t i=0 ; i<3*nucleus->num ; i++) {
+    nucleus->coord[i] = coord[i];
+  }
+  
+  trio_exit_code rc = trio_hdf5_write_nucleus((trio_hdf5_t*) file, nucleus);
+  assert (rc == TRIO_SUCCESS);
+
+  trio_hdf5_free_nucleus(nucleus);
+  
+  return TRIO_SUCCESS;
 }
+
     #+end_src
 *** TODO Read/Write the charge attribute
 

src/trio_hdf5.c

@@ -7,6 +7,62 @@
 
 #include "trio_hdf5.h"
 
+  #define NUCLEUS_GROUP_NAME  "nucleus"
+  #define NUCLEUS_NUM_NAME    "nucleus_num"
+  #define NUCLEUS_CHARGE_NAME "nucleus_charge"
+  #define NUCLEUS_COORD_NAME  "nucleus_coord"
+
+/* 
+ * Currently H5LTread_dataset_ is used instead of this function 
+ * but keep it for later if we decide to get rid of the H5LT API 
+ */
+dset_t* trio_hdf5_read_dset_low(const trio_hdf5_t* file, const char *dset_name, void *buf) {
+
+  assert (file != NULL);
+  assert (dset_name != NULL);
+  assert (buf != NULL);
+  /*
+   * Low-level implementation. Involves dealing with all HDF5 handles and dimensions
+   */
+  dset_t* dset = (dset_t*) malloc(sizeof(dset_t));
+  assert (dset != NULL);
+
+  dset->dset_id = H5Dopen(file->nucleus_group, 
+                          dset_name, 
+                          H5P_DEFAULT);
+  assert (dset->dset_id > 0);
+  /*
+   * Get dataspace, datatype and dimensions
+   * dspace and dtype handles created below have to be closed when not used
+   */    
+  dset->dspace_id = H5Dget_space(dset->dset_id);
+  assert (dset->dspace_id > 0);
+  
+  dset->dtype_id = H5Dget_type(dset->dset_id);
+  assert (dset->dtype_id > 0);
+  
+  /* Check dimensions. Usefull, but then additional parameters 
+   * ranks and dims[] have to be passed to the function
+  int rrank;
+  const int rank = 1;
+  hsize_t dims[1] = {0};
+  rrank = H5Sget_simple_extent_dims(nucleus->h5_charge->dspace_id, 
+                                    dims, NULL);
+  assert (rrank == rank);
+  for (int i=0; i<rank; i++){
+    assert (dims[i] > 0);
+  }
+  */
+  herr_t status;
+  status = H5Dread(dset->dset_id, dset->dtype_id,
+               H5S_ALL, H5S_ALL, H5P_DEFAULT,
+	       buf);
+  assert (status >= 0);
+
+  return dset;
+
+}
+
 trio_exit_code trio_hdf5_init(trio_t* file) {
 
   trio_hdf5_t* f = (trio_hdf5_t*) file;
@@ -26,18 +82,13 @@ trio_exit_code trio_hdf5_init(trio_t* file) {
   }  
 
   /* Create groups in the hdf5 file */
-  const char* nucleus_group_name = "nucleus";
-  //const char* electron_group_name = "electron";
- 
   if (f_ishere == 0){
-    f->nucleus_group = H5Gcreate(f->file_id, nucleus_group_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+    f->nucleus_group = H5Gcreate(f->file_id, NUCLEUS_GROUP_NAME, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-    //f->electron_group = H5Gcreate(f->file_id, electron_group_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+    //f->electron_group = H5Gcreate(f->file_id, ELECTRON_GROUP_NAME, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
   } else { 
-    f->nucleus_group = H5Gopen(f->file_id, nucleus_group_name, H5P_DEFAULT);
+    f->nucleus_group = H5Gopen(f->file_id, NUCLEUS_GROUP_NAME, H5P_DEFAULT);
-    //f->electron_group = H5Gopen(f->file_id, electron_group_name, H5P_DEFAULT);   
+    //f->electron_group = H5Gopen(f->file_id, ELECTRON_GROUP_NAME, H5P_DEFAULT);   
   }
-  /* not sure if assert statement here makes sence 
-     H5Gcreate will raise its own H5 error if somethings is wrong*/
   assert (f->nucleus_group > 0L);
   //assert (f->electron_group > 0L);
 
@@ -55,6 +106,10 @@ trio_exit_code trio_hdf5_finalize(trio_t* file) {
   H5Gclose(f->electron_group);
   f->electron_group = 0; 
 */
+
+  H5Fclose(f->file_id);
+  f->file_id = 0;
+
   return TRIO_SUCCESS;
 }
 
@@ -70,14 +125,16 @@ h5nucleus_t* trio_hdf5_read_nucleus(const trio_hdf5_t* file) {
   nucleus->h5_coord  = NULL;
   nucleus->h5_charge = NULL;
 
-  /* Try to open the file. If HDF5 cannot open, return */ 
+  /* Check that the file was opened/created correctly, return */ 
   if (file->file_id < 0) return nucleus;
 
+  /* Quit if the dimensioning attribute is missing in the file */
+  if (H5Aexists(file->nucleus_group, NUCLEUS_NUM_NAME) == 0) return nucleus;
+
   herr_t status;  
   /* Read the nucleus_num attribute of nucleus group */
-  const char *num_name = "nucleus_num";
   hid_t num_id;
-  num_id = H5Aopen(file->nucleus_group, num_name, H5P_DEFAULT);
+  num_id = H5Aopen(file->nucleus_group, NUCLEUS_NUM_NAME, H5P_DEFAULT);
   assert (num_id > 0);
 
   status = H5Aread(num_id, H5T_NATIVE_ULLONG, &(nucleus->num));
@@ -87,84 +144,150 @@ h5nucleus_t* trio_hdf5_read_nucleus(const trio_hdf5_t* file) {
   nucleus->charge = (double*) calloc(nucleus->num, sizeof(double));
   assert (nucleus->charge != NULL);
 
-  /*
+  /* High-level H5LT API. No need to deal with dataspaces and datatypes */
-   * High-level H5LT API. No need to deal with dataspaces and datatypes
-   */
-  /*  
   status = H5LTread_dataset_double(file->nucleus_group, 
-                                "nucleus_charge", 
+                                NUCLEUS_CHARGE_NAME, 
-				nucleus->charge)
-  */
-  /*
-   * Low-level implementation. Involves dealing with all HDF5 handles and dimensions
-   */
-  nucleus->h5_charge = (dset_t*) malloc(sizeof(dset_t));
-  assert (nucleus->h5_charge != NULL);
-
-  nucleus->h5_charge->dset_id = H5Dopen(file->nucleus_group, 
-                                        "nucleus_charge", 
-                                        H5P_DEFAULT);
-  assert (nucleus->h5_charge->dset_id > 0);
-  /*
-   * Get dataspace, datatype and dimensions
-   * dspace and dtype handles created below have to be closed when not used
-   */    
-  nucleus->h5_charge->dspace_id = H5Dget_space(nucleus->h5_charge->dset_id);
-  assert (nucleus->h5_charge->dspace_id > 0);
-  
-  nucleus->h5_charge->dtype_id = H5Dget_type(nucleus->h5_charge->dset_id);
-  assert (nucleus->h5_charge->dtype_id > 0);
-  
-  int rrank;
-  const int rank = 1;
-  hsize_t dims[1] = {0};
-  rrank = H5Sget_simple_extent_dims(nucleus->h5_charge->dspace_id, 
-                                    dims, NULL);
-  assert (rrank == rank);
-  for (int i=0; i<rank; i++){
-    assert (dims[i] > 0);
-  }
-
-  status = H5Dread(nucleus->h5_charge->dset_id, nucleus->h5_charge->dtype_id,
-               H5S_ALL, H5S_ALL, H5P_DEFAULT,
 				nucleus->charge);
-  assert (status >= 0);
 
   /* Allocate and read nucleus_coord array */
   nucleus->coord = (double*) calloc(3 * nucleus->num, sizeof(double));
   assert (nucleus->coord != NULL);
 
+  /* High-level H5LT API. No need to deal with dataspaces and datatypes */
   status = H5LTread_dataset_double(file->nucleus_group, 
-                                    "nucleus_coord", 
+                                    NUCLEUS_COORD_NAME, 
 				    nucleus->coord);
   assert (status >= 0);
-  /* Print arrays */
-  /*
-  for (size_t i=0 ; i<nucleus->num ; i++) {
-    printf("%lf \n", nucleus->charge[i]);
-  }
 
-  for (size_t i=0 ; i<3*nucleus->num ; i++) {
+  /* Low-level read. Do not forget to close the associated IDs (dset,dtype,dspace) 
-      printf("%lf \n", nucleus->coord[i]);
+   *  when not used anymore, see below. Note how this function is similar to H5LTread_dataset_double
-  }
    */
+  /*
+  nucleus->h5_coord = trio_hdf5_read_dset_low(file, NUCLEUS_COORD_NAME, 
+  		  		nucleus->coord); 
+
+  H5Sclose(nucleus->h5_coord->dspace_id);
+  H5Tclose(nucleus->h5_coord->dtype_id);
+  H5Dclose(nucleus->h5_coord->dset_id);
+  */
+
   H5Aclose(num_id);
-  H5Sclose(nucleus->h5_charge->dspace_id);
-  H5Tclose(nucleus->h5_charge->dtype_id);
-  H5Dclose(nucleus->h5_charge->dset_id);
-  H5Fclose(file->file_id);
 
   return nucleus;
 }
 
 
 trio_exit_code trio_hdf5_write_nucleus(const trio_hdf5_t* file, h5nucleus_t* nucleus) {
-  assert (nucleus != NULL);
 
   assert (file != NULL);
+  assert (nucleus != NULL);
 
-  // TODO
+  herr_t status;
-  return TRIO_FAILURE;
+  hid_t dspace, dtype; 
+  hid_t attr_id;
+    
+  dtype = H5Tcopy(H5T_NATIVE_ULLONG);
+  /* Write the dimensioning variables */
+  if (H5Aexists(file->nucleus_group, NUCLEUS_NUM_NAME) == 0) {
+    dspace = H5Screate(H5S_SCALAR);
+    attr_id = H5Acreate(file->nucleus_group, NUCLEUS_NUM_NAME, dtype, dspace,
+                     H5P_DEFAULT, H5P_DEFAULT);
+    assert (attr_id > 0);
+
+    /* High-level routine does not work for some reason
+     * status = H5LTset_attribute_ulong (file->nucleus_group, "nucleus", NUCLEUS_NUM_NAME, 
+     *                                  &(nucleus->num), 1);
+     */
+  } else {
+    attr_id = H5Aopen(file->nucleus_group, NUCLEUS_NUM_NAME, H5P_DEFAULT);
+    assert (attr_id > 0);
+  }
+
+  status = H5Awrite(attr_id, dtype, &(nucleus->num));
+  assert (status >= 0); 
+
+  H5Aclose(attr_id);
+  
+  /* Write arrays */
+  hid_t dset_id;
+  int charge_rank = 1;
+  const hsize_t charge_dims[1] = {nucleus->num};
+
+  if ( H5LTfind_dataset(file->nucleus_group, NUCLEUS_CHARGE_NAME) != 1) { 
+    
+    status = H5LTmake_dataset_double (file->nucleus_group, NUCLEUS_CHARGE_NAME, 
+                                    charge_rank, charge_dims, nucleus->charge);
+    assert (status >= 0);
+
+  } else {   
+
+    dset_id = H5Dopen(file->nucleus_group, NUCLEUS_CHARGE_NAME, H5P_DEFAULT);
+    assert (dset_id > 0);
+    
+    dspace = H5Dget_space(dset_id);
+    assert (dspace > 0);
+
+    dtype = H5Dget_type(dset_id);
+    assert (dtype > 0);
+
+    
+    int rrank;
+    hsize_t dims[1] = {0};
+    rrank = H5Sget_simple_extent_dims(dspace,
+                                    dims, NULL);
+    assert (rrank == charge_rank);
+    // disabling asserts like this allows to overwrite _num variable
+    for (int i=0; i<rrank; i++){
+      assert (dims[i] == charge_dims[i]);
+    }
+
+    status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, nucleus->charge);
+    assert (status >= 0);
+
+    H5Sclose(dspace);
+    H5Tclose(dtype);
+    H5Dclose(dset_id);
+
+  }
+
+  int coord_rank = 2;
+  const hsize_t coord_dims[2] = {nucleus->num, 3};
+  if ( H5LTfind_dataset(file->nucleus_group, NUCLEUS_COORD_NAME) != 1) { 
+    status = H5LTmake_dataset_double (file->nucleus_group, NUCLEUS_COORD_NAME, 
+                                    coord_rank, coord_dims, nucleus->coord);
+    assert (status >= 0);
+
+  } else {   
+
+    dset_id = H5Dopen(file->nucleus_group, NUCLEUS_COORD_NAME, H5P_DEFAULT);
+    assert (dset_id > 0);
+    
+    dspace = H5Dget_space(dset_id);
+    assert (dspace > 0);
+
+    dtype = H5Dget_type(dset_id);
+    assert (dtype > 0);
+
+    
+    int rrank;
+    hsize_t dims[2] = {0, 0};
+    rrank = H5Sget_simple_extent_dims(dspace,
+                                    dims, NULL);
+    assert (rrank == coord_rank);
+    for (int i=0; i<rrank; i++){
+      assert (dims[i] == coord_dims[i]);
+    }
+
+    status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, nucleus->coord);
+    assert (status >= 0);
+
+    H5Sclose(dspace);
+    H5Tclose(dtype);
+    H5Dclose(dset_id);
+
+  }
+
+  return TRIO_SUCCESS;
 }
 
 trio_exit_code trio_hdf5_free_nucleus(h5nucleus_t* nucleus) {
@@ -188,14 +311,77 @@ trio_exit_code trio_hdf5_free_nucleus(h5nucleus_t* nucleus) {
   return TRIO_SUCCESS;
 }
 
+trio_exit_code trio_hdf5_read_nucleus_num(const trio_t* file, uint64_t* num) {
+
+  assert (file != NULL);
+  assert (num  != NULL);
+
+  h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file);
+  
+  if (nucleus == NULL) {
+    return TRIO_FAILURE;
+  }
+
+  /**/ *num = nucleus->num;
+
+  trio_hdf5_free_nucleus(nucleus);
+  return TRIO_SUCCESS;
+}
+
+ 
+trio_exit_code trio_hdf5_write_nucleus_num(const trio_t* file, const uint64_t num) {
+
+  assert (file != NULL);
+  assert (num > 0L);
+  
+  h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file);
+
+  assert (nucleus != NULL);
+  
+  if (nucleus->num != num) {
+
+    if (nucleus->num != 0) {
+	printf("%ld -> %ld %s \n", num, nucleus->num, 
+	       "This variable alreasy exists. Overwriting it is not supported");
+        trio_hdf5_free_nucleus(nucleus);
+    	return TRIO_FAILURE;
+    }
+
+    nucleus->num = num;
+
+    if (nucleus->charge != NULL) free(nucleus->charge);
+    nucleus->charge = NULL;
+
+    nucleus->charge = (double*) calloc(num, sizeof(double));
+    assert (nucleus->charge != NULL);
+    
+    if (nucleus->coord  != NULL) free(nucleus->coord );
+    nucleus->coord = NULL;
+
+    nucleus->coord = (double*) calloc(3*num, sizeof(double));
+    assert (nucleus->coord != NULL);
+
+  } else {
+    nucleus->num = num;
+  }
+  
+  trio_exit_code rc = trio_hdf5_write_nucleus((trio_hdf5_t*) file, nucleus);
+  assert (rc == TRIO_SUCCESS);
+
+  trio_hdf5_free_nucleus(nucleus);
+  
+  return TRIO_SUCCESS;
+}
+
 trio_exit_code trio_hdf5_read_nucleus_coord(const trio_t* file, double* coord) {
 
   assert (file != NULL);
+  assert (coord != NULL);
+
   h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file);
 
   if (nucleus == NULL) return TRIO_FAILURE;
-
+  assert (nucleus->coord != NULL);
-  assert (coord != NULL);
   
   for (size_t i=0 ; i<3*nucleus->num ; i++) {
     coord[i] = nucleus->coord[i];
@@ -207,8 +393,23 @@ trio_exit_code trio_hdf5_read_nucleus_coord(const trio_t* file, double* coord) {
 
  
 trio_exit_code trio_hdf5_write_nucleus_coord(const trio_t* file, const double* coord) {
+
   assert (file != NULL);
   assert (coord != NULL);
-  // TODO
+  
-  return TRIO_FAILURE;
+  h5nucleus_t* nucleus = trio_hdf5_read_nucleus((trio_hdf5_t*) file);
+
+  if (nucleus == NULL) return TRIO_FAILURE;
+  assert (nucleus->coord != NULL); 
+ 
+  for (size_t i=0 ; i<3*nucleus->num ; i++) {
+    nucleus->coord[i] = coord[i];
+  }
+  
+  trio_exit_code rc = trio_hdf5_write_nucleus((trio_hdf5_t*) file, nucleus);
+  assert (rc == TRIO_SUCCESS);
+
+  trio_hdf5_free_nucleus(nucleus);
+  
+  return TRIO_SUCCESS;
 }

src/trio_hdf5.h

@@ -62,8 +62,10 @@ trio_exit_code trio_hdf5_init(trio_t* file);
 
 trio_exit_code trio_hdf5_finalize(trio_t* file);
 
+trio_exit_code trio_hdf5_read_nucleus_num(const trio_t* file, uint64_t* num);
+trio_exit_code trio_hdf5_write_nucleus_num(const trio_t* file, const uint64_t num);
+
 trio_exit_code trio_hdf5_read_nucleus_coord(const trio_t* file, double* coord);
-// TODO
+trio_exit_code trio_hdf5_write_nucleus_coord(const trio_t* file, const double* coord);
-//trio_exit_code trio_hdf5_write_nucleus_coord(const trio_t* file, const double* coord);
 
 #endif