Python unit.sqrt函数代码示例

def test_propertyestimator():
    """Test PhysicalPropertyEstimator on synthetic data.

    """
    # Get synthetic dataset and parameters
    from openforcefield.tests import testsystems
    dataset = testsystems.TestDataset()
    parameters = testsystems.TestParameterSet()

    # Create a PropertyEstimator
    from openforcefield.propertyestimator import PropertyEstimator
    estimator = PropertyEstimator()

    # Compute properties
    computed_properties = estimator.computeProperties(dataset, parameters)

    # Assert that computed properties are within statistical error
    for (measured_property, computed_property) in zip(dataset, computed_properties):
        error_value = (computed_property.value - measured_property.value)
        error_uncertainty = unit.sqrt(computed_property.uncertainty**2 + measured_property.uncertainty**2)
        relative_error = error_value / error_uncertainty
        if (relative_error > SIGMA_CUTOFF):
            msg  = 'Computed property %s differs from measured property by more than SIGMA_CUTOFF (%f):\n' % (computed_property.name, SIGMA_CUTOFF)
            msg += 'Measured: %12.3f +- %12.3f %s' % (measured_property.value / measured_property.unit, measured_property.uncertainty / measured_property.unit, str(measured_property.unit))
            msg += 'Computed: %12.3f +- %12.3f %s' % (computed_property.value / measured_property.unit, computed_property.uncertainty / measured_property.unit, str(measured_property.unit))
            msg += 'ERROR   : %12.3f +- %12.3f %s (%12.3f SIGMA)' % (error_value / measured_property.unit, error_uncertainty / measured_property.unit, str(measured_property.unit), relative_error)
            raise Exception(msg)

    def createDisulfideBonds(self, positions):
        """Identify disulfide bonds based on proximity and add them to the
        Topology.

        Parameters
        ----------
        positions : list
            The list of atomic positions based on which to identify bonded atoms
        """
        def isCyx(res):
            names = [atom.name for atom in res._atoms]
            return 'SG' in names and 'HG' not in names

        cyx = [res for res in self.residues() if res.name == 'CYS' and isCyx(res)]
        atomNames = [[atom.name for atom in res._atoms] for res in cyx]
        for i in range(len(cyx)):
            sg1 = cyx[i]._atoms[atomNames[i].index('SG')]
            pos1 = positions[sg1.index]
            for j in range(i):
                sg2 = cyx[j]._atoms[atomNames[j].index('SG')]
                pos2 = positions[sg2.index]
                delta = [x-y for (x,y) in zip(pos1, pos2)]
                distance = sqrt(delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2])
                if distance < 0.3*nanometers:
                    self.addBond(sg1, sg2)

    def _automatic_parameter_selection(self, positions, receptor_atoms, ligand_atoms):
        """
        Determine parameters and restrained atoms automatically, rejecting choices where standard state correction will be incorrectly computed.

        Parameters
        ----------
        positions : simtk.unit.Quantity of natoms x 3 with units compatible with nanometers
            Reference positions to use for imposing restraints
        receptor_atoms : list of int
            A complete list of receptor atoms
        ligand_atoms : list of int
            A complete list of ligand atoms
        """
        NSIGMA = 4
        temperature = 300 * unit.kelvin
        kT = kB * temperature
        attempt = 0
        MAX_ATTEMPTS = 100
        reject = True
        logger.debug('Automatically selecting restraint atoms and parameters:')
        while reject and attempt < MAX_ATTEMPTS:
            logger.debug('Attempt %d / %d at automatically selecting atoms and restraint parameters...' % (attempt, MAX_ATTEMPTS))

            # Select atoms to be used in restraint.
            self._restraint_atoms = self._select_restraint_atoms(positions, receptor_atoms, ligand_atoms)

            # Determine restraint parameters
            self._determine_restraint_parameters()

            # Terminate if we satisfy criteria
            reject = False
            for name in ['A', 'B']:
                theta0 = self._parameters['theta_' + name + '0']
                K = self._parameters['K_theta' + name]
                sigma = unit.sqrt(NSIGMA * kT / (K/2.))
                if (theta0 < sigma) or (theta0 > (np.pi*unit.radians - sigma)):
                    logger.debug('Reject because theta_' + name + '0 is too close to 0 or pi for standard state correction to be accurate.')
                    reject = True

            r0 = self._parameters['r_aA0']
            K = self._parameters['K_r']
            sigma = unit.sqrt(NSIGMA * kT / (K/2.))
            if (r0 < sigma):
                logger.debug('Reject because r_aA0 is too close to 0 for standard state correction to be accurate.')
                reject = True

            attempt += 1

    def end_to_end_CA_distance(self, topology, positions):
        residues = list(topology.residues())
        # get the index of the first and last alpha carbons
        i1 = [a.index for a in residues[0].atoms() if a.name == 'CA'][0]
        i2 = [a.index for a in residues[-1].atoms() if a.name == 'CA'][0]

        # get the current distanc be between the two alpha carbons
        return i1, i2, sqrt(sum((positions[i1] - positions[i2])**2))

