Merge branch 'master' of https://github.com/SimTk/openmm

wrappers/python/simtk/openmm/app/internal/amber_file_parser.py

@@ -41,9 +41,9 @@ USE OR OTHER DEALINGS IN THE SOFTWARE.
 
 import os
 import os.path
-import copy
 import re
 import math
+import warnings
 
 try:
     import numpy
@@ -428,15 +428,25 @@ class PrmtopLoader(object):
         charges=self.getCharges()
         nonbondTerms = self.getNonbondTerms()
         for ii in range(0,len(dihedralPointers),5):
-             if int(dihedralPointers[ii+2])>0 and int(dihedralPointers[ii+3])>0:
-                 iAtom = int(dihedralPointers[ii])//3
-                 lAtom = int(dihedralPointers[ii+3])//3
-                 chargeProd = charges[iAtom]*charges[lAtom]
-                 (rVdwI, epsilonI) = nonbondTerms[iAtom]
-                 (rVdwL, epsilonL) = nonbondTerms[lAtom]
-                 rMin = (rVdwI+rVdwL)
-                 epsilon = math.sqrt(epsilonI*epsilonL)
-                 returnList.append((iAtom, lAtom, chargeProd, rMin, epsilon))
+            if int(dihedralPointers[ii+2])>0 and int(dihedralPointers[ii+3])>0:
+                iAtom = int(dihedralPointers[ii])//3
+                lAtom = int(dihedralPointers[ii+3])//3
+                iidx = int(dihedralPointers[ii+4]) - 1
+                chargeProd = charges[iAtom]*charges[lAtom]
+                (rVdwI, epsilonI) = nonbondTerms[iAtom]
+                (rVdwL, epsilonL) = nonbondTerms[lAtom]
+                rMin = (rVdwI+rVdwL)
+                epsilon = math.sqrt(epsilonI*epsilonL)
+                try:
+                    iScee = float(self._raw_data["SCEE_SCALE_FACTOR"][iidx])
+                except KeyError:
+                    iScee = 1.2
+                try:
+                    iScnb = float(self._raw_data["SCNB_SCALE_FACTOR"][iidx])
+                except KeyError:
+                    iScnb = 2.0
+
+                returnList.append((iAtom, lAtom, chargeProd, rMin, epsilon, iScee, iScnb))
         return returnList
 
     def getExcludedAtoms(self):
@@ -502,11 +512,15 @@ class PrmtopLoader(object):
                 units.Quantity(y, units.angstrom),
                 units.Quantity(z, units.angstrom))
 
+    @property
+    def has_scee_scnb(self):
+        return ("SCEE_SCALE_FACTOR" in self._raw_data and "SCNB_SCALE_FACTOR" in self._raw_data)
+
 #=============================================================================================
 # AMBER System builder (based on, but not identical to, systemManager from 'zander')
 #=============================================================================================
 
-def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmodel=None, soluteDielectric=1.0, solventDielectric=78.5, nonbondedCutoff=None, nonbondedMethod='NoCutoff', scee=1.2, scnb=2.0, mm=None, verbose=False, EwaldErrorTolerance=None, flexibleConstraints=True, rigidWater=True):
+def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmodel=None, soluteDielectric=1.0, solventDielectric=78.5, nonbondedCutoff=None, nonbondedMethod='NoCutoff', scee=None, scnb=None, mm=None, verbose=False, EwaldErrorTolerance=None, flexibleConstraints=True, rigidWater=True):
     """
     Create an OpenMM System from an Amber prmtop file.
 
@@ -520,8 +534,8 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
       soluteDielectric (float) - The solute dielectric constant to use in the implicit solvent model (default: 1.0)
       solventDielectric (float) - The solvent dielectric constant to use in the implicit solvent model (default: 78.5)
       nonbondedCutoff (float) - if specified, will set nonbondedCutoff (default: None)
-      scnb (float) - 1-4 Lennard-Jones scaling factor (default: 1.2)
-      scee (float) - 1-4 electrostatics scaling factor (default: 2.0)
+      scnb (float) - 1-4 Lennard-Jones scaling factor (default: taken from prmtop or 1.2 if not present there)
+      scee (float) - 1-4 electrostatics scaling factor (default: taken from prmtop or 2.0 if not present there)
       mm - if specified, this module will be used in place of pyopenmm (default: None)
       verbose (boolean) - if True, print out information on progress (default: False)
       flexibleConstraints (boolean) - if True, flexible bonds will be added in addition ot constrained bonds
@@ -571,6 +585,10 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
     if prmtop.getIfBox()>1:
         raise Exception("only standard periodic boxes are currently supported")
 
+    if prmtop.has_scee_scnb and (scee is not None or scnb is not None):
+        warnings.warn("1-4 scaling parameters in topology file are being ignored. "
+            "This is not recommended unless you know what you are doing.")
+
     # Use pyopenmm implementation of OpenMM by default.
     if mm is None:
         mm = simtk.openmm
@@ -615,7 +633,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
     if shake == 'h-angles':
         numConstrainedBonds = system.getNumConstraints()
         atomConstraints = [[]]*system.getNumParticles()
-        for i in range(system.getNumConstraints()):
+        for i in range(numConstrainedBonds):
             c = system.getConstraintParameters(i)
             distance = c[2].value_in_unit(units.nanometer)
             atomConstraints[c[0]].append((c[1], distance))
@@ -712,9 +730,12 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
     # Add 1-4 Interactions
     excludedAtomPairs = set()
     sigmaScale = 2**(-1./6.)
-    for (iAtom, lAtom, chargeProd, rMin, epsilon) in prmtop.get14Interactions():
-        chargeProd /= scee
-        epsilon /= scnb
+    _scee, _scnb = scee, scnb
+    for (iAtom, lAtom, chargeProd, rMin, epsilon, iScee, iScnb) in prmtop.get14Interactions():
+        if scee is None: _scee = iScee
+        if scnb is None: _scnb = iScnb
+        chargeProd /= _scee
+        epsilon /= _scnb
         sigma = rMin * sigmaScale
         force.addException(iAtom, lAtom, chargeProd, sigma, epsilon)
         excludedAtomPairs.add(min((iAtom, lAtom), (lAtom, iAtom)))

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

@@ -313,6 +313,56 @@ class PDBFile(object):
          - topology (Topology) The Topology defining the molecular system being written
          - file (file=stdout) A file to write the file to
         """
+        # Identify bonds that should be listed as CONECT records.
+        
+        standardResidues = ['ALA', 'ASN', 'CYS', 'GLU', 'HIS', 'LEU', 'MET', 'PRO', 'THR', 'TYR',
+                            'ARG', 'ASP', 'GLN', 'GLY', 'ILE', 'LYS', 'PHE', 'SER', 'TRP', 'VAL',
+                            'A', 'G', 'C', 'U', 'I', 'DA', 'DG', 'DC', 'DT', 'DI', 'HOH']
+        conectBonds = []
+        for atom1, atom2 in topology.bonds():
+            if atom1.residue.name not in standardResidues or atom2.residue.name not in standardResidues:
+                conectBonds.append((atom1, atom2))
+            elif atom1.name == 'SG' and atom2.name == 'SG' and atom1.residue.name == 'CYS' and atom2.residue.name == 'CYS':
+                conectBonds.append((atom1, atom2))
+        if len(conectBonds) > 0:
+            
+            # Work out the index used in the PDB file for each atom.
+            
+            atomIndex = {}
+            nextAtomIndex = 0
+            prevChain = None
+            for chain in topology.chains():
+                for atom in chain.atoms():
+                    if atom.residue.chain != prevChain:
+                        nextAtomIndex += 1
+                        prevChain = atom.residue.chain
+                    atomIndex[atom] = nextAtomIndex
+                    nextAtomIndex += 1
+            
+            # Record which other atoms each atom is bonded to.
+            
+            atomBonds = {}
+            for atom1, atom2 in conectBonds:
+                index1 = atomIndex[atom1]
+                index2 = atomIndex[atom2]
+                if index1 not in atomBonds:
+                    atomBonds[index1] = []
+                if index2 not in atomBonds:
+                    atomBonds[index2] = []
+                atomBonds[index1].append(index2)
+                atomBonds[index2].append(index1)
+            
+            # Write the CONECT records.
+            
+            for index1 in sorted(atomBonds):
+                bonded = atomBonds[index1]
+                while len(bonded) > 4:
+                    print >>file, "CONECT%5d%5d%5d%5d" % (index1, bonded[0], bonded[1], bonded[2])
+                    del bonded[:4]
+                line = "CONECT%5d" % index1
+                for index2 in bonded:
+                    line = "%s%5d" % (line, index2)
+                print >>file, line
         print >>file, "END"
 
 

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

@@ -79,5 +79,6 @@ class PDBReporter(object):
         self._nextModel += 1
 
     def __del__(self):
-        PDBFile.writeFooter(self._topology, self._out)
+        if self._topology is not None:
+            PDBFile.writeFooter(self._topology, self._out)
         self._out.close()