def testCollectionQuantities(self):
     """ Tests the use of collections as Quantity values """
     s = [1, 2, 3] * u.centimeters
     self.assertEqual(str(s), '[1, 2, 3] cm')
     self.assertTrue(u.is_quantity(s))
     s2 = s / u.millimeters
     self.assertEqual(s2, [10.0, 20.0, 30.0])
     self.assertEqual(s2, s.value_in_unit(u.millimeters))
     # Test 2-D list
     s = [[1, 2, 3], [4, 5, 6]]
     s *= u.centimeters
     self.assertTrue(u.is_quantity(s))
     s2 = s / u.millimeters
     self.assertEqual(s2, [[10.0, 20.0, 30.0], [40.0, 50.0, 60.0]])
     self.assertEqual(s.value_in_unit(u.millimeters), s2)
     # Test tuples
     s = (1, 2, 3) * u.centimeters
     self.assertTrue(u.is_quantity(s))
     self.assertEqual(str(s), '(1, 2, 3) cm')
     s2 = s / u.millimeters
     self.assertEqual(s2, (10, 20, 30))
     self.assertIsInstance(s2, tuple)
     self.assertEqual(s.value_in_unit(u.millimeters), s2)
     self.assertIsInstance(s.value_in_unit(u.millimeters), tuple)
     x = [1, 2, 3] * u.centimeters
     x *= u.meters
     self.assertEqual(x, [100, 200, 300] * u.centimeters**2)

 def testCollectionQuantityOperations(self):
     """ Tests that Quantity collections behave correctly """
     # Tests that __getitem__ returns a unit
     s = [1, 2, 3, 4] * u.angstroms
     self.assertTrue(u.is_quantity(s[0]))
     for i, val in enumerate(s):
         self.assertTrue(u.is_quantity(val))
         self.assertEqual(val, (i+1) * u.angstroms)
     # Tests that __setitem__ fails when an incompatible type is added
     def fail(s): s[0] = 5
     self.assertRaises(AttributeError, lambda: fail(s))
     def fail(s): s[0] = 5 * u.joules
     self.assertRaises(TypeError, lambda: fail(s))
     def fail(s): s[0] /= 10 * u.meters
     self.assertRaises(AttributeError, lambda: fail(s))
     # Tests that __setitem__ converts to the unit of the container
     s[0] = 1 * u.nanometers
     self.assertEqual(s[0]._value, 10)
     # Tests standard unit conversions
     x = [1, 2, 3] * u.centimeters
     self.assertEqual(x / u.millimeters, [10, 20, 30])
     # Test the construction of a container in which each element is a
     # Quantity, passed to the Quantity constructor
     x = u.Quantity([1*u.angstrom, 2*u.nanometer, 3*u.angstrom])
     self.assertEqual(x._value, [1, 20, 3])
     self.assertEqual(x.unit, u.angstrom)
     x = u.Quantity((1, 2, 3))
     self.assertTrue(u.is_quantity(x))
     self.assertTrue(x.unit.is_dimensionless())
     x = u.Quantity(([1*u.angstrom, 2*u.nanometer, 3*u.angstrom],
                     [1*u.angstrom, 4*u.nanometer, 3*u.angstrom]))
     self.assertEqual(x._value, ([1, 20, 3], [1, 40, 3]))
     self.assertEqual(x.unit, u.angstrom)
     self.assertTrue(u.is_quantity(u.Quantity([])))

 def testScalarQuantityConstructor(self):
     """ Tests creating a Quantity using the Quantity constructor """
     self.assertTrue(u.is_quantity(u.Quantity(5, u.centimeters)))
     self.assertTrue(u.is_quantity(u.Quantity(5, u.centimeters**-1)))
     x = u.Quantity(value=5.0, unit=100.0*u.meters)
     self.assertTrue(u.is_quantity(x))
     self.assertEqual(x, 500*u.meters)

    def getEffectiveEnergy(self, totalEnergy, groupEnergy):
        """Given the actual potential energy of the system, return the value of the effective potential.

        Parameters
        ----------
        totalEnergy : energy
            the actual potential energy of the whole system
        groupEnergy : energy
            the actual potential energy of the boosted force group

        Returns
        -------
        energy
            the value of the effective potential
        """
        alphaTotal = self.getAlphaTotal()
        ETotal = self.getETotal()
        alphaGroup = self.getAlphaGroup()
        EGroup = self.getEGroup()
        if not is_quantity(totalEnergy):
            totalEnergy = totalEnergy*kilojoules_per_mole # Assume kJ/mole
        if not is_quantity(groupEnergy):
            groupEnergy = groupEnergy*kilojoules_per_mole # Assume kJ/mole
        dE = 0.0*kilojoules_per_mole
        if (totalEnergy < ETotal):
            dE = dE + (ETotal-totalEnergy)*(ETotal-totalEnergy)/(alphaTotal+ETotal-totalEnergy)
        if (groupEnergy < EGroup):
            dE = dE + (EGroup-groupEnergy)*(EGroup-groupEnergy)/(alphaGroup+EGroup-groupEnergy)
        return totalEnergy+dE

 def testNumpyDeepCopy(self):
     """ Check that deepcopy on numpy array does not strip units """
     x = u.Quantity(np.zeros((2, 3)), u.nanometer)
     y = copy.deepcopy(x)
     self.assertTrue(np.all(x == y))
     self.assertTrue(u.is_quantity(x))
     self.assertTrue(u.is_quantity(y))

    def writeModel(self, positions, unitCellDimensions=None, periodicBoxVectors=None):
        """Write out a model to the DCD file.

        The periodic box can be specified either by the unit cell dimensions (for a rectangular box), or the full set of box
        vectors (for an arbitrary triclinic box).  If neither is specified, the box vectors specified in the Topology will be
        used.  Regardless of the value specified, no dimensions will be written if the Topology does not represent a periodic system.

        Parameters:
         - positions (list) The list of atomic positions to write
         - unitCellDimensions (Vec3=None) The dimensions of the crystallographic unit cell.
         - periodicBoxVectors (tuple of Vec3=None) The vectors defining the periodic box.
        """
        if len(list(self._topology.atoms())) != len(positions):
            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.

        self._modelCount += 1
        file.seek(8, os.SEEK_SET)
        file.write(struct.pack('<i', self._modelCount))
        file.seek(20, os.SEEK_SET)
        file.write(struct.pack('<i', self._firstStep+self._modelCount*self._interval))

        # Write the data.

        file.seek(0, os.SEEK_END)
        boxVectors = self._topology.getPeriodicBoxVectors()
        if boxVectors is not None:
            if periodicBoxVectors is not None:
                boxVectors = periodicBoxVectors
            elif unitCellDimensions is not None:
                if is_quantity(unitCellDimensions):
                    unitCellDimensions = unitCellDimensions.value_in_unit(nanometers)
                boxVectors = (Vec3(unitCellDimensions[0], 0, 0), Vec3(0, unitCellDimensions[1], 0), Vec3(0, 0, unitCellDimensions[2]))*nanometers
            (a_length, b_length, c_length, alpha, beta, gamma) = computeLengthsAndAngles(boxVectors)
            a_length = a_length * 10.  # computeLengthsAndAngles returns unitless nanometers, but need angstroms here.
            b_length = b_length * 10.  # computeLengthsAndAngles returns unitless nanometers, but need angstroms here.
            c_length = c_length * 10.  # computeLengthsAndAngles returns unitless nanometers, but need angstroms here.
            angle1 = math.sin(math.pi/2-gamma)
            angle2 = math.sin(math.pi/2-beta)
            angle3 = math.sin(math.pi/2-alpha)
            file.write(struct.pack('<i6di', 48, a_length, angle1, b_length, angle2, angle3, c_length, 48))
        length = struct.pack('<i', 4*len(positions))
        for i in range(3):
            file.write(length)
            data = array.array('f', (10*x[i] for x in positions))
            data.tofile(file)
            file.write(length)

 def testAngleQuantities(self):
     """ Tests angle measurements """
     self.assertEqual(1.0*u.radians / u.degrees, 180 / math.pi)
     self.assertTrue(u.is_quantity(1.0*u.radians))
     self.assertTrue(u.is_quantity(1.0*u.degrees))
     self.assertEqual((1.0*u.radians).in_units_of(u.degrees),
                      (180 / math.pi)*u.degrees)
     self.assertEqual(90*u.degrees/u.radians, math.pi/2)
     q = 90 * u.degrees + 0.3 * u.radians
     self.assertEqual(q._value, 90 + 180*0.3/math.pi)
     self.assertEqual(q.unit, u.degrees)

    def _initializeConstants(self, simulation):
        """Initialize a set of constants required for the reports

        Parameters
        - simulation (Simulation) The simulation to generate a report for
        """
        system = simulation.system
        if self._temperature:
            # Compute the number of degrees of freedom.
            dof = 0
            for i in range(system.getNumParticles()):
                if system.getParticleMass(i) > 0*unit.dalton:
                    dof += 3
            dof -= system.getNumConstraints()
            if any(type(system.getForce(i)) == mm.CMMotionRemover for i in range(system.getNumForces())):
                dof -= 3
            self._dof = dof
        if self._density:
            if self._totalMass is None:
                # Compute the total system mass.
                self._totalMass = 0*unit.dalton
                for i in range(system.getNumParticles()):
                    self._totalMass += system.getParticleMass(i)
            elif not unit.is_quantity(self._totalMass):
                self._totalMass = self._totalMass*unit.dalton

 def testDimensionless(self):
     """ Tests the properties of unit.dimensionless """
     x = 5 * u.dimensionless
     y = u.Quantity(5, u.dimensionless)
     self.assertTrue(u.is_quantity(x))
     self.assertTrue(u.is_quantity(y))
     self.assertNotEqual(x, 5)
     self.assertNotEqual(y, 5)
     self.assertEqual(x, y)
     self.assertEqual(x.value_in_unit_system(u.si_unit_system), 5)
     self.assertEqual(x.value_in_unit_system(u.cgs_unit_system), 5)
     self.assertEqual(x.value_in_unit_system(u.md_unit_system), 5)
     x = u.Quantity(1.0, u.dimensionless)
     y = u.Quantity(1.0, u.dimensionless)
     self.assertIsNot(x, y)
     self.assertEqual(x, y)

    def __init__(self, file, topology, dt, firstStep=0, interval=1):
        """Create a DCD file and write out the header.

        Parameters:
         - file (file) A file to write to
         - topology (Topology) The Topology defining the molecular system being written
         - dt (time) The time step used in the trajectory
         - firstStep (int=0) The index of the first step in the trajectory
         - interval (int=1) The frequency (measured in time steps) at which states are written to the trajectory
        """
        self._file = file
        self._topology = topology
        self._firstStep = firstStep
        self._interval = interval
        self._modelCount = 0
        if is_quantity(dt):
            dt = dt.value_in_unit(picoseconds)
        dt /= 0.04888821
        boxFlag = 0
        if topology.getUnitCellDimensions() is not None:
            boxFlag = 1
        header = struct.pack('<i4c9if', 84, b'C', b'O', b'R', b'D', 0, firstStep, interval, 0, 0, 0, 0, 0, 0, dt)
        header += struct.pack('<13i', boxFlag, 0, 0, 0, 0, 0, 0, 0, 0, 24, 84, 164, 2)
        header += struct.pack('<80s', b'Created by OpenMM')
        header += struct.pack('<80s', b'Created '+time.asctime(time.localtime(time.time())).encode('ascii'))
        header += struct.pack('<4i', 164, 4, len(list(topology.atoms())), 4)
        file.write(header)

def _strip_optunit(thing, unit):
    """
    Strips optional units, converting to specified unit type. If no unit
    present, it just returns the number
    """
    if u.is_quantity(thing):
        return thing.value_in_unit(unit)
    return thing

    def runForClockTime(self, time, checkpointFile=None, stateFile=None, checkpointInterval=None):
        """Advance the simulation by integrating time steps until a fixed amount of clock time has elapsed.

        This is useful when you have a limited amount of computer time available, and want to run the longest simulation
        possible in that time.  This method will continue taking time steps until the specified clock time has elapsed,
        then return.  It also can automatically write out a checkpoint and/or state file before returning, so you can
        later resume the simulation.  Another option allows it to write checkpoints or states at regular intervals, so
        you can resume even if the simulation is interrupted before the time limit is reached.

        Parameters
        ----------
        time : time
            the amount of time to run for.  If no units are specified, it is
            assumed to be a number of hours.
        checkpointFile : string or file=None
            if specified, a checkpoint file will be written at the end of the
            simulation (and optionally at regular intervals before then) by
            passing this to saveCheckpoint().
        stateFile : string or file=None
            if specified, a state file will be written at the end of the
            simulation (and optionally at regular intervals before then) by
            passing this to saveState().
        checkpointInterval : time=None
            if specified, checkpoints and/or states will be written at regular
            intervals during the simulation, in addition to writing a final
            version at the end.  If no units are specified, this is assumed to
            be in hours.
        """
        if unit.is_quantity(time):
            time = time.value_in_unit(unit.hours)
        if unit.is_quantity(checkpointInterval):
            checkpointInterval = checkpointInterval.value_in_unit(unit.hours)
        endTime = datetime.now()+timedelta(hours=time)
        while (datetime.now() < endTime):
            if checkpointInterval is None:
                nextTime = endTime
            else:
                nextTime = datetime.now()+timedelta(hours=checkpointInterval)
                if nextTime > endTime:
                    nextTime = endTime
            self._simulate(endTime=nextTime)
            if checkpointFile is not None:
                self.saveCheckpoint(checkpointFile)
            if stateFile is not None:
                self.saveState(stateFile)

 def testNumpyDivision(self):
     """ Tests that division of numpy Quantities works correctly """
     x = u.Quantity(np.asarray([1., 2.]), u.nanometers)
     y = u.Quantity(np.asarray([3., 4.]), u.picoseconds)
     xy = x / y
     self.assertTrue(u.is_quantity(xy))
     self.assertEqual(xy.unit, u.nanometers/u.picoseconds)
     self.assertEqual(xy[0].value_in_unit(u.nanometers/u.picoseconds), 1/3)
     self.assertEqual(xy[1].value_in_unit(u.nanometers/u.picoseconds), 0.5)

 def testNumpyFunctions(self):
     """ Tests various numpy attributes that they result in Quantities """
     a = u.Quantity(np.arange(10), u.seconds)
     self.assertEqual(a.max(), 9*u.seconds)
     self.assertEqual(a.min(), 0*u.seconds)
     self.assertEqual(a.mean(), 4.5*u.seconds)
     self.assertAlmostEqualQuantities(a.std(), 2.8722813232690143*u.seconds)
     b = a.reshape((5, 2))
     self.assertTrue(u.is_quantity(b))

    def setUnitCellDimensions(self, dimensions):
        """Set the dimensions of the crystallographic unit cell.

        This method is an alternative to setPeriodicBoxVectors() for the case of a rectangular box.  It sets
        the box vectors to be orthogonal to each other and to have the specified lengths."""
        if dimensions is None:
            self._periodicBoxVectors = None
        else:
            if is_quantity(dimensions):
                dimensions = dimensions.value_in_unit(nanometers)
            self._periodicBoxVectors = (Vec3(dimensions[0], 0, 0), Vec3(0, dimensions[1], 0), Vec3(0, 0, dimensions[2]))*nanometers

    def writeModel(topology, positions, file=sys.stdout, modelIndex=None):
        """Write out a model to a PDB file.

        Parameters:
         - topology (Topology) The Topology defining the model to write
         - positions (list) The list of atomic positions to write
         - file (file=stdout) A file to write the model to
         - modelIndex (int=None) If not None, the model will be surrounded by MODEL/ENDMDL records with this index
        """
        if len(list(topology.atoms())) != len(positions):
            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:
            print >>file, "MODEL     %4d" % modelIndex
        for (chainIndex, chain) in enumerate(topology.chains()):
            chainName = chr(ord('A')+chainIndex%26)
            residues = list(chain.residues())
            for (resIndex, res) in enumerate(residues):
                if len(res.name) > 3:
                    resName = res.name[:3]
                else:
                    resName = res.name
                for atom in res.atoms():
                    if len(atom.name) < 4 and atom.name[:1].isalpha() and (atom.element is None or len(atom.element.symbol) < 2):
                        atomName = ' '+atom.name
                    elif len(atom.name) > 4:
                        atomName = atom.name[:4]
                    else:
                        atomName = atom.name
                    coords = positions[posIndex]
                    if atom.element is not None:
                        symbol = atom.element.symbol
                    else:
                        symbol = ' '
                    line = "ATOM  %5d %-4s %3s %s%4d    %s%s%s  1.00  0.00          %2s  " % (
                        atomIndex%100000, atomName, resName, chainName,
                        (resIndex+1)%10000, _format_83(coords[0]),
                        _format_83(coords[1]), _format_83(coords[2]), symbol)
                    assert len(line) == 80, 'Fixed width overflow detected'
                    print >>file, line
                    posIndex += 1
                    atomIndex += 1
                if resIndex == len(residues)-1:
                    print >>file, "TER   %5d      %3s %s%4d" % (atomIndex, resName, chainName, resIndex+1)
                    atomIndex += 1
        if modelIndex is not None:
            print >>file, "ENDMDL"

 def setPeriodicBoxVectors(self, vectors):
     """Set the vectors defining the periodic box."""
     if vectors is not None:
         if not is_quantity(vectors[0][0]):
             vectors = vectors*nanometers
         if vectors[0][1] != 0*nanometers or vectors[0][2] != 0*nanometers:
             raise ValueError("First periodic box vector must be parallel to x.");
         if vectors[1][2] != 0*nanometers:
             raise ValueError("Second periodic box vector must be in the x-y plane.");
         if vectors[0][0] <= 0*nanometers or vectors[1][1] <= 0*nanometers or vectors[2][2] <= 0*nanometers or vectors[0][0] < 2*abs(vectors[1][0]) or vectors[0][0] < 2*abs(vectors[2][0]) or vectors[1][1] < 2*abs(vectors[2][1]):
             raise ValueError("Periodic box vectors must be in reduced form.");
     self._periodicBoxVectors = deepcopy(vectors)

    def writeModel(topology, positions, file=sys.stdout, modelIndex=1, keepIds=False):
        """Write out a model to a PDBx/mmCIF file.

        Parameters
        ----------
        topology : Topology
            The Topology defining the model to write
        positions : list
            The list of atomic positions to write
        file : file=stdout
            A file to write the model to
        modelIndex : int=1
            The model number of this frame
        keepIds : bool=False
            If True, keep the residue and chain IDs specified in the Topology
            rather than generating new ones.  Warning: It is up to the caller to
            make sure these are valid IDs that satisfy the requirements of the
            PDBx/mmCIF format.  Otherwise, the output file will be invalid.
        """
        if len(list(topology.atoms())) != len(positions):
            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
        for (chainIndex, chain) in enumerate(topology.chains()):
            if keepIds:
                chainName = chain.id
            else:
                chainName = chr(ord('A')+chainIndex%26)
            residues = list(chain.residues())
            for (resIndex, res) in enumerate(residues):
                if keepIds:
                    resId = res.id
                    resIC = (res.insertionCode if len(res.insertionCode) > 0 else '.')
                else:
                    resId = resIndex + 1
                    resIC = '.'
                for atom in res.atoms():
                    coords = positions[posIndex]
                    if atom.element is not None:
                        symbol = atom.element.symbol
                    else:
                        symbol = '?'
                    line = "ATOM  %5d %-3s %-4s . %-4s %s ? %5s %s %10.4f %10.4f %10.4f  0.0  0.0  ?  ?  ?  ?  ?  .  %5s %4s %s %4s %5d"
                    print(line % (atomIndex, symbol, atom.name, res.name, chainName, resId, resIC, coords[0], coords[1], coords[2],
                                  resId, res.name, chainName, atom.name, modelIndex), file=file)
                    posIndex += 1
                    atomIndex += 1

    def __init__(self, topology, positions):
        """Create a new Modeller object

        Parameters:
         - topology (Topology) the initial Topology of the model
         - positions (list) the initial atomic positions
        """
        ## The Topology describing the structure of the system
        self.topology = topology
        if not is_quantity(positions):
            positions = positions*nanometer
        ## The list of atom positions
        self.positions = positions

 def getEffectiveEnergy(self, groupEnergy):
     """Given the actual group energy of the system, return the value of the effective potential.
     
     Parameters:
       - groupEnergy (energy): the actual potential energy of the boosted force group
     Returns: the value of the effective potential
     """
     alphaGroup = self.getAlphaGroup()
     EGroup = self.getEGroup()
     if not is_quantity(groupEnergy):
         groupEnergy = groupEnergy*kilojoules_per_mole # Assume kJ/mole
     dE = 0.0*kilojoules_per_mole
     if (groupEnergy < EGroup):
         dE = dE + (EGroup-groupEnergy)*(EGroup-groupEnergy)/(alphaGroup+EGroup-groupEnergy)
     return groupEnergy+dE