def check_harmonic_oscillator_ncmc(ncmc_nsteps=50, ncmc_integrator="VV"):
    """
    Test NCMC switching of a 3D harmonic oscillator.
    In this test, the oscillator center is dragged in space, and we check the computed free energy difference with BAR, which should be 0.
    """
    # Parameters for 3D harmonic oscillator
    mass = 39.948 * unit.amu # mass of particle (argon)
    sigma = 5.0 * unit.angstrom # standard deviation of harmonic oscillator
    collision_rate = 5.0/unit.picosecond # collision rate
    temperature = 300.0 * unit.kelvin # temperature
    platform_name = 'Reference' # platform anme
    NSIGMA_MAX = 6.0 # number of standard errors away from analytical solution tolerated before Exception is thrown

    # Compute derived quantities.
    kT = kB * temperature # thermal energy
    beta = 1.0 / kT # inverse energy
    K = kT / sigma**2 # spring constant
    tau = 2 * math.pi * unit.sqrt(mass/K) # time constant
    timestep = tau / 20.0
    platform = openmm.Platform.getPlatformByName(platform_name)

    # Create a 3D harmonic oscillator with context parameter controlling center of oscillator.
    system = openmm.System()
    system.addParticle(mass)
    energy_expression = '(K/2.0) * ((x-x0)^2 + y^2 + z^2);'
    force = openmm.CustomExternalForce(energy_expression)
    force.addGlobalParameter('K', K.in_unit_system(unit.md_unit_system))
    force.addGlobalParameter('x0', 0.0)
    force.addParticle(0, [])
    system.addForce(force)

    # Set the positions at the origin.
    positions = unit.Quantity(np.zeros([1, 3], np.float32), unit.angstroms)
    functions = { 'x0' : 'lambda' } # drag spring center x0

    from perses.annihilation import NCMCVVAlchemicalIntegrator, NCMCGHMCAlchemicalIntegrator
    if ncmc_integrator=="VV":
        ncmc_insert = NCMCVVAlchemicalIntegrator(temperature, system, functions, direction='insert', nsteps=ncmc_nsteps, timestep=timestep) # 'insert' drags lambda from 0 -> 1
        ncmc_delete = NCMCVVAlchemicalIntegrator(temperature, system, functions, direction='delete', nsteps=ncmc_nsteps, timestep=timestep) # 'insert' drags lambda from 0 -> 1
    elif ncmc_integrator=="GHMC":
        ncmc_insert = NCMCGHMCAlchemicalIntegrator(temperature, system, functions, direction='insert', collision_rate=9.1/unit.picoseconds, nsteps=ncmc_nsteps, timestep=timestep) # 'insert' drags lambda from 0 -> 1
        ncmc_delete = NCMCGHMCAlchemicalIntegrator(temperature, system, functions, direction='delete', collision_rate=9.1/unit.picoseconds, nsteps=ncmc_nsteps, timestep=timestep) # 'insert' drags lambda from 0 -> 1
    else:
        raise Exception("%s not recognized as integrator name. Options are VV and GHMC" % ncmc_integrator)

    # Run NCMC switching trials where the spring center is switched with lambda: 0 -> 1 over a finite number of steps.
    w_f = collect_switching_data(system, positions, functions, temperature, collision_rate, timestep, platform, ncmc_integrator=ncmc_insert, ncmc_nsteps=ncmc_nsteps, direction='insert')
    w_r = collect_switching_data(system, positions, functions, temperature, collision_rate, timestep, platform, ncmc_integrator=ncmc_delete, ncmc_nsteps=ncmc_nsteps, direction='delete')

    from pymbar import BAR
    [df, ddf] = BAR(w_f, w_r, method='self-consistent-iteration')
    print('%8.3f +- %.3f kT' % (df, ddf))
    if (abs(df) > NSIGMA_MAX * ddf):
        msg = 'Delta F (%d steps switching) = %f +- %f kT; should be within %f sigma of 0' % (ncmc_nsteps, df, ddf, NSIGMA_MAX)
        msg += '\n'
        msg += 'w_f = %s\n' % str(w_f)
        msg += 'w_r = %s\n' % str(w_r)
        raise Exception(msg)

