Python AbsoluteAlchemicalFactory.perturbContext方法代码示例

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import perturbContext [as 别名]
    def _compute_energies(self):
        """
        Compute energies of all replicas at all states.

        TODO

        * We have to re-order Context initialization if we have variable box volume
        * Parallel implementation

        """

        # Create and cache Integrator and Context if needed.
        if not hasattr(self, '_context'):
            self._cache_context()

        logger.debug("Computing energies...")
        start_time = time.time()

        # Retrieve context.
        context = self._context

        if self.mpicomm:
            # MPI version.

            # Compute energies for this node's share of states.
            for state_index in range(self.mpicomm.rank, self.nstates, self.mpicomm.size):
                # Set alchemical state.
                AbsoluteAlchemicalFactory.perturbContext(context, self.states[state_index].alchemical_state)
                for replica_index in range(self.nstates):
                    self.u_kl[replica_index,state_index] = self.states[state_index].reduced_potential(self.replica_positions[replica_index], box_vectors=self.replica_box_vectors[replica_index], context=context)

            # Send final energies to all nodes.
            energies_gather = self.mpicomm.allgather(self.u_kl[:,self.mpicomm.rank:self.nstates:self.mpicomm.size])
            for state_index in range(self.nstates):
                source = state_index % self.mpicomm.size # node with trajectory data
                index = state_index // self.mpicomm.size # index within trajectory batch
                self.u_kl[:,state_index] = energies_gather[source][:,index]

        else:
            # Serial version.
            for state_index in range(self.nstates):
                # Set alchemical state.
                AbsoluteAlchemicalFactory.perturbContext(context, self.states[state_index].alchemical_state)
                for replica_index in range(self.nstates):
                    self.u_kl[replica_index,state_index] = self.states[state_index].reduced_potential(self.replica_positions[replica_index], box_vectors=self.replica_box_vectors[replica_index], context=context)

        end_time = time.time()
        elapsed_time = end_time - start_time
        time_per_energy= elapsed_time / float(self.nstates)**2
        logger.debug("Time to compute all energies %.3f s (%.3f per energy calculation)." % (elapsed_time, time_per_energy))

        return

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import perturbContext [as 别名]
    def _minimize_replica(self, replica_index):
        """
        Minimize the specified replica.

        """
        # Create and cache Integrator and Context if needed.
        if not hasattr(self, '_context'):
            self._cache_context()

        context = self._context

        # Retrieve thermodynamic state.
        state_index = self.replica_states[replica_index] # index of thermodynamic state that current replica is assigned to
        state = self.states[state_index] # thermodynamic state

        # Set alchemical state.
        AbsoluteAlchemicalFactory.perturbContext(context, state.alchemical_state)

        # Set box vectors.
        box_vectors = self.replica_box_vectors[replica_index]
        context.setPeriodicBoxVectors(box_vectors[0,:], box_vectors[1,:], box_vectors[2,:])

        # Set positions.
        positions = self.replica_positions[replica_index]
        context.setPositions(positions)

        # Report initial energy
        logger.debug("Replica %5d/%5d: initial energy %8.3f kT", replica_index, self.nstates, state.reduced_potential(positions, box_vectors=box_vectors, context=context))

        # Minimize energy.
        self.mm.LocalEnergyMinimizer.minimize(context, self.minimize_tolerance, self.minimize_max_iterations)

        # Store final positions
        positions = context.getState(getPositions=True, enforcePeriodicBox=state.system.usesPeriodicBoundaryConditions()).getPositions(asNumpy=True)
        self.replica_positions[replica_index] = positions

        logger.debug("Replica %5d/%5d: final   energy %8.3f kT", replica_index, self.nstates, state.reduced_potential(positions, box_vectors=box_vectors, context=context))

        return

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import perturbContext [as 别名]
    def _create_phase(self, phase, reference_system, positions, atom_indices, thermodynamic_state, protocols=None, options=None, mpicomm=None):
        """
        Create a repex object for a specified phase.

        Parameters
        ----------
        phase : str
           The phase being initialized (one of ['complex', 'solvent', 'vacuum'])
        reference_system : simtk.openmm.System
           The reference system object from which alchemical intermediates are to be construcfted.
        positions : list of simtk.unit.Qunatity objects containing (natoms x 3) positions (as np or lists)
           The list of positions to be used to seed replicas in a round-robin way.
        atom_indices : dict
           atom_indices[phase][component] is the set of atom indices associated with component, where component is ['ligand', 'receptor']
        thermodynamic_state : ThermodynamicState
           Thermodynamic state from which reference temperature and pressure are to be taken.
        protocols : dict of list of AlchemicalState, optional, default=None
           If specified, the alchemical protocol protocols[phase] will be used for phase 'phase' instead of the default.
        options : dict of str, optional, default=None
           If specified, these options will override default repex simulation options.

        """


        # Combine simulation options with defaults to create repex options.
        repex_options = dict(self.default_options.items() + options.items())

        # Make sure positions argument is a list of coordinate snapshots.
        if hasattr(positions, 'unit'):
            # Wrap in list.
            positions = [positions]

        # Check the dimensions of positions.
        for index in range(len(positions)):
            # Make sure it is recast as a np array.
            positions[index] = unit.Quantity(np.array(positions[index] / positions[index].unit), positions[index].unit)

            [natoms, ndim] = (positions[index] / positions[index].unit).shape
            if natoms != reference_system.getNumParticles():
                raise Exception("positions argument must be a list of simtk.unit.Quantity of (natoms,3) lists or np array with units compatible with nanometers.")

        # Create metadata storage.
        metadata = dict()

        # Make a deep copy of the reference system so we don't accidentally modify it.
        reference_system = copy.deepcopy(reference_system)

        # TODO: Use more general approach to determine whether system is periodic.
        is_periodic = self._is_periodic(reference_system)

        # Make sure pressure is None if not periodic.
        if not is_periodic: thermodynamic_state.pressure = None

        # Compute standard state corrections for complex phase.
        metadata['standard_state_correction'] = 0.0
        # TODO: Do we need to include a standard state correction for other phases in periodic boxes?
        if phase == 'complex-implicit':
            # Impose restraints for complex system in implicit solvent to keep ligand from drifting too far away from receptor.
            logger.debug("Creating receptor-ligand restraints...")
            reference_positions = positions[0]
            if self.restraint_type == 'harmonic':
                restraints = HarmonicReceptorLigandRestraint(thermodynamic_state, reference_system, reference_positions, atom_indices['receptor'], atom_indices['ligand'])
            elif self.restraint_type == 'flat-bottom':
                restraints = FlatBottomReceptorLigandRestraint(thermodynamic_state, reference_system, reference_positions, atom_indices['receptor'], atom_indices['ligand'])
            else:
                raise Exception("restraint_type of '%s' is not supported." % self.restraint_type)

            force = restraints.getRestraintForce() # Get Force object incorporating restraints
            reference_system.addForce(force)
            metadata['standard_state_correction'] = restraints.getStandardStateCorrection() # standard state correction in kT
        elif phase == 'complex-explicit':
            # For periodic systems, we do not use a restraint, but must add a standard state correction for the box volume.
            # TODO: What if the box volume fluctuates during the simulation?
            box_vectors = reference_system.getDefaultPeriodicBoxVectors()
            box_volume = thermodynamic_state._volume(box_vectors)
            STANDARD_STATE_VOLUME = 1660.53928 * unit.angstrom**3
            metadata['standard_state_correction'] = np.log(STANDARD_STATE_VOLUME / box_volume) # TODO: Check sign.

        # Use default alchemical protocols if not specified.
        if not protocols:
            protocols = self.default_protocols

        # Create alchemically-modified states using alchemical factory.
        logger.debug("Creating alchemically-modified states...")
        #factory = AbsoluteAlchemicalFactory(reference_system, ligand_atoms=atom_indices['ligand'], test_positions=positions[0], platform=repex_options['platform'])
        factory = AbsoluteAlchemicalFactory(reference_system, ligand_atoms=atom_indices['ligand'])
        alchemical_states = protocols[phase]
        alchemical_system = factory.alchemically_modified_system
        thermodynamic_state.system = alchemical_system

        # Check systems for finite energies.
        finite_energy_check = False
        if finite_energy_check:
            logger.debug("Checking energies are finite for all alchemical systems.")
            integrator = openmm.VerletIntegrator(1.0 * unit.femtosecond)
            context = openmm.Context(alchemical_system, integrator)
            context.setPositions(positions[0])
            for alchemical_state in alchemical_states:
                AbsoluteAlchemicalFactory.perturbContext(context, alchemical_state)
                potential = context.getState(getEnergy=True).getPotentialEnergy()
