Made energy minimization more robust

createSoftForcefield.py

@@ -143,7 +143,8 @@ print ' </PeriodicTorsionForce>'
 
 print ' <Script>'
 print """import simtk.openmm as mm
-nb = mm.CustomNonbondedForce('1/((r/0.2)^4+1)')
+nb = mm.CustomNonbondedForce('C/((r/0.2)^4+1)')
+nb.addGlobalParameter('C', 1.0)
 sys.addForce(nb)
 for i in range(sys.getNumParticles()):
     nb.addParticle([])

pdbfixer.py

@@ -397,7 +397,23 @@ class PDBFixer(object):
             context = mm.Context(system, integrator)
             context.setPositions(newPositions)
             mm.LocalEnergyMinimizer.minimize(context)
+            state = context.getState(getPositions=True)
+            nearest = self._findNearestDistance(context, newTopology, newAtoms)
+            if nearest < 0.15:
+                
+                # Some atoms are very close together.  Run some dynamics while slowly increasing the strength of the
+                # repulsive interaction to try to improve the result.
+                
+                for i in range(10):
+                    context.setParameter('C', 0.15*(i+1))
                     integrator.step(1000)
+                    d = self._findNearestDistance(context, newTopology, newAtoms)
+                    if d > nearest:
+                        nearest = d
+                        state = context.getState(getPositions=True)
+                        if nearest >= 0.15:
+                            break
+                context.setState(state)
                 mm.LocalEnergyMinimizer.minimize(context)
                 state = context.getState(getPositions=True)
             
@@ -472,6 +488,17 @@ class PDBFixer(object):
                 forcefield._templateSignatures[signature] = [template]
         return forcefield
     
+    def _findNearestDistance(self, context, topology, newAtoms):
+        positions = context.getState(getPositions=True).getPositions(asNumpy=True).value_in_unit(unit.nanometer)
+        atomResidue = [atom.residue for atom in topology.atoms()]
+        nearest = sys.float_info.max
+        for atom in newAtoms:
+            p = positions-positions[atom.index]
+            dist = math.sqrt(min(np.dot(p[i], p[i]) for i in range(len(atomResidue)) if atomResidue[i] != atom.residue))
+            if dist < nearest:
+                nearest = dist
+        return nearest
+
 
 if __name__=='__main__':
     if len(sys.argv) < 2:

soft.xml

@@ -4492,7 +4492,8 @@
  </PeriodicTorsionForce>
  <Script>
 import simtk.openmm as mm
-nb = mm.CustomNonbondedForce('1/((r/0.2)^4+1)')
+nb = mm.CustomNonbondedForce('C/((r/0.2)^4+1)')
+nb.addGlobalParameter('C', 1.0)
 sys.addForce(nb)
 for i in range(sys.getNumParticles()):
     nb.addParticle([])