def _identifyRestrainedAtoms(self, coordinates, minimum_mass=10.0*units.amu, natoms_to_pick=5, minimum_distance=5.0*units.angstroms, maximum_distance=15.0*units.angstroms, verbose=False):
        """
        Identify atoms to restraint using a simple automated algorithm.

        ARGUMENTS
    
        restraint_filename        restraint filename
        system                    system object to which restraints are to be added
        positions                 reference configuration from which to extract equilibrium distances
        receptor_atoms            atom(s) in receptor (zero-indexed)
        ligand_atoms              atom(s) in ligand (zero-indexed)
        
        OPTIONAL ARGUMENTS
        
        minimum_mass (simtk.unit.Quantity compatible with simtk.unit.amu) - minimum atomic mass for heavy atoms
        natoms_to_pick (int) - number of receptor atoms to pick
        minimum_distance (simtk.unit.Quantity compatible with simtk.unit.angstrom) - minimum distance from ligand atom to receptor atoms to use as restraints
        maximum_distance (simtk.unit.Quantity compatible with simtk.unit.angstrom) - maximum distance from ligand atom to receptor atoms to use as restraints

        verbose (boolean) - if True, more debug information is printed during standard execution (default: False)

        RETURNS

        restrained_receptor_atoms (set integers)
        restrained_ligand_atoms (set of integers)
        
        NOTES
        
        The strategy used here is to first choose a ligand heavy atom near the centroid of the initial ligand geometry.  
        All protein heavy atoms within a given distance range of this designated ligand atom (e.g. between 5 and 15 A) are then identified,
        and a subset (e.g. 5 atoms) are selected that are far apart from each other.

        TESTS
        
        >>> # Create a test system.
        >>> import simtk.pyopenmm.extras.testsystems as testsystems
        >>> [system, coordinates] = testsystems.LysozymeImplicit()
        >>> # Initialize receptor-ligand restraints.
        >>> receptor_atoms = range(0,2603)
        >>> ligand_atoms = range(2603,2621)
        >>> sigma = 1.0 * units.nanometer
        >>> temperature = 298.0 * units.kelvin
        >>> restraints = ReceptorLigandRestraints(system, ligand_atoms, receptor_atoms, temperature)
        >>> # Identify atoms to be restrained.
        >>> [restrained_receptor_atoms, restrained_ligand_atoms] = restraints._identifyRestrainedAtoms(coordinates)        
        >>> print restrained_receptor_atoms
        >>> print restrained_ligand_atoms

        """
    
        # PARAMETERS
        minimum_mass = 10.0 * units.amu # minimum atomic mass for heavy atoms
        natoms_to_pick = 5 # number of receptor atoms to pick
        minimum_distance = 5.0 * units.angstroms # minimum distance from ligand atom to receptor atoms to use as restraints
        maximum_distance = 15.0 * units.angstroms # maximum distance from ligand atom to receptor atoms to use as restraints

        # Compute centroid of ligand.
        centroid = positions[ligand_atoms,:].mean(axis=0)
        if verbose: print "Ligand centroid is at (%8.3f, %8.3f, %8.3f) nm\n" % (centroid[0] / units.nanometers, centroid[1] / units.nanometers, centroid[2] / units.nanometers)

        restrained_ligand_atoms = set()
        restrained_receptor_atoms = set()
    
        # Pick a ligand heavy atom closest to the ligand center of geometry.    
        closest_distance = None
        for ligand_atom in ligand_atoms:
            delta = positions[ligand_atom,:] - centroid
            distance = units.sqrt(delta**2)
            if ((closest_distance is None) or (distance < closest_distance)) and (system.getParticleMass(ligand_atom) >= minimum_mass):
                closest_distance = distance
                closest_atom = ligand_atom
        restrained_ligand_atoms.add(closest_atom)
        if verbose: print "Identified ligand atom %d as atom closest to centroid (distance %.3f nm)\n" % (closest_atom, closest_distance / units.nanometers)
    
        # First pass: Identify all receptor heavy atoms within specified distance range of chosen ligand atom.
        atom_pool = list() # pool of all receptor atoms that fit criteria
        ligand_atom_position = positions[closest_atom,:] # ligand atom position
        for receptor_atom in receptor_atoms:
            delta = ligand_atom_position - positions[receptor_atom,:];
            distance = units.sqrt(delta**2)
            if ((distance >= minimum_distance) and (distance < maximum_distance) and (system.getParticleMass(receptor_atom) >= minimum_mass)):
                atom_pool.append(receptor_atom);

        # Sanity check: Did we identify the expected number of receptor atoms?
        if (len(atom_pool) < natoms_to_pick):
            raise Exception("identifyRestrainedAtoms: Not enough receptor atoms match criterion.  Expected at least %d, found only %lu.\n" % (natoms_to_pick, len(atom_pool)))

        # Identify a given number of receptor heavy atoms that are close to the designated ligand atom but well-separated from each other. These will be restrained.
        # TODO: Make this less ugly, using simplified routines to find closest and farthest atoms in sets.
        restrained_receptor_atoms = [atom_index.pop(0)] # list of atoms in receptor where restraints will terminate
        while (len(restrained_receptor_atoms) < natoms_to_pick):
            # Find atom in 'atom_pool' that is farthest away from all atoms in 'restrained_receptor_atoms'.
            farthest_distance = None
            farthest_atom = None
            for receptor_atom in atom_pool:
                closest_distance = None
                closest_atom = None
                for restrained_atom in restrained_receptor_atoms:
                    delta = positions[receptor_atom,:] - positions[restrained_atom,:]
                    distance = units.sqrt(delta**2)
                    if ((closest_distance is None) or (distance < closest_distance)):
                        closest_distance = distance
                        closest_atom = restrained_atom
                if (fathest_distance is None) or (closest_distance > farthest_distance):
                    fathest_distance = closest_distance
                    farthest_atom = receptor_atom
            # Add the farthest receptor atom to the set of restrained atoms.
            restrained_receptor_atoms.append(farthest_atom)
            atom_pool.remove(farthest_atom)

        if verbose: 
            print "Identified atoms to restrain:"
            print "  receptor: %s\n" % str(restrained_receptor_atoms)
            print "  ligand: %s\n" % str(restrined_ligand_atoms)

        return [restrained_receptor_atoms, restrained_ligand_atoms]

    def createRestraintsAutomatically(reference_positions):
        """
        Create restraints automatically based on the given reference positions.
        
        ARGUMENTS

        reference_positions (simtk.unit.Quantity array compatible with simtk.unit.nanometers) - positions to use in determining restraints

        """
        [restrained_receptor_atoms, restrained_ligand_atoms] = self.identifyRestrainedAtoms(reference_positions)
        return self.imposeRestraints(restrained_receptor_atoms, restrained_ligand_atoms)
            
#=============================================================================================
# Harmonic protein-ligand restraint.
#=============================================================================================

class HarmonicReceptorLigandRestraints(object):
    """

    """
    def __init__(self, system, ligand_atoms, receptor_atoms, temperature):
        ReceptorLigandRestraints.__init__(self, system, ligand_atoms, receptor_atoms, temperature)
    
    def getStandardStateCorrection():
        """
        Compute the standard state correction for a single harmonic distance restraint with nonzero equilibrium length.

        ARGUMENTS
        beta (simtk.unit.Quantity) - inverse temperature (1/energy)
        K (simtk.unit.Quantity) - spring constant (length/distance)
        r0 (simtk.unit.Quantity) - equilibrium spring distance (length)

        OPTIONAL ARGUMENTS
        verbose (bool) - if True, will print extra debug information (default: False)

        RETURN VALUES
        DeltaG (float) - computed standard-state correction in dimensionless units (kT);
        Equivalent to the free energy of releasing restraints and confining into a box of standard
        state size.

        EXAMPLE

        Compute standard state correction for spring with K = 10.0 kcal/mol/nm, r0 = 1.0 nm, at 298 K.

        >>> temperature = 298.0 * units.kelvin                              # temperature
        >>> kB = units.BOLTZMANN_CONSTANT_kB                                # Boltzmann constant
        >>> kT = kB * temperature                                           # thermal energy
        >>> beta = 1.0 / kT                                                 # inverse temperature
        >>> K    = 10.0 * units.kilocalories_per_mole / units.nanometers**2 # spring constant
        >>> r0   = 1.0 * units.nanometers                                   # equilibrium distance
        >>> DeltaG = _standard_state_correction(beta, K, r0)                # compute standard state correction in dimensionless units (kT)
        >>> print 'DeltaG = %8.3f kT' % DeltaG
        DeltaG =    4.620 kT

        """

        # Compute standard deviation of equilibrium distribution.
        # \sigma^2 = 1 / (\beta K)
        sigma = 1.0 / units.sqrt(beta * K) # standard deviation of equilibrium distribution in 1D
        if verbose: print "sigma = %f.3 nm" % (sigma / units.nanometers)

        # Use numerical quadrature to integrate partition function (equal to effective volume of sampled shell).
        r_min = 0.0 * units.nanometers # initial value for integration
        r_max = r0 + 10.0 * sigma # maximum radius to integrate to
        integrand = lambda r : 4.0 * math.pi * r**2 * math.exp(-beta * (K/2.0) * (r*units.nanometers - r0)**2)
        (shell_volume, shell_volume_error) = scipy.integrate.quad(lambda r : integrand(r), r_min / units.nanometers, r_max / units.nanometers) * units.nanometers**3 # integrate shell volume
        if verbose: print "shell_volume = %f nm^3" % (shell_volume / units.nanometers**3)
    
        # Compute standard-state volume for a single molecule in a box of size (1 L) / (avogadros number)
        liter = 1000.0 * units.centimeters**3 # one liter
        box_volume = liter / units.AVOGADRO_NUMBER_NA # standard state volume
        if verbose: print "box_volume = %f nm^3" % (box_volume / units.nanometers**3)
    
        # Compute standard state correction for releasing shell restraints into standard-state box (in units of kT).
        DeltaG = - math.log(box_volume / shell_volume)
        if verbose: print "Standard state correction: %.3f kT" % DeltaG
    
        # Return standard state correction (in kT).
        return DeltaG