Merge pull request #1 from GMPavanLab/addTests

Adding test

.gitignore

@@ -8,6 +8,8 @@
 *.xyz
 test.py
 log.lammps
+launch.json
+
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[cod]

Changelog.md

@@ -0,0 +1,8 @@
+# Changes From 0.1
+
+- Tests
+- Changed default imports for HDF5er
+- Adding override option to HDF5er.MDA2HDF5
+- Added 3rd neighbours in the References
+- Added attributes to HDF5er.MDA2HDF5
+- In the referenceMaker: added the possibility tyo choose lmax and nmax for SOAP

Readme.md

@@ -1,4 +1,4 @@
-# SOAPify and ReferenceMaker
+# SOAPify, HDF5er and ReferenceMaker
 
 This sequence of commands will setup the environment:
 
@@ -12,11 +12,15 @@ pip install -r requirements
 
 This package contains a toolbox to calculate the [SOAP fingerprints](https://doi.org/10.1103/PhysRevB.87.184115) of a system of atoms.
 
+## HDF5er
+
+This package contains a small toolbox to create [hdf5 files](https://www.hdfgroup.org/) with [h5py](https://www.h5py.org/) from trajectory and topology files. The format we use **do not** align with [h5md](https://www.nongnu.org/h5md/h5md.html).
+
 ## ReferenceMaker
 
 The ReferenceMaker package contains a set of function that can create a reference file to be used with the SOAPify package.
 
-ReferenceMaker function can be called with custommade scripts, but the user can create a list of SOAP references with the following:
+ReferenceMaker function can be called with custom made scripts, but the user can create a list of SOAP references with the following:
 
 ```bash
 python3 -m ReferenceMaker

setup.py

@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="SOAPify-GMPavanLab",
-    version="0.0.1",
+    version="0.0.2",
     author="Daniele Rapetti",
     author_email="[email protected]",
     description="A package for creating and studying SOAP fingerprints",
@@ -22,5 +22,5 @@
     ],
     package_dir={"": "src"},
     packages=setuptools.find_packages(where="src"),
-    python_requires=">=3.6",
+    python_requires=">=3.8",  # due to MDA
 )

src/HDF5er/MDA2HDF5.py

@@ -2,11 +2,14 @@
 import h5py
 from .HDF5erUtils import exportChunk2HDF5
 
+__all__ = ["Universe2HDF5", "MDA2HDF5"]
+
 
 def Universe2HDF5(
     MDAUniverseOrSelection: "MDAnalysis.Universe | MDAnalysis.AtomGroup",
     trajFolder: h5py.Group,
     trajChunkSize: int = 100,
+    trajslice: slice = slice(None),
 ):
     """Uploads an mda.Universe or an mda.AtomGroup to a h5py.Group in an hdf5 file
 
@@ -17,11 +20,7 @@ def Universe2HDF5(
     """
 
     atoms = MDAUniverseOrSelection.atoms
-    universe = (
-        MDAUniverseOrSelection
-        if MDAUniverseOrSelection is MDAnalysis.Universe
-        else MDAUniverseOrSelection.universe
-    )
+    universe = MDAUniverseOrSelection.universe
     nat = len(atoms)
 
     if "Types" not in list(trajFolder.keys()):
@@ -45,7 +44,7 @@ def Universe2HDF5(
     first = 0
     boxes = []
     atomicframes = []
-    for frame in universe.trajectory:
+    for frame in universe.trajectory[trajslice]:
         boxes.append(universe.dimensions)
         atomicframes.append(atoms.positions)
         frameNum += 1
@@ -66,6 +65,9 @@ def MDA2HDF5(
     targetHDF5File: str,
     groupName: str,
     trajChunkSize: int = 100,
+    override: bool = False,
+    attrs: dict = None,
+    trajslice: slice = slice(None),
 ):
     """Opens or creates the given HDF5 file, request the user's chosen group, then uploads an mda.Universe or an mda.AtomGroup to a h5py.Group in an hdf5 file
 
@@ -76,7 +78,16 @@ def MDA2HDF5(
         targetHDF5File (str): the name of HDF5 file
         groupName (str): the name of the group in wich save the trajectory data within the `targetHDF5File`
         trajChunkSize (int, optional): The desired dimension of the chunks of data that are stored in the hdf5 file. Defaults to 100.
+        override (bool, optional): If true the hdf5 file will be completely overwritten. Defaults to False.
     """
-    with h5py.File(targetHDF5File, "a") as newTraj:
+    with h5py.File(targetHDF5File, "w" if override else "a") as newTraj:
         trajGroup = newTraj.require_group(f"Trajectories/{groupName}")
-        Universe2HDF5(MDAUniverseOrSelection, trajGroup, trajChunkSize=trajChunkSize)
+        Universe2HDF5(
+            MDAUniverseOrSelection,
+            trajGroup,
+            trajChunkSize=trajChunkSize,
+            trajslice=trajslice,
+        )
+        if attrs:
+            for key in attrs.keys():
+                trajGroup.attrs.create(key, attrs[key])

src/HDF5er/ase2HDF5.py

@@ -4,6 +4,10 @@
 import numpy as np
 import h5py
 
+# this is "legacy code", use with caution
+
+__all__ = ["xyz2hdf5Converter", "HDF52AseAtomsChunckedwithSymbols"]
+
 # TODO: convert use exportChunk2HDF5
 def xyz2hdf5Converter(xyzName: str, boxfilename: str, group: h5py.Group):
     """Generate an HDF5 trajectory from an xyz file and a box file
@@ -76,6 +80,7 @@ def xyz2hdf5Converter(xyzName: str, boxfilename: str, group: h5py.Group):
             group["Trajectory"][first:frameNum] = atomicframes
 
 
+# TODO: using slices is not the best compromise here
 def HDF52AseAtomsChunckedwithSymbols(
     groupTraj: h5py.Group,
     chunkTraj: "tuple[slice]",
@@ -99,13 +104,13 @@ def HDF52AseAtomsChunckedwithSymbols(
     for frame, box in zip(
         groupTraj["Trajectory"][chunkTraj], groupTraj["Box"][chunkBox]
     ):
-        theBox = [[box[0], 0, 0], [0, box[1], 0], [0, 0, box[2]]]
+        # theBox = [[box[0], 0, 0], [0, box[1], 0], [0, 0, box[2]]]
         # celldisp = -box[0:3] / 2
         atoms.append(
             aseAtoms(
                 symbols=symbols,
                 positions=frame,
-                cell=theBox,
+                cell=box,
                 pbc=True,
                 # celldisp=celldisp,
             )

src/ReferenceMaker/__main__.py

@@ -72,6 +72,8 @@
 pair_coeff = f"pair_coeff	1 1 {chosen['pair_coeff']}"
 atomMass = chosen["atomMass"]
 diameter = chosen["diameter"]
+SOAPlmax = 8
+SOAPnmax = 8
 referencesFileName = f"{kind}References.hdf5"
 latticebcc = 2 * diameter / numpy.sqrt(3.0)
 latticefcc = diameter * numpy.sqrt(2.0)
@@ -83,6 +85,7 @@
         (2.82843, "LattUn"),  # 1 lattice unit
         (4.0, "4R"),  # 4r
         (5.65685, "2LattUn"),  # 2 lattice unit
+        (3.45659722, "3rdNeighFCC"),  # 2 lattice unit
     ]
 ]
 for d in [
@@ -147,5 +150,11 @@
         ]
         lmp.commands_list(commands_list)
 
-referenceSaponificator(rcuts=rcuts, referencesFileName=referencesFileName, kind=kind)
+referenceSaponificator(
+    rcuts=rcuts,
+    referencesFileName=referencesFileName,
+    kind=kind,
+    SOAPlmax=SOAPlmax,
+    SOAPnmax=SOAPnmax,
+)
 distanceVisualizer(rcuts=rcuts, referencesFile=referencesFileName, kind=kind)

src/ReferenceMaker/distanceVisualizer.py

@@ -48,7 +48,7 @@ def distanceVisualizer(rcuts: "list[radiusInfo]", referencesFile: str, kind: str
 
         df = pd.DataFrame(distanceMatrix, index=legend, columns=legend)
         outFname = f"distmat_{kind}Classification_{rcut.name}"
-        # df.to_csv(f"{outFname}.csv")
+        df.to_csv(f"{outFname}.csv", float_format="%14.8g")
         center = df.values[np.tril_indices(df.shape[0], k=1)].mean()
 
         g = sns.clustermap(

src/ReferenceMaker/referenceSaponificator.py

@@ -115,6 +115,8 @@ def RefCreator(
     structureFile: str,
     masks: "list[dict[str,numpy.ndarray]]",
     PBC: bool = False,
+    SOAPlmax: int = 8,
+    SOAPnmax: int = 8,
 ):
     """Creates a series of SOAP fingerprints and store them in a subgroup in the requested group
 
@@ -140,19 +142,30 @@ def RefCreator(
     outgroup = targetFile.require_group(f"{groupname}/{rcut.name}")
     outgroup.attrs["rcutName"] = rcut.name
     outgroup.attrs["rcut"] = rcut.rcut
-    outgroup.attrs["lmax"] = 8
-    outgroup.attrs["nmax"] = 8
+    outgroup.attrs["lmax"] = SOAPlmax
+    outgroup.attrs["nmax"] = SOAPnmax
     for mask in masks:
         structure = structureFile
         export2hdf5(
-            fingerprintMaker(structure, rcut.rcut, mask["mask"], SOAPpbc=PBC),
+            fingerprintMaker(
+                structure,
+                rcut.rcut,
+                mask["mask"],
+                SOAPpbc=PBC,
+                SOAPlmax=SOAPlmax,
+                SOAPnmax=SOAPnmax,
+            ),
             outgroup,
             mask["name"],
         )
 
 
 def referenceSaponificator(
-    rcuts: "list[radiusInfo]", referencesFileName: str, kind: str
+    rcuts: "list[radiusInfo]",
+    referencesFileName: str,
+    kind: str,
+    SOAPlmax: int = 8,
+    SOAPnmax: int = 8,
 ):
     """generates the SOAP fingerprints for the reference systems.
 
@@ -195,8 +208,8 @@ def referenceSaponificator(
             outgroup = referenceFile.require_group(f"Bulk/{rcut.name}")
             outgroup.attrs["rcutName"] = rcut.name
             outgroup.attrs["rcut"] = rcut.rcut
-            outgroup.attrs["lmax"] = 8
-            outgroup.attrs["nmax"] = 8
+            outgroup.attrs["lmax"] = SOAPlmax
+            outgroup.attrs["nmax"] = SOAPnmax
             for structure in [
                 f"{kind}_bcc.data",
                 f"{kind}_fcc.data",
@@ -206,7 +219,11 @@ def referenceSaponificator(
                 structureName = structure.rsplit(sep=".", maxsplit=1)[0]
                 structureName = structureName.rsplit(sep="_", maxsplit=1)[1]
                 export2hdf5(
-                    fingerprintMaker(structure, rcut.rcut), outgroup, structureName
+                    fingerprintMaker(
+                        structure, rcut.rcut, SOAPlmax=SOAPlmax, SOAPnmax=SOAPnmax
+                    ),
+                    outgroup,
+                    structureName,
                 )
 
             RefCreator(
@@ -221,6 +238,8 @@ def referenceSaponificator(
                     {"name": "fiveFoldedAxis_ico", "mask": maskFiveFoldedAxis_ico},
                 ],
                 PBC=False,
+                SOAPlmax=SOAPlmax,
+                SOAPnmax=SOAPnmax,
             )
             RefCreator(
                 rcut,
@@ -234,6 +253,8 @@ def referenceSaponificator(
                     {"name": "edges_th", "mask": maskEdges_th},
                 ],
                 PBC=False,
+                SOAPlmax=SOAPlmax,
+                SOAPnmax=SOAPnmax,
             )
 
             RefCreator(
@@ -247,6 +268,8 @@ def referenceSaponificator(
                     {"name": "fiveFoldedAxis_dh", "mask": maskFiveFoldedAxis_dh},
                 ],
                 PBC=False,
+                SOAPlmax=SOAPlmax,
+                SOAPnmax=SOAPnmax,
             )