Repository: mosdef-hub/gmso

Issue:
Wrong expression in convert to foyer xml
The SymPy expressions are not consistent with foyer forcefield file.  E.g.,:

            "expression": "epsilon * ((sigma/r)**12 - (sigma/r)**6)",

https://github.com/mosdef-hub/gmso/blob/cc0f42d8cfda9b4236c0f1f4b724dc99036427f3/gmso/external/convert_foyer_xml.py#L249


These needs the factor of 4 out front.  Similarly harmonic potentials need the 1/2:

            "expression": "k * (r-r_eq)**2",

https://github.com/mosdef-hub/gmso/blob/cc0f42d8cfda9b4236c0f1f4b724dc99036427f3/gmso/external/convert_foyer_xml.py#L293

The output from, e.g., the LAMMPS writer is still correct, because it assumes the expression consistent with the original foyer XML file and does not directly examine the SymPy expression.

Discussion:
daico007: Pending unit tests 
codecov[bot]: # [Codecov](https://codecov.io/gh/mosdef-hub/gmso/pull/629?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mosdef-hub) Report
> Merging [#629](https://codecov.io/gh/mosdef-hub/gmso/pull/629?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mosdef-hub) (08ed956) into [master](https://codecov.io/gh/mosdef-hub/gmso/commit/cc0f42d8cfda9b4236c0f1f4b724dc99036427f3?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mosdef-hub) (cc0f42d) will **decrease** coverage by `0.03%`.
> The diff coverage is `83.33%`.

[![Impacted file tree graph](https://codecov.io/gh/mosdef-hub/gmso/pull/629/graphs/tree.svg?width=650&height=150&src=pr&token=rqPGwmXDzu&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mosdef-hub)](https://codecov.io/gh/mosdef-hub/gmso/pull/629?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mosdef-hub)

```diff
@@            Coverage Diff             @@
##           master     #629      +/-   ##
==========================================
- Coverage   90.78%   90.75%   -0.04%     
==========================================
  Files          56       56              
  Lines        4625     4629       +4     
==========================================
+ Hits         4199     4201       +2     
- Misses        426      428       +2     
```


| [Impacted Files](https://codecov.io/gh/mosdef-hub/gmso/pull/629?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mosdef-hub) | Coverage Δ | |
|---|---|---|
| [gmso/external/convert\_foyer\_xml.py](https://codecov.io/gh/mosdef-hub/gmso/pull/629/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mosdef-hub#diff-Z21zby9leHRlcm5hbC9jb252ZXJ0X2ZveWVyX3htbC5weQ==) | `92.23% <83.33%> (-0.33%)` | :arrow_down: |
| [gmso/core/bond\_type.py](https://codecov.io/gh/mosdef-hub/gmso/pull/629/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mosdef-hub#diff-Z21zby9jb3JlL2JvbmRfdHlwZS5weQ==) | `96.29% <0.00%> (-3.71%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/mosdef-hub/gmso/pull/629?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mosdef-hub).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mosdef-hub)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/mosdef-hub/gmso/pull/629?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mosdef-hub). Last update [cc0f42d...08ed956](https://codecov.io/gh/mosdef-hub/gmso/pull/629?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mosdef-hub). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mosdef-hub).

chrisiacovella: Regarding the function expressions: Do we want to also change the Urey-Bradley terms to have a prefactor of 1/2 for consistency with how we define harmonic potentials in the Foyer XML?  

The standard CHARMM forcefield definition does not include the 1/2 term for UB terms (the factor of 1/2 is dropped from all harmonic terms...so NOT having 1/2 would be consistent with CHARMM, but not consistent with how you'd have to define the bonds and angles anyway).   
daico007: Good catch, I will udpate the U-B expression 
bc118: The current equation compare functions I have (not in PR yet) should scale any different input scaling factor the user inputs.  The only standard required should really be the equation, with any scaling factor.  
chrisiacovella: So this is converting from the old foyer xml files where there is an implicitly defined functional form (not explicitly stated).   But yeah, more generally, the writers need to be doing these checks for equivalence and doing the appropriate scaling (I'm working on that right now for the lammps writer). 
bc118: I would use my equation compare functions which will be in a PR soon to scale the K constants for the potentials.  They are here for now (https://github.com/bc118/mbuild/blob/gmso_gomc/mbuild/utils/gmso_equation_compare.py)

Patch:
diff --git a/gmso/external/convert_foyer_xml.py b/gmso/external/convert_foyer_xml.py
index 8185b77..fa31c52 100644
--- a/gmso/external/convert_foyer_xml.py
+++ b/gmso/external/convert_foyer_xml.py
@@ -65,9 +65,11 @@ def from_foyer_xml(
     ff_root = foyer_xml_tree.getroot()
     name = ff_root.attrib.get("name")
     version = ff_root.attrib.get("version")
+    combining_rule = ff_root.attrib.get("combining_rule")
     f_kwargs = {
         "name": name,
         "version": version,
+        "combining_rule": combining_rule,
         "coulomb14scale": [],
         "lj14scale": [],
         "atom_types": [],
@@ -143,6 +145,10 @@ def _write_gmso_xml(gmso_xml, **kwargs):
         forcefield.attrib["version"] = "0.0.1"
 
     ffMeta = _create_sub_element(forcefield, "FFMetaData")
+    if kwargs.get("combining_rule"):
+        ffMeta.attrib["combining_rule"] = kwargs.get("combining_rule")
+    else:
+        ffMeta.attrib["combining_rule"] = "geometric"
     if kwargs["coulomb14scale"]:
         ffMeta.attrib["electrostatics14Scale"] = kwargs["coulomb14scale"]
 
@@ -246,7 +252,7 @@ def _write_nbforces(forcefield, ff_kwargs):
         forcefield,
         "AtomTypes",
         attrib_dict={
-            "expression": "epsilon * ((sigma/r)**12 - (sigma/r)**6)",
+            "expression": "4 * epsilon * ((sigma/r)**12 - (sigma/r)**6)",
         },
     )
     parameters_units = {"epsilon": "kJ/mol", "sigma": "nm"}
@@ -290,7 +296,7 @@ def _write_harmonic_bonds(forcefield, ff_kwargs):
         forcefield,
         "BondTypes",
         attrib_dict={
-            "expression": "k * (r-r_eq)**2",
+            "expression": "1/2 * k * (r-r_eq)**2",
         },
     )
 
@@ -323,7 +329,7 @@ def _write_harmonic_angles(forcefield, ff_kwargs):
         forcefield,
         "AngleTypes",
         attrib_dict={
-            "expression": "k * (theta - theta_eq)**2",
+            "expression": "1/2 * k * (theta - theta_eq)**2",
         },
     )
 
@@ -357,7 +363,7 @@ def _write_ub_angles(forcefield, ff_kwargs):
         forcefield,
         "AngleTypes",
         attrib_dict={
-            "expression": "k * (w - w_0) ** 2",
+            "expression": "1/2 * k * (w - w_0) ** 2",
         },
     )
 
diff --git a/gmso/formats/lammpsdata.py b/gmso/formats/lammpsdata.py
index eb4bae2..1256ca7 100644
--- a/gmso/formats/lammpsdata.py
+++ b/gmso/formats/lammpsdata.py
@@ -6,7 +6,7 @@ import warnings
 
 import numpy as np
 import unyt as u
-from sympy import sympify
+from sympy import simplify, sympify
 from unyt.array import allclose_units
 
 from gmso.core.angle import Angle
@@ -183,48 +183,90 @@ def write_lammpsdata(topology, filename, atom_style="full"):
             # Pair coefficients
             data.write("\nPair Coeffs # lj\n\n")
             for idx, param in enumerate(topology.atom_types):
-                data.write(
-                    "{}\t{:.5f}\t{:.5f}\n".format(
-                        idx + 1,
-                        param.parameters["epsilon"]
-                        .in_units(u.Unit("kcal/mol"))
-                        .value,
-                        param.parameters["sigma"].in_units(u.angstrom).value,
-                    )
+                # expected expression for lammps for standard LJ
+                lj_expression = "4.0 * epsilon * ((sigma/r)**12 - (sigma/r)**6)"
+                scaling_factor = simplify(lj_expression) / simplify(
+                    param.expression
                 )
 
-            if topology.bonds:
-                data.write("\nBond Coeffs\n\n")
-                for idx, bond_type in enumerate(topology.bond_types):
+                if scaling_factor.is_real:
                     data.write(
                         "{}\t{:.5f}\t{:.5f}\n".format(
                             idx + 1,
-                            bond_type.parameters["k"]
-                            .in_units(u.Unit("kcal/mol/angstrom**2"))
+                            param.parameters["epsilon"]
+                            .in_units(u.Unit("kcal/mol"))
                             .value
-                            / 2,
-                            bond_type.parameters["r_eq"]
-                            .in_units(u.Unit("angstrom"))
+                            / float(scaling_factor),
+                            param.parameters["sigma"]
+                            .in_units(u.angstrom)
                             .value,
                         )
                     )
+                else:
+                    raise ValueError(
+                        'Pair Style "{}" is invalid or is not currently supported'.format(
+                            param.expression
+                        )
+                    )
+            if topology.bonds:
+                data.write("\nBond Coeffs\n\n")
+                for idx, bond_type in enumerate(topology.bond_types):
+
+                    # expected harmonic potential expression for lammps
+                    bond_expression = "k * (r-r_eq)**2"
+
+                    scaling_factor = simplify(bond_expression) / simplify(
+                        bond_type.expression
+                    )
+
+                    if scaling_factor.is_real:
+                        data.write(
+                            "{}\t{:.5f}\t{:.5f}\n".format(
+                                idx + 1,
+                                bond_type.parameters["k"]
+                                .in_units(u.Unit("kcal/mol/angstrom**2"))
+                                .value
+                                / float(scaling_factor),
+                                bond_type.parameters["r_eq"]
+                                .in_units(u.Unit("angstrom"))
+                                .value,
+                            )
+                        )
+                    else:
+                        raise ValueError(
+                            'Bond Style "{}" is invalid or is not currently supported'.format(
+                                bond_type.expression
+                            )
+                        )
 
             if topology.angles:
                 data.write("\nAngle Coeffs\n\n")
                 for idx, angle_type in enumerate(topology.angle_types):
-                    data.write(
-                        "{}\t{:.5f}\t{:.5f}\n".format(
-                            idx + 1,
-                            angle_type.parameters["k"]
-                            .in_units(u.Unit("kcal/mol/radian**2"))
-                            .value
-                            / 2,
-                            angle_type.parameters["theta_eq"]
-                            .in_units(u.Unit("degree"))
-                            .value,
-                        )
+                    # expected lammps harmonic angle expression
+                    angle_expression = "k * (theta - theta_eq)**2"
+                    scaling_factor = simplify(angle_expression) / simplify(
+                        angle_type.expression
                     )
 
+                    if scaling_factor.is_real:
+                        data.write(
+                            "{}\t{:.5f}\t{:.5f}\n".format(
+                                idx + 1,
+                                angle_type.parameters["k"]
+                                .in_units(u.Unit("kcal/mol/radian**2"))
+                                .value
+                                / float(scaling_factor),
+                                angle_type.parameters["theta_eq"]
+                                .in_units(u.Unit("degree"))
+                                .value,
+                            )
+                        )
+                    else:
+                        raise ValueError(
+                            'Angle Style "{}" is invalid or is not currently supported'.format(
+                                angle_type.expression
+                            )
+                        )
             # TODO: Write out multiple dihedral styles
             if topology.dihedrals:
                 data.write("\nDihedral Coeffs\n\n")
@@ -403,7 +445,7 @@ def get_units(unit_style):
     # Need separate angle units for harmonic force constant and angle
     unit_style_dict = {
         "real": {
-            "mass": u.g,
+            "mass": u.g / u.mol,
             "distance": u.angstrom,
             "energy": u.kcal / u.mol,
             "angle_k": u.radian,
@@ -444,7 +486,7 @@ def _get_connection(filename, topology, unit_style, connection_type):
                 * 2
             )
             c_type.parameters["r_eq"] = float(line.split()[2]) * (
-                get_units(unit_style)["distance"] ** 2
+                get_units(unit_style)["distance"]
             )
         elif connection_type == "angle":
             c_type = AngleType(name=line.split()[0])
diff --git a/gmso/utils/schema/ff-gmso.xsd b/gmso/utils/schema/ff-gmso.xsd
index 43a843f..f4bc4e9 100644
--- a/gmso/utils/schema/ff-gmso.xsd
+++ b/gmso/utils/schema/ff-gmso.xsd
@@ -211,6 +211,7 @@
             <xsd:sequence>
                 <xsd:element name="Units" type="Units" minOccurs="0" maxOccurs="1"/>
             </xsd:sequence>
+            <xsd:attribute name="combining_rule" type="xsd:string" use="optional"/>
             <xsd:attribute name="electrostatics14Scale" type="xsd:float" use="optional"/>
             <xsd:attribute name="nonBonded14Scale" type="xsd:float" use="optional"/>
             <xsd:attribute name="combiningRule" type="xsd:string" use="optional"/>


Test patch:
diff --git a/gmso/tests/test_convert_foyer_xml.py b/gmso/tests/test_convert_foyer_xml.py
index ff1fbc8..2f7d593 100644
--- a/gmso/tests/test_convert_foyer_xml.py
+++ b/gmso/tests/test_convert_foyer_xml.py
@@ -85,13 +85,16 @@ class TestXMLConversion(BaseTest):
         assert foyer_fullerene.atom_types["C"].description == "carbon"
         assert foyer_fullerene.atom_types["C"].definition == "[C;r5;r6]"
         assert foyer_fullerene.atom_types["C"].expression == sympify(
-            "epsilon*(-sigma**6/r**6 + sigma**12/r**12)"
+            "4*epsilon*(-sigma**6/r**6 + sigma**12/r**12)"
         )
 
     def test_foyer_bonds(self, foyer_fullerene):
         assert len(foyer_fullerene.bond_types) == 1
         assert "C~C" in foyer_fullerene.bond_types
 
+        assert foyer_fullerene.bond_types["C~C"].expression == sympify(
+            "1/2 * k * (r-r_eq)**2"
+        )
         assert (
             sympify("r")
             in foyer_fullerene.bond_types["C~C"].independent_variables
@@ -108,6 +111,9 @@ class TestXMLConversion(BaseTest):
         assert len(foyer_fullerene.angle_types) == 1
         assert "C~C~C" in foyer_fullerene.angle_types
 
+        assert foyer_fullerene.angle_types["C~C~C"].expression == sympify(
+            "1/2 * k * (theta - theta_eq)**2"
+        )
         assert (
             sympify("theta")
             in foyer_fullerene.angle_types["C~C~C"].independent_variables
diff --git a/gmso/tests/test_lammps.py b/gmso/tests/test_lammps.py
index be7ae23..a74788a 100644
--- a/gmso/tests/test_lammps.py
+++ b/gmso/tests/test_lammps.py
@@ -30,7 +30,7 @@ class TestLammpsWriter(BaseTest):
     def test_read_lammps(self, filename=get_path("data.lammps")):
         top = gmso.Topology.load(filename)
 
-    def test_read_sites(self, filename=get_path("data.lammps")):
+    def test_read_box(self, filename=get_path("data.lammps")):
         read = gmso.Topology.load(filename)
 
         assert read.box == Box(lengths=[1, 1, 1])
@@ -46,7 +46,7 @@ class TestLammpsWriter(BaseTest):
         masses = [i.mass for i in read.atom_types]
 
         assert_allclose_units(
-            masses, u.unyt_array(1.0079, u.g), rtol=1e-5, atol=1e-8
+            masses, u.unyt_array(1.0079, u.g / u.mol), rtol=1e-5, atol=1e-8
         )
 
     def test_read_charge(self, filename=get_path("data.lammps")):
@@ -106,3 +106,84 @@ class TestLammpsWriter(BaseTest):
             rtol=1e-5,
             atol=1e-8,
         )
+
+    def test_read_n_bonds(self, typed_ethane):
+        typed_ethane.save("ethane.lammps")
+        read = gmso.Topology.load("ethane.lammps")
+
+        assert read.n_bonds == 7
+
+    def test_read_n_angles(self, typed_ethane):
+        typed_ethane.save("ethane.lammps")
+        read = gmso.Topology.load("ethane.lammps")
+
+        assert read.n_angles == 12
+
+    def test_read_bond_params(self, typed_ethane):
+        typed_ethane.save("ethane.lammps")
+        read = gmso.Topology.load("ethane.lammps")
+        bond_params = [i.parameters for i in read.bond_types]
+
+        assert_allclose_units(
+            bond_params[0]["k"],
+            u.unyt_array(680, (u.kcal / u.mol / u.angstrom / u.angstrom)),
+            rtol=1e-5,
+            atol=1e-8,
+        )
+        assert_allclose_units(
+            bond_params[0]["r_eq"],
+            u.unyt_array(1.09, (u.angstrom)),
+            rtol=1e-5,
+            atol=1e-8,
+        )
+        assert_allclose_units(
+            bond_params[1]["k"],
+            u.unyt_array(536, (u.kcal / u.mol / u.angstrom / u.angstrom)),
+            rtol=1e-5,
+            atol=1e-8,
+        )
+        assert_allclose_units(
+            bond_params[1]["r_eq"],
+            u.unyt_array(1.529, (u.angstrom)),
+            rtol=1e-5,
+            atol=1e-8,
+        )
+
+    def test_read_angle_params(self, typed_ethane):
+        typed_ethane.save("ethane.lammps")
+        read = gmso.Topology.load("ethane.lammps")
+        angle_params = [i.parameters for i in read.angle_types]
+
+        assert_allclose_units(
+            angle_params[0]["k"],
+            u.unyt_array(75, (u.kcal / u.mol / u.radian / u.radian)),
+            rtol=1e-5,
+            atol=1e-8,
+        )
+        assert_allclose_units(
+            angle_params[0]["theta_eq"],
+            u.unyt_array(110.7, (u.degree)),
+            rtol=1e-5,
+            atol=1e-8,
+        )
+        assert_allclose_units(
+            angle_params[1]["k"],
+            u.unyt_array(66, (u.kcal / u.mol / u.radian / u.radian)),
+            rtol=1e-5,
+            atol=1e-8,
+        )
+        assert_allclose_units(
+            angle_params[1]["theta_eq"],
+            u.unyt_array(107.8, (u.degree)),
+            rtol=1e-5,
+            atol=1e-8,
+        )
+
+
+"""
+    def test_read_n_diherals(self, typed_ethane):
+        typed_ethane.save("ethane.lammps")
+        read = gmso.Topology.load("ethane.lammps")
+
+        assert read.n_dihedrals == 9
+"""