#.........这里部分代码省略.........

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import perturbContext [as 别名]
    def _propagate_replica(self, replica_index):
        """
        Attempt a Monte Carlo rotation/translation move followed by dynamics.

        """

        # Create and cache Integrator and Context if needed.
        if not hasattr(self, '_context'):
            self._cache_context()

        # Retrieve state.
        state_index = self.replica_states[replica_index] # index of thermodynamic state that current replica is assigned to
        state = self.states[state_index] # thermodynamic state

        # Retrieve cached integrator and context.
        integrator = self._integrator
        context = self._context

        # Set thermodynamic parameters for this state.
        integrator.setTemperature(state.temperature)
        integrator.setRandomNumberSeed(int(np.random.randint(0, MAX_SEED)))
        if state.temperature and state.pressure:
            forces = { state.system.getForce(index).__class__.__name__ : state.system.getForce(index) for index in range(state.system.getNumForces()) }

            if 'MonteCarloAnisotropicBarostat' in forces:
                raise Exception('MonteCarloAnisotropicBarostat is unsupported.')

            if 'MonteCarloBarostat' in forces:
                barostat = forces['MonteCarloBarostat']
                # Set temperature and pressure.
                barostat.setTemperature(state.temperature)
                barostat.setDefaultPressure(state.pressure)
                context.setParameter(barostat.Pressure(), state.pressure) # must be set in context
                barostat.setRandomNumberSeed(int(np.random.randint(0, MAX_SEED)))

        # Set alchemical state.
        AbsoluteAlchemicalFactory.perturbContext(context, state.alchemical_state)

        #
        # Attempt a Monte Carlo rotation/translation move.
        #

        # Attempt gaussian trial displacement with stddev 'self.displacement_sigma'.
        # TODO: Can combine these displacements and/or use cached potential energies to speed up this phase.
        # TODO: Break MC displacement and rotation into member functions and write separate unit tests.
        if self.mc_displacement and (self.mc_atoms is not None):
            initial_time = time.time()
            # Store original positions and energy.
            original_positions = self.replica_positions[replica_index]
            u_old = state.reduced_potential(original_positions, context=context)
            # Make symmetric Gaussian trial displacement of ligand.
            perturbed_positions = self.propose_displacement(self.displacement_sigma, original_positions, self.mc_atoms)
            u_new = state.reduced_potential(perturbed_positions, context=context)
            # Accept or reject with Metropolis criteria.
            du = u_new - u_old
            if (du <= 0.0) or (np.random.rand() < np.exp(-du)):
                self.displacement_trials_accepted += 1
                self.replica_positions[replica_index] = perturbed_positions
            #print "translation du = %f (%d)" % (du, self.displacement_trials_accepted)
            # Print timing information.
            final_time = time.time()
            elapsed_time = final_time - initial_time
            self.displacement_trial_time += elapsed_time

        # Attempt random rotation of ligand.
        if self.mc_rotation and (self.mc_atoms is not None):
            initial_time = time.time()
            # Store original positions and energy.
            original_positions = self.replica_positions[replica_index]
            u_old = state.reduced_potential(original_positions, context=context)
            # Compute new potential.
            perturbed_positions = self.propose_rotation(original_positions, self.mc_atoms)
            u_new = state.reduced_potential(perturbed_positions, context=context)
            du = u_new - u_old
            if (du <= 0.0) or (np.random.rand() < np.exp(-du)):
                self.rotation_trials_accepted += 1
                self.replica_positions[replica_index] = perturbed_positions
            #print "rotation du = %f (%d)" % (du, self.rotation_trials_accepted)
            # Accumulate timing information.
            final_time = time.time()
            elapsed_time = final_time - initial_time
            self.rotation_trial_time += elapsed_time

        #
        # Propagate with dynamics.
        #

        start_time = time.time()

        # Set box vectors.
        box_vectors = self.replica_box_vectors[replica_index]
        context.setPeriodicBoxVectors(box_vectors[0,:], box_vectors[1,:], box_vectors[2,:])
        # Set positions.
        positions = self.replica_positions[replica_index]
        context.setPositions(positions)
        setpositions_end_time = time.time()
        # Assign Maxwell-Boltzmann velocities.
        context.setVelocitiesToTemperature(state.temperature, int(np.random.randint(0, MAX_SEED)))
        setvelocities_end_time = time.time()
        # Run dynamics.
