Update the CHARMM API to be consistent with the rest of OpenMM.

wrappers/python/simtk/openmm/app/charmmparameterset.py

@@ -220,7 +220,7 @@ class CharmmParameterSet(object):
                     atomic_number = get_by_symbol(elem).atomic_number
                 except (IndexError, KeyError):
                     # Figure it out from the mass
-                    atomic_number = element_by_mass(mass).atomic_number
+                    atomic_number = Element.getByMass(mass).atomic_number
                 atype = AtomType(name=name, number=idx, mass=mass,
                                  atomic_number=atomic_number)
                 self.atom_types_str[atype.name] = atype
@@ -441,7 +441,7 @@ class CharmmParameterSet(object):
                 atomic_number = get_by_symbol(elem).atomic_number
             except (IndexError, KeyError):
                 # Figure it out from the mass
-                atomic_number = element_by_mass(mass).atomic_number
+                atomic_number = Element.getByMass(mass).atomic_number
             atype = AtomType(name=name, number=idx, mass=mass,
                              atomic_number=atomic_number)
             self.atom_types_str[atype.name] = atype
@@ -488,7 +488,7 @@ class CharmmParameterSet(object):
             time.
 
         Example:
-        >>> params = CharmmParameterSet.loadSet(pfile='charmm.prm').condense()
+        >>> params = CharmmParameterSet('charmm.prm').condense()
         """
         # First scan through all of the bond types
         self._condense_types(self.bond_types)
@@ -529,18 +529,3 @@ class CharmmParameterSet(object):
                     typedict[key2] = typedict[key1]
 
 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-def element_by_mass(mass):
-    """ Determines what element the given atom is based on its mass """
-
-    diff = mass
-    best_guess = 'EP'
-
-    for key in Element._elements_by_atomic_number:
-        element = Element._elements_by_atomic_number[key]
-        massdiff = abs(element.mass.value_in_unit(u.daltons) - mass)
-        if massdiff < diff:
-            best_guess = element
-            diff = massdiff
-
-    return best_guess

wrappers/python/simtk/openmm/app/charmmpsffile.py

@@ -53,7 +53,7 @@ def _catchindexerror(func):
 class CharmmPsfFile(object):
     """
     A chemical structure instantiated from CHARMM files. You can instantiate a
-    ProteinStructure from a PSF file using the load_from_psf constructor
+    CharmmPsfFile from a PSF file using the load_from_psf constructor
 
     Example:
     >>> cs = CharmmPsfFile("testfiles/test.psf")
@@ -80,7 +80,7 @@ class CharmmPsfFile(object):
     etc.)
 
     Example:
-    >>> cs = ProteinStructure.load_from_psf("testfiles/test.psf")
+    >>> cs = CharmmPsfFile("testfiles/test.psf")
     >>> len(cs.atom_list)
     33
     >>> len(cs.bond_list)
@@ -102,8 +102,8 @@ class CharmmPsfFile(object):
     @_catchindexerror
     def __init__(self, psf_name):
         """
-        Opens and parses a PSF file, then instantiates a ProteinStructure
-        instance from the data.
+        Opens and parses a PSF file, then instantiates a CharmmPsfFile instance
+        from the data.
             
         Parameters:
             psf_name (str) : Name of the PSF file (it must exist)
@@ -112,7 +112,7 @@ class CharmmPsfFile(object):
             IOError : If file "psf_name" does not exist
             CharmmPSFError: If any parsing errors are encountered
         """
-        conv = ProteinStructure._convert
+        conv = CharmmPsfFile._convert
         # Make sure the file exists
         if not os.path.exists(psf_name):
             raise IOError('Could not find PSF file %s' % psf_name)
@@ -164,7 +164,7 @@ class CharmmPsfFile(object):
         # Eat the next line
         psf.readline()
         # Now get the number of bonds
-        nbond, holder = ProteinStructure._parse_psf_section(psf, int)
+        nbond, holder = CharmmPsfFile._parse_psf_section(psf, int)
         bond_list = TrackedList()
         if len(holder) != nbond * 2:
             raise CharmmPSFError('Got %d indexes for %d bonds' %
@@ -175,7 +175,7 @@ class CharmmPsfFile(object):
             bond_list.append(Bond(atom_list[id1], atom_list[id2]))
         bond_list.changed = False
         # Now get the number of angles and the angle list
-        ntheta, holder = ProteinStructure._parse_psf_section(psf, int)
+        ntheta, holder = CharmmPsfFile._parse_psf_section(psf, int)
         angle_list = TrackedList()
         if len(holder) != ntheta * 3:
             raise CharmmPSFError('Got %d indexes for %d angles' %
@@ -189,7 +189,7 @@ class CharmmPsfFile(object):
             )
         angle_list.changed = False
         # Now get the number of torsions and the torsion list
-        nphi, holder = ProteinStructure._parse_psf_section(psf, int)
+        nphi, holder = CharmmPsfFile._parse_psf_section(psf, int)
         dihedral_list = TrackedList()
         if len(holder) != nphi * 4:
             raise CharmmPSFError('Got %d indexes for %d torsions' %
@@ -205,7 +205,7 @@ class CharmmPsfFile(object):
             )
         dihedral_list.changed = False
         # Now get the number of improper torsions
-        nimphi, holder = ProteinStructure._parse_psf_section(psf, int)
+        nimphi, holder = CharmmPsfFile._parse_psf_section(psf, int)
         improper_list = TrackedList()
         if len(holder) != nimphi * 4:
             raise CharmmPSFError('Got %d indexes for %d impropers' %
@@ -221,7 +221,7 @@ class CharmmPsfFile(object):
             )
         improper_list.changed = False
         # Now handle the donors (what is this used for??)
-        ndon, holder = ProteinStructure._parse_psf_section(psf, int)
+        ndon, holder = CharmmPsfFile._parse_psf_section(psf, int)
         donor_list = TrackedList()
         if len(holder) != ndon * 2:
             raise CharmmPSFError('Got %d indexes for %d donors' %
@@ -232,7 +232,7 @@ class CharmmPsfFile(object):
             donor_list.append(AcceptorDonor(atom_list[id1], atom_list[id2]))
         donor_list.changed = False
         # Now handle the acceptors (what is this used for??)
-        nacc, holder = ProteinStructure._parse_psf_section(psf, int)
+        nacc, holder = CharmmPsfFile._parse_psf_section(psf, int)
         acceptor_list = TrackedList()
         if len(holder) != nacc * 2:
             raise CharmmPSFError('Got %d indexes for %d acceptors' %
@@ -244,9 +244,9 @@ class CharmmPsfFile(object):
         acceptor_list.changed = False
         # Now get the NNB section. Not sure what this section is for or what it
         # does...
-        nnb, holder = ProteinStructure._parse_psf_section(psf, int)
+        nnb, holder = CharmmPsfFile._parse_psf_section(psf, int)
         # Now get the group sections
-        pointers, holder = ProteinStructure._parse_psf_section(psf, int)
+        pointers, holder = CharmmPsfFile._parse_psf_section(psf, int)
         group_list = TrackedList()
         try:
             ngrp, nst2 = pointers
@@ -273,7 +273,7 @@ class CharmmPsfFile(object):
         # and the number of entries in the group. If the # of entries is
         # NATOM, assume we have MOLNT section. Warn if the MOLNT section is
         # 'wrong'...
-        pointer, holder = ProteinStructure._parse_psf_section(psf, int)
+        pointer, holder = CharmmPsfFile._parse_psf_section(psf, int)
 
         # Assign all of the atoms to molecules recursively
         set_molecules(atom_list)
@@ -292,7 +292,7 @@ class CharmmPsfFile(object):
                                      'section.')
             psf.readline() # blank
             # Now we get to the cross-term section
-            ncrterm, holder = ProteinStructure._parse_psf_section(psf, int)
+            ncrterm, holder = CharmmPsfFile._parse_psf_section(psf, int)
         else:
             ncrterm = pointer
         # At this point, ncrterm and holder are both set to the CMAP list for
@@ -370,7 +370,7 @@ class CharmmPsfFile(object):
             - data (list) : A list of all data in the parsed section converted
                     to `dtype'
         """
-        conv = ProteinStructure._convert
+        conv = CharmmPsfFile._convert
         words = psf.readline().split()
         if len(words) == 1:
             pointers = conv(words[0], int, 'pointer')
@@ -400,7 +400,7 @@ class CharmmPsfFile(object):
             - vmd (bool) : If True, it will write out a PSF in the format that
                     VMD prints it in (i.e., no NUMLP/NUMLPH or MOLNT sections)
         Example:
-            >>> cs = ProteinStructure.load_from_psf('testfiles/test.psf')
+            >>> cs = CharmmPsfFile.load_from_psf('testfiles/test.psf')
             >>> cs.write_psf('testfiles/test2.psf')
         """
         # See if this is an extended format
@@ -769,8 +769,13 @@ class CharmmPsfFile(object):
             if atom.residue.idx != last_residue:
                 last_residue = atom.residue.idx
                 residue = topology.addResidue(atom.residue.resname, chain)
-            atomic_num = atom.type.atomic_number
-            elem = element.Element.getByAtomicNumber(atomic_num)
+            if atom.type is not None:
+                # This is the most reliable way of determining the element
+                atomic_num = atom.type.atomic_number
+                elem = element.Element.getByAtomicNumber(atomic_num)
+            else:
+                # If we didn't load a parameter set yet, try based on mass
+                elem = element.Element.getByMass(atom.mass)
             topology.addAtom(atom.name, elem, residue)
 
         # Add all of the bonds
@@ -933,7 +938,7 @@ class CharmmPsfFile(object):
         # back up the dihedral list
         dihedral_list = self.dihedral_list
         # Load the parameter set
-        self.load_parameters(params)
+        self.load_parameters(params.condense())
         hasbox = self.topology.getUnitCellDimensions() is not None
         # Set the cutoff distance in nanometers
         cutoff = None

wrappers/python/simtk/openmm/app/element.py

@@ -34,7 +34,7 @@ __author__ = "Christopher M. Bruns"
 __version__ = "1.0"
 
 
-from simtk.unit import daltons
+from simtk.unit import daltons, is_quantity
 import copy_reg
 
 class Element(object):
@@ -93,6 +93,23 @@ class Element(object):
     def getByAtomicNumber(atomic_number):
         return Element._elements_by_atomic_number[atomic_number]
 
+    @staticmethod
+    def getByMass(mass):
+        # Assume masses are in daltons if they are not units
+        if not is_quantity(mass):
+            mass = mass * daltons
+        diff = mass
+        best_guess = 'EP'
+
+        for key in Element._elements_by_atomic_number:
+            element = Element._elements_by_atomic_number[key]
+            massdiff = abs(element.mass - mass)
+            if massdiff < diff:
+                best_guess = element
+                diff = massdiff
+
+        return best_guess
+
     @property
     def atomic_number(self):
         return self._atomic_number

wrappers/python/tests/TestCharmmFiles.py

@@ -13,20 +13,13 @@ class TestCharmmFiles(unittest.TestCase):
         """Set up the tests by loading the input files."""
 
         # alanine tripeptide; no waters
-        self.psf_c = CharmmPSF.load_from_psf('systems/ala_ala_ala.psf')
-        self.psf_x = CharmmPSF.load_from_psf('systems/ala_ala_ala.xpsf')
-        self.psf_v = CharmmPSF.load_from_psf('systems/ala_ala_ala.vpsf')
-        params = CharmmParameterSet.load_set(
-                            tfile='systems/charmm22.rtf',
-                            pfile='systems/charmm22.par',
-        ).condense()
+        self.psf_c = CharmmPsfFile('systems/ala_ala_ala.psf')
+        self.psf_x = CharmmPsfFile('systems/ala_ala_ala.xpsf')
+        self.psf_v = CharmmPsfFile('systems/ala_ala_ala.vpsf')
+        self.params = CharmmParameterSet(
+                            'systems/charmm22.rtf', 'systems/charmm22.par')
         self.pdb = PDBFile('systems/ala_ala_ala.pdb')
 
-        # Load parameters
-        self.psf_c.load_parameters(params)
-        self.psf_x.load_parameters(params)
-        self.psf_v.load_parameters(params)
-
     def test_NonbondedMethod(self):
         """Test both non-periodic methods for the systems"""
 
@@ -34,7 +27,7 @@ class TestCharmmFiles(unittest.TestCase):
                      CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic}
         for top in (self.psf_c, self.psf_x, self.psf_v):
             for method in methodMap:
-                system = top.createSystem(nonbondedMethod=method)
+                system = top.createSystem(self.params, nonbondedMethod=method)
                 forces = system.getForces()
                 self.assertTrue(any(isinstance(f, NonbondedForce) and 
                                     f.getNonbondedMethod()==methodMap[method] 
@@ -45,7 +38,7 @@ class TestCharmmFiles(unittest.TestCase):
 
         for top in (self.psf_c, self.psf_x, self.psf_v):
             for method in [CutoffNonPeriodic]:
-                system = top.createSystem(nonbondedMethod=method, 
+                system = top.createSystem(self.params, nonbondedMethod=method,
                                           nonbondedCutoff=2*nanometer, 
                                           constraints=HBonds)
                 cutoff_distance = 0.0*nanometer
@@ -59,21 +52,21 @@ class TestCharmmFiles(unittest.TestCase):
         """Test both options (True and False) for the removeCMMotion parameter."""
 
         for b in [True, False]:
-            system = self.psf_c.createSystem(removeCMMotion=b)
+            system = self.psf_c.createSystem(self.params, removeCMMotion=b)
             self.assertEqual(any(isinstance(f, CMMotionRemover) for f in system.getForces()), b)
 
     def test_ImplicitSolvent(self):
         """Test implicit solvent using the implicitSolvent parameter.
 
         """
-        system = self.psf_v.createSystem(implicitSolvent=OBC2)
+        system = self.psf_v.createSystem(self.params, implicitSolvent=OBC2)
         self.assertTrue(any(isinstance(f, CustomGBForce) for f in system.getForces()))
 
     def test_ImplicitSolventParameters(self):
         """Test that solventDielectric and soluteDielectric are passed correctly.
 
         """
-        system = self.psf_x.createSystem(implicitSolvent=GBn,
+        system = self.psf_x.createSystem(self.params, implicitSolvent=GBn,
                                          solventDielectric=50.0, 
                                          soluteDielectric = 0.9)
         found_matching_solvent_dielectric=False
@@ -97,8 +90,8 @@ class TestCharmmFiles(unittest.TestCase):
         
         topology = self.psf_v.topology
         hydrogenMass = 4*amu
-        system1 = self.psf_v.createSystem()
-        system2 = self.psf_v.createSystem(hydrogenMass=hydrogenMass)
+        system1 = self.psf_v.createSystem(self.params)
+        system2 = self.psf_v.createSystem(self.params, hydrogenMass=hydrogenMass)
         for atom in topology.atoms():
             if atom.element == elem.hydrogen:
                 self.assertNotEqual(hydrogenMass, system1.getParticleMass(atom.index))