def computeHarmonicOscillatorExpectations(K, mass, temperature):
    """
    Compute mean and variance of potential and kinetic energies for a 3D harmonic oscillator.
    
    NOTES

    Numerical quadrature is used to compute the mean and standard deviation of the potential energy.
    Mean and standard deviation of the kinetic energy, as well as the absolute free energy, is computed analytically.
    
    ARGUMENTS
    
    K (simtk.unit.Quantity) - spring constant
    mass (simtk.unit.Quantity) - mass of particle
    temperature (simtk.unit.Quantity) - temperature
    
    RETURNS
    
    values (dict)
    
    """

    values = dict()

    # Compute thermal energy and inverse temperature from specified temperature.
    kB = units.BOLTZMANN_CONSTANT_kB * units.AVOGADRO_CONSTANT_NA
    kT = kB * temperature # thermal energy
    beta = 1.0 / kT # inverse temperature
   
    # Compute standard deviation along one dimension.
    sigma = 1.0 / units.sqrt(beta * K) 

    # Define limits of integration along r.
    r_min = 0.0 * units.nanometers # initial value for integration
    r_max = 10.0 * sigma      # maximum radius to integrate to

    # Compute mean and std dev of potential energy.
    V = lambda r : (K/2.0) * (r*units.nanometers)**2 / units.kilojoules_per_mole # potential in kJ/mol, where r in nm
    q = lambda r : 4.0 * math.pi * r**2 * math.exp(-beta * (K/2.0) * (r*units.nanometers)**2) # q(r), where r in nm
    (IqV2, dIqV2) = scipy.integrate.quad(lambda r : q(r) * V(r)**2, r_min / units.nanometers, r_max / units.nanometers)
    (IqV, dIqV)   = scipy.integrate.quad(lambda r : q(r) * V(r), r_min / units.nanometers, r_max / units.nanometers)
    (Iq, dIq)     = scipy.integrate.quad(lambda r : q(r), r_min / units.nanometers, r_max / units.nanometers)
    values['potential'] = dict()
    values['potential']['mean'] = (IqV / Iq) * units.kilojoules_per_mole
    values['potential']['stddev'] = (IqV2 / Iq) * units.kilojoules_per_mole   
    
    # Compute mean and std dev of kinetic energy.
    values['kinetic'] = dict()
    values['kinetic']['mean'] = (3./2.) * kT
    values['kinetic']['stddev'] = math.sqrt(3./2.) * kT

    # Compute dimensionless free energy.
    # f = - \ln \int_{-\infty}^{+\infty} \exp[-\beta K x^2 / 2] 
    #   = - \ln \int_{-\infty}^{+\infty} \exp[-x^2 / 2 \sigma^2] 
    #   = - \ln [\sqrt{2 \pi} \sigma]
    values['f'] = - numpy.log(2 * numpy.pi * (sigma / units.angstroms)**2) * (3.0/2.0)

    return values   

 def testUnitMathModule(self):
     """ Tests the unit_math functions on Quantity objects """
     self.assertEqual(u.sqrt(1.0*u.kilogram*u.joule),
                      1.0*u.kilogram*u.meter/u.second)
     self.assertEqual(u.sqrt(1.0*u.kilogram*u.calorie),
                      math.sqrt(4.184)*u.kilogram*u.meter/u.second)
     self.assertEqual(u.sqrt(9), 3) # Test on a scalar
     self.assertEqual(u.sin(90*u.degrees), 1)
     self.assertEqual(u.sin(math.pi/2*u.radians), 1)
     self.assertEqual(u.sin(math.pi/2), 1)
     self.assertEqual(u.cos(180*u.degrees), -1)
     self.assertEqual(u.cos(math.pi*u.radians), -1)
     self.assertEqual(u.cos(math.pi), -1)
     self.assertAlmostEqual(u.tan(45*u.degrees), 1)
     self.assertAlmostEqual(u.tan(math.pi/4*u.radians), 1)
     self.assertAlmostEqual(u.tan(math.pi/4), 1)
     acos = u.acos(1.0)
     asin = u.asin(1.0)
     atan = u.atan(1.0)
     self.assertTrue(u.is_quantity(acos))
     self.assertTrue(u.is_quantity(asin))
     self.assertTrue(u.is_quantity(atan))
     self.assertEqual(acos.unit, u.radians)
     self.assertEqual(asin.unit, u.radians)
     self.assertEqual(atan.unit, u.radians)
     self.assertEqual(acos.value_in_unit(u.degrees), 0)
     self.assertEqual(acos / u.radians, 0)
     self.assertEqual(asin.value_in_unit(u.degrees), 90)
     self.assertEqual(asin / u.radians, math.pi/2)
     self.assertAlmostEqual(atan.value_in_unit(u.degrees), 45)
     self.assertAlmostEqual(atan / u.radians, math.pi/4)
     # Check some sequence maths
     seq = [1, 2, 3, 4] * u.meters
     self.assertEqual(u.sum(seq), 10*u.meters)
     self.assertEqual(u.dot(seq, seq), (1+4+9+16)*u.meters**2)
     self.assertEqual(u.norm(seq), math.sqrt(30)*u.meters)

