Parmed Conversion Improvement (#658)

* Include ImproperType reading in "from_parmed"
This PR is branched from PR #644. The parmed conversion to GMSO
topology was not able to read in various Improper connections to
GMSO. It would aggregate them with the Dihedral connections, from
when GMSO did not have explicit treatment of impropers, similar to
ParmEd. This PR aims to check a ParmEd object for structure.impropers
or structure.dihedrals with the dihedral.impropers flag set to
True and convert those to `gmso.Impropers`. Likewise, it looks to
read in periodic or harmonic impropers and generate the correct
`gmso.ImproperType` parametric potential for that type, if the flag
`refer_types` is set to True in the `from_parmed` call.

* Add pmd_improper_types_map function to convert_parmed
* Generate a parametric potential with harmonic type
expression from pmd.impropers
* Generate a parametric potential with periodic type
expression from pmd.dihedrals and dihedral.improper=True
* Validate testing with foyer opls_validation when using
functions in PR #644

Harmonic type expression for impropers is taken from
https://manual.gromacs.org/current/reference-manual/functions/bonded-interactions.html#improper-dihedrals-harmonic-type
This uses Xi as the independent variable. If someone testing this PR can
validate that this is the formalism we want to use when reading in
ImproperTypes from ParmEd improper_types, that would be a good sanity
check.
The periodic type expression is taken from
https://manual.gromacs.org/current/reference-manual/functions/bonded-interactions.html#equation-eqnperiodicpropdihedral
with phi as the independent variable.

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

* Add test for number of impropers from NNdimethylformamide molecule

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

* WIP- include gmso in impropertype testing

* WIP- Add support for NoneTypes in site charge/mass

* Address Review comments and create impropertype expression from HarmonicImproperPotential in library

* WIP- Add testing for impropers

* Modification for improper potential forms from the PotentialTemplateLibrary

* WIP- Use id instead of object for storing mapping; impropertypes test

* WIP- combined test for dihedral/impropertypes

* Fix test case and id mapping in conversion

* Add tests for harmonic impropers

* Add proper properties for improper_types

* WIP-Additional test cases

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Umesh Timalsina <umesh.timalsina@vanderbilt.edu>

gmso/external/convert_parmed.py

@@ -13,27 +13,30 @@
     element_by_name,
     element_by_symbol,
 )
+from gmso.lib.potential_templates import PotentialTemplateLibrary
 from gmso.utils.io import has_parmed, import_
 
 if has_parmed:
     pmd = import_("parmed")
 
+lib = PotentialTemplateLibrary()
+
 
 def from_parmed(structure, refer_type=True):
     """Convert a parmed.Structure to a gmso.Topology.
 
     Convert a parametrized or un-parametrized parmed.Structure object to a topology.Topology.
     Specifically, this method maps Structure to Topology and Atom to Site.
-    At this point, this method can only convert AtomType, BondType and AngleType.
-    Conversion of DihedralType will be implement in the near future.
+    This method can only convert AtomType, BondType AngleType, DihedralType, and
+    ImproperType.
 
     Parameters
     ----------
     structure : parmed.Structure
         parmed.Structure instance that need to be converted.
     refer_type : bool, optional, default=True
         Whether or not to transfer AtomType, BondType, AngleType,
-        and DihedralType information
+        DihedralType, and ImproperType information
 
     Returns
     -------
@@ -72,11 +75,23 @@ def from_parmed(structure, refer_type=True):
             structure, angle_types_member_map=angle_types_map
         )
         # Consolidate parmed dihedraltypes and relate to topology dihedraltypes
-        dihedral_types_map = _get_types_map(structure, "dihedrals")
+        # TODO: CCC seperate structure dihedrals.improper = False
+        dihedral_types_map = _get_types_map(
+            structure, "dihedrals", impropers=False
+        )
         dihedral_types_map.update(_get_types_map(structure, "rb_torsions"))
         pmd_top_dihedraltypes = _dihedral_types_from_pmd(
             structure, dihedral_types_member_map=dihedral_types_map
         )
+        # Consolidate parmed dihedral/impropertypes and relate to topology impropertypes
+        # TODO: CCC seperate structure dihedrals.improper = True
+        improper_types_map = _get_types_map(structure, "impropers")
+        improper_types_map.update(
+            _get_types_map(structure, "dihedrals"), impropers=True
+        )
+        pmd_top_impropertypes = _improper_types_from_pmd(
+            structure, improper_types_member_map=improper_types_map
+        )
 
     subtops = list()
     for residue in structure.residues:
@@ -158,41 +173,73 @@ def from_parmed(structure, refer_type=True):
         top.add_connection(top_connection, update_types=False)
 
     for dihedral in structure.dihedrals:
-        # Generate dihedral parameters for DihedralType that gets passed
-        # to Dihedral
+        # Generate parameters for ImproperType or DihedralType that gets passed
+        # to corresponding Dihedral or Improper
         # These all follow periodic torsions functions
-        # (even if they are improper dihedrals)
         # Which are the default expression in top.DihedralType
         # These periodic torsion dihedrals get stored in top.dihedrals
+        # and periodic torsion impropers get stored in top.impropers
+
         if dihedral.improper:
             warnings.warn(
                 "ParmEd improper dihedral {} ".format(dihedral)
                 + "following periodic torsion "
                 + "expression detected, currently accounted for as "
-                + "topology.Dihedral with a periodic torsion expression"
-            )
-        if refer_type and isinstance(dihedral.type, pmd.DihedralType):
-            top_connection = gmso.Dihedral(
-                connection_members=[
-                    site_map[dihedral.atom1],
-                    site_map[dihedral.atom2],
-                    site_map[dihedral.atom3],
-                    site_map[dihedral.atom4],
-                ],
-                dihedral_type=pmd_top_dihedraltypes[dihedral.type],
+                + "topology.Improper with a periodic improper expression"
             )
-        # No bond parameters, make Connection with no connection_type
+
+            if refer_type and isinstance(dihedral.type, pmd.DihedralType):
+                # TODO: Improper atom order is not always clear in a Parmed object.
+                # This reader assumes the order of impropers is central atom first,
+                # so that is where the central atom is located. This decision comes
+                # from .top files in utils/files/NN-dimethylformamide.top, which
+                # clearly places the periodic impropers with central atom listed first,
+                # and that is where the atom is placed in the parmed.dihedrals object.
+                top_connection = gmso.Improper(
+                    connection_members=[
+                        site_map[dihedral.atom1],
+                        site_map[dihedral.atom2],
+                        site_map[dihedral.atom3],
+                        site_map[dihedral.atom4],
+                    ],
+                    improper_type=pmd_top_impropertypes[id(dihedral.type)],
+                )
+            # No bond parameters, make Connection with no connection_type
+            else:
+                top_connection = gmso.Improper(
+                    connection_members=[
+                        site_map[dihedral.atom1],
+                        site_map[dihedral.atom2],
+                        site_map[dihedral.atom3],
+                        site_map[dihedral.atom4],
+                    ],
+                    improper_type=None,
+                )
+            top.add_connection(top_connection, update_types=False)
+
         else:
-            top_connection = gmso.Dihedral(
-                connection_members=[
-                    site_map[dihedral.atom1],
-                    site_map[dihedral.atom2],
-                    site_map[dihedral.atom3],
-                    site_map[dihedral.atom4],
-                ],
-                dihedral_type=None,
-            )
-        top.add_connection(top_connection, update_types=False)
+            if refer_type and isinstance(dihedral.type, pmd.DihedralType):
+                top_connection = gmso.Dihedral(
+                    connection_members=[
+                        site_map[dihedral.atom1],
+                        site_map[dihedral.atom2],
+                        site_map[dihedral.atom3],
+                        site_map[dihedral.atom4],
+                    ],
+                    dihedral_type=pmd_top_dihedraltypes[id(dihedral.type)],
+                )
+            # No bond parameters, make Connection with no connection_type
+            else:
+                top_connection = gmso.Dihedral(
+                    connection_members=[
+                        site_map[dihedral.atom1],
+                        site_map[dihedral.atom2],
+                        site_map[dihedral.atom3],
+                        site_map[dihedral.atom4],
+                    ],
+                    dihedral_type=None,
+                )
+            top.add_connection(top_connection, update_types=False)
 
     for rb_torsion in structure.rb_torsions:
         # Generate dihedral parameters for DihedralType that gets passed
@@ -214,7 +261,7 @@ def from_parmed(structure, refer_type=True):
                     site_map[rb_torsion.atom3],
                     site_map[rb_torsion.atom4],
                 ],
-                dihedral_type=pmd_top_dihedraltypes[rb_torsion.type],
+                dihedral_type=pmd_top_dihedraltypes[id(rb_torsion.type)],
             )
         # No bond parameters, make Connection with no connection_type
         else:
@@ -229,6 +276,36 @@ def from_parmed(structure, refer_type=True):
             )
         top.add_connection(top_connection, update_types=False)
 
+    for improper in structure.impropers:
+        if refer_type and isinstance(improper.type, pmd.ImproperType):
+            # TODO: Improper atom order is not always clear in a Parmed object.
+            # This reader assumes the order of impropers is central atom first,
+            # so that is where the central atom is located. This decision comes
+            # from .top files in utils/files/NN-dimethylformamide.top, which
+            # clearly places the periodic impropers with central atom listed first,
+            # and that is where the atom is placed in the parmed.dihedrals object.
+            top_connection = gmso.Improper(
+                connection_members=[
+                    site_map[improper.atom1],
+                    site_map[improper.atom2],
+                    site_map[improper.atom3],
+                    site_map[improper.atom4],
+                ],
+                improper_type=pmd_top_impropertypes[improper.type],
+            )
+            # No bond parameters, make Connection with no connection_type
+        else:
+            top_connection = gmso.Improper(
+                connection_members=[
+                    site_map[improper.atom1],
+                    site_map[improper.atom2],
+                    site_map[improper.atom3],
+                    site_map[improper.atom4],
+                ],
+                improper_type=None,
+            )
+        top.add_connection(top_connection, update_types=False)
+
     top.update_topology()
 
     top.combining_rule = structure.combining_rule
@@ -400,7 +477,7 @@ def _dihedral_types_from_pmd(structure, dihedral_types_member_map=None):
         top_dihedraltype = gmso.DihedralType(
             potential_expression=expr, member_types=member_types
         )
-        pmd_top_dihedraltypes[dihedraltype] = top_dihedraltype
+        pmd_top_dihedraltypes[id(dihedraltype)] = top_dihedraltype
 
     for dihedraltype in structure.rb_torsion_types:
         dihedral_params = {
@@ -422,10 +499,66 @@ def _dihedral_types_from_pmd(structure, dihedral_types_member_map=None):
             independent_variables="phi",
             member_types=member_types,
         )
-        pmd_top_dihedraltypes[dihedraltype] = top_dihedraltype
+        pmd_top_dihedraltypes[id(dihedraltype)] = top_dihedraltype
     return pmd_top_dihedraltypes
 
 
+def _improper_types_from_pmd(structure, improper_types_member_map=None):
+    """Convert ParmEd improper types to GMSO ImproperType.
+
+    This function take in a Parmed Structure, iterate through its
+    improper_types and dihedral_types with the `improper=True` flag,
+    create a corresponding GMSO.ImproperType, and finally return
+    a dictionary containing all pairs of pmd.ImproperType
+    (or pmd.DihedralType) and GMSO.ImproperType
+
+    Parameters
+    ----------
+    structure: pmd.Structure
+        Parmed Structure that needed to be converted.
+    improper_types_member_map: optional, dict, default=None
+        The member types (atomtype string) for each atom associated with the improper_types the structure
+
+    Returns
+    -------
+    pmd_top_impropertypes : dict
+        A dictionary linking a pmd.ImproperType or pmd.DihedralType
+        object to its corresponding GMSO.ImproperType object.
+    """
+    pmd_top_impropertypes = dict()
+    improper_types_member_map = _assert_dict(
+        improper_types_member_map, "improper_types_member_map"
+    )
+
+    for dihedraltype in structure.dihedral_types:
+        improper_params = {
+            "k": (dihedraltype.phi_k * u.Unit("kcal / mol")),
+            "phi_eq": (dihedraltype.phase * u.degree),
+            "n": dihedraltype.per * u.dimensionless,
+        }
+        expr = lib["PeriodicImproperPotential"]
+        member_types = improper_types_member_map.get(id(dihedraltype))
+        top_impropertype = gmso.ImproperType.from_template(
+            potential_template=expr, parameters=improper_params
+        )
+        pmd_top_impropertypes[id(dihedraltype)] = top_impropertype
+        top_impropertype.member_types = member_types
+
+    for impropertype in structure.improper_types:
+        improper_params = {
+            "k": (impropertype.psi_k * u.Unit("kcal/mol")),
+            "phi_eq": (impropertype.psi_eq * u.Unit("kcal/mol")),
+        }
+        expr = lib["HarmonicImproperPotential"]
+        member_types = improper_types_member_map.get(id(impropertype))
+        top_impropertype = gmso.ImproperType.from_template(
+            potential_template=expr, parameters=improper_params
+        )
+        top_impropertype.member_types = member_types
+        pmd_top_impropertypes[impropertype] = top_impropertype
+    return pmd_top_impropertypes
+
+
 def to_parmed(top, refer_type=True):
     """Convert a gmso.topology.Topology to a parmed.Structure.
 
@@ -460,6 +593,8 @@ def to_parmed(top, refer_type=True):
             top.box.lengths.to("angstrom").value,
             top.box.angles.to("degree").value,
         )
+        if top.box
+        else None
     )
 
     # Maps
@@ -479,8 +614,10 @@ def to_parmed(top, refer_type=True):
         pmd_atom = pmd.Atom(
             atomic_number=atomic_number,
             name=site.name,
-            mass=site.mass.to(u.amu).value,
-            charge=site.charge.to(u.elementary_charge).value,
+            mass=site.mass.to(u.amu).value if site.mass else None,
+            charge=site.charge.to(u.elementary_charge).value
+            if site.charge
+            else None,
         )
         pmd_atom.xx, pmd_atom.xy, pmd_atom.xz = site.position.to(
             "angstrom"
@@ -758,35 +895,56 @@ def _dihedral_types_from_gmso(top, structure, dihedral_map):
     structure.rb_torsions.claim()
 
 
-def _get_types_map(structure, attr):
+def _get_types_map(structure, attr, impropers=False):
     """Build `member_types` map for atoms, bonds, angles and dihedrals."""
-    assert attr in {"atoms", "bonds", "angles", "dihedrals", "rb_torsions"}
+    assert attr in {
+        "atoms",
+        "bonds",
+        "angles",
+        "dihedrals",
+        "rb_torsions",
+        "impropers",
+    }
     type_map = {}
     for member in getattr(structure, attr):
-        conn_type_id, member_types = _get_member_types_map_for(member)
+        conn_type_id, member_types = _get_member_types_map_for(
+            member, impropers
+        )
         if conn_type_id not in type_map and all(member_types):
             type_map[conn_type_id] = member_types
     return type_map
 
 
-def _get_member_types_map_for(member):
+def _get_member_types_map_for(member, impropers=False):
     if isinstance(member, pmd.Atom):
         return id(member.atom_type), member.type
-    if isinstance(member, pmd.Bond):
+    elif isinstance(member, pmd.Bond):
         return id(member.type), (member.atom1.type, member.atom2.type)
-    if isinstance(member, pmd.Angle):
+    elif isinstance(member, pmd.Angle):
         return id(member.type), (
             member.atom1.type,
             member.atom2.type,
             member.atom3.type,
         )
-    if isinstance(member, pmd.Dihedral):
-        return id(member.type), (
-            member.atom1.type,
-            member.atom2.type,
-            member.atom3.type,
-            member.atom4.type,
-        )
+    elif not impropers:  # return dihedrals
+        if isinstance(member, pmd.Dihedral) and not member.improper:
+            return id(member.type), (
+                member.atom1.type,
+                member.atom2.type,
+                member.atom3.type,
+                member.atom4.type,
+            )
+    elif impropers:  # return impropers
+        if (isinstance(member, pmd.Dihedral) and member.improper) or isinstance(
+            member, pmd.Improper
+        ):
+            return id(member.type), (
+                member.atom1.type,
+                member.atom2.type,
+                member.atom3.type,
+                member.atom4.type,
+            )
+    return None, (None, None)
 
 
 def _assert_dict(input_dict, param):