11import glob
22import itertools as it
33import os
4-
4+ from scipy . constants import N_A , k , calorie
55import parmed as pmd
66from pkg_resources import resource_filename
77import pytest
@@ -96,6 +96,13 @@ def mol2file(request):
9696 )
9797
9898
99+ def parmed2trappe (val ):
100+ return val * calorie / (k * N_A / 1000. )
101+
102+ def angle_parmed2trappe (val ):
103+ return 2 * parmed2trappe (val )
104+
105+
99106class TestTrappeReferences (object ):
100107 @staticmethod
101108 def check_atoms (ref , struc ):
@@ -106,8 +113,9 @@ def check_atoms(ref, struc):
106113
107114 assert_almost_equal (foyer_atom .charge , float (charge ))
108115 assert_almost_equal (foyer_atom .sigma , float (sigma ))
109- # TODO: figure out units here
110- #assert_almost_equal(foyer_atom.epsilon, float(epsilon))
116+ assert_almost_equal (parmed2trappe (foyer_atom .epsilon ),
117+ float (epsilon ),
118+ decimal = 3 )
111119
112120 @staticmethod
113121 def check_bonds (ref , struc ):
@@ -122,18 +130,43 @@ def check_bonds(ref, struc):
122130 # lookup foyer bond and compare length
123131 for ref_bond in ref_bonds :
124132 _ , desc , _ , length = ref_bond .strip ().split (',' )
125- i , j = map (int , ref_bond . split ( ',' )[ 1 ] .strip ("\" \' " ).split ('-' ))
133+ i , j = map (int , desc .strip ("\" \' " ).split ('-' ))
126134 foyer_bond = bond_dict [i , j ]
127135
128- assert_almost_equal (float ( length ), foyer_bond .type .req )
136+ assert_almost_equal (foyer_bond .type .req , float ( length ) )
129137 # trappe is fixed bond length, no point in K comparison
130138
131139 @staticmethod
132140 def check_angles (ref , struc ):
133141 ref_angles = ref ['bend' ]
134-
135142 assert len (ref_angles ) == len (struc .angles )
136143
144+ angle_dict = {}
145+ for a in struc .angles :
146+ i , j , k = a .atom1 .idx , a .atom2 .idx , a .atom3 .idx
147+ if k < i :
148+ i , j , k = k , j , i
149+ angle_dict [i + 1 , j + 1 , k + 1 ] = a
150+
151+ for ref_angle in ref_angles :
152+ _ , desc , _ , angle , k_angle = ref_angle .strip ().split ("," )
153+ i , j , k = map (int , desc .strip ("\" \' " ).split ('-' ))
154+
155+ foyer_angle = angle_dict [i , j , k ]
156+
157+ assert_almost_equal (foyer_angle .type .theteq ,
158+ float (angle ),
159+ decimal = 4 )
160+ assert_almost_equal (angle_parmed2trappe (foyer_angle .type .k ),
161+ float (k_angle ),
162+ decimal = 3 )
163+
164+ @staticmethod
165+ def check_dihedrals (ref , struc ):
166+ ref_dihedrals = ref ['torsion' ]
167+
168+ assert len (ref_dihedrals ) == len (struc .dihedrals )
169+
137170 @pytest .mark .parametrize ('trappe_ref,mol2file' , [
138171 ('methyl_acrylate/trappe_parameters_96.csv' , 'methyl_acrylate' ),
139172 ], indirect = True )
@@ -144,6 +177,8 @@ def test_trappe_reference(self, trappe_ref, mol2file):
144177
145178 self .check_angles (trappe_ref , mol2file )
146179
180+ self .check_dihedrals (trappe_ref , mol2file )
181+
147182
148183if __name__ == '__main__' :
149184 TestTraPPE ().find_correctly_implemented ()
0 commit comments