#.........这里部分代码省略.........

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import perturbContext [as 别名]
    def _propagate_replica(self, replica_index):
        """
        Attempt a Monte Carlo rotation/translation move followed by dynamics.

        """

        # Create and cache Integrator and Context if needed.
        if not hasattr(self, '_context'):
            self._cache_context()

        # Retrieve state.
        state_index = self.replica_states[replica_index] # index of thermodynamic state that current replica is assigned to
        state = self.states[state_index] # thermodynamic state

        # Retrieve cached integrator and context.
        integrator = self._integrator
        context = self._context

        # Set thermodynamic parameters for this state.
        integrator.setTemperature(state.temperature)
        integrator.setRandomNumberSeed(int(np.random.randint(0, MAX_SEED)))
        if state.temperature and state.pressure:
            forces = { state.system.getForce(index).__class__.__name__ : state.system.getForce(index) for index in range(state.system.getNumForces()) }

            if 'MonteCarloAnisotropicBarostat' in forces:
                raise Exception('MonteCarloAnisotropicBarostat is unsupported.')

            if 'MonteCarloBarostat' in forces:
                barostat = forces['MonteCarloBarostat']
                # Set temperature and pressure.
                try:
                    barostat.setDefaultTemperature(state.temperature)
                except AttributeError:  # versions previous to OpenMM0.8
                    barostat.setTemperature(state.temperature)
                barostat.setDefaultPressure(state.pressure)
                context.setParameter(barostat.Pressure(), state.pressure) # must be set in context
                barostat.setRandomNumberSeed(int(np.random.randint(0, MAX_SEED)))

        # Set alchemical state.
        AbsoluteAlchemicalFactory.perturbContext(context, state.alchemical_state)

        # Set box vectors.
        box_vectors = self.replica_box_vectors[replica_index]
        context.setPeriodicBoxVectors(box_vectors[0,:], box_vectors[1,:], box_vectors[2,:])

        # Check if initial potential energy is NaN.
        reduced_potential = state.reduced_potential(self.replica_positions[replica_index], box_vectors=box_vectors, context=context)
        if np.isnan(reduced_potential):
            raise Exception('Initial potential for replica %d state %d is NaN before Monte Carlo displacement/rotation' % (replica_index, state_index))

        #
        # Attempt a Monte Carlo rotation/translation move.
        #

        # Attempt gaussian trial displacement with stddev 'self.displacement_sigma'.
        # TODO: Can combine these displacements and/or use cached potential energies to speed up this phase.
        # TODO: Break MC displacement and rotation into member functions and write separate unit tests.
        if self.mc_displacement and (self.mc_atoms is not None):
            initial_time = time.time()
            # Store original positions and energy.
            original_positions = self.replica_positions[replica_index]
            u_old = state.reduced_potential(original_positions, box_vectors=box_vectors, context=context)
            # Make symmetric Gaussian trial displacement of ligand.
            perturbed_positions = self.propose_displacement(self.displacement_sigma, original_positions, self.mc_atoms)
            u_new = state.reduced_potential(perturbed_positions, box_vectors=box_vectors, context=context)
            # Accept or reject with Metropolis criteria.
            du = u_new - u_old
            if (not np.isnan(u_new)) and ((du <= 0.0) or (np.random.rand() < np.exp(-du))):
                self.displacement_trials_accepted += 1
                self.replica_positions[replica_index] = perturbed_positions
            #print "translation du = %f (%d)" % (du, self.displacement_trials_accepted)
            # Print timing information.
            final_time = time.time()
            elapsed_time = final_time - initial_time
            self.displacement_trial_time += elapsed_time

        # Attempt random rotation of ligand.
        if self.mc_rotation and (self.mc_atoms is not None):
            initial_time = time.time()
            # Store original positions and energy.
            original_positions = self.replica_positions[replica_index]
            u_old = state.reduced_potential(original_positions, box_vectors=box_vectors, context=context)
            # Compute new potential.
            perturbed_positions = self.propose_rotation(original_positions, self.mc_atoms)
            u_new = state.reduced_potential(perturbed_positions, box_vectors=box_vectors, context=context)
            du = u_new - u_old
            if (not np.isnan(u_new)) and ((du <= 0.0) or (np.random.rand() < np.exp(-du))):
                self.rotation_trials_accepted += 1
                self.replica_positions[replica_index] = perturbed_positions
            #print "rotation du = %f (%d)" % (du, self.rotation_trials_accepted)
            # Accumulate timing information.
            final_time = time.time()
            elapsed_time = final_time - initial_time
            self.rotation_trial_time += elapsed_time

        #
        # Propagate with dynamics.
        #

        start_time = time.time()