def computeHarmonicOscillatorExpectations(K, mass, temperature):
   """
   Compute mean and variance of potential and kinetic energies for harmonic oscillator.

   Numerical quadrature is used.

   ARGUMENTS

   K - spring constant
   mass - mass of particle
   temperature - temperature

   RETURNS

   values

   """

   values = dict()

   # Compute thermal energy and inverse temperature from specified temperature.
   kB = units.BOLTZMANN_CONSTANT_kB * units.AVOGADRO_CONSTANT_NA
   kT = kB * temperature # thermal energy
   beta = 1.0 / kT # inverse temperature
   
   # Compute standard deviation along one dimension.
   sigma = 1.0 / units.sqrt(beta * K) 

   # Define limits of integration along r.
   r_min = 0.0 * units.nanometers # initial value for integration
   r_max = 10.0 * sigma      # maximum radius to integrate to

   # Compute mean and std dev of potential energy.
   V = lambda r : (K/2.0) * (r*units.nanometers)**2 / units.kilojoules_per_mole # potential in kJ/mol, where r in nm
   q = lambda r : 4.0 * math.pi * r**2 * math.exp(-beta * (K/2.0) * (r*units.nanometers)**2) # q(r), where r in nm
   (IqV2, dIqV2) = scipy.integrate.quad(lambda r : q(r) * V(r)**2, r_min / units.nanometers, r_max / units.nanometers)
   (IqV, dIqV)   = scipy.integrate.quad(lambda r : q(r) * V(r), r_min / units.nanometers, r_max / units.nanometers)
   (Iq, dIq)     = scipy.integrate.quad(lambda r : q(r), r_min / units.nanometers, r_max / units.nanometers)
   values['potential'] = dict()
   values['potential']['mean'] = (IqV / Iq) * units.kilojoules_per_mole
   values['potential']['stddev'] = (IqV2 / Iq) * units.kilojoules_per_mole   
   
   # Compute mean and std dev of kinetic energy.
   values['kinetic'] = dict()
   values['kinetic']['mean'] = (3./2.) * kT
   values['kinetic']['stddev'] = math.sqrt(3./2.) * kT

   return values

    def __init__(self, **kwargs):
        super(HarmonicOscillatorSimulatedTempering, self).__init__(**kwargs)
        self.description = 'Harmonic oscillator simulated tempering simulation'

        # Create topology, positions, and system.
        from openmmtools.testsystems import HarmonicOscillator
        K = 1.0 * unit.kilocalories_per_mole / unit.angstroms**2 # 3D harmonic oscillator spring constant
        mass = 39.948 * unit.amu # 3D harmonic oscillator particle mass
        period = 2.0 * np.pi * unit.sqrt(mass / K) # harmonic oscillator period
        timestep = 0.01 * period
        testsystem = HarmonicOscillator(K=K, mass=mass)
        self.topology = testsystem.topology
        self.positions = testsystem.positions
        self.system = testsystem.system

        # Create thermodynamic states.
        Tmin = 100 * unit.kelvin
        Tmax = 1000 * unit.kelvin
        ntemps = 8 # number of temperatures
        from sams import ThermodynamicState
        temperatures = unit.Quantity(np.logspace(np.log10(Tmin / unit.kelvin), np.log10(Tmax / unit.kelvin), ntemps), unit.kelvin)
        self.thermodynamic_states = [ ThermodynamicState(system=self.system, temperature=temperature) for temperature in temperatures ]

        # Compute analytical logZ for each thermodynamic state.
        self.logZ = np.zeros([ntemps], np.float64)
        for (index, thermodynamic_state) in enumerate(self.thermodynamic_states):
            beta = thermodynamic_state.beta
            self.logZ[index] = - 1.5 * np.log(beta * K * unit.angstrom**2)
        self.logZ[:] -= self.logZ[0]

        # Create SAMS samplers
        from sams.samplers import SamplerState, MCMCSampler, ExpandedEnsembleSampler, SAMSSampler
        thermodynamic_state_index = 0 # initial thermodynamic state index
        thermodynamic_state = self.thermodynamic_states[thermodynamic_state_index]
        sampler_state = SamplerState(positions=self.positions)
        self.mcmc_sampler = MCMCSampler(sampler_state=sampler_state, thermodynamic_state=thermodynamic_state, ncfile=self.ncfile)
        self.mcmc_sampler.pdbfile = open('output.pdb', 'w')
        self.mcmc_sampler.topology = self.topology
        self.mcmc_sampler.timestep = timestep
        self.mcmc_sampler.collision_rate = 1.0 / (100 * timestep)
        self.mcmc_sampler.nsteps = 1000
        self.mcmc_sampler.verbose = True
        self.exen_sampler = ExpandedEnsembleSampler(self.mcmc_sampler, self.thermodynamic_states)
        self.exen_sampler.verbose = True
        self.sams_sampler = SAMSSampler(self.exen_sampler, update_stages='two-stage', update_method='optimal')
        self.sams_sampler.verbose = True

 def get14Interactions(self):
     """Return list of atom pairs, chargeProduct, rMin and epsilon for each 1-4 interaction"""
     dihedralPointers = self._raw_data["DIHEDRALS_INC_HYDROGEN"] \
                       +self._raw_data["DIHEDRALS_WITHOUT_HYDROGEN"]
     returnList=[]
     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 = units.sqrt(epsilonI*epsilonL)
              returnList.append((iAtom, lAtom, chargeProd, rMin, epsilon))
     return returnList

def generateMaxwellBoltzmannVelocities(system, temperature):
   """Generate Maxwell-Boltzmann velocities.
   
   ARGUMENTS
   
   system (simtk.openmm.System) - the system for which velocities are to be assigned
   temperature (simtk.unit.Quantity of temperature) - the temperature at which velocities are to be assigned
   
   RETURNS
   
   velocities (simtk.unit.Quantity of numpy Nx3 array, units length/time) - particle velocities
   
   TODO

   This could be sped up by introducing vector operations.
   
   """
   
   # Get number of atoms
   natoms = system.getNumParticles()
   
   # Create storage for velocities.        
   velocities = units.Quantity(numpy.zeros([natoms, 3], numpy.float32), units.nanometer / units.picosecond) # velocities[i,k] is the kth component of the velocity of atom i
   
   # Compute thermal energy and inverse temperature from specified temperature.
   kB = units.BOLTZMANN_CONSTANT_kB * units.AVOGADRO_CONSTANT_NA
   kT = kB * temperature # thermal energy
   beta = 1.0 / kT # inverse temperature
   
   # Assign velocities from the Maxwell-Boltzmann distribution.
   for atom_index in range(natoms):
      mass = system.getParticleMass(atom_index) # atomic mass
      sigma = units.sqrt(kT / mass) # standard deviation of velocity distribution for each coordinate for this atom
      for k in range(3):
         velocities[atom_index,k] = sigma * numpy.random.normal()

   # Return velocities
   return velocities

    def _assign_Maxwell_Boltzmann_velocities(self, system, temperature):
        """Generate Maxwell-Boltzmann velocities.

        @param system the system for which velocities are to be assigned
        @type simtk.chem.openmm.System or System
        
        @param temperature the temperature at which velocities are to be assigned
        @type Quantity with units of temperature

        @return velocities drawn from the Maxwell-Boltzmann distribution at the appropriate temperature
        @returntype (natoms x 3) numpy array wrapped in Quantity with units of velocity

        TODO

        This could be sped up by introducing vector operations.

        """

        # Get number of atoms
        natoms = system.getNumParticles()

        # Create storage for velocities.        
        velocities = units.Quantity(numpy.zeros([natoms, 3], numpy.float32), units.nanometer / units.picosecond) # velocities[i,k] is the kth component of the velocity of atom i
  
        # Compute thermal energy and inverse temperature from specified temperature.
        kT = kB * temperature # thermal energy
        beta = 1.0 / kT # inverse temperature
  
        # Assign velocities from the Maxwell-Boltzmann distribution.
        for atom_index in range(natoms):
            mass = system.getParticleMass(atom_index) # atomic mass
            sigma = units.sqrt(kT / mass) # standard deviation of velocity distribution for each coordinate for this atom
            for k in range(3):
                velocities[atom_index,k] = sigma * numpy.random.normal()

        # Return velocities
        return velocities

# Roughly corresponds to conditions from http://www.cstl.nist.gov/srs/LJ_PURE/mc.htm
nparticles = 500
mass = 39.9 * unit.amu
sigma = 3.4 * unit.angstrom
epsilon = 0.238 * unit.kilocalories_per_mole
#reduced_density = 0.860     # reduced_density = density * (sigma**3)
reduced_density = 0.960     # reduced_density = density * (sigma**3)
reduced_temperature = 0.850 # reduced_temperature = kB * temperature / epsilon
reduced_pressure = 1.2660   # reduced_pressure = pressure * (sigma**3) / epsilon

platform_name = 'CUDA'    # OpenMM platform name to use for simulation
platform = openmm.Platform.getPlatformByName(platform_name)
data_directory = 'data'     # Directory in which data is to be stored

r0 = 2.0**(1./6.) * sigma   # minimum potential distance for Lennard-Jones interaction
characteristic_timescale = unit.sqrt((mass * r0**2) / (72 * epsilon)) # characteristic timescale for bound Lennard-Jones interaction
                                                            # http://borisv.lk.net/matsc597c-1997/simulations/Lecture5/node3.html
timestep = 0.01 * characteristic_timescale # integrator timestep

# From http://www.cstl.nist.gov/srs/LJ_PURE/md.htm
#characteristic_timescale = unit.sqrt(mass * sigma**2 / epsilon)
#timestep = 0.05 * characteristic_timescale

print "characteristic timescale = %.3f ps" % (characteristic_timescale / unit.picoseconds)
print "timestep = %.12f ps" % (timestep / unit.picoseconds)

collision_rate = 5.0 / unit.picoseconds # collision rate for Langevin thermostat
barostat_frequency = 25 # number of steps between barostat updates

# Set parameters for number of simulation replicates, number of iterations per simulation, and number of steps per iteration.
nreplicates = 100

#    ntrials = integrator.getGlobalVariable(global_variables['ntrials'])
#    print "accepted %d / %d (%.3f %%)" % (naccept, ntrials, float(naccept)/float(ntrials)*100.0)

    # Accumulate statistics.
    x_n[iteration] = state.getPositions(asNumpy=True)[0,0] / units.angstroms
    potential_n[iteration] = final_potential_energy / kT
    kinetic_n[iteration] = final_kinetic_energy / kT
    temperature_n[iteration] = instantaneous_temperature / units.kelvin
    delta_n[iteration] = delta_total_energy / kT


# Compute expected statistics for harmonic oscillator.
K = 100.0 * units.kilocalories_per_mole / units.angstroms**2
beta = 1.0 / kT
x_mean_exact = 0.0 # mean, in angstroms
x_std_exact = 1.0 / units.sqrt(beta * K) / units.angstroms # std dev, in angstroms

# Analyze statistics.
g = statisticalInefficiency(potential_n) 
Neff = niterations / g # number of effective samples

x_mean = x_n.mean()
dx_mean = x_n.std() / numpy.sqrt(Neff)
x_mean_error = x_mean - x_mean_exact

x_var = x_n.var()
dx_var = x_var * numpy.sqrt(2. / (Neff-1))

x_std = x_n.std()
dx_std = 0.5 * dx_var / x_std 
x_std_error = x_std - x_std_exact

        test_success = True
        for platform_index in range(openmm.Platform.getNumPlatforms()):
            try:

                platform = openmm.Platform.getPlatform(platform_index)
                platform_name = platform.getName()
                [platform_potential, platform_force] = compute_potential_and_force(system, positions, platform)

                # Compute error in potential.
                potential_error = platform_potential - reference_potential

                # Compute per-atom RMS (magnitude) and RMS error in force.
                force_unit = units.kilocalories_per_mole / units.nanometers
                natoms = system.getNumParticles()
                force_mse = (((reference_force - platform_force) / force_unit)**2).sum() / natoms * force_unit**2
                force_rmse = units.sqrt(force_mse)

                force_ms = ((platform_force / force_unit)**2).sum() / natoms * force_unit**2
                force_rms = units.sqrt(force_ms)

                print "%32s %16.6f kcal/mol %16.6f kcal/mol %16.6f kcal/mol %16.6f kcal/mol" % (platform_name, platform_potential / units.kilocalories_per_mole, potential_error / units.kilocalories_per_mole, force_rms / force_unit, force_rmse / force_unit)

                # Mark whether tolerance is exceeded or not.
                if abs(potential_error) > ENERGY_TOLERANCE:
                    test_success = False
                    print "%32s WARNING: Potential energy error (%.6f kcal/mol) exceeds tolerance (%.6f kcal/mol).  Test failed." % ("", potential_error/units.kilocalories_per_mole, ENERGY_TOLERANCE/units.kilocalories_per_mole)
                if abs(force_rmse) > FORCE_RMSE_TOLERANCE:
                    test_success = False
                    print "%32s WARNING: Force RMS error (%.6f kcal/mol) exceeds tolerance (%.6f kcal/mol).  Test failed." % ("", force_rmse/force_unit, FORCE_RMSE_TOLERANCE/force_unit)
                    if debug:
                        for atom_index in range(natoms):