Support list of parameter values in Forcefield XML (#681)

* Support list of parameter values in Forcefield XML * fix indentation for aesthetic purpose Co-authored-by: Co Quach <daico007@gmail.com>
mosdef-hub · Jul 28, 2022 · cc23eb8 · cc23eb8
1 parent c336588
commit cc23eb8
Show file tree

Hide file tree

Showing 4 changed files with 84 additions and 2 deletions.
diff --git a/gmso/tests/files/sequence_of_parameters_ff.xml b/gmso/tests/files/sequence_of_parameters_ff.xml
@@ -0,0 +1,53 @@
+<ForceField version="0.0.0" name="Flexible carbon">
+    <FFMetaData>
+        <Units energy="kJ/mol" distance="nm" mass="amu" charge="coulomb"/>
+    </FFMetaData>
+    <AtomTypes expression="4*epsilon*((sigma/r)**12 - (sigma/r)**6)">
+        <ParametersUnitDef parameter="sigma" unit="nm"/>
+	<ParametersUnitDef parameter="epsilon" unit="kJ/mol"/>
+	<AtomType name="C" element="C" charge="0.0" mass="12.011" definition="C" description="Generic sp2 carbon">
+            <Parameters>
+                <Parameter name="sigma" value="0.339966950842"/>
+                <Parameter name="epsilon" value="0.359824"/>
+            </Parameters>
+        </AtomType>
+    </AtomTypes>
+    <BondTypes expression="0.5 * k * (r-r_eq)**2">
+        <ParametersUnitDef parameter="r_eq" unit="nm"/>
+	<ParametersUnitDef parameter="k" unit="kJ/(mol*nm**2)"/>
+        <BondType name="BondType1" type1='C' type2='C'>
+            <Parameters>
+		<Parameter name='r_eq' value="0.1324"/>
+                <Parameter name='k' value="493460.96"/>
+            </Parameters>
+        </BondType>
+    </BondTypes>
+    <AngleTypes expression="0.5 * k * (theta-theta_eq)**2">
+        <ParametersUnitDef parameter="theta_eq" unit="radian"/>
+	<ParametersUnitDef parameter="k" unit="kJ/(mol*radian**2)"/>
+        <AngleType name="AngleType1" type1='C' type2='C' type3="C">
+            <Parameters>
+                <Parameter name='theta_eq' value="2.12598556185"/>
+		<Parameter name='k' value="584.42112"/>
+            </Parameters>
+        </AngleType>
+    </AngleTypes>
+    <DihedralTypes expression="k * (1 + cos(n * theta - theta_0))">
+ 	<ParametersUnitDef parameter="k" unit="kJ/mol"/>
+        <ParametersUnitDef parameter="n" unit="dimensionless"/>
+	<ParametersUnitDef parameter="theta_0" unit="radian"/>
+        <DihedralType name="DihedralType1" type1="" type2="C" type3="C" type4="">
+            <Parameters>
+                <Parameter name='k'>
+                    <Value>38</Value>
+                    <Value>45</Value>
+                </Parameter>
+                <Parameter name='n' value="2"/>
+                <Parameter name='theta_0'>
+                    <Value>25</Value>
+                    <Value>32</Value>
+                </Parameter>
+            </Parameters>
+        </DihedralType>
+    </DihedralTypes>
+</ForceField>
diff --git a/gmso/tests/test_forcefield.py b/gmso/tests/test_forcefield.py
@@ -634,3 +634,11 @@ def test_write_xml(self, opls_ethane_foyer):
     def test_write_not_xml(self, opls_ethane_foyer):
         with pytest.raises(ForceFieldError):
             opls_ethane_foyer.xml("bad_path")
+
+    def test_valid_sequence(self):
+        for j in range(10):
+            ff = ForceField(get_path("sequence_of_parameters_ff.xml"), "r")
+            dih_with_list = ff.dihedral_types["*~C~C~*"]
+            params = dih_with_list.get_parameters()
+            assert u.allclose_units(params["theta_0"], [25, 32] * u.radian)
+            assert u.allclose_units(params["k"], [38, 45] * u.kJ / u.mol)
diff --git a/gmso/utils/ff_utils.py b/gmso/utils/ff_utils.py
@@ -2,6 +2,7 @@
 import os
 import re
 
+import numpy as np
 import unyt as u
 from lxml import etree
 from sympy import sympify
@@ -70,7 +71,24 @@ def _parse_params_values(parent_tag, units_dict, child_tag, expression=None):
             )
         param_name = param.attrib["name"]
         param_unit = units_dict[param_name]
-        param_value = u.unyt_quantity(float(param.attrib["value"]), param_unit)
+        if param.attrib.get("value"):
+            param_value = u.unyt_quantity(
+                float(param.attrib["value"]), param_unit
+            )
+        else:
+            children = param.getchildren()
+            if len(children) == 0:
+                raise ForceFieldParseError(
+                    f"Neither a single value nor a sequence of values "
+                    f"is specified for parameter {param_name}, please specify "
+                    f"either a single value as an attribute value or a sequence "
+                    f"of values."
+                )
+            value_array = np.array(
+                [value.text for value in children], dtype=float
+            )
+            param_value = u.unyt_array(value_array, param_unit)
+
         params_dict[param_name] = param_value
     param_ref_dict = units_dict
     if child_tag == "DihedralType":

diff --git a/gmso/utils/schema/ff-gmso.xsd b/gmso/utils/schema/ff-gmso.xsd
@@ -193,8 +193,11 @@
             <xsd:sequence>
                 <xsd:element name="Parameter" minOccurs="0" maxOccurs="unbounded">
                     <xsd:complexType>
+                        <xsd:sequence>
+                            <xsd:element name="Value" minOccurs="0" maxOccurs="unbounded" type="xsd:float"/>
+                        </xsd:sequence>
                         <xsd:attribute name="name" use="required" type="xsd:string"/>
-                        <xsd:attribute name="value" use="required" type="xsd:float"/>
+                        <xsd:attribute name="value" use="optional" type="xsd:float"/>
                     </xsd:complexType>
                 </xsd:element>
             </xsd:sequence>