Reporters check for NaNs or infinities in output

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

@@ -35,8 +35,9 @@ import array
 import os
 import time
 import struct
-from simtk.unit import picoseconds, nanometers, angstroms, is_quantity
-    
+import math
+from simtk.unit import picoseconds, nanometers, angstroms, is_quantity, norm
+
 class DCDFile(object):
     """DCDFile provides methods for creating DCD files.
     
@@ -89,6 +90,10 @@ class DCDFile(object):
             raise ValueError('The number of positions must match the number of atoms') 
         if is_quantity(positions):
             positions = positions.value_in_unit(nanometers)
+        if any(math.isnan(norm(pos)) for pos in positions):
+            raise ValueError('Particle position is NaN')
+        if any(math.isinf(norm(pos)) for pos in positions):
+            raise ValueError('Particle position is infinite')
         file = self._file
         
         # Update the header.

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

@@ -33,12 +33,13 @@ __version__ = "1.0"
 
 import os
 import sys
+import math
 import xml.etree.ElementTree as etree
 from copy import copy
 from simtk.openmm import Vec3
 from simtk.openmm.app.internal.pdbstructure import PdbStructure
 from simtk.openmm.app import Topology
-from simtk.unit import nanometers, angstroms, is_quantity
+from simtk.unit import nanometers, angstroms, is_quantity, norm
 import element as elem
 try:
     import numpy
@@ -240,6 +241,10 @@ class PDBFile(object):
             raise ValueError('The number of positions must match the number of atoms') 
         if is_quantity(positions):
             positions = positions.value_in_unit(angstroms)
+        if any(math.isnan(norm(pos)) for pos in positions):
+            raise ValueError('Particle position is NaN')
+        if any(math.isinf(norm(pos)) for pos in positions):
+            raise ValueError('Particle position is infinite')
         atomIndex = 1
         posIndex = 0
         if modelIndex is not None:

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

@@ -34,6 +34,7 @@ __version__ = "1.0"
 import simtk.openmm as mm
 import simtk.unit as unit
 from simtk.openmm.app import PDBFile
+import math
     
 class StateDataReporter(object):
     """StateDataReporter outputs information about a simulation, such as energy and temperature, to a file.
@@ -136,6 +137,15 @@ class StateDataReporter(object):
             print >>self._out, '#"%s"' % ('"'+self._separator+'"').join(headers)
             self._hasInitialized = True
 
+        # Check for errors.
+        
+        if self._needEnergy:
+            energy = (state.getKineticEnergy()+state.getPotentialEnergy()).value_in_unit(unit.kilojoules_per_mole)
+            if math.isnan(energy):
+                raise ValueError('Energy is NaN')
+            if math.isinf(energy):
+                raise ValueError('Energy is infinite')
+                
         # Write the values.
         
         values = []