#.........这里部分代码省略.........

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import perturbContext [as 别名]

#.........这里部分代码省略.........
            metadata['standard_state_correction'] = - np.log(V0 / box_volume)
        elif is_complex_implicit:
            # For implicit solvent/vacuum complex systems, we require a restraint
            # to keep the ligand from drifting too far away from receptor.
            raise ValueError('A receptor-ligand system in implicit solvent or '
                             'vacuum requires a restraint.')

        # Create alchemically-modified states using alchemical factory.
        logger.debug("Creating alchemically-modified states...")
        try:
            alchemical_indices = atom_indices['ligand_counterions'] + atom_indices['ligand']
        except KeyError:
            alchemical_indices = atom_indices['ligand']
        factory = AbsoluteAlchemicalFactory(reference_system, ligand_atoms=alchemical_indices,
                                            **self._alchemy_parameters)
        alchemical_system = factory.alchemically_modified_system
        thermodynamic_state.system = alchemical_system

        # Create the expanded cutoff decoupled state
        if fully_interacting_expanded_state is None:
            noninteracting_expanded_state = None
        else:
            # Create the system for noninteracting
            expanded_factory = AbsoluteAlchemicalFactory(fully_interacting_expanded_state.system,
                                                         ligand_atoms=alchemical_indices,
                                                         **self._alchemy_parameters)
            noninteracting_expanded_system = expanded_factory.alchemically_modified_system
            # Set all USED alchemical interactions to the decoupled state
            alchemical_state = alchemical_states[-1]
            AbsoluteAlchemicalFactory.perturbSystem(noninteracting_expanded_system, alchemical_state)

            # Construct thermodynamic states
            noninteracting_expanded_state = copy.deepcopy(thermodynamic_state)
            noninteracting_expanded_state.system = noninteracting_expanded_system

        # Check systems for finite energies.
        # TODO: Refactor this into another function.
        finite_energy_check = False
        if finite_energy_check:
            logger.debug("Checking energies are finite for all alchemical systems.")
            integrator = openmm.VerletIntegrator(1.0 * unit.femtosecond)
            context = openmm.Context(alchemical_system, integrator)
            context.setPositions(positions[0])
            for index, alchemical_state in enumerate(alchemical_states):
                AbsoluteAlchemicalFactory.perturbContext(context, alchemical_state)
                potential = context.getState(getEnergy=True).getPotentialEnergy()
                if np.isnan(potential / unit.kilocalories_per_mole):
                    raise Exception("Energy for system %d is NaN." % index)
            del context, integrator
            logger.debug("All energies are finite.")

        # Randomize ligand position if requested, but only for implicit solvent systems.
        if self._randomize_ligand and is_complex_implicit:
            logger.debug("Randomizing ligand positions and excluding overlapping configurations...")
            randomized_positions = list()
            nstates = len(alchemical_states)
            for state_index in range(nstates):
                positions_index = np.random.randint(0, len(positions))
                current_positions = positions[positions_index]
                new_positions = ModifiedHamiltonianExchange.randomize_ligand_position(current_positions,
                                                                                      atom_indices['receptor'], atom_indices['ligand'],
                                                                                      self._randomize_ligand_sigma_multiplier * restraints.getReceptorRadiusOfGyration(),
                                                                                      self._randomize_ligand_close_cutoff)
                randomized_positions.append(new_positions)
            positions = randomized_positions
        if self._randomize_ligand and is_complex_explicit:
            logger.warning("Ligand randomization requested, but will not be performed for explicit solvent simulations.")

        # Identify whether any atoms will be displaced via MC, unless option is turned off.
        mc_atoms = None
        if self._mc_displacement_sigma:
            mc_atoms = list()
            if 'ligand' in atom_indices:
                mc_atoms = atom_indices['ligand']

        # Set up simulation.
        # TODO: Support MPI initialization?
        logger.debug("Creating replica exchange object...")
        store_filename = os.path.join(self._store_directory, alchemical_phase.name + '.nc')
        self._store_filenames[alchemical_phase.name] = store_filename
        simulation = ModifiedHamiltonianExchange(store_filename, platform=self._platform)
        simulation.create(thermodynamic_state, alchemical_states, positions,
                          displacement_sigma=self._mc_displacement_sigma, mc_atoms=mc_atoms,
                          options=repex_parameters, metadata=metadata,
                          fully_interacting_expanded_state = fully_interacting_expanded_state,
                          noninteracting_expanded_state = noninteracting_expanded_state)

        # Initialize simulation.
        # TODO: Use the right scheme for initializing the simulation without running.
        #logger.debug("Initializing simulation...")
        #simulation.run(0)

        # Clean up simulation.
        del simulation

        # Add to list of phases that have been set up.
        self._phases.append(alchemical_phase.name)
        self._phases.sort()  # sorting ensures all MPI processes run the same phase

        return

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import perturbContext [as 别名]

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

        # Compute standard state corrections for complex phase.
        metadata['standard_state_correction'] = 0.0
        # TODO: Do we need to include a standard state correction for other phases in periodic boxes?
        if is_complex_implicit:
            # Impose restraints for complex system in implicit solvent to keep ligand from drifting too far away from receptor.
            logger.debug("Creating receptor-ligand restraints...")
            reference_positions = positions[0]
            if self._restraint_type == 'harmonic':
                restraints = HarmonicReceptorLigandRestraint(thermodynamic_state, reference_system, reference_positions, atom_indices['receptor'], atom_indices['ligand'])
            elif self._restraint_type == 'flat-bottom':
                restraints = FlatBottomReceptorLigandRestraint(thermodynamic_state, reference_system, reference_positions, atom_indices['receptor'], atom_indices['ligand'])
            else:
                raise Exception("restraint_type of '%s' is not supported." % self._restraint_type)

            force = restraints.getRestraintForce() # Get Force object incorporating restraints
            reference_system.addForce(force)
            metadata['standard_state_correction'] = restraints.getStandardStateCorrection() # standard state correction in kT
        elif is_complex_explicit:
            # For periodic systems, we do not use a restraint, but must add a standard state correction for the box volume.
            # TODO: What if the box volume fluctuates during the simulation?
            box_vectors = reference_system.getDefaultPeriodicBoxVectors()
            box_volume = thermodynamic_state._volume(box_vectors)
            STANDARD_STATE_VOLUME = 1660.53928 * unit.angstrom**3
            metadata['standard_state_correction'] = - np.log(STANDARD_STATE_VOLUME / box_volume)

        # Create alchemically-modified states using alchemical factory.
        logger.debug("Creating alchemically-modified states...")
        try:
            alchemical_indices = atom_indices['ligand_counterions'] + atom_indices['ligand']
        except KeyError:
            alchemical_indices = atom_indices['ligand']
        factory = AbsoluteAlchemicalFactory(reference_system, ligand_atoms=alchemical_indices,
                                            **self._alchemy_parameters)
        alchemical_system = factory.alchemically_modified_system
        thermodynamic_state.system = alchemical_system

        # Check systems for finite energies.
        # TODO: Refactor this into another function.
        finite_energy_check = False
        if finite_energy_check:
            logger.debug("Checking energies are finite for all alchemical systems.")
            integrator = openmm.VerletIntegrator(1.0 * unit.femtosecond)
            context = openmm.Context(alchemical_system, integrator)
            context.setPositions(positions[0])
            for alchemical_state in alchemical_states:
                AbsoluteAlchemicalFactory.perturbContext(context, alchemical_state)
                potential = context.getState(getEnergy=True).getPotentialEnergy()
                if np.isnan(potential / unit.kilocalories_per_mole):
                    raise Exception("Energy for system %d is NaN." % index)
            del context, integrator
            logger.debug("All energies are finite.")

        # Randomize ligand position if requested, but only for implicit solvent systems.
        if self._randomize_ligand and is_complex_implicit:
            logger.debug("Randomizing ligand positions and excluding overlapping configurations...")
            randomized_positions = list()
            nstates = len(alchemical_states)
            for state_index in range(nstates):
                positions_index = np.random.randint(0, len(positions))
                current_positions = positions[positions_index]
                new_positions = ModifiedHamiltonianExchange.randomize_ligand_position(current_positions,
                                                                                      atom_indices['receptor'], atom_indices['ligand'],
                                                                                      self._randomize_ligand_sigma_multiplier * restraints.getReceptorRadiusOfGyration(),
                                                                                      self._randomize_ligand_close_cutoff)
                randomized_positions.append(new_positions)
            positions = randomized_positions
        if self._randomize_ligand and is_complex_explicit:
            logger.warning("Ligand randomization requested, but will not be performed for explicit solvent simulations.")

        # Identify whether any atoms will be displaced via MC, unless option is turned off.
        mc_atoms = None
        if self._mc_displacement_sigma:
            mc_atoms = list()
            if 'ligand' in atom_indices:
                mc_atoms = atom_indices['ligand']

        # Set up simulation.
        # TODO: Support MPI initialization?
        logger.debug("Creating replica exchange object...")
        store_filename = os.path.join(self._store_directory, alchemical_phase.name + '.nc')
        self._store_filenames[alchemical_phase.name] = store_filename
        simulation = ModifiedHamiltonianExchange(store_filename)
        simulation.create(thermodynamic_state, alchemical_states, positions,
                          displacement_sigma=self._mc_displacement_sigma, mc_atoms=mc_atoms,
                          options=repex_parameters, metadata=metadata,
                          fully_interacting_state=fully_interacting_state)

        # Initialize simulation.
        # TODO: Use the right scheme for initializing the simulation without running.
        #logger.debug("Initializing simulation...")
        #simulation.run(0)

        # Clean up simulation.
        del simulation

        # Add to list of phases that have been set up.
        self._phases.append(alchemical_phase.name)

        return