Fix parsing of some CHARMM parameter files defining Amber force field. Also fix

assignment of 1-4 electrostatic scaling factor when it is on the NONBONDED line
(rather than the line after). Add tests for the above. Ignore .coverage file.
Improve test coverage of charmm/parameters.py

.gitignore

@@ -4,3 +4,4 @@
 *~
 build
 a.out
+.coverage

parmed/charmm/parameters.py

@@ -33,6 +33,15 @@ class _EmptyStringIterator(object):
     def __getitem__(self, idx):
         return ''
 
+def _typeconv(name):
+    if isinstance(name, integer_types):
+        return name
+    if name.upper() == name:
+        return name
+    # Lowercase letters present -- decorate the type name with LTU --
+    # Lower To Upper
+    return '%sLTU' % name.upper()
+
 class CharmmParameterSet(ParameterSet):
     """
     Stores a parameter set defined by CHARMM files. It stores the equivalent of
@@ -171,14 +180,6 @@ class CharmmParameterSet(ParameterSet):
             upper-case
         """
         new_params = cls()
-        def typeconv(name):
-            if isinstance(name, integer_types):
-                return name
-            if name.upper() == name:
-                return name
-            # Lowercase letters present -- decorate the type name with LTU --
-            # Lower To Upper
-            return '%sLTU' % name.upper()
         if copy:
             do_copy = lambda x: _copy(x)
         else:
@@ -187,9 +188,9 @@ class CharmmParameterSet(ParameterSet):
         id_typemap = dict()
         def copy_paramtype(key, typ, dict):
             if isinstance(key, string_types):
-                key = typeconv(key)
+                key = _typeconv(key)
             elif isinstance(key, tuple):
-                key = tuple(typeconv(k) for k in key)
+                key = tuple(_typeconv(k) for k in key)
             # NoUreyBradley should never be copied
             if typ is NoUreyBradley:
                 dict[key] = NoUreyBradley
@@ -201,13 +202,13 @@ class CharmmParameterSet(ParameterSet):
                 dict[key] = newtype
 
         for key, atom_type in iteritems(params.atom_types_tuple):
-            atom_type.name = typeconv(atom_type.name)
+            atom_type.name = _typeconv(atom_type.name)
             copy_paramtype(key, atom_type, new_params.atom_types_tuple)
         for typename, atom_type in iteritems(params.atom_types):
-            atom_type.name = typeconv(atom_type.name)
+            atom_type.name = _typeconv(atom_type.name)
             copy_paramtype(typename, atom_type, new_params.atom_types)
         for idx, atom_type in iteritems(params.atom_types_int):
-            atom_type.name = typeconv(atom_type.name)
+            atom_type.name = _typeconv(atom_type.name)
             copy_paramtype(idx, atom_type, new_params.atom_types_int)
 
         for key, typ in iteritems(params.bond_types):
@@ -224,10 +225,9 @@ class CharmmParameterSet(ParameterSet):
             copy_paramtype(key, typ, new_params.improper_types)
         for key, typ in iteritems(params.cmap_types):
             if len(key) == 8:
+                key = (key[0], key[1], key[2], key[3], key[7])
                 copy_paramtype(key, typ, new_params.cmap_types)
             elif len(key) == 5:
-                key = (key[0], key[1], key[2], key[3],
-                       key[1], key[2], key[3], key[4])
                 copy_paramtype(key, typ, new_params.cmap_types)
         for key, typ in iteritems(params.nbfix_types):
             copy_paramtype(key, typ, new_params.nbfix_types)
@@ -302,7 +302,7 @@ class CharmmParameterSet(ParameterSet):
             inst.read_topology_file(tfile)
         if pfile is not None:
             inst.read_parameter_file(pfile)
-        if isinstance(sfiles, str):
+        if isinstance(sfiles, string_types):
             # The API docstring requests a list, but allow for users to pass a
             # string with a single filename instead
             inst.read_stream_file(sfiles)
@@ -370,13 +370,13 @@ class CharmmParameterSet(ParameterSet):
             if line.startswith('BOND'):
                 section = 'BONDS'
                 continue
-            if line.startswith('ANGLE'):
+            if line.startswith('ANGLE') or line.startswith('THETA'):
                 section = 'ANGLES'
                 continue
-            if line.startswith('DIHEDRAL'):
+            if line.startswith('DIHEDRAL') or line.startswith('PHI'):
                 section = 'DIHEDRALS'
                 continue
-            if line.startswith('IMPROPER'):
+            if line.startswith('IMPROPER') or line.startswith('IMPHI'):
                 section = 'IMPROPER'
                 continue
             if line.startswith('CMAP'):
@@ -402,6 +402,11 @@ class CharmmParameterSet(ParameterSet):
                             # as the one we specified before
                             if abs(self.dihedral_types[0][0].scee-scee) > TINY:
                                 raise CharmmError('Inconsistent 1-4 scalings')
+                        else:
+                            scee = 1 / scee
+                            for key, dtl in iteritems(self.dihedral_types):
+                                for dt in dtl:
+                                    dt.scee = scee
                     elif word.upper().startswith('GEOM'):
                         if (self._declared_nbrules and
                                 self.combining_rule != 'geometric'):

test/test_parmed_charmm.py

@@ -3,14 +3,17 @@ Tests for the parmed/charmm subpackage
 """
 from __future__ import division, print_function
 
+from collections import OrderedDict
 import numpy as np
+import os
+import parmed as pmd
 from parmed.utils.six import iteritems, string_types
 from parmed.utils.six.moves import StringIO
 from parmed.charmm import charmmcrds, parameters, psf
 from parmed.charmm._charmmfile import CharmmFile, CharmmStreamFile
 from parmed import exceptions, topologyobjects as to, load_file, ParameterSet
 import parmed.unit as u
-import os
+import random
 import unittest
 import utils
 from utils import HAS_GROMACS
@@ -423,15 +426,60 @@ class TestCharmmPsf(unittest.TestCase):
 
 class TestCharmmParameters(utils.FileIOTestCase):
     """ Test CHARMM Parameter file parsing """
-    
+
+    def testPrivateFunctions(self):
+        """ Tests private helper functions for CharmmParameterSet """
+        # EmptyStringIterator
+        si = parameters._EmptyStringIterator()
+        it = iter(si)
+        for i in range(random.randint(100, 1000)):
+            self.assertEqual(next(it), '')
+        self.assertEqual(si[random.randint(0, 10000)], '')
+        # _typeconv
+        randint = random.randint(0, 100000)
+        self.assertEqual(parameters._typeconv(randint), randint)
+        self.assertEqual(parameters._typeconv('NOCHNG'), 'NOCHNG')
+        self.assertEqual(parameters._typeconv('NoCh'), 'NOCHLTU')
+
+    def testE14FAC(self):
+        """ Test reading CHARMM parameter files with 1-4 EEL scaling """
+        params = parameters.CharmmParameterSet(
+                get_fn('parm14sb_all.prm'),
+        )
+        for i, tortype in iteritems(params.dihedral_types):
+            for typ in tortype:
+                self.assertAlmostEqual(typ.scee, 1.2)
+        params = parameters.CharmmParameterSet(
+                get_fn('parm14sb_all_2.prm')
+        )
+        for i, tortype in iteritems(params.dihedral_types):
+            for typ in tortype:
+                self.assertAlmostEqual(typ.scee, 1.2)
+
     def testSingleParameterset(self):
         """ Test reading a single parameter set """
+        # Make sure we error if trying to load parameters before topology
         self.assertRaises(RuntimeError, lambda: parameters.CharmmParameterSet(
                                                 get_fn('par_all22_prot.inp')))
-        params = parameters.CharmmParameterSet(
+        # Test error handling for loading files with unsupported extensions
+        self.assertRaises(ValueError, lambda:
+                parameters.CharmmParameterSet(get_fn('trx.prmtop'))
+        )
+        self.assertRaises(ValueError, lambda: parameters.CharmmParameterSet('x.inp'))
+        self._check_single_paramset(
+                parameters.CharmmParameterSet(
                                 get_fn('top_all22_prot.inp'),
                                 get_fn('par_all22_prot.inp'),
+                )
+        )
+        self._check_single_paramset(
+                parameters.CharmmParameterSet.load_set(
+                                tfile=get_fn('top_all22_prot.inp'),
+                                pfile=get_fn('par_all22_prot.inp'),
+                )
         )
+
+    def _check_single_paramset(self, params):
         for i, tup in enumerate(params.atom_types_tuple):
             name, num = tup
             self.assertTrue(params.atom_types_tuple[tup] is
@@ -502,6 +550,11 @@ class TestCharmmParameters(utils.FileIOTestCase):
     def testParamFileOnly(self):
         """ Test reading only a parameter file with no RTF (CHARMM36) """
         parameters.CharmmParameterSet(get_fn('par_all36_carb.prm')).condense()
+        # Make sure read_parameter_file can accept a list of lines *without*
+        # comments
+        with CharmmFile(get_fn('par_all36_carb.prm'), 'r') as f:
+            params = parameters.CharmmParameterSet()
+            params.read_parameter_file(f.readlines())
 
     def testCollection(self):
         """ Test reading a large number of parameter files """
@@ -543,10 +596,18 @@ class TestCharmmParameters(utils.FileIOTestCase):
                                 get_fn('test.par', written=True)
         )
         params3 = parameters.CharmmParameterSet(get_fn('test.str', written=True))
+        params4 = parameters.CharmmParameterSet.load_set(
+                sfiles=get_fn('test.str', written=True)
+        )
+        params5 = parameters.CharmmParameterSet.load_set(
+                sfiles=[get_fn('test.str', written=True)]
+        )
 
         # Check that all of the params are equal
         self._compare_paramsets(params, params2, copy=True)
         self._compare_paramsets(params, params3, copy=True)
+        self._compare_paramsets(params, params4, copy=True)
+        self._compare_paramsets(params, params5, copy=True)
 
     def testCGenFF(self):
         """ Test parsing stream files generated by CGenFF """
@@ -640,6 +701,56 @@ class TestCharmmParameters(utils.FileIOTestCase):
         for name in chparams3.atom_types:
             self.assertTrue(name.endswith('LTU'))
 
+        # Load a parameter set with NoUreyBradley to make sure it's retained as
+        # a singleton. Also build a list of atom type tuples
+        for i, (typstr, typ) in enumerate(iteritems(params1.atom_types)):
+            params1.atom_types_tuple[(typstr, i+1)] = typ
+            params1.atom_types_int[i+1] = typ
+        for key in params1.angle_types:
+            params1.urey_bradley_types[key] = to.NoUreyBradley
+        chparams = parameters.CharmmParameterSet.from_parameterset(params1)
+        for _, item in iteritems(chparams.urey_bradley_types):
+            self.assertIs(item, to.NoUreyBradley)
+        for (typstr1, typ1), (typstr2, typ2) in zip(iteritems(params1.atom_types),
+                                                    iteritems(params1.atom_types)):
+            self.assertEqual(typstr1, typstr2)
+            self.assertEqual(typ1, typ2)
+        for (typstr1, typ1), (typstr2, typ2) in zip(iteritems(params1.atom_types_int),
+                                                    iteritems(params1.atom_types_int)):
+            self.assertEqual(typstr1, typstr2)
+            self.assertEqual(typ1, typ2)
+        for (typstr1, typ1), (typstr2, typ2) in zip(iteritems(params1.atom_types_tuple),
+                                                    iteritems(params1.atom_types_tuple)):
+            self.assertEqual(typstr1, typstr2)
+            self.assertEqual(typ1, typ2)
+
+        # Convert from the GROMACS topology file parameter set to
+        # CharmmParameterSet when loaded from the CHARMM force field
+        gmx = pmd.gromacs.GromacsTopologyFile(
+                os.path.join(pmd.gromacs.GROMACS_TOPDIR,
+                             'charmm27.ff', 'forcefield.itp')
+        )
+        from_gmx = parameters.CharmmParameterSet.from_parameterset(gmx.parameterset)
+        gmx = pmd.gromacs.GromacsTopologyFile(
+                os.path.join(pmd.gromacs.GROMACS_TOPDIR,
+                             'charmm27.ff', 'forcefield.itp')
+        )
+        # Make sure it can handle 8-key CMAPs
+        types = OrderedDict()
+        for key, typ in iteritems(gmx.parameterset.cmap_types):
+            assert len(key) == 5, 'Unexpected cmap key length'
+            types[(key[0], key[1], key[2], key[3],
+                   key[1], key[2], key[3], key[4])] = typ
+        gmx.parameterset.cmap_types = types
+        gmx.parameterset.nbfix_types[('X', 'Y')] = (2.0, 3.0)
+        from_gmx2 = parameters.CharmmParameterSet.from_parameterset(gmx.parameterset)
+        for (key1, typ1), (key2, typ2) in zip(iteritems(from_gmx.cmap_types),
+                                              iteritems(from_gmx2.cmap_types)):
+            self.assertEqual(key1, key2)
+            self.assertEqual(typ1, typ2)
+        self.assertEqual(len(from_gmx2.nbfix_types), 1)
+        self.assertEqual(from_gmx2.nbfix_types[('X', 'Y')], (2.0, 3.0))
+
     def _check_uppercase_types(self, params):
         for aname, atom_type in iteritems(params.atom_types):
             self.assertEqual(aname, aname.upper())