diff --git a/foyer/forcefields/trappe-ua.xml b/foyer/forcefields/trappe-ua.xml index cb700626..d94368c4 100644 --- a/foyer/forcefields/trappe-ua.xml +++ b/foyer/forcefields/trappe-ua.xml @@ -2,217 +2,253 @@ - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - + + + + + + + + + + - - - - + + + + + + + + + + + + + + + - - - + - - - - - - - - - - + + - - - - - - - - - - + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - + + + - - - - - - - - - - - - - - - - + + + + + + + + - - - - - - - - - - - - - - - - + + + + + + + + - - - - - - - - - - + + + + + + + + + + - - - - + + - - - - + - - - - + - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + + + + + diff --git a/foyer/tests/implemented_trappe_tests.txt b/foyer/tests/implemented_trappe_tests.txt index 0a666ee1..71dc0716 100644 --- a/foyer/tests/implemented_trappe_tests.txt +++ b/foyer/tests/implemented_trappe_tests.txt @@ -8,3 +8,9 @@ methane methanol neopentane propane +dimethylether +aldehyde_1 +aldehyde_2 +aldehyde_3 +aldehyde_4 +methyl_acrylate diff --git a/foyer/tests/test_trappe.py b/foyer/tests/test_trappe.py index eae113ed..4fdce9c6 100755 --- a/foyer/tests/test_trappe.py +++ b/foyer/tests/test_trappe.py @@ -1,7 +1,7 @@ import glob import itertools as it import os - +from scipy.constants import N_A, k, calorie import parmed as pmd from pkg_resources import resource_filename import pytest @@ -9,6 +9,8 @@ from foyer import Forcefield from foyer.tests.utils import atomtype +from numpy.testing import assert_almost_equal + TRAPPE_UA = Forcefield(name='trappe-ua') TRAPPE_TESTFILES_DIR = resource_filename('foyer', 'trappe_validation') @@ -45,12 +47,199 @@ def find_correctly_implemented(self): @pytest.mark.parametrize('mol_name', correctly_implemented) def test_atomtyping(self, mol_name, testfiles_dir=TRAPPE_TESTFILES_DIR): - files = glob.glob(os.path.join(testfiles_dir, mol_name, '*')) + files = glob.glob(os.path.join(testfiles_dir, mol_name, '*.mol2')) for mol_file in files: _, ext = os.path.splitext(mol_file) mol2_path = os.path.join(testfiles_dir, mol_name, mol_file) structure = pmd.load_file(mol2_path, structure=True) atomtype(structure, TRAPPE_UA, non_atomistic=True) +@pytest.fixture +def trappe_ref(request): + # load raw trappe csv file and process slightly + container = { + '(pseudo)atom': [], + 'stretch': [], + 'bend': [], + 'torsion': [], + } + # split csv file into sections + csv_file = os.path.join(TRAPPE_TESTFILES_DIR, + request.param) + with open(csv_file, 'r') as inf: + for line in inf: + if line.startswith('#'): + section_clue = line.split(",")[1] + section = container[section_clue] + line = next(inf) + while not line == '\n': + # slurp until newline hit + section.append(line) + try: + line = next(inf) + except StopIteration: + break + + return container + + +@pytest.fixture +def mol2file(request): + # load a mol2 file and apply foyer's trappe-ua forcefield + fname = request.param + molfile = glob.glob(os.path.join(TRAPPE_TESTFILES_DIR, fname, '*'))[0] + structure = pmd.load_file(molfile, structure=True) + + return TRAPPE_UA.apply(structure, + assert_angle_params=False, + assert_dihedral_params=False, + ) + + +def parmed2trappe(val): + # convert kcal/mol (parmed units) to K/k_B (trappe units) + return val * calorie / (k * N_A / 1000.) + +def angle_parmed2trappe(val): + # trappe angle units need doubling + return 2 * parmed2trappe(val) + +def assert_dihedrals_match(foyer_d, reference_params): + # Check that a dihedral from foyer matches the text in the csv file + # foyer_d - parmed rb_torsion + # reference_params - dihedral parameters from trappe csv file + def conv1(t): + # convert RB to original Trappe dihedral form: + # V = c0 + c1(1 + cos phi) + c1(1 - cos 2phi) + c3(1 + cos 3phi) + 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 second Trappe dihedral form: + # V = c0 + c1 cos phi + c2 cos 2phi + c3 cos 3phi + c4 cos 4phi + 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, has an extra column + foyer_params = conv2(foyer_d.type) + elif len(reference_params) == 4: + # using original dihedral form + foyer_params = conv1(foyer_d.type) + + assert_almost_equal(foyer_params, reference_params, decimal=3) + +class TestTrappeReferences(object): + @staticmethod + def check_atoms(ref, struc): + ref_atoms = sorted(ref['(pseudo)atom'], + key=lambda x: int(x.split(',')[0])) + for ref_atom, foyer_atom in zip(ref_atoms, struc.atoms): + _, _, _, epsilon, sigma, charge = ref_atom.strip().split(',') + + assert_almost_equal(foyer_atom.charge, float(charge)) + assert_almost_equal(foyer_atom.sigma, float(sigma)) + assert_almost_equal(parmed2trappe(foyer_atom.epsilon), + float(epsilon), + decimal=3) + + @staticmethod + def check_bonds(ref, struc): + ref_bonds = ref['stretch'] + # organise foyer's bonds into dict of indices + bond_dict = {} + for b in struc.bonds: + i, j = sorted((b.atom1.idx, b.atom2.idx)) + bond_dict[i+1, j+1] = b + + # iterate reference file bonds + # lookup foyer bond and compare length + for ref_bond in ref_bonds: + _, desc, _, length = ref_bond.strip().split(',') + i, j = map(int, desc.strip("\"\'").split('-')) + foyer_bond = bond_dict[i, j] + + assert_almost_equal(foyer_bond.type.req, float(length)) + # trappe is fixed bond length, no point in K comparison + + @staticmethod + def check_angles(ref, struc): + ref_angles = ref['bend'] + assert len(ref_angles) == len(struc.angles) + + angle_dict = {} + for a in struc.angles: + i, j, k = a.atom1.idx, a.atom2.idx, a.atom3.idx + if k < i: + i, j, k = k, j, i + angle_dict[i+1, j+1, k+1] = a + + for ref_angle in ref_angles: + _, desc, _, angle, k_angle = ref_angle.strip().split(',') + i, j, k = map(int, desc.strip("\"\'").split('-')) + + foyer_angle = angle_dict[i, j, k] + + assert_almost_equal(foyer_angle.type.theteq, + float(angle), + decimal=4) + assert_almost_equal(angle_parmed2trappe(foyer_angle.type.k), + float(k_angle), + decimal=3) + + @staticmethod + def check_dihedrals(ref, struc): + ref_dihedrals = ref['torsion'] + + 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('-')) + if l < i: + i, j, k, l = l, k, j, i + + foyer_d = dih_dict[i, j, k, l] + + assert_dihedrals_match(foyer_d, params) + + @pytest.mark.parametrize('trappe_ref,mol2file', [ + ('butadiene/trappe_parameters_94.csv', 'butadiene'), + ('isoprene/trappe_parameters_95.csv', 'isoprene'), + ('methyl_acrylate/trappe_parameters_96.csv', 'methyl_acrylate'), + ], indirect=True) + def test_trappe_reference(self, trappe_ref, mol2file): + self.check_atoms(trappe_ref, mol2file) + + self.check_bonds(trappe_ref, mol2file) + + self.check_angles(trappe_ref, mol2file) + + self.check_dihedrals(trappe_ref, mol2file) + + if __name__ == '__main__': TestTraPPE().find_correctly_implemented() diff --git a/foyer/trappe_validation/aldehyde_1/aldehyde_1.mol2 b/foyer/trappe_validation/aldehyde_1/aldehyde_1.mol2 new file mode 100644 index 00000000..2d85b21c --- /dev/null +++ b/foyer/trappe_validation/aldehyde_1/aldehyde_1.mol2 @@ -0,0 +1,16 @@ +@MOLECULE +RES +3 2 1 0 1 +SMALL +USER_CHARGES +@CRYSIN + 5.0000 5.0000 5.0000 90.0000 90.0000 90.0000 1 1 +@ATOM + 1 _CH3 0.0000 0.0000 0.0000 CH3_sp3-Cald 1 RES -0.043000 + 2 _HC 0.0000 0.0000 0.0000 CH_ald 1 RES 0.525000 + 3 O 0.0000 0.0000 0.0000 O_ald 1 RES -0.482000 +@BOND + 1 1 2 1 + 2 3 2 1 +@SUBSTRUCTURE + 1 RES 1 RESIDUE 0 1 ROOT 0 diff --git a/foyer/trappe_validation/aldehyde_2/aldehyde_2.mol2 b/foyer/trappe_validation/aldehyde_2/aldehyde_2.mol2 new file mode 100644 index 00000000..8ba38e62 --- /dev/null +++ b/foyer/trappe_validation/aldehyde_2/aldehyde_2.mol2 @@ -0,0 +1,18 @@ +@MOLECULE +RES +4 3 1 0 1 +SMALL +USER_CHARGES +@CRYSIN + 5.0000 5.0000 5.0000 90.0000 90.0000 90.0000 1 1 +@ATOM + 1 _CH3 0.0000 0.0000 0.0000 CH3_sp3 1 RES 0.000000 + 2 _CH2 0.0000 0.0000 0.0000 CH2_sp3-Cald 1 RES -0.043000 + 3 _HC 0.0000 0.0000 0.0000 CH_ald 1 RES 0.525000 + 4 O 0.0000 0.0000 0.0000 O_ald 1 RES -0.482000 +@BOND + 1 1 2 1 + 2 2 3 1 + 3 3 4 1 +@SUBSTRUCTURE + 1 RES 1 RESIDUE 0 1 ROOT 0 diff --git a/foyer/trappe_validation/aldehyde_3/aldehyde_3.mol2 b/foyer/trappe_validation/aldehyde_3/aldehyde_3.mol2 new file mode 100644 index 00000000..26d844ce --- /dev/null +++ b/foyer/trappe_validation/aldehyde_3/aldehyde_3.mol2 @@ -0,0 +1,20 @@ +@MOLECULE +RES +5 4 1 0 1 +SMALL +USER_CHARGES +@CRYSIN + 5.0000 5.0000 5.0000 90.0000 90.0000 90.0000 1 1 +@ATOM + 1 _CH3 0.0000 0.0000 0.0000 CH3_sp3 1 RES 0.000000 + 2 _CH3 0.0000 0.0000 0.0000 CH3_sp3 1 RES 0.000000 + 3 _HC 0.0000 0.0000 0.0000 CH_sp3-Cald 1 RES -0.043000 + 4 _HC 0.0000 0.0000 0.0000 CH_ald 1 RES 0.525000 + 5 O 0.0000 0.0000 0.0000 O_ald 1 RES -0.482000 +@BOND + 1 2 3 1 + 2 3 1 1 + 3 4 3 1 + 4 4 5 1 +@SUBSTRUCTURE + 1 RES 1 RESIDUE 0 1 ROOT 0 diff --git a/foyer/trappe_validation/aldehyde_4/aldehyde_4.mol2 b/foyer/trappe_validation/aldehyde_4/aldehyde_4.mol2 new file mode 100644 index 00000000..386cc424 --- /dev/null +++ b/foyer/trappe_validation/aldehyde_4/aldehyde_4.mol2 @@ -0,0 +1,22 @@ +@MOLECULE +RES +6 5 1 0 1 +SMALL +USER_CHARGES +@CRYSIN + 5.0000 5.0000 5.0000 90.0000 90.0000 90.0000 1 1 +@ATOM + 1 _CH3 0.0000 0.0000 0.0000 CH3_sp3 1 RES 0.000000 + 2 _CH3 0.0000 0.0000 0.0000 CH3_sp3 1 RES 0.000000 + 3 _CH3 0.0000 0.0000 0.0000 CH3_sp3 1 RES 0.000000 + 4 C 0.0000 0.0000 0.0000 C_sp3-Cald 1 RES -0.043000 + 5 _HC 0.0000 0.0000 0.0000 CH_ald 1 RES 0.525000 + 6 O 0.0000 0.0000 0.0000 O_ald 1 RES -0.482000 +@BOND + 1 1 4 1 + 2 2 4 1 + 3 3 4 1 + 4 4 5 1 + 5 5 6 1 +@SUBSTRUCTURE + 1 RES 1 RESIDUE 0 1 ROOT 0 diff --git a/foyer/trappe_validation/butadiene/butadiene.mol2 b/foyer/trappe_validation/butadiene/butadiene.mol2 new file mode 100644 index 00000000..d1ce61ec --- /dev/null +++ b/foyer/trappe_validation/butadiene/butadiene.mol2 @@ -0,0 +1,18 @@ +@MOLECULE +RES +4 3 1 0 1 +SMALL +NO_CHARGES +@CRYSIN + 9.5000 5.0000 5.0000 90.0000 90.0000 90.0000 1 1 +@ATOM + 1 _CH2 0.0000 0.0000 0.0000 CH2_sp2 1 RES + 2 _HC 1.5000 0.0000 0.0000 CH_sp2-Csp2 1 RES + 3 _HC 3.0000 0.0000 0.0000 CH_sp2-Csp2 1 RES + 4 _CH2 4.5000 0.0000 0.0000 CH2_sp2 1 RES +@BOND + 1 1 2 1 + 2 2 3 1 + 3 3 4 1 +@SUBSTRUCTURE + 1 RES 1 RESIDUE 0 1 ROOT 0 diff --git a/foyer/trappe_validation/butadiene/trappe_parameters_94.csv b/foyer/trappe_validation/butadiene/trappe_parameters_94.csv new file mode 100644 index 00000000..3bcb5dc3 --- /dev/null +++ b/foyer/trappe_validation/butadiene/trappe_parameters_94.csv @@ -0,0 +1,19 @@ +"1,3-butadiene" + +#,(pseudo)atom,type,"epsilon/kB [K]","sigma [Ang.]","charge [e]" +1,CH2,[CH2]=CHx,85.0,3.675,0.000 +4,CH2,[CH2]=CHx,85.0,3.675,0.000 +2,CH,CHx=[CH](sp2)-CHy(sp2),52.0,3.710,0.000 +3,CH,CHx=[CH](sp2)-CHy(sp2),52.0,3.710,0.000 + +#,stretch,type,"length [Ang.]" +1,"'1 - 2'",CHx=CHy,1.3300 +2,"'2 - 3'",CHx-CHy,1.5400 +3,"'3 - 4'",CHx=CHy,1.3300 + +#,bend,type,"theta [degrees]","k_theta/kB [K/rad^2]" +1,"'1 - 2 - 3'",CHx=(CH)-CHx,119.70,70420 +2,"'2 - 3 - 4'",CHx=(CH)-CHx,119.70,70420 + +#,torsion,type,"c0/kB [K]","c1/kB [K]","c2/kB [K]","c3/kB [K]","c4/kB [K]" +1,"'1 - 2 - 3 - 4'",CH2=(CH)-(CH)=CH2,2034.58,531.57,-1239.35,460.04,196.38 diff --git a/foyer/trappe_validation/dimethylether/dme.mol2 b/foyer/trappe_validation/dimethylether/dme.mol2 new file mode 100644 index 00000000..95ac71b1 --- /dev/null +++ b/foyer/trappe_validation/dimethylether/dme.mol2 @@ -0,0 +1,14 @@ +@MOLECULE +DME +3 2 1 0 1 +SMALL +USER_CHARGES +@ATOM + 1 _CH3 12.5360 1.0140 22.0490 CH3_eth 1 DME 0.265000 + 2 O 13.5570 1.3980 22.8620 O_eth 1 DME 0.000000 + 3 _CH3 13.5710 2.7160 23.1980 CH3_eth 1 DME 0.265000 +@BOND + 1 1 2 1 + 2 2 3 1 +@SUBSTRUCTURE + 1 DME 1 RESIDUE 0 1 ROOT 0 diff --git a/foyer/trappe_validation/isoprene/isoprene.mol2 b/foyer/trappe_validation/isoprene/isoprene.mol2 new file mode 100644 index 00000000..3695102c --- /dev/null +++ b/foyer/trappe_validation/isoprene/isoprene.mol2 @@ -0,0 +1,20 @@ +@MOLECULE +RES +5 4 1 0 1 +SMALL +NO_CHARGES +@CRYSIN + 9.5000 6.5000 5.0000 90.0000 90.0000 90.0000 1 1 +@ATOM + 1 _CH2 3.0000 1.5000 0.0000 CH2_sp2 1 RES + 2 C 3.0000 0.0000 0.0000 C_sp2-Csp2 1 RES + 3 _HC 1.5000 0.0000 0.0000 CH_sp2-Csp2 1 RES + 4 _CH2 0.0000 0.0000 0.0000 CH2_sp2 1 RES + 5 _CH3 4.5000 0.0000 0.0000 CH3_sp3 1 RES +@BOND + 1 2 1 1 + 2 2 5 1 + 3 3 2 1 + 4 4 3 1 +@SUBSTRUCTURE + 1 RES 1 RESIDUE 0 1 ROOT 0 diff --git a/foyer/trappe_validation/isoprene/trappe_parameters_95.csv b/foyer/trappe_validation/isoprene/trappe_parameters_95.csv new file mode 100644 index 00000000..3771b338 --- /dev/null +++ b/foyer/trappe_validation/isoprene/trappe_parameters_95.csv @@ -0,0 +1,23 @@ +isoprene + +#,(pseudo)atom,type,"epsilon/kB [K]","sigma [Ang.]","charge [e]" +5,CH3,[CH3]-CHx,98.0,3.750,0.000 +1,CH2,[CH2]=CHx,85.0,3.675,0.000 +4,CH2,[CH2]=CHx,85.0,3.675,0.000 +3,CH,CHx=[CH](sp2)-CHy(sp2),52.0,3.710,0.000 +2,C,CHx=[C](sp2)-CHy(sp2),22.0,3.850,0.000 + +#,stretch,type,"length [Ang.]" +1,"'1 - 2'",CHx=CHy,1.3300 +2,"'2 - 3'",CHx-CHy,1.5400 +3,"'2 - 5'",CHx-CHy,1.5400 +4,"'3 - 4'",CHx=CHy,1.3300 + +#,bend,type,"theta [degrees]","k_theta/kB [K/rad^2]" +1,"'1 - 2 - 3'",CHx=(C)-CHy,119.70,70420 +2,"'1 - 2 - 5'",CHx=(C)-CHy,119.70,70420 +3,"'2 - 3 - 4'",CHx=(CH)-CHx,119.70,70420 +4,"'3 - 2 - 5'",CHx-(C)[=CHz]-CHy,119.70,70420 + +#,torsion,type,"c0/kB [K]","c1/kB [K]","c2/kB [K]","c3/kB [K]","c4/kB [K]" +2,"'5 - 2 - 3 - 4'",CH2=(CH)-(C)[=CHx]-CH3,1861.27,-349.97,-1048.70,-580.54,117.92 diff --git a/foyer/trappe_validation/methyl_acrylate/methyl_acrylate.mol2 b/foyer/trappe_validation/methyl_acrylate/methyl_acrylate.mol2 new file mode 100644 index 00000000..d3e6d361 --- /dev/null +++ b/foyer/trappe_validation/methyl_acrylate/methyl_acrylate.mol2 @@ -0,0 +1,22 @@ +@MOLECULE +RES +6 5 1 0 1 +SMALL +USER_CHARGES +@CRYSIN + 11.0000 6.5000 5.0000 90.0000 90.0000 90.0000 1 1 +@ATOM + 1 _CH2 6.0000 0.0000 0.0000 CH2_sp2 1 RES 0.000000 + 2 _HC 4.5000 0.0000 0.0000 CH_sp2-Csp2 1 RES 0.000000 + 3 C 3.0000 0.0000 0.0000 C_acryl 1 RES 0.400000 + 4 O 1.5000 0.0000 0.0000 O_acryl_sp3 1 RES -0.250000 + 5 O 3.0000 1.5000 0.0000 O_acryl_sp2 1 RES -0.400000 + 6 _CH3 0.0000 0.0000 0.0000 CH3_eth 1 RES 0.250000 +@BOND + 1 2 3 1 + 2 2 1 1 + 3 3 4 1 + 4 5 3 1 + 5 6 4 1 +@SUBSTRUCTURE + 1 RES 1 RESIDUE 0 1 ROOT 0 diff --git a/foyer/trappe_validation/methyl_acrylate/trappe_parameters_96.csv b/foyer/trappe_validation/methyl_acrylate/trappe_parameters_96.csv new file mode 100644 index 00000000..c27e5552 --- /dev/null +++ b/foyer/trappe_validation/methyl_acrylate/trappe_parameters_96.csv @@ -0,0 +1,29 @@ +"methyl acrylate" + +#,(pseudo)atom,type,"epsilon/kB [K]","sigma [Ang.]","charge [e]" +6,CH3,[CH3]-O-CHx,98.0,3.750,0.250 +1,CH2,[CH2]=CHx,85.0,3.675,0.000 +3,C,CHx-O-[C]=O,40.0,3.820,0.400 +5,O,CHx-O-C=[O],79.0,3.050,-0.400 +4,O,CHx-[O]-C=O,55.0,2.800,-0.250 +2,CH,CHx=[CH](sp2)-CHy(sp2),52.0,3.710,0.000 + +#,stretch,type,"length [Ang.]" +1,"'1 - 2'",CHx=CHy,1.3300 +2,"'2 - 3'",CHx(sp2)-(CHy=O),1.5200 +3,"'3 - 4'",O-(C=O),1.3440 +4,"'3 - 5'",OC=O,1.2000 +5,"'4 - 6'",CHx-O,1.4100 + +#,bend,type,"theta [degrees]","k_theta/kB [K/rad^2]" +1,"'1 - 2 - 3'",CHx=(CH)-CHx,119.70,70420 +2,"'2 - 3 - 4'",O-(C)[=O]-CHx,110.00,70600 +3,"'2 - 3 - 5'",O=(C)[-OCHx]-CHx,125.00,62500 +4,"'3 - 4 - 6'",C[=O]-(O)-CHx,115.00,62500 +5,"'4 - 3 - 5'",O=(C)-OCHx,125.00,62500 + +#,torsion,type,"c0/kB [K]","c1/kB [K]","c2/kB [K]","c3/kB [K]","c4/kB [K]" +1,"'1 - 2 - 3 - 4'",CH2=(CH)-(C)[=O]-O,823.03,47.91,-773.13,1.99,0.00 +2,"'1 - 2 - 3 - 5'",CH2=(CH)-(C)[-O]=O,823.03,-47.91,-773.13,-1.99,0.00 +3,"'2 - 3 - 4 - 6'",CHx-(O)-(C)(sp2)-C(sp2),1820.74,417.41,-1373.14,30.19,0.00 +4,"'5 - 3 - 4 - 6'",CHx-(O)-(C)=O,1820.74,-417.41,-1373.14,-30.19,0.00