fix(file IO), fix(global initialization), fix(mass output), feat(channel-bound atoms) (#7)

* Minor file IO bug fix

Fixed a minor bug related to file IO.
Description:
Not providing a .rad file for read_from_CIF, read_from_V1 and read_from_ARC resulted in a type error.

Code to reproduce:
from pyzeo.netstorage import AtomNetwork
atmnet = AtomNetwork.read_from_CIF("EDI.cif")

Example error output:
File "src/pyzeo/extension.pyx", line 512, in pyzeo.extension.AtomNetwork.read_from_CIF
TypeError: expected bytes, NoneType found

Comment:
A unit test "test_read_file.py" was added for validation.

* Minor global state pollution bug fix

Fixed a minor bug related to global state pollution of atomic property tables.
Description:
The global atomic property tables are not reevaluted in subsequent
calculations if only using default zeo++ provided values.

Example:
Reading an AtomNetwork from a file using default Zeo++ atomic radii does not correctly reinitialize the default atomic radii if an AtomNetwork was read in beforehand using custom atomic radii.

Code to reproduce:
from pyzeo.netstorage import AtomNetwork
from pyzeo.extension import lookupRadius
atmnet = AtomNetwork.read_from_CSSR("MgO_vac1.cssr", rad_file="MgO.rad")
atmnet = AtomNetwork.read_from_CSSR("EDI.cssr")

element = "O"
encoded_str = element.encode("utf-8")

element_radius = lookupRadius(encoded_str)

assert element_radius == 1.52, \
    f"Wrong element radius. Expected 1.52, got {element_radius}"

Example error output:                                                                                             AssertionError: Wrong element radius. Expected 1.52, got 1.84

Comment:
A test "test_rad_table_init.py" was added for validation.

* Bug fix for mass-related properties

Mass-related properties in surface area and volume calculations are now displayed correctly.

Description:
Previously mass-related properties (Density, surface area and volume
referenced to mass) were not displayed correctly.

Example code:
from pyzeo.netstorage import AtomNetwork
from pyzeo.area_volume import volume
from pyzeo.high_accuracy import high_accuracy_atomnet

atmnet = AtomNetwork.read_from_CSSR("EDI.cssr")
ha_atmnet = atmnet.copy()
high_accuracy_atomnet(ha_atmnet, "DEF")

vol_str = volume(
    atmnet,
    1.2,
    1.2,
    50000,
    high_accuracy=True,
    high_accuracy_atmnet=ha_atmnet
    )

decoded_str = vol_str.decode("utf-8")
print(decoded_str)

Previous example output:
@ No filename Unitcell_volume: 307.484   Density: 0   AV_A^3: 22.6493 AV_Volume_fraction: 0.07366 AV_cm^3/g: inf NAV_A^3: 0 NAV_Volume_fraction: 0 NAV_cm^3/g: -nan

New example output:
@ No filename Unitcell_volume: 307.484   Density: 1.62239   AV_A^3: 22.6493 AV_Volume_fraction: 0.07366 AV_cm^3/g: 0.0454022 NAV_A^3: 0 NAV_Volume_fraction: 0 NAV_cm^3/g: 0

Comment:
The tests "test_areavol_mass.py" and "test_surface_area_mass.py" were added for validation of the mass-related outputs in the volume and surface area calculations. The test "test_mass_table_init.py" was added to test the correct initialization of the atomic mass table.

* Added feature: Retrieve channel-bound atom IDs

The IDs of atoms that bound a channel (set of accessible Voronoi nodes) that is found in the structure via Voronoi decomposition, can now be returned.

New python methods:
-VoronoiNetwork.find_channels(channel_radius: float) -> List[Channel], List[bool]
    Find channels in a Voronoi network.

    Identifies channels (sets of accessible Voronoi nodes) within the associated Voronoi network for a given probe radius.

    Args:
       channel_radius : float
         Radius of probe used to determine the accessibility of void space.

    Returns:
      1) channels: List of Channel objects representing accessible channels
      2) access_info: List of Booleans, where access_info[i] indicates
               if Voronoi node i is accessible for the probe

-Channel.find_bounding_atoms(atmnet: AtomNetwork, vornet: VoronoiNetwork) -> List[int]
    Get the IDs of all atoms that bound this channel.
    An atom is considered to bound a channel if a Voronoi node of the channel is a member of the atom's Voronoi cell.

    Args:
      atmnet: AtomNetwork object
        The atom network structure of the associated channel
      vornet: VoronoiNetwork object
        The Voronoi network structure of the associated channel

    Returns:
      List of indices of the atoms that bound this channel

New python properties:
-AtomNetwork.atoms: List[Atom]
-Atom.type: str

Newly declared Cython attributes:
-VoronoiNetwork.bvcells
-VoronoiNetwork.has_bvcells

Basic usage example:

    from pyzeo.netstorage import AtomNetwork

    atmnet = AtomNetwork.read_from_CSSR("EDI.cssr")
    probe_radius = 1.5

    vornet, _, _ = atmnet.perform_voronoi_decomposition()
    channels, access_info = vornet.find_channels(probe_radius)
    num_accessible = access_info.count(True)

    print(f"Found {len(channels)} channels in the structure")
    print(f"Accessibility info for {len(access_info)} nodes")
    print(f"Number of accessible nodes {num_accessible}")

    for channel in channels:
      atom_ids = channel.find_bounding_atoms(atmnet, vornet)

* Restructured tests and small test bug fix

Author: Tom-Schaechtel

src/channel.cc

@@ -801,6 +801,7 @@ CHANNEL::CHANNEL(PORE *p){
     unitCells = p->unitCells;
     ucNodes = p->ucNodes;
     dimensionality = p->dimensionality;
+    idMappings = p->idMappings;
     basis[0][0]=p->basis[0][0];basis[0][1]=p->basis[0][1];basis[0][2]=p->basis[0][2];
     basis[1][0]=p->basis[1][0];basis[1][1]=p->basis[1][1];basis[1][2]=p->basis[1][2];
     basis[2][0]=p->basis[2][0];basis[2][1]=p->basis[2][1];basis[2][2]=p->basis[2][2];

src/networkinfo.cc

@@ -33,6 +33,7 @@ stripAtomNameInternalFlag = value;
 
 /** Fills the radius table with several default values. */
 void initializeRadTable(){
+radTable.clear();
 //radTable.insert(pair <string,double> ("Symbol",  vdW_Radius));
 radTable.insert(pair <string,double> ("H",  1.09));
 radTable.insert(pair <string,double> ("D",  1.09));
@@ -154,6 +155,7 @@ radTable.insert(pair <string,double> ("Ds",  2));
  * threshold +0.4A
  * */
 void initializeCovRadTable(){
+covRadTable.clear();
 //covRadTable.insert(pair <string,double> ("Symbol",  Covalent_Radius));
 covRadTable.insert(pair <string,double> ("H",  0.23));
 covRadTable.insert(pair <string,double> ("D",  0.23));
@@ -271,6 +273,7 @@ covRadTable.insert(pair <string,double> ("Ds",  1.5));
 /** Fills the mass table with several default values
  ** in units of g/mole. */
 void initializeMassTable(){
+massTable.clear();
 //massTable.insert(pair <string,double> ("Symbol",  Atomic_Weight));
 massTable.insert(pair <string,double> ("H",  1.008));
 massTable.insert(pair <string,double> ("D",  2.0));
@@ -388,6 +391,7 @@ massTable.insert(pair <string,double> ("Ds",  271));
 
 /** Fills atom metal/nonmetal information **/
 void  initializeAtomCharacterTable(){
+ atomicCharacterTable.clear();
   // value set to true if atom is metal
 
 //  atomicCharacterTable.insert(pair <string,bool> ("Symbol",true));
@@ -508,6 +512,7 @@ void  initializeAtomCharacterTable(){
 /** Fills the atomic number table with atomic number of all elements
  ** */
 void initializeAtomicNumberTable(){
+  atomicNumberTable.clear();
   atomicNumberTable.insert(pair <string,int> ("H", 1));
   atomicNumberTable.insert(pair <string,int> ("D", 1));
   atomicNumberTable.insert(pair <string,int> ("He", 2));

src/pyzeo/extension.pxd

@@ -12,6 +12,7 @@ from libcpp.string cimport string
 from libcpp.vector cimport vector
 from libcpp.map cimport map as cmap
 from libcpp.set cimport set as cset
+from libcpp cimport bool as cpp_bool
 
 #=============================================================================
 # geometry
@@ -78,12 +79,14 @@ cdef extern from "../networkinfo.h":
 cdef extern from "../channel.h":
     cdef cppclass CHANNEL:
         CHANNEL() except +
+        cmap[int, int] idMappings
+        void findBoundingAtoms(ATOM_NETWORK*, vector[BASIC_VCELL], vector[int])
 
 cdef extern from "../channel.h" namespace "CHANNEL":
-    cdef c_findChannelsInDijkstraNet "findChannels"(DIJKSTRA_NETWORK*, 
-            vector[bint] *, vector[CHANNEL] *)
-    cdef c_findChannelsInVorNet "findChannels"(VORONOI_NETWORK*, double, 
-            vector[bint] *, vector[CHANNEL] *)
+    cdef void c_findChannelsInDijkstraNet "CHANNEL::findChannels"(DIJKSTRA_NETWORK*, 
+            vector[cpp_bool] *, vector[CHANNEL] *)
+    cdef void c_findChannelsInVorNet "CHANNEL::findChannels"(VORONOI_NETWORK*, double, 
+            vector[cpp_bool] *, vector[CHANNEL] *)
 
 cdef class Channel:
     cdef CHANNEL* thisptr
@@ -141,7 +144,7 @@ cdef extern from "../networkstorage.h":
         ATOM() except +
         double x, y, z
         double radius
-        #string type
+        string type
         #int specialID
         double mass
         double charge 
@@ -200,7 +203,7 @@ cdef extern from "../network.h":
 
     cdef void loadRadii(ATOM_NETWORK*)
 
-    cdef void loadMass(bool, ATOM_NETWORK*)
+    cdef void loadMass(bint, ATOM_NETWORK*)
 
 cdef extern from "../area_and_volume.h":
     cdef void visVoro(char* name, double probeRad, int skel_a, int skel_b, int skel_c,
@@ -230,8 +233,13 @@ cdef class VoronoiNode:
 cdef class VoronoiNetwork:
     """ 
     Cython wrapper class for Zeo++ VORONOI_NETWORK class.
+    Contains a pointer to VORONOI_NETWORK, a vector
+    of BASIC_VCELL objects and a flag denoting whether
+    the BASIC_VCELL vector has already been populated.
     """
     cdef VORONOI_NETWORK* thisptr
+    cdef vector[BASIC_VCELL] bvcells
+    cdef bint has_bvcells
 
 #=============================================================================
 # netstorage
@@ -377,4 +385,4 @@ cdef extern from "../holograms.h":
 #=============================================================================
 # string_add
 cdef extern from "../string_additions.h":
-    cdef int strCmpList(vector[string] list, string str)
+    cdef int strCmpList(vector[string] list, string str)