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):

    nsteps = 250 # number of steps per interation
    deviceid = 0
    print 'nsteps: ', nsteps

    # Create system.     
    [system, coordinates] = wcadimer.WCADimer()

    # Form vectors of masses and sqrt(kT/m) for force propagation and velocity randomization.
    print "Creating masses array..."
    nparticles = system.getNumParticles()
    masses = units.Quantity(numpy.zeros([nparticles,3], numpy.float64), units.amu)
    for particle_index in range(nparticles):
        masses[particle_index,:] = system.getParticleMass(particle_index)
    sqrt_kT_over_m = units.Quantity(numpy.zeros([nparticles,3], numpy.float64), units.nanometers / units.picosecond)
    for particle_index in range(nparticles):
        sqrt_kT_over_m[particle_index,:] = units.sqrt(kT / masses[particle_index,0]) # standard deviation of velocity distribution for each coordinate for this atom

    # List all available platforms
    print "Available platforms:"
    for platform_index in range(openmm.Platform.getNumPlatforms()):
        platform = openmm.Platform.getPlatform(platform_index)
        print "%5d %s" % (platform_index, platform.getName())
    print ""

    # Select platform.
    #platform = openmm.Platform.getPlatformByName("CPU")
    platform = openmm.Platform.getPlatformByName("CUDA")
    min_platform = openmm.Platform.getPlatformByName("Reference")

    for prop in platform.getPropertyNames():
        print prop, platform.getPropertyDefaultValue(prop)

kB = units.BOLTZMANN_CONSTANT_kB * units.AVOGADRO_CONSTANT_NA

#=============================================================================================
# DEFAULT PARAMETERS
#=============================================================================================

natoms = 216                             # number of particles

# WCA fluid parameters (argon).
mass     = 39.9 * units.amu              # reference mass
sigma    = 3.4 * units.angstrom          # reference lengthscale
epsilon  = 120.0 * units.kelvin * kB     # reference energy

r_WCA    = 2.**(1./6.) * sigma           # interaction truncation range
tau = units.sqrt(sigma**2 * mass / epsilon) # characteristic timescale

# Simulation conditions.
temperature = 0.824 / (kB / epsilon)     # temperature
kT = kB * temperature                    # thermal energy
beta = 1.0 / kT                          # inverse temperature    
density  = 0.96 / sigma**3               # default density
stable_timestep = 0.001 * tau            # stable timestep
collision_rate = 1 / tau                 # collision rate for Langevin interator

# Dimer potential parameters.
h = 5.0 * kT                             # barrier height
r0 = r_WCA                               # compact state distance
w = 0.5 * r_WCA                          # extended state distance is r0 + 2*w

#=============================================================================================

def runTest(path, systemName, platformName, eps, tolerance, precision, queue):
    """This function runs the test for a single system."""
    print "\nsystemName =", systemName

    # load the system from the xml file
    system = loadXMLFile(path, systemName)
    
    # set Ewald error tolerance to a sufficiently small value so the requested accuracy will be achievable
    for f in system.getForces():
        try:
            f.setEwaldErrorTolerance(tolerance/2)
        except:
            pass

    # read in the particle positions from the .pos file
    positions = loadPosFile(path, systemName)

    numParticles = len(positions)
    print "numParticles =", numParticles, "(from the .pos file)"

    integrator = mm.LangevinIntegrator(300*unit.kelvin, 1/unit.picosecond, 0.002*unit.picoseconds)

    print "mm.Platform.getNumPlatforms =", mm.Platform.getNumPlatforms() 
    platform = mm.Platform.getPlatformByName(platformName)
    print "platform.getName() =", platform.getName()

    print "Building \'context\' on \'platform\'."
    context = mm.Context(system, integrator, platform, properties)
    context.setPositions(positions)
    state0 = context.getState(getForces=True)
    forces0 = state0.getForces()


    # make sure the size of the force vector is equal to the number of particles
    assert (len(forces0)==numParticles)

    # calculate the norm of the forces
    force0NormSum = 0.0*unit.kilojoules**2/unit.mole**2/unit.nanometer**2
    for f in forces0:
        force0NormSum += unit.dot(f,f)
    force0Norm = unit.sqrt(force0NormSum)
    print "force0Norm =", force0Norm

    epsilon = eps*unit.nanometer
    step = epsilon/force0Norm
    print "step =", step

    # perturb the coordinates along the direction of forces0 and evaluate the energy
    context.setPositions([p-2*f*step for p,f in zip(positions, forces0)])
    pe1 = context.getState(getEnergy=True).getPotentialEnergy()
    context.setPositions([p-f*step for p,f in zip(positions, forces0)])
    pe2 = context.getState(getEnergy=True).getPotentialEnergy()
    context.setPositions([p+f*step for p,f in zip(positions, forces0)])
    pe3 = context.getState(getEnergy=True).getPotentialEnergy()
    context.setPositions([p+2*f*step for p,f in zip(positions, forces0)])
    pe4 = context.getState(getEnergy=True).getPotentialEnergy()

    # use a finite difference approximation to calculate the expected force0Norm
    expectedForceNorm = (-pe1+8*pe2-8*pe3+pe4)/(12*epsilon)
    relativeDifference = abs(expectedForceNorm-force0Norm)/force0Norm

    print "pe1 =", pe1
    print "pe2 =", pe2
    print "pe3 =", pe3
    print "pe4 =", pe4
    print "expectedForceNorm =", expectedForceNorm
    print "relativeDifference =", relativeDifference

    # check that the energy is within the desired range
    if relativeDifference > tolerance:
        print "*** ERROR EXCEEDS TOLERANCE ***"
        queue.put(False)
    else:
        queue.put(True)
    print "Test of", systemName, "complete."
    del(context)

    def createMergedTopology(self, alchemical_lambda, mm=None, verbose=False):
        """
        Create a merged topology file with the specified alchemical lambda value for interpolating between molecules A and B.

        ARGUMENTS

        alchemical_lambda (float) - the alchemical lambda in interval [0,1] for interpolating between molecule A (alchemical_lambda = 0) and molecule B (alchemical_lambda = 1),

        OPTIONAL ARGUMENTS

        mm (implements simtk.openmm interface) - OpenMM API implementation to use (default: simtk.openmm)

        TODO

        EXAMPLES

        NOTES

        Merged molecule will contain atom groups in this order:

        S_AB : [atoms in A and B]
        S_A  : [atoms in A and not B]
        S_B  : [atoms in B and not A]

        Masses in S_AB are the geometric product of their masses in A and B.

        """

        # Record timing statistics.
        if verbose: print "Creating merged topology corresponding to alchemical lamdba of %f..." % alchemical_lambda

        # Get local references to systems A and B.
        system_A = self.system_A
        system_B = self.system_B
        corresponding_atoms = self.corresponding_atoms

        #
        # Construct atom sets and correspondence lists for A, B, and merged topology.
        # 

        # Determine number of atoms in each system.
        natoms_A = system_A.getNumParticles()  # number of atoms in molecule A
        natoms_B = system_B.getNumParticles()  # number of atoms in molecule B
        natoms_AandB = len(corresponding_atoms)  # number of atoms in both A and B
        natoms_AnotB = natoms_A - natoms_AandB          # number of atoms in A and not B
        natoms_BnotA = natoms_B - natoms_AandB          # number of atoms in B and not A
        natoms_merged = natoms_AandB + natoms_AnotB + natoms_BnotA  # number of atoms in merged topology

        # Determine sets of atoms shared and not shared.
        atomset_A_AandB = set([ index_A for (index_A, index_B) in corresponding_atoms ]) # atoms in molecule A and B (A molecule numbering)
        atomset_A_AnotB = set(range(natoms_A)) - atomset_A_AandB       # atoms in molecule A and not B (A molecule numbering)
        atomset_A_BnotA = set()

        atomset_B_AandB = set([ index_B for (index_A, index_B) in corresponding_atoms ]) # atoms in molecule B and A (B molecule numbering)
        atomset_B_BnotA = set(range(natoms_B)) - atomset_B_AandB       # atoms in molecule B and not A (B molecule numbering)
        atomset_B_AnotB = set()

        atomset_merged_AandB = set(range(natoms_AandB))                                                               # atoms present in A and B (merged molecule numbering)
        atomset_merged_AnotB = set(range(natoms_AandB, natoms_AandB + natoms_AnotB))                                  # atoms present in A and not B (merged molcule numbering)
        atomset_merged_BnotA = set(range(natoms_AandB + natoms_AnotB, natoms_AandB + natoms_AnotB + natomsBnotA))     # atoms present in B and not A (merged molecule numbering)

        # Construct lists of corresponding atom indices.
        atom_index = dict()

        #
        # Construct merged OpenMM system.
        #

        import simtk.unit as units
        
        # Select OpenMM API implementation to use.
        if not mm:
            import simtk.openmm 
            mm = simtk.openmm

        # Create new System object.
        system = mm.System()

        # Populate merged sytem with atoms.
        # Masses of atoms in both A and B are geometric mean; otherwise standard mass.
        # Add particles in A and B.
        for (index_A, index_B) in corresponding_atoms:
            # Create particle with geometric mean of masses.
            mass_A = system_A.getParticleMass(index_A)
            mass_B = system_B.getParticleMass(index_B)
            mass = units.sqrt(mass_A * mass_B)
            system.addParticle(mass)
        for index_A in atomlist_A_AnotB:
            mass_A = system_A.getParticleMass(index_A)
            system.addParticle(mass_A)
        for index_B in atomlist_B_BnotA:
            mass_B = system_B.getParticleMass(index_B)
            system.addParticle(mass_B)
        
        # Define helper function.
        def find_force(system, classname):
            """
            Find the specified Force object in an OpenMM System by classname.

            ARGUMENTS
#.........这里部分代码省略.........