added round trip tests for dihedrals in trappe

Author: richardjgowers

foyer/forcefields/trappe-ua.xml

@@ -189,6 +189,21 @@
   <!-- CHX-CHY-[C]=[O] -->
   <Proper class1="CHx" class2="CHx" type3="CH_ald" class4="O_sp2" c0="9.32159" c1="-1.18722" c2="-0.96182" c3="9.7522" c4="0.0" c5="0.0"/>
   <Proper class1="CHx" class2="CHx" type3="C_ket" class4="O_sp2" c0="9.32159" c1="-1.18722" c2="-0.96182" c3="9.7522" c4="0.0" c5="0.0"/>
+  <!-- acrylates -->
+  <!-- maerzke T4.1 C-O-C=O -->
+  <Proper class1="CHx" type2="O_acryl_sp3" type3="C_acryl" type4="O_acryl_sp2" c0="-26.55539" c1="2.7175" c2="-22.83383" c3="1.00405" c4="0.0" c5="0.0"/>
+  <!-- maerzke T4.2 C-O-C(sp2)-C(sp2) -->
+  <Proper class1="CHx" type2="O_acryl_sp3" type3="C_acryl" class4="CHx_sp2" c0="-26.55539" c1="-2.7175" c2="-22.83383" c3="-1.00405" c4="0.0" c5="0.0"/>
+  <!-- maerzke T4.3 C(sp2)=C(sp2)-C=O -->
+  <Proper class1="CHx_sp2" class2="CHx_sp2" type3="C_acryl" type4="O_acryl_sp2" c0="-13.27121" c1="0.34871" c2="-12.85632" c3="0.06618" c4="0.0" c5="0.0"/>
+  <!-- maerzke T4.4 C(sp2)=C-C-O -->
+  <Proper class1="CHx_sp2" class2="CHx_sp2" type3="C_acryl" type4="O_acryl_sp3" c0="-13.27121" c1="-0.34871" c2="-12.85632" c3="-0.06618" c4="0.0" c5="0.0"/>
+  <!-- maerzke T4.5 O=C-C(sp2)-C -->
+  <Proper type1="O_acryl_sp2" type2="C_acryl" class3="CHx_sp2" class4="CHx" c0="-0.49271" c1="1.84298" c2="0.84109" c3="-0.80218" c4="1.89237" c5="0.0"/>
+  <!-- maerzke T4.6 O-C-C(sp2)-C -->
+  <Proper type1="O_acryl_sp3" type2="C_acryl" class3="CHx_sp2" class4="CHx" c0="-0.49271" c1="-1.84298" c2="0.84109" c3="0.80218" c4="1.89237" c5="0.0"/>
+  <!-- maerzke T4.8 C=C-C=C -->
+  <Proper type1="CH2_sp2" type2="CH_sp2-Csp2" type3="CH_sp2-Csp2" type4="CH2_sp2" c0="-28.85375" c1="7.05523" c2="-33.6714" c3="-15.29994" c4="13.06235" c5="0.0"/>
   <!-- 
   No parameters for CHX-[CH2]-[CH]-OH, CHX-[CH2]-[C]-OH, CHX-[CH]-[CH2]-OH,
   CHX-[CH]-[CH]-OH, CHX-[CH]-[C]-OH, CHX-[C]-[CH2]-OH, CHX-[C]-[CH]-OH,
@@ -197,7 +212,7 @@
   <!-- alkenes -->
   <Proper class1="CHx" type2="CH2_sp3" type3="CH_sp2" class4="CHx_sp2" c0="2.27051" c1="-7.75806" c2="1.82536" c3="9.38669" c4="0.0" c5="0.0"/>
   <Proper class1="CHx" type2="CH2_sp3-Cald" type3="CH_sp2" class4="CHx_sp2" c0="2.27051" c1="-7.75806" c2="1.82536" c3="9.38669" c4="0.0" c5="0.0"/>
-  <Proper type1="CH2_sp2" type2="CH_sp2-Csp2" type3="CH_sp2-Csp2" type4="CH2_sp2" c0="4.55208" c1="7.05523" c2="20.60905" c3="-15.29994" c4="0.0" c5="0.0"/>
+
   <Proper type1="CH2_sp2" type2="CH_sp2-Csp2" type3="C_sp2-Csp2" class4="CHx" c0="-9.69998" c1="-11.57082" c2="17.43875" c3="19.30751" c4="0.0" c5="0.0"/>
  </RBTorsionForce>
  <NonbondedForce coulomb14scale="0" lj14scale="0">

foyer/tests/test_trappe.py

@@ -102,6 +102,41 @@ def parmed2trappe(val):
 def angle_parmed2trappe(val):
     return 2 * parmed2trappe(val)
 
+def assert_dihedrals_match(foyer_d, reference_params):
+    def conv1(t):
+        # convert RB to Trappe#1
+        raise ValueError
+        R0, R1, R2, R3 = t.c0, t.c1, t.c2, t.c3
+        c0 = R0 - R1 + R2 - R3
+        c1 = R1 + 3/4. * R3
+        c2 = -0.5 * R2
+        c3 = 1/4. * R3
+
+        return c0, c1, c2, c3
+
+    def conv2(t):
+        # convert RB to Trappe#2
+        R0, R1, R2, R3, R4 = t.c0, t.c1, t.c2, t.c3, t.c4
+        print(list(map(lambda x: x * calorie, [R0, R1, R2, R3, R4])))
+        c0 = - R0 + 0.5 * R2 + 3/8. * R4
+        c1 = - R1 - 3/4. * R3
+        c2 = 0.5 * R2 + 0.5 * R4
+        c3 = - 1/4. * R3
+        c4 = 1/8. * R4
+
+        return [parmed2trappe(v) for v in [c0, c1, c2, c3, c4]]
+
+    reference_params = list(map(float, reference_params))
+    # check parameters in *foyer_d* match the raw text from trappe
+    # trappe's params are in one of two functional forms.....
+    if len(reference_params) == 5:
+        # using later form
+        foyer_params = conv2(foyer_d.type)
+    else:
+        # using original dihedral form
+        foyer_params = conv1(foyer_d.type)
+
+    assert_almost_equal(foyer_params, reference_params, decimal=3)
 
 class TestTrappeReferences(object):
     @staticmethod
@@ -149,7 +184,7 @@ def check_angles(ref, struc):
             angle_dict[i+1, j+1, k+1] = a
 
         for ref_angle in ref_angles:
-            _, desc, _, angle, k_angle = ref_angle.strip().split(",")
+            _, desc, _, angle, k_angle = ref_angle.strip().split(',')
             i, j, k = map(int, desc.strip("\"\'").split('-'))
 
             foyer_angle = angle_dict[i, j, k]
@@ -165,7 +200,23 @@ def check_angles(ref, struc):
     def check_dihedrals(ref, struc):
         ref_dihedrals = ref['torsion']
 
-        assert len(ref_dihedrals) == len(struc.dihedrals)
+        assert len(ref_dihedrals) == len(struc.rb_torsions)
+
+        dih_dict = {}
+        for d in struc.rb_torsions:
+            i, j, k, l = d.atom1.idx, d.atom2.idx, d.atom3.idx, d.atom4.idx
+
+            if l < i:
+                i, j, k, l = l, k, j, i
+            dih_dict[i+1, j+1, k+1, l+1] = d
+
+        for ref_d in ref_dihedrals:
+            _, desc, _, *params = ref_d.strip().split(',')
+            i, j, k, l = map(int, desc.strip("\"\'").split('-'))
+
+            foyer_d = dih_dict[i, j, k, l]
+
+            assert_dihedrals_match(foyer_d, params)
 
     @pytest.mark.parametrize('trappe_ref,mol2file', [
         ('methyl_acrylate/trappe_parameters_96.csv', 'methyl_acrylate'),