Modify scaling factors to a numpy array (#669)

* Modify scaling factors to a numpy array

* Misc Changes

1. Split property into two parts. global_scaling_factors and per_molecule_scaling_factors
2. Handle Nan/NoneType cases, non float cases
3. Add additional functions for handling scaling factors
4. Change formats/json.py to incorporate new properties
5. Change function name to prefix `get_`, formats/mcf.py change

* WIP- Refactor global_scaling_factors to scaling_factors

* Additional tests; Simplificaiton of get_scaling_factors method

* Atomtype parameterization methods for parameterizing a topology in GMSO. (#644)

* Atomtype parameterization methods for parameterizing a topology in GMSO.

This PR will add the atomtyping module to GMSO for passing a GMSO.Forcefield
and a GMSO.Topology and match atomtype using foyer as the backend. Then,
the corresponding connection types will be found in the Forcefield and
applied to the connections in the topology.

* Create parameterize.py, which has the apply function which can take a
topology, and a gmso forcefield to apply to it. This can use subgraph
isomorphism to identify molecular structures in the topology through
the bondgraph and bin those into unique molecules that match a specified
forcefield. This apply function can also do the standard atomtyping of
the entire topology in one step.
* Create isomorph.py which uses networkx graphs to identify disconnected
components and isomorphism to identify repeated structures.

* Move module imports for apply into atomtyping to prevent circular imports

* Add a quick fix which will do atomtyping if no residue flag is passed

* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci

* correctly update sites when adding subtop

* Changes to doc strings for clarity. Add a subtop_label variable to generalize where the molecule definition is pulled from

* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci

* More modular architecture

* WIP- Add testing, minor refactoring of the APIs.

* WIP- Better error handling

* Misc Changes

1. Update Topology after parametrization
2. Add dependency for forcefield_utilities in env files
3. Add tests for trappe forcefield
4. Patch parmed (should be moved to a cal's PR #658

* WIP- Skip opls for now, full TrappE tests

* Avoid accidental overwriting of typemap when using isomorphism

* WIP- Remove unused import

* Use enumerate for atom index while converting to TopologyGraph

* Fix argument order

* WIP- Add test for subtopology parameterization

* Make opls/trappe global fixtures, Add tests for isomorphism

* Further testing isomorphism

* REVERT - skip OPLS tests

* Copy scaling factors and combining rules after parametrization

* Proper OPLS tests

* WIP- Refactor the test module

* WIP- Remove unused import

* WIP- Add test for parameterization with impropers

* WIP- Additional impropers test; Separate module for testing impropers

* Minor refacotors; additional edge cases coverage/tests

* Docstring minor fix

* Remove rel_to_module as is obsolete in forcefield_utilities

* Change trappe_ua to trappe-ua for correct loading

* fix typo, add note about specific use case

* pip install forcefield-utilites until new release

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Co Quach <[email protected]>
Co-authored-by: Co Quach <[email protected]>
Co-authored-by: Umesh Timalsina <[email protected]>

* Properly apply different scaling factors while parameterizing

Co-authored-by: CalCraven <[email protected]>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Co Quach <[email protected]>
Co-authored-by: Co Quach <[email protected]>

gmso/core/topology.py

@@ -25,6 +25,8 @@
     identify_connections as _identify_connections,
 )
 
+scaling_interaction_idxes = {"12": 0, "13": 1, "14": 2}
+
 
 class Topology(object):
     """A topology.
@@ -153,14 +155,14 @@ def __init__(self, name="Topology", box=None):
         self._subtops = IndexedSet()
         self._combining_rule = "lorentz"
         self._pairpotential_types = IndexedSet()
-        self._scaling_factors = {
-            "nonBonded12Scale": 0.0,
-            "nonBonded13Scale": 0.0,
-            "nonBonded14Scale": 0.5,
-            "electrostatics12Scale": 0.0,
-            "electrostatics13Scale": 0.0,
-            "electrostatics14Scale": 0.5,
-        }
+        self._scaling_factors = np.array(
+            [
+                [0.0, 0.0, 0.5],  # lj scales
+                [0.0, 0.0, 0.5],  # electrostatics scale
+            ],
+            dtype=float,
+        )
+        self._molecule_scaling_factors = {}
         self.is_updated = True
         self._potentials_count = {
             "atom_types": 0,
@@ -217,32 +219,11 @@ def combining_rule(self, rule):
 
     @property
     def scaling_factors(self):
-        """Return the scaling factors for the topology."""
-        return self._scaling_factors
-
-    @scaling_factors.setter
-    def scaling_factors(self, scaling_factors):
-        """Set the scaling factors for the topology."""
-        expected_items = [
-            "nonBonded12Scale",
-            "nonBonded13Scale",
-            "nonBonded14Scale",
-            "electrostatics12Scale",
-            "electrostatics13Scale",
-            "electrostatics14Scale",
-        ]
-        if not isinstance(scaling_factors, dict):
-            raise GMSOError("Scaling factors should be a dictionary")
-        for item in expected_items:
-            if item not in scaling_factors.keys():
-                raise GMSOError(
-                    f"Expected {expected_items} as keys in the scaling factors"
-                )
-        for val in scaling_factors.values():
-            if val < 0.0 or val > 1.0:
-                raise GMSOError("Scaling factors should be between 0.0 and 1.0")
+        return self._scaling_factors.copy()
 
-        self._scaling_factors = scaling_factors
+    @property
+    def molecule_scaling_factors(self):
+        return {k: v.copy() for k, v in self._molecule_scaling_factors.items()}
 
     @property
     def positions(self):
@@ -637,6 +618,137 @@ def pairpotential_type_expressions(self):
             set([ptype.expression for ptype in self._pairpotential_types])
         )
 
+    def get_lj_scale(self, *, molecule_id=None, interaction=None):
+        """Return the selected lj_scales defined for this topology."""
+        return self._get_scaling_factor(molecule_id, interaction, "lj_scale", 0)
+
+    def set_lj_scale(self, value, *, molecule_id=None, interaction=None):
+        """Set the correct lj_scaling factors for this topology."""
+        self._set_scaling_factor(value, molecule_id, interaction, "lj_scale", 0)
+
+    def get_scaling_factors(self, *, molecule_id=None):
+        """Get the scaling factor of this topology or for a particular molecule"""
+        return np.vstack(
+            [
+                self.get_lj_scale(molecule_id=molecule_id),
+                self.get_electrostatics_scale(molecule_id=molecule_id),
+            ]
+        )
+
+    def set_scaling_factors(self, lj, electrostatics, *, molecule_id=None):
+        """Set both lj and electrostatics scaling factors."""
+        self.set_lj_scale(
+            lj,
+            molecule_id=molecule_id,
+            interaction=None,
+        )
+
+        self.set_electrostatics_scale(
+            electrostatics,
+            molecule_id=molecule_id,
+        )
+
+    def get_electrostatics_scale(self, *, molecule_id=None, interaction=None):
+        """Return the selected electrostatics_scale defined for this topology.
+
+        Parameters
+        ----------
+        molecule_id: str, default=None
+            The molecule id that this scaling factor applies to, if None
+            this will return the Topology's global scaling factors
+
+        interaction: str, one of {'12', '13', '14'}, default=None
+            The interaction for which to return the scaling factor for, if None
+            a 3 tuple
+
+        Raises
+        ------
+        GMSOError
+            If the specified parameters can't return a scaling factor
+        """
+        return self._get_scaling_factor(
+            molecule_id, interaction, "electrostatics_scale", 1
+        )
+
+    def set_electrostatics_scale(
+        self, value, *, molecule_id=None, interaction=None
+    ):
+        """Set the correct lj_scaling factors for this topology.
+
+        Parameters
+        ----------
+        value: float, numpy.ndarray, list, or tuple of floats
+            The value to set for this scale
+
+        molecule_id: str, default=None
+            The molecule id that this scaling factor applies to, if None
+            this will return the Topology's global scaling factors
+
+        interaction: str, one of {'12', '13', '14'}, default=None
+            The interaction for which to return the scaling factor for, if None
+            a 3 tuple
+
+        Raises
+        ------
+        GMSOError
+            If the specified parameters can't return a scaling factor
+        """
+        self._set_scaling_factor(
+            value, molecule_id, interaction, "electrostatics_scale", 1
+        )
+
+    def _get_scaling_factor(self, molecule_id, interaction, name, index):
+        """Get the scaling factor according to molecule_id, interaction, and name."""
+        if molecule_id is None:
+            all_scales = self._scaling_factors
+        else:
+            if molecule_id not in self._molecule_scaling_factors:
+                raise GMSOError(
+                    f"Scaling factors for molecule `{molecule_id}` is not defined "
+                    f"in the topology. Please use appropriate molecule_id"
+                )
+            all_scales = self._molecule_scaling_factors[molecule_id]
+
+        if interaction is None:
+            return all_scales[index].copy()
+        else:
+            if interaction not in scaling_interaction_idxes:
+                raise GMSOError(f"Unknown `{name}` interaction `{interaction}`")
+            return all_scales[index][scaling_interaction_idxes[interaction]]
+
+    def _set_scaling_factor(self, value, molecule_id, interaction, name, index):
+        """Set the scaling factor according to molecule_id, interaction, and name."""
+        org_value = value
+        value = np.array(value, dtype=float).reshape(-1)
+
+        if any(np.isnan(value)):
+            raise ValueError(
+                f"Cannot assign a nan/NoneType to `{name}`. "
+                f"Provided value: {org_value}"
+            )
+
+        if value.shape != (1,) and value.shape != (3,):
+            raise ValueError(
+                f"Cannot determine the appropriate shape for {org_value} to "
+                f"assign it to `{name}`"
+            )
+
+        if molecule_id is None:
+            all_scales = self._scaling_factors
+        else:
+            if molecule_id not in self._molecule_scaling_factors:
+                self._molecule_scaling_factors[
+                    molecule_id
+                ] = self._scaling_factors.copy()
+            all_scales = self._molecule_scaling_factors[molecule_id]
+
+        if interaction is None:
+            all_scales[index] = value
+        else:
+            if interaction not in scaling_interaction_idxes:
+                raise GMSOError(f"Unknown `{name}` interaction `{interaction}`")
+            all_scales[index][scaling_interaction_idxes[interaction]] = value
+
     def add_site(self, site, update_types=False):
         """Add a site to the topology.
 

gmso/formats/json.py

@@ -63,7 +63,12 @@ def _to_json(top, types=False, update=True):
 
     json_dict = {
         "name": top._name,
-        "scaling_factors": top.scaling_factors,
+        "scaling_factors": {
+            "scaling_factors": top.scaling_factors.tolist(),
+            "molecule_scaling_factors": {
+                k: v.tolist() for k, v in top.molecule_scaling_factors.items()
+            },
+        },
         "subtopologies": [],
         "box": top.box.json_dict() if top.box else None,
         "atoms": [],
@@ -145,6 +150,14 @@ def _to_json(top, types=False, update=True):
     return json_dict
 
 
+def _set_scaling_factors(top, scaling_factors):
+    """Set the global/permolecule scaling factors."""
+    global_scaling_factor = scaling_factors["scaling_factors"]
+    top.set_scaling_factors(global_scaling_factor[0], global_scaling_factor[1])
+    for k, v in scaling_factors["molecule_scaling_factors"].items():
+        top.set_scaling_factors(v[0], v[1], molecule_id=k)
+
+
 def _from_json(json_dict):
     """Convert a json_dict into a topology.
 
@@ -170,7 +183,7 @@ def _from_json(json_dict):
     top = Topology(
         name=json_dict["name"],
     )
-    top.scaling_factors = json_dict["scaling_factors"]
+    _set_scaling_factors(top, json_dict["scaling_factors"])
     id_to_type_map = {}
     for atom_dict in json_dict["atoms"]:
         atom_type_id = atom_dict.pop("atom_type", None)

gmso/formats/mcf.py

@@ -625,7 +625,8 @@ def _write_intrascaling_information(mcf, top):
         The 1-4 scaling parameter for Coulombic interactions
 
     """
-    sf = top.scaling_factors
+    nbonded_sf = top.get_lj_scale()
+    electstatic_sf = top.get_electrostatics_scale()
     header = (
         "\n!Intra Scaling\n"
         "!vdw_scaling    1-2 1-3 1-4 1-N\n"
@@ -634,20 +635,8 @@ def _write_intrascaling_information(mcf, top):
     )
 
     mcf.write(header)
-    mcf.write(
-        "{:.4f} {:.4f} {:.4f} 1.0000\n".format(
-            sf["nonBonded12Scale"],
-            sf["nonBonded13Scale"],
-            sf["nonBonded14Scale"],
-        )
-    )
-    mcf.write(
-        "{:.4f} {:.4f} {:.4f} 1.0000\n".format(
-            sf["electrostatics12Scale"],
-            sf["electrostatics13Scale"],
-            sf["electrostatics14Scale"],
-        )
-    )
+    mcf.write("{:.4f} {:.4f} {:.4f} 1.0000\n".format(*nbonded_sf))
+    mcf.write("{:.4f} {:.4f} {:.4f} 1.0000\n".format(*electstatic_sf))
 
 
 def _check_compatibility(top):

gmso/parameterization/topology_parameterizer.py

@@ -181,41 +181,58 @@ def _parameterize(self, subtop_or_top, typemap, is_subtop=False):
             subtop_or_top, forcefield, is_subtop=is_subtop
         )
 
-    def _verify_forcefields_metadata(self):
-        """Verify all the provided forcefields have the same scaling factors and combining rule."""
+    def _set_combining_rule(self):
+        """Verify all the provided forcefields have the same combining rule and set it for the Topology."""
         if isinstance(self.forcefields, dict):
-            ffs = list(self.forcefields.values())
-            init_scaling_factors = ffs[0].scaling_factors
-            init_combining_rule = ffs[0].combining_rule
-            for ff in ffs[1:]:
-                if ff.scaling_factors != init_scaling_factors:
-                    raise ParameterizationError(
-                        "Scaling factors of the provided forcefields do not"
-                        "match, please provide forcefields with same scaling"
-                        "factors that apply to a Topology"
-                    )
+            all_comb_rules = set(
+                ff.combining_rule for ff in self.forcefields.values()
+            )
+        else:
+            all_comb_rules = {self.forcefields.combining_rule}
 
-                if ff.combining_rule != init_combining_rule:
-                    raise ParameterizationError(
-                        "Combining rules of the provided forcefields do not"
-                        "match, please provide forcefields with same scaling"
-                        "factors that apply to a Topology"
-                    )
-            return init_scaling_factors, init_combining_rule
+        if not len(all_comb_rules) == 1:
+            raise ParameterizationError(
+                "Combining rules of the provided forcefields do not"
+                "match, please provide forcefields with same scaling"
+                "factors that apply to a Topology"
+            )
+        self.topology.combining_rule = all_comb_rules.pop()
+
+    def _set_scaling_factors(self):
+        """Set either per-molecule or global scaling factors for the topology based on the forcefields provided."""
+        # ToDo: Set other scaling factors by extending the forcefield schema
+        # ToDo: What to do when all the scaling factors matchup? Should we promote them to be global?
+        if isinstance(self.forcefields, Dict):
+            for subtop_id, ff in self.forcefields.items():
+                self.topology.set_lj_scale(
+                    ff.scaling_factors["nonBonded14Scale"],
+                    interaction="14",
+                    molecule_id=subtop_id,
+                )
+                self.topology.set_electrostatics_scale(
+                    ff.scaling_factors["electrostatics14Scale"],
+                    interaction="14",
+                    molecule_id=subtop_id,
+                )
         else:
-            return (
-                self.forcefields.scaling_factors,
-                self.forcefields.combining_rule,
+            self.topology.set_lj_scale(
+                self.forcefields.scaling_factors["nonBonded14Scale"],
+                interaction="14",
+            )
+            self.topology.set_electrostatics_scale(
+                self.forcefields.scaling_factors["electrostatics14Scale"],
+                interaction="14",
             )
 
     def run_parameterization(self):
         """Run parameterization of the topology with give forcefield(s) and configuration."""
-        scaling_factors, combining_rule = self._verify_forcefields_metadata()
         if self.topology.is_typed():
             raise ParameterizationError(
                 "Cannot parameterize a typed topology. Please provide a topology without any types"
             )
 
+        self._set_combining_rule()  # Fail Early if no match
+
         if self.config.identify_connections:
             """ToDo: This mutates the topology and is agnostic to downstream
             errors. So, here we should use index only option"""
@@ -262,8 +279,7 @@ def run_parameterization(self):
                 is_subtop=False,  # This will be removed from the future iterations
             )
 
-        self.topology.scaling_factors.update(scaling_factors)
-        self.topology.combining_rule = combining_rule
+        self._set_scaling_factors()  # Set global or per molecule scaling factors
         self.topology.update_topology()
 
     @staticmethod

gmso/tests/base_test.py

@@ -414,10 +414,13 @@ def test_topology_equivalence(top1, top2):
                 return False, "Unequal number of impropers"
             if top1.name != top2.name:
                 return False, "Dissimilar names"
-
-            if top1.scaling_factors != top2.scaling_factors:
-                return False, f"Mismatch in scaling factors"
-
+            if not np.allclose(top1.scaling_factors, top2.scaling_factors):
+                return False, "Mismatch in scaling factors"
+            for k, v in top1.molecule_scaling_factors.items():
+                if k not in top2.scaling_factors:
+                    return False, "Mismatch in scaling factors"
+                elif not np.allclose(v, top2.molecule_scaling_factors[k]):
+                    return False, "Mismatch in scaling factors"
             if not have_equivalent_boxes(top1, top2):
                 return (
                     False,