Skip to content

Commit 19e642d

Browse files
richardjgowersrichardjgowers
committed
added more Trappe parameters
added atom types for - alkenes - ethers - acrylates
1 parent ef39f62 commit 19e642d

File tree

14 files changed

+621
-169
lines changed

14 files changed

+621
-169
lines changed

‎foyer/forcefields/trappe-ua.xml

+203-167
Large diffs are not rendered by default.

‎foyer/tests/implemented_trappe_tests.txt

+6
Original file line numberDiff line numberDiff line change
@@ -8,3 +8,9 @@ methane
88
methanol
99
neopentane
1010
propane
11+
dimethylether
12+
aldehyde_1
13+
aldehyde_2
14+
aldehyde_3
15+
aldehyde_4
16+
methyl_acrylate

‎foyer/tests/test_trappe.py

+191-2
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,16 @@
11
import glob
22
import itertools as it
33
import os
4-
4+
from scipy.constants import N_A, k, calorie
55
import parmed as pmd
66
from pkg_resources import resource_filename
77
import pytest
88

99
from foyer import Forcefield
1010
from foyer.tests.utils import atomtype
1111

12+
from numpy.testing import assert_almost_equal
13+
1214
TRAPPE_UA = Forcefield(name='trappe-ua')
1315

1416
TRAPPE_TESTFILES_DIR = resource_filename('foyer', 'trappe_validation')
@@ -45,12 +47,199 @@ def find_correctly_implemented(self):
4547

4648
@pytest.mark.parametrize('mol_name', correctly_implemented)
4749
def test_atomtyping(self, mol_name, testfiles_dir=TRAPPE_TESTFILES_DIR):
48-
files = glob.glob(os.path.join(testfiles_dir, mol_name, '*'))
50+
files = glob.glob(os.path.join(testfiles_dir, mol_name, '*.mol2'))
4951
for mol_file in files:
5052
_, ext = os.path.splitext(mol_file)
5153
mol2_path = os.path.join(testfiles_dir, mol_name, mol_file)
5254
structure = pmd.load_file(mol2_path, structure=True)
5355
atomtype(structure, TRAPPE_UA, non_atomistic=True)
5456

57+
@pytest.fixture
58+
def trappe_ref(request):
59+
# load raw trappe csv file and process slightly
60+
container = {
61+
'(pseudo)atom': [],
62+
'stretch': [],
63+
'bend': [],
64+
'torsion': [],
65+
}
66+
# split csv file into sections
67+
csv_file = os.path.join(TRAPPE_TESTFILES_DIR,
68+
request.param)
69+
with open(csv_file, 'r') as inf:
70+
for line in inf:
71+
if line.startswith('#'):
72+
section_clue = line.split(",")[1]
73+
section = container[section_clue]
74+
line = next(inf)
75+
while not line == '\n':
76+
# slurp until newline hit
77+
section.append(line)
78+
try:
79+
line = next(inf)
80+
except StopIteration:
81+
break
82+
83+
return container
84+
85+
86+
@pytest.fixture
87+
def mol2file(request):
88+
# load a mol2 file and apply foyer's trappe-ua forcefield
89+
fname = request.param
90+
molfile = glob.glob(os.path.join(TRAPPE_TESTFILES_DIR, fname, '*'))[0]
91+
structure = pmd.load_file(molfile, structure=True)
92+
93+
return TRAPPE_UA.apply(structure,
94+
assert_angle_params=False,
95+
assert_dihedral_params=False,
96+
)
97+
98+
99+
def parmed2trappe(val):
100+
# convert kcal/mol (parmed units) to K/k_B (trappe units)
101+
return val * calorie / (k * N_A / 1000.)
102+
103+
def angle_parmed2trappe(val):
104+
# trappe angle units need doubling
105+
return 2 * parmed2trappe(val)
106+
107+
def assert_dihedrals_match(foyer_d, reference_params):
108+
# Check that a dihedral from foyer matches the text in the csv file
109+
# foyer_d - parmed rb_torsion
110+
# reference_params - dihedral parameters from trappe csv file
111+
def conv1(t):
112+
# convert RB to original Trappe dihedral form:
113+
# V = c0 + c1(1 + cos phi) + c1(1 - cos 2phi) + c3(1 + cos 3phi)
114+
R0, R1, R2, R3 = t.c0, t.c1, t.c2, t.c3
115+
c0 = - R0 + R1 + R2 + R3
116+
c1 = - R1 - 3/4. * R3
117+
c2 = - 0.5 * R2
118+
c3 = - 1/4. * R3
119+
120+
return c0, c1, c2, c3
121+
122+
def conv2(t):
123+
# convert RB to second Trappe dihedral form:
124+
# V = c0 + c1 cos phi + c2 cos 2phi + c3 cos 3phi + c4 cos 4phi
125+
R0, R1, R2, R3, R4 = t.c0, t.c1, t.c2, t.c3, t.c4
126+
print(list(map(lambda x: x * calorie, [R0, R1, R2, R3, R4])))
127+
c0 = - R0 + 0.5 * R2 + 3/8. * R4
128+
c1 = - R1 - 3/4. * R3
129+
c2 = 0.5 * R2 + 0.5 * R4
130+
c3 = - 1/4. * R3
131+
c4 = 1/8. * R4
132+
133+
return [parmed2trappe(v) for v in [c0, c1, c2, c3, c4]]
134+
135+
reference_params = list(map(float, reference_params))
136+
# check parameters in *foyer_d* match the raw text from trappe
137+
# trappe's params are in one of two functional forms.....
138+
if len(reference_params) == 5:
139+
# using later form, has an extra column
140+
foyer_params = conv2(foyer_d.type)
141+
elif len(reference_params) == 4:
142+
# using original dihedral form
143+
foyer_params = conv1(foyer_d.type)
144+
145+
assert_almost_equal(foyer_params, reference_params, decimal=3)
146+
147+
class TestTrappeReferences(object):
148+
@staticmethod
149+
def check_atoms(ref, struc):
150+
ref_atoms = sorted(ref['(pseudo)atom'],
151+
key=lambda x: int(x.split(',')[0]))
152+
for ref_atom, foyer_atom in zip(ref_atoms, struc.atoms):
153+
_, _, _, epsilon, sigma, charge = ref_atom.strip().split(',')
154+
155+
assert_almost_equal(foyer_atom.charge, float(charge))
156+
assert_almost_equal(foyer_atom.sigma, float(sigma))
157+
assert_almost_equal(parmed2trappe(foyer_atom.epsilon),
158+
float(epsilon),
159+
decimal=3)
160+
161+
@staticmethod
162+
def check_bonds(ref, struc):
163+
ref_bonds = ref['stretch']
164+
# organise foyer's bonds into dict of indices
165+
bond_dict = {}
166+
for b in struc.bonds:
167+
i, j = sorted((b.atom1.idx, b.atom2.idx))
168+
bond_dict[i+1, j+1] = b
169+
170+
# iterate reference file bonds
171+
# lookup foyer bond and compare length
172+
for ref_bond in ref_bonds:
173+
_, desc, _, length = ref_bond.strip().split(',')
174+
i, j = map(int, desc.strip("\"\'").split('-'))
175+
foyer_bond = bond_dict[i, j]
176+
177+
assert_almost_equal(foyer_bond.type.req, float(length))
178+
# trappe is fixed bond length, no point in K comparison
179+
180+
@staticmethod
181+
def check_angles(ref, struc):
182+
ref_angles = ref['bend']
183+
assert len(ref_angles) == len(struc.angles)
184+
185+
angle_dict = {}
186+
for a in struc.angles:
187+
i, j, k = a.atom1.idx, a.atom2.idx, a.atom3.idx
188+
if k < i:
189+
i, j, k = k, j, i
190+
angle_dict[i+1, j+1, k+1] = a
191+
192+
for ref_angle in ref_angles:
193+
_, desc, _, angle, k_angle = ref_angle.strip().split(',')
194+
i, j, k = map(int, desc.strip("\"\'").split('-'))
195+
196+
foyer_angle = angle_dict[i, j, k]
197+
198+
assert_almost_equal(foyer_angle.type.theteq,
199+
float(angle),
200+
decimal=4)
201+
assert_almost_equal(angle_parmed2trappe(foyer_angle.type.k),
202+
float(k_angle),
203+
decimal=3)
204+
205+
@staticmethod
206+
def check_dihedrals(ref, struc):
207+
ref_dihedrals = ref['torsion']
208+
209+
assert len(ref_dihedrals) == len(struc.rb_torsions)
210+
211+
dih_dict = {}
212+
for d in struc.rb_torsions:
213+
i, j, k, l = d.atom1.idx, d.atom2.idx, d.atom3.idx, d.atom4.idx
214+
215+
if l < i:
216+
i, j, k, l = l, k, j, i
217+
dih_dict[i+1, j+1, k+1, l+1] = d
218+
219+
for ref_d in ref_dihedrals:
220+
_, desc, _, *params = ref_d.strip().split(',')
221+
i, j, k, l = map(int, desc.strip("\"\'").split('-'))
222+
if l < i:
223+
i, j, k, l = l, k, j, i
224+
225+
foyer_d = dih_dict[i, j, k, l]
226+
227+
assert_dihedrals_match(foyer_d, params)
228+
229+
@pytest.mark.parametrize('trappe_ref,mol2file', [
230+
('butadiene/trappe_parameters_94.csv', 'butadiene'),
231+
('isoprene/trappe_parameters_95.csv', 'isoprene'),
232+
('methyl_acrylate/trappe_parameters_96.csv', 'methyl_acrylate'),
233+
], indirect=True)
234+
def test_trappe_reference(self, trappe_ref, mol2file):
235+
self.check_atoms(trappe_ref, mol2file)
236+
237+
self.check_bonds(trappe_ref, mol2file)
238+
239+
self.check_angles(trappe_ref, mol2file)
240+
241+
self.check_dihedrals(trappe_ref, mol2file)
242+
243+
55244
if __name__ == '__main__':
56245
TestTraPPE().find_correctly_implemented()

‎foyer/trappe_validation/aldehyde_1/aldehyde_1.mol2

+16
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
@<TRIPOS>MOLECULE
2+
RES
3+
3 2 1 0 1
4+
SMALL
5+
USER_CHARGES
6+
@<TRIPOS>CRYSIN
7+
5.0000 5.0000 5.0000 90.0000 90.0000 90.0000 1 1
8+
@<TRIPOS>ATOM
9+
1 _CH3 0.0000 0.0000 0.0000 CH3_sp3-Cald 1 RES -0.043000
10+
2 _HC 0.0000 0.0000 0.0000 CH_ald 1 RES 0.525000
11+
3 O 0.0000 0.0000 0.0000 O_ald 1 RES -0.482000
12+
@<TRIPOS>BOND
13+
1 1 2 1
14+
2 3 2 1
15+
@<TRIPOS>SUBSTRUCTURE
16+
1 RES 1 RESIDUE 0 1 ROOT 0

‎foyer/trappe_validation/aldehyde_2/aldehyde_2.mol2

+18
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
@<TRIPOS>MOLECULE
2+
RES
3+
4 3 1 0 1
4+
SMALL
5+
USER_CHARGES
6+
@<TRIPOS>CRYSIN
7+
5.0000 5.0000 5.0000 90.0000 90.0000 90.0000 1 1
8+
@<TRIPOS>ATOM
9+
1 _CH3 0.0000 0.0000 0.0000 CH3_sp3 1 RES 0.000000
10+
2 _CH2 0.0000 0.0000 0.0000 CH2_sp3-Cald 1 RES -0.043000
11+
3 _HC 0.0000 0.0000 0.0000 CH_ald 1 RES 0.525000
12+
4 O 0.0000 0.0000 0.0000 O_ald 1 RES -0.482000
13+
@<TRIPOS>BOND
14+
1 1 2 1
15+
2 2 3 1
16+
3 3 4 1
17+
@<TRIPOS>SUBSTRUCTURE
18+
1 RES 1 RESIDUE 0 1 ROOT 0

‎foyer/trappe_validation/aldehyde_3/aldehyde_3.mol2

+20
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
@<TRIPOS>MOLECULE
2+
RES
3+
5 4 1 0 1
4+
SMALL
5+
USER_CHARGES
6+
@<TRIPOS>CRYSIN
7+
5.0000 5.0000 5.0000 90.0000 90.0000 90.0000 1 1
8+
@<TRIPOS>ATOM
9+
1 _CH3 0.0000 0.0000 0.0000 CH3_sp3 1 RES 0.000000
10+
2 _CH3 0.0000 0.0000 0.0000 CH3_sp3 1 RES 0.000000
11+
3 _HC 0.0000 0.0000 0.0000 CH_sp3-Cald 1 RES -0.043000
12+
4 _HC 0.0000 0.0000 0.0000 CH_ald 1 RES 0.525000
13+
5 O 0.0000 0.0000 0.0000 O_ald 1 RES -0.482000
14+
@<TRIPOS>BOND
15+
1 2 3 1
16+
2 3 1 1
17+
3 4 3 1
18+
4 4 5 1
19+
@<TRIPOS>SUBSTRUCTURE
20+
1 RES 1 RESIDUE 0 1 ROOT 0

‎foyer/trappe_validation/aldehyde_4/aldehyde_4.mol2

+22
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
@<TRIPOS>MOLECULE
2+
RES
3+
6 5 1 0 1
4+
SMALL
5+
USER_CHARGES
6+
@<TRIPOS>CRYSIN
7+
5.0000 5.0000 5.0000 90.0000 90.0000 90.0000 1 1
8+
@<TRIPOS>ATOM
9+
1 _CH3 0.0000 0.0000 0.0000 CH3_sp3 1 RES 0.000000
10+
2 _CH3 0.0000 0.0000 0.0000 CH3_sp3 1 RES 0.000000
11+
3 _CH3 0.0000 0.0000 0.0000 CH3_sp3 1 RES 0.000000
12+
4 C 0.0000 0.0000 0.0000 C_sp3-Cald 1 RES -0.043000
13+
5 _HC 0.0000 0.0000 0.0000 CH_ald 1 RES 0.525000
14+
6 O 0.0000 0.0000 0.0000 O_ald 1 RES -0.482000
15+
@<TRIPOS>BOND
16+
1 1 4 1
17+
2 2 4 1
18+
3 3 4 1
19+
4 4 5 1
20+
5 5 6 1
21+
@<TRIPOS>SUBSTRUCTURE
22+
1 RES 1 RESIDUE 0 1 ROOT 0

‎foyer/trappe_validation/butadiene/butadiene.mol2

+18
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
@<TRIPOS>MOLECULE
2+
RES
3+
4 3 1 0 1
4+
SMALL
5+
NO_CHARGES
6+
@<TRIPOS>CRYSIN
7+
9.5000 5.0000 5.0000 90.0000 90.0000 90.0000 1 1
8+
@<TRIPOS>ATOM
9+
1 _CH2 0.0000 0.0000 0.0000 CH2_sp2 1 RES
10+
2 _HC 1.5000 0.0000 0.0000 CH_sp2-Csp2 1 RES
11+
3 _HC 3.0000 0.0000 0.0000 CH_sp2-Csp2 1 RES
12+
4 _CH2 4.5000 0.0000 0.0000 CH2_sp2 1 RES
13+
@<TRIPOS>BOND
14+
1 1 2 1
15+
2 2 3 1
16+
3 3 4 1
17+
@<TRIPOS>SUBSTRUCTURE
18+
1 RES 1 RESIDUE 0 1 ROOT 0

‎foyer/trappe_validation/butadiene/trappe_parameters_94.csv

+19
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
"1,3-butadiene"
2+
3+
#,(pseudo)atom,type,"epsilon/kB [K]","sigma [Ang.]","charge [e]"
4+
1,CH2,[CH2]=CHx,85.0,3.675,0.000
5+
4,CH2,[CH2]=CHx,85.0,3.675,0.000
6+
2,CH,CHx=[CH](sp2)-CHy(sp2),52.0,3.710,0.000
7+
3,CH,CHx=[CH](sp2)-CHy(sp2),52.0,3.710,0.000
8+
9+
#,stretch,type,"length [Ang.]"
10+
1,"'1 - 2'",CHx=CHy,1.3300
11+
2,"'2 - 3'",CHx-CHy,1.5400
12+
3,"'3 - 4'",CHx=CHy,1.3300
13+
14+
#,bend,type,"theta [degrees]","k_theta/kB [K/rad^2]"
15+
1,"'1 - 2 - 3'",CHx=(CH)-CHx,119.70,70420
16+
2,"'2 - 3 - 4'",CHx=(CH)-CHx,119.70,70420
17+
18+
#,torsion,type,"c0/kB [K]","c1/kB [K]","c2/kB [K]","c3/kB [K]","c4/kB [K]"
19+
1,"'1 - 2 - 3 - 4'",CH2=(CH)-(CH)=CH2,2034.58,531.57,-1239.35,460.04,196.38

‎foyer/trappe_validation/dimethylether/dme.mol2

+14
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
@<TRIPOS>MOLECULE
2+
DME
3+
3 2 1 0 1
4+
SMALL
5+
USER_CHARGES
6+
@<TRIPOS>ATOM
7+
1 _CH3 12.5360 1.0140 22.0490 CH3_eth 1 DME 0.265000
8+
2 O 13.5570 1.3980 22.8620 O_eth 1 DME 0.000000
9+
3 _CH3 13.5710 2.7160 23.1980 CH3_eth 1 DME 0.265000
10+
@<TRIPOS>BOND
11+
1 1 2 1
12+
2 2 3 1
13+
@<TRIPOS>SUBSTRUCTURE
14+
1 DME 1 RESIDUE 0 1 ROOT 0

‎foyer/trappe_validation/isoprene/isoprene.mol2

+20
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
@<TRIPOS>MOLECULE
2+
RES
3+
5 4 1 0 1
4+
SMALL
5+
NO_CHARGES
6+
@<TRIPOS>CRYSIN
7+
9.5000 6.5000 5.0000 90.0000 90.0000 90.0000 1 1
8+
@<TRIPOS>ATOM
9+
1 _CH2 3.0000 1.5000 0.0000 CH2_sp2 1 RES
10+
2 C 3.0000 0.0000 0.0000 C_sp2-Csp2 1 RES
11+
3 _HC 1.5000 0.0000 0.0000 CH_sp2-Csp2 1 RES
12+
4 _CH2 0.0000 0.0000 0.0000 CH2_sp2 1 RES
13+
5 _CH3 4.5000 0.0000 0.0000 CH3_sp3 1 RES
14+
@<TRIPOS>BOND
15+
1 2 1 1
16+
2 2 5 1
17+
3 3 2 1
18+
4 4 3 1
19+
@<TRIPOS>SUBSTRUCTURE
20+
1 RES 1 RESIDUE 0 1 ROOT 0

0 commit comments

Comments
 (0)