WIP - Amber writer

eex/__init__.py

@@ -10,6 +10,7 @@
 from . import translators
 from . import utility
 from . import form_converters
+from . import metadata
 from .units import ureg
 
 

eex/datalayer.py

@@ -1162,7 +1162,7 @@ def add_terms(self, order, df):
         # Finally store the dataframe
         return self.store.add_table("term" + str(order), df)
 
-    def remove_terms(self, order, index=None, propogate=False):
+    def remove_terms(self, order, index=None, propagate=False):
         """
         Removes terms using a index notation.
 
@@ -1172,7 +1172,7 @@ def remove_terms(self, order, index=None, propogate=False):
             The order (number of atoms) involved in the expression i.e. 2, "two"
         index: list
             The indices of the terms to be removed. If index is None, all terms of that order will be removed.
-        propogate: bool
+        propagate: bool
             This flag indicates if higher order terms should be removed with the removed term.
 
         Returns
@@ -1186,7 +1186,7 @@ def remove_terms(self, order, index=None, propogate=False):
         order_list = [order]
 
         # If this action should be propagated, orders will be added to order list.
-        if propogate == True:
+        if propagate == True:
             order_list.extend([x for x in [3,4] if x > order])
 
         # Figure out atom numbers for this removal.
@@ -1203,14 +1203,14 @@ def remove_terms(self, order, index=None, propogate=False):
             atoms = atom_df[cols]
 
         # Loop through order to be removed
-        for ord in order_list:
+        for current_order in order_list:
 
             # Get column names for order ord
-            cols = metadata.get_term_metadata(ord, "index_columns")
+            cols = metadata.get_term_metadata(current_order, "index_columns")
             cols = [x for x in cols if 'atom' in x]
 
             # Get terms for order ord
-            terms = self.get_terms(ord)
+            terms = self.get_terms(current_order)
 
             # If atoms list is empty, remove_index = None (ie, all removed). Otherwise, only remove interactions
             # for specified atoms.
@@ -1227,23 +1227,23 @@ def remove_terms(self, order, index=None, propogate=False):
                     remove_index.extend(matching_ind)
 
             # Use FL remove function.
-            self.store.remove_table("term" + str(ord), remove_index)
+            self.store.remove_table("term" + str(current_order), remove_index)
 
-            df = self.get_terms(ord)
+            df = self.get_terms(current_order)
 
             # Redo term count
-            self._term_count[ord] = {}
-            self._term_count[ord]["total"] = 0
+            self._term_count[current_order] = {}
+            self._term_count[current_order]["total"] = 0
 
             if not df.empty:
                 uvals, ucnts = np.unique(df["term_index"], return_counts=True)
                 for uval, cnt in zip(uvals, ucnts):
-                    if uval not in self._term_count[ord]:
-                        self._term_count[ord][uval] = cnt
+                    if uval not in self._term_count[current_order]:
+                        self._term_count[current_order][uval] = cnt
                     else:
-                        self._term_count[ord][uval] += cnt
+                        self._term_count[current_order][uval] += cnt
 
-                    self._term_count[ord]["total"] += cnt
+                    self._term_count[current_order]["total"] += cnt
 
         return True
 

eex/filelayer.py

@@ -75,7 +75,11 @@ def add_table(self, key, data):
             do_append = False
             self.created_tables.append(key)
 
-        data.to_hdf(self.store, key, format="t", append=do_append)
+        if do_append:
+            self.store.append(key, data)
+        else:
+            data.to_hdf(self.store, key, format="t", append=do_append)
+
         return True
 
     def read_table(self, key, rows=None, where=None, chunksize=None):
@@ -128,6 +132,11 @@ def remove_table(self, key, index=None):
             for i in ret:
                 self.store.remove(key, start=i[0], stop=i[-1])
 
+            # If everything is removed by index, we should drop this key from the list of created tables
+            if self.read_table(key).empty:
+                self.created_tables.remove(key)
+
+
     def close(self):
         """
         Closes the FL file.

eex/tests/conftest.py

@@ -71,14 +71,23 @@ def build_butane(ff=True, nb=True, scale=True):
         dl.add_dihedrals(dihedral_df)
 
         if ff:
+            # Add harmonic bond
             dl.add_term_parameter(3, "harmonic", {'K': 62.100, 'theta0': 114}, uid=0,
                                   utype={'K': 'kcal * mol ** -1 * radian ** -2',
                                          'theta0': 'degree'})
 
+            # Add harmonic angle
             dl.add_term_parameter(2, "harmonic", {'K': 300.9, 'R0': 1.540}, uid=0,
                                   utype={'K': "kcal * mol **-1 * angstrom ** -2",
                                          'R0': "angstrom"})
 
+            # Add opls dihedral
+            dl.add_term_parameter(4, "opls", {'K_1': 1.41103414, 'K_2': -0.27101489,
+                                              'K_3': 3.14502869, 'K_4': 0}, uid=0, utype={'K_1': 'kcal * mol ** -1',
+                                                                                          'K_2': 'kcal * mol ** -1',
+                                                                                          'K_3': 'kcal * mol ** -1',
+                                                                                          'K_4': 'kcal * mol ** -1'})
+
         if nb:
             dl.add_nb_parameter(atom_type=1, nb_name="LJ",
                                 nb_model="epsilon/sigma", nb_parameters={'sigma': 3.75, 'epsilon': 0.1947460018},

eex/tests/scratch/test_output

@@ -3,56 +3,50 @@ LAMMPS data file generated by MolSSI EEX
  4 atoms
  2 atom types
  3 bonds
- 2 bond types
+ 1 bond types
  2 angles
  1 angle types
- 4 dihedrals
- 4 dihedral types
+ 1 dihedrals
+ 1 dihedral types
  0 impropers
  0 improper types
- 0.000000 100.000000 xlo xhi
- 0.000000 100.000000 ylo yhi
- 0.000000 100.000000 zlo zhi
+-47.690000 52.310000 xlo xhi
+-50.000000 50.000000 ylo yhi
+-50.000000 50.000000 zlo zhi
 
+Pair Coeffs
 
-PairIJ Coeffs
-
- 1  1 0.19474595 3.75000000
- 1  2 0.13342411 3.85000000
- 2  2 0.09141137 3.95000000
+ 1 0.09141129 3.95000000
+ 2 0.19474600 3.75000000
 
 Bond Coeffs
 
  1 300.90000000 1.54000000
- 2 300.90000000 1.54000000
 
 Angle Coeffs
 
- 1 62.10011250 114.00004886
+ 1 62.10011249 114.00000000
 
 Dihedral Coeffs
 
- 1 0.00000000 0.00000000 0.00000000
- 2 0.70551689 1.00000000 0.00000000
- 3 -0.13550741 2.00000000 -180.00007714
- 4 1.57251395 3.00000000 0.00000000
+ 1 1.41103414 -0.27101489 3.14502869 0.00000000
 
  Masses
 
-1 15.0452
-2 14.02658
+2 15.03452
+1 14.02658
 
  Atoms
 
-1  1  1 0.00000000  0.00000000 -0.45970000 -1.53020000
-2  1  2 0.00000000  0.00000000  0.00000000  0.00000000
-3  1  2 0.00000000  0.00000000  1.59800000  0.00000000
-4  1  1 0.00000000 -1.47400000  1.57300000 -0.61670000
+1  1  2 0.00000000  0.00000000 -0.45970000 -1.53020000
+2  1  1 0.00000000  0.00000000  0.00000000  0.00000000
+3  1  1 0.00000000  0.00000000  1.59800000  0.00000000
+4  1  2 0.00000000 -1.47400000  1.57300000 -0.61670000
 
  Bonds
 
 1 1 1 2
-2 2 2 3
+2 1 2 3
 3 1 3 4
 
  Angles
@@ -63,7 +57,4 @@ Dihedral Coeffs
  Dihedrals
 
 1 1 1 2 3 4
-2 2 1 2 3 4
-3 3 1 2 3 4
-4 4 1 2 3 4
 

eex/tests/test_amber.py

@@ -206,6 +206,21 @@ def test_amber_compatibility_no_scaling(butane_dl):
         eex.translators.amber.write_amber_file(dl, oname)
 
 
+def test_amber_compatibility_dihedral_conversion(butane_dl):
+
+    dl = butane_dl()
+
+    oname = eex_find_files.get_scratch_directory("dl_compatibility.prmtop")
+
+    eex.translators.amber.amber_write._check_dl_compatibility(dl)
+
+    print(dl.get_term_count(4))
+
+    print(dl.get_terms(4))
+
+    # eex.translators.amber.write_amber_file(dl, oname)
+
+
 
 
 

eex/tests/test_compare_alkanes.py

@@ -74,7 +74,7 @@ def test_alkane(lammps_bench, program):
 
 # write test that loads in amber --> lammps --> EEX compare datalayers
 
-@pytest.mark.parametrize("program1", ["amber"]) #build dl
+@pytest.mark.parametrize("program1", ["amber", "lammps"]) #build dl
 @pytest.mark.parametrize("program2", ["amber", "lammps"]) #write dl
 @pytest.mark.parametrize("molecule", _alkane_molecules)
 def test_translation(program1, program2, molecule):

eex/tests/test_datalayer.py

@@ -877,7 +877,7 @@ def test_remove_terms_by_index(butane_dl):
     # Check that bonds are what we expect
     assert(sorted(bonds1.loc[2].values) == sorted(bonds2.values[0]))
 
-    # Here, since propogate was set to false. Angles and dihedrals remain unchanged.
+    # Here, since propagate was set to false. Angles and dihedrals remain unchanged.
     assert(dl.get_term_count(3)['total'] == 2)
 
     assert (dl.get_term_count(4)['total'] == 1)
@@ -902,7 +902,7 @@ def test_remove_terms_by_index_nonconsecutive(butane_dl):
     # Check that bonds are what we expect
     assert(sorted(bonds1.loc[1].values) == sorted(bonds2.values[0]))
 
-    # Here, since propogate was set to false. Angles and dihedrals remain unchanged.
+    # Here, since propagate was set to false. Angles and dihedrals remain unchanged.
     assert(dl.get_term_count(3)['total'] == 2)
 
     assert (dl.get_term_count(4)['total'] == 1)
@@ -921,7 +921,7 @@ def test_remove_terms_propagate(butane_dl):
     assert(dl.get_term_count(3)['total'] == 2)
     assert (dl.get_term_count(4)['total'] == 1)
 
-    dl.remove_terms(2, propogate=True)
+    dl.remove_terms(2, propagate=True)
 
     # Assert all have been removed.
     bonds = dl.get_terms(2)
@@ -933,7 +933,7 @@ def test_remove_terms_propagate(butane_dl):
 
     return True
 
-def test_remove_terms_by_index_propogate(butane_dl):
+def test_remove_terms_by_index_propagate(butane_dl):
 
     dl = butane_dl()
 
@@ -947,16 +947,16 @@ def test_remove_terms_by_index_propogate(butane_dl):
 
     assert(not bonds.empty)
 
-    # Remove one bond - choose to propogate this so dihedral and angle should also be removed.
-    dl.remove_terms(2, index=[0], propogate=True)
+    # Remove one bond - choose to propagate this so dihedral and angle should also be removed.
+    dl.remove_terms(2, index=[0], propagate=True)
 
     assert(dl.get_term_count(2)['total'] == 2)
     assert (dl.get_term_count(3)['total'] == 1)
     assert(dl.get_term_count(4)['total'] == 0)
 
     return True
 
-def test_remove_terms_by_index_nonconsecutive_propogate(butane_dl):
+def test_remove_terms_by_index_nonconsecutive_propagate(butane_dl):
     dl = butane_dl()
 
     bonds = dl.get_terms(2)
@@ -970,11 +970,51 @@ def test_remove_terms_by_index_nonconsecutive_propogate(butane_dl):
     assert (dl.get_term_count(3)['total'] == 2)
     assert (dl.get_term_count(4)['total'] == 1)
 
-    # Remove two bonds - choose to propogate this so dihedral and both angles should also be removed.
-    dl.remove_terms(2, index=[0, 2], propogate=True)
+    # Remove two bonds - choose to propagate this so dihedral and both angles should also be removed.
+    dl.remove_terms(2, index=[0, 2], propagate=True)
 
     assert(dl.get_term_count(2)['total'] == 1)
     assert (dl.get_term_count(3)['total'] == 0)
     assert(dl.get_term_count(4)['total'] == 0)
 
     return True
+
+def test_remove_and_add_terms(butane_dl):
+    dl = butane_dl()
+
+    bonds = dl.get_terms(2)
+
+    assert(not bonds.empty)
+
+    dl.remove_terms(2)
+
+    bonds_new = dl.get_terms(2)
+
+    assert(bonds_new.empty)
+    assert(dl.get_term_count(2)['total'] == 0)
+
+    dl.add_bonds(bonds)
+
+    added_bonds = dl.get_terms(2)
+
+    assert (not added_bonds.empty)
+
+def test_remove_and_add_terms_index(butane_dl):
+    dl = butane_dl()
+
+    bonds = dl.get_terms(2)
+
+    assert(not bonds.empty)
+
+    dl.remove_terms(2, index=[0])
+
+    assert(dl.get_term_count(2)['total'] == 2)
+
+    dl.add_bonds(bonds.iloc[[0]])
+
+    added_bonds = dl.get_terms(2)
+
+    assert (not added_bonds.empty)
+
+    for col in added_bonds.columns:
+        assert(set(added_bonds[col].values) == set(bonds[col].values))