Merge pull request #45 from openforcefield/parmed-2

Update ParmEd converter for new components

.github/workflows/integration_tests.yaml

@@ -25,7 +25,7 @@ jobs:
     steps:
     - uses: actions/checkout@v2
 
-    - uses: goanpeca/setup-miniconda@v1
+    - uses: conda-incubator/setup-miniconda@v2
       with:
         python-version: ${{ matrix.python-version }}
         activate-environment: test

openff/system/components/potentials.py

@@ -10,7 +10,7 @@
 class Potential(BaseModel):
     """Base class for storing applied parameters"""
 
-    parameters: Dict[str, Optional[unit.Quantity]] = dict()
+    parameters: Dict[str, Optional[Union[unit.Quantity, List, int]]] = dict()
 
     class Config:
         arbitrary_types_allowed = True

openff/system/components/smirnoff.py

@@ -126,11 +126,13 @@ def store_potentials(self, parameter_handler: ProperTorsionHandler) -> None:
             # ParameterHandler.get_parameter returns a list, although this
             # should only ever be length 1
             parameter_type = parameter_handler.get_parameter({"smirks": smirks})[0]
+            n_terms = len(parameter_type.k)
             potential = Potential(
                 parameters={
-                    "k": simtk_to_pint(parameter_type.k),
-                    "periodicity": simtk_to_pint(parameter_type.periodicity),
-                    "phase": simtk_to_pint(parameter_type.phase),
+                    "k": [simtk_to_pint(val) for val in parameter_type.k],
+                    "periodicity": parameter_type.periodicity,
+                    "phase": [simtk_to_pint(val) for val in parameter_type.phase],
+                    "n_terms": n_terms,
                 },
             )
             self.potentials[smirks] = potential
@@ -143,6 +145,9 @@ class SMIRNOFFvdWHandler(PotentialHandler):
     independent_variables: Set[str] = {"r"}
     slot_map: Dict[str, str] = dict()
     potentials: Dict[str, Potential] = dict()
+    scale_13: float = 0.0
+    scale_14: float = 0.5
+    scale_15: float = 1.0
 
     def store_matches(
         self,
@@ -191,6 +196,9 @@ class SMIRNOFFElectrostaticsHandler(BaseModel):
     expression: str = "coul"
     independent_variables: Set[str] = {"r"}
     charge_map: Dict[tuple, unit.Quantity] = dict()
+    scale_13: float = 0.0
+    scale_14: float = 0.8333333333
+    scale_15: float = 1.0
 
     def store_charges(
         self,

openff/system/components/system.py

@@ -1,8 +1,14 @@
-from typing import Dict
+from typing import Dict, Optional
 
-from pydantic import BaseModel
+import numpy as np
+from openforcefield.topology.topology import Topology
+from pydantic import BaseModel, validator
+from simtk import unit as omm_unit
 
 from openff.system.components.potentials import PotentialHandler
+from openff.system.interop.parmed import to_parmed
+from openff.system.types import UnitArray
+from openff.system.utils import simtk_to_pint
 
 
 class System(BaseModel):
@@ -13,7 +19,26 @@ class System(BaseModel):
     """
 
     handlers: Dict[str, PotentialHandler] = dict()
+    topology: Optional[Topology] = None
+    box: Optional[UnitArray] = None
+    positions: Optional[UnitArray] = None
 
     class Config:
         arbitrary_types_allowed = True
-        validate_assignment = True
+
+    @validator("box")
+    def validate_box(cls, val):
+        if val is None:
+            return val
+        if type(val) == omm_unit.Quantity:
+            val = simtk_to_pint(val)
+        if val.shape == (3, 3):
+            return val
+        elif val.shape == (3,):
+            val = val * np.eye(3)
+            return val
+        else:
+            raise ValueError  # InvalidBoxError
+
+
+System.to_parmed = to_parmed

openff/system/exceptions.py

@@ -49,3 +49,9 @@ class InvalidBoxError(TypeError):
     """
     Generic exception for errors reading box data
     """
+
+
+class InterMolEnergyComparisonError(AssertionError):
+    """
+    Exception for when energies derived from InterMol do not match
+    """

openff/system/interop/parmed.py

@@ -1,47 +1,72 @@
-import re
 from typing import Any
 
 import numpy as np
 import parmed as pmd
 
-from .. import unit
+from openff.system import unit
+
+kcal_mol = unit.Unit("kilocalories / mol")
 
 
 def to_parmed(off_system: Any) -> pmd.Structure:
     """Convert an OpenFF System to a ParmEd Structure"""
     structure = pmd.Structure()
-    _convert_box(off_system, structure)
+    _convert_box(off_system.box, structure)
+
+    if "Electrostatics" in off_system.handlers.keys():
+        has_electrostatics = True
+        electrostatics_handler = off_system.handlers["Electrostatics"]
+    else:
+        has_electrostatics = False
 
     for topology_molecule in off_system.topology.topology_molecules:
         for atom in topology_molecule.atoms:
+            atomic_number = atom.atomic_number
+            element = pmd.periodic_table.Element[atomic_number]
+            mass = pmd.periodic_table.Mass[element]
             structure.add_atom(
-                pmd.Atom(atomic_number=atom.atomic_number), resname="FOO", resnum=0
+                pmd.Atom(
+                    atomic_number=atomic_number,
+                    mass=mass,
+                ),
+                resname="FOO",
+                resnum=0,
             )
 
-    if "Bonds" in off_system.term_collection.terms:
-        bond_term = off_system.term_collection.terms["Bonds"]
-        for bond, smirks in bond_term.smirks_map.items():
-            idx_1, idx_2 = bond
-            pot = bond_term.potentials[bond_term.smirks_map[bond]]
-            k = pot.parameters["k"].m / 2
-            length = pot.parameters["length"].m
-            bond_type = pmd.BondType(k=k, req=length)
+    if "Bonds" in off_system.handlers.keys():
+        bond_handler = off_system.handlers["Bonds"]
+        bond_map = dict()
+        for bond_slot, smirks in bond_handler.slot_map.items():
+            idx_1, idx_2 = eval(bond_slot)
+            try:
+                bond_type = bond_map[smirks]
+            except KeyError:
+                pot = bond_handler.potentials[smirks]
+                k = pot.parameters["k"].m / 2
+                length = pot.parameters["length"].m
+                bond_type = pmd.BondType(k=k, req=length)
+                bond_map[smirks] = bond_type
+                del pot, k, length
+            if bond_type not in structure.bond_types:
+                structure.bond_types.append(bond_type)
             structure.bonds.append(
                 pmd.Bond(
                     atom1=structure.atoms[idx_1],
                     atom2=structure.atoms[idx_2],
                     type=bond_type,
                 )
             )
+            del bond_type, idx_1, idx_2, smirks, bond_slot
 
-    if "Angles" in off_system.term_collection.terms:
-        angle_term = off_system.term_collection.terms["Angles"]
-        for angle, smirks in angle_term.smirks_map.items():
-            idx_1, idx_2, idx_3 = angle
-            pot = angle_term.potentials[angle_term.smirks_map[angle]]
+    if "Angles" in off_system.handlers.keys():
+        angle_term = off_system.handlers["Angles"]
+        for angle, smirks in angle_term.slot_map.items():
+            idx_1, idx_2, idx_3 = eval(angle)
+            pot = angle_term.potentials[smirks]
             # TODO: Look at cost of redundant conversions, to ensure correct units of .m
-            k = pot.parameters["k"].magnitude  # kcal/mol/rad**2
+            k = pot.parameters["k"].magnitude / 2  # kcal/mol/rad**2
             theta = pot.parameters["angle"].magnitude  # degree
+            # TODO: Look up if AngleType already exists in struct
             angle_type = pmd.AngleType(k=k, theteq=theta)
             structure.angles.append(
                 pmd.Angle(
@@ -51,66 +76,110 @@ def to_parmed(off_system: Any) -> pmd.Structure:
                     type=angle_type,
                 )
             )
-
-    if "ProperTorsions" in off_system.term_collection.terms:
-        proper_term = off_system.term_collection.terms["ProperTorsions"]
-        for proper, smirks in proper_term.smirks_map.items():
-            idx_1, idx_2, idx_3, idx_4 = proper
-            pot = proper_term.potentials[proper_term.smirks_map[proper]]
-            # TODO: Better way of storing periodic data in generally, probably need to improve Potential
-            n = re.search(r"\d", "".join(pot.parameters.keys())).group()
-            k = pot.parameters["k" + n].m  # kcal/mol
-            periodicity = pot.parameters["periodicity" + n].m  # dimless
-            phase = pot.parameters["phase" + n].m  # degree
-
-            dihedral_type = pmd.DihedralType(per=periodicity, phi_k=k, phase=phase)
-            structure.dihedrals.append(
-                pmd.Dihedral(
-                    atom1=structure.atoms[idx_1],
-                    atom2=structure.atoms[idx_2],
-                    atom3=structure.atoms[idx_3],
-                    atom4=structure.atoms[idx_4],
-                    type=dihedral_type,
+            structure.angle_types.append(angle_type)
+
+    # ParmEd treats 1-4 scaling factors at the level of each DihedralType,
+    # whereas SMIRNOFF captures them at the level of the non-bonded handler,
+    # so they need to be stored here for processing dihedrals
+    vdw_14 = off_system.handlers["vdW"].scale_14
+    if has_electrostatics:
+        coul_14 = off_system.handlers["Electrostatics"].scale_14
+    else:
+        coul_14 = 1.0
+    vdw_handler = off_system.handlers["vdW"]
+    if "ProperTorsions" in off_system.handlers.keys():
+        proper_term = off_system.handlers["ProperTorsions"]
+        for proper, smirks in proper_term.slot_map.items():
+            idx_1, idx_2, idx_3, idx_4 = eval(proper)
+            pot = proper_term.potentials[smirks]
+            for n in range(pot.parameters["n_terms"]):
+                k = pot.parameters["k"][n].to(kcal_mol).magnitude
+                periodicity = pot.parameters["periodicity"][n]
+                phase = pot.parameters["phase"][n].magnitude
+                dihedral_type = pmd.DihedralType(
+                    phi_k=k,
+                    per=periodicity,
+                    phase=phase,
+                    scnb=1 / vdw_14,
+                    scee=1 / coul_14,
                 )
-            )
-
-    if "ImroperTorsions" in off_system.term_collection.terms:
-        improper_term = off_system.term_collection.terms["ImproperTorsions"]
-        for improper, smirks in improper_term.smirks_map.items():
-            idx_1, idx_2, idx_3, idx_4 = improper
-            pot = improper_term.potentials[improper_term.smirks_map[improper]]
-            # TODO: Better way of storing periodic data in generally, probably need to improve Potential
-            n = re.search(r"\d", "".join(pot.parameters.keys())).group()
-            k = pot.parameters["k" + n].m  # kcal/mol
-            periodicity = pot.parameters["periodicity" + n].m  # dimless
-            phase = pot.parameters["phase" + n].m  # degree
-
-            dihedral_type = pmd.DihedralType(per=periodicity, phi_k=k, phase=phase)
-            structure.dihedrals.append(
-                pmd.Dihedral(
-                    atom1=structure.atoms[idx_1],
-                    atom2=structure.atoms[idx_2],
-                    atom3=structure.atoms[idx_3],
-                    atom4=structure.atoms[idx_4],
-                    type=dihedral_type,
+                structure.dihedrals.append(
+                    pmd.Dihedral(
+                        atom1=structure.atoms[idx_1],
+                        atom2=structure.atoms[idx_2],
+                        atom3=structure.atoms[idx_3],
+                        atom4=structure.atoms[idx_4],
+                        type=dihedral_type,
+                    )
                 )
-            )
-
-    vdw_term = off_system.term_collection.terms["vdW"]
+                structure.dihedral_types.append(dihedral_type)
+                vdw1 = vdw_handler.potentials[vdw_handler.slot_map[str((idx_1,))]]
+                vdw4 = vdw_handler.potentials[vdw_handler.slot_map[str((idx_4,))]]
+                sig1, eps1 = _lj_params_from_potential(vdw1)
+                sig4, eps4 = _lj_params_from_potential(vdw4)
+                sig = (sig1 + sig4) * 0.5
+                eps = (eps1 * eps4) ** 0.5
+                nbtype = pmd.NonbondedExceptionType(
+                    rmin=sig * 2 ** (1 / 6), epsilon=eps * vdw_14, chgscale=coul_14
+                )
+                structure.adjusts.append(
+                    pmd.NonbondedException(
+                        structure.atoms[idx_1], structure.atoms[idx_4], type=nbtype
+                    )
+                )
+                structure.adjust_types.append(nbtype)
+
+    #    if False:  # "ImroperTorsions" in off_system.term_collection.terms:
+    #        improper_term = off_system.term_collection.terms["ImproperTorsions"]
+    #        for improper, smirks in improper_term.smirks_map.items():
+    #            idx_1, idx_2, idx_3, idx_4 = improper
+    #            pot = improper_term.potentials[improper_term.smirks_map[improper]]
+    #            # TODO: Better way of storing periodic data in generally, probably need to improve Potential
+    #            n = re.search(r"\d", "".join(pot.parameters.keys())).group()
+    #            k = pot.parameters["k" + n].m  # kcal/mol
+    #            periodicity = pot.parameters["periodicity" + n].m  # dimless
+    #            phase = pot.parameters["phase" + n].m  # degree
+    #
+    #            dihedral_type = pmd.DihedralType(per=periodicity, phi_k=k, phase=phase)
+    #            structure.dihedrals.append(
+    #                pmd.Dihedral(
+    #                    atom1=structure.atoms[idx_1],
+    #                    atom2=structure.atoms[idx_2],
+    #                    atom3=structure.atoms[idx_3],
+    #                    atom4=structure.atoms[idx_4],
+    #                    type=dihedral_type,
+    #                )
+    #            )
+
+    vdw_handler = off_system.handlers["vdW"]
     for pmd_idx, pmd_atom in enumerate(structure.atoms):
-        potential = vdw_term.potentials[vdw_term.smirks_map[(pmd_idx,)]]
+        smirks = vdw_handler.slot_map[str((pmd_idx,))]
+        potential = vdw_handler.potentials[smirks]
+        element = pmd.periodic_table.Element[pmd_atom.element]
         sigma, epsilon = _lj_params_from_potential(potential)
-        pmd_atom.sigma = sigma
-        pmd_atom.epsilon = epsilon
 
-    electrostatics_term = off_system.term_collection.terms["Electrostatics"]
-    for pmd_idx, pmd_atom in enumerate(structure.atoms):
-        partial_charge = (
-            electrostatics_term.potentials[str(pmd_idx)]
-            .to(unit.elementary_charge)
-            .magnitude
+        atom_type = pmd.AtomType(
+            name=element + str(pmd_idx + 1),
+            number=pmd_idx,
+            atomic_number=pmd_atom.atomic_number,
+            mass=pmd.periodic_table.Mass[element],
         )
-        pmd_atom.charge = partial_charge
+
+        atom_type.set_lj_params(eps=epsilon, rmin=sigma * 2 ** (1 / 6) / 2)
+        pmd_atom.atom_type = atom_type
+        pmd_atom.type = atom_type.name
+        pmd_atom.name = pmd_atom.type
+
+    for pmd_idx, pmd_atom in enumerate(structure.atoms):
+        if has_electrostatics:
+            partial_charge = electrostatics_handler.charge_map[(pmd_idx,)]
+            partial_charge_unitless = partial_charge.to(
+                unit.elementary_charge
+            ).magnitude
+            pmd_atom.charge = float(partial_charge_unitless)
+            pmd_atom.atom_type.charge = float(partial_charge_unitless)
+        else:
+            pmd_atom.charge = 0
 
     # Assign dummy residue names, GROMACS will not accept empty strings
     for res in structure.residues:
@@ -121,15 +190,18 @@ def to_parmed(off_system: Any) -> pmd.Structure:
     return structure
 
 
-def _convert_box(off_system: Any, structure: pmd.Structure) -> None:
+def _convert_box(box: unit.Quantity, structure: pmd.Structure) -> None:
     # TODO: Convert box vectors to box lengths + angles
-    lengths = off_system.box.to(unit("angstrom")).diagonal().magnitude
-    angles = 3 * [90]
-    structure.box = np.hstack([lengths, angles])
+    if box is None:
+        structure.box = [0, 0, 0, 90, 90, 90]
+    else:
+        lengths = box.to(unit("angstrom")).diagonal().magnitude
+        angles = 3 * [90]
+        structure.box = np.hstack([lengths, angles])
 
 
 def _lj_params_from_potential(potential):
-    sigma = potential.parameters["sigma"].to(unit.angstrom)
-    epsilon = potential.parameters["epsilon"].to(unit.Unit("kilocalorie/mol"))
+    sigma = potential.parameters["sigma"].to(unit.angstrom).magnitude
+    epsilon = potential.parameters["epsilon"].to(kcal_mol).magnitude
 
-    return sigma.magnitude, epsilon.magnitude
+    return sigma, epsilon