Python alchemy.AbsoluteAlchemicalFactory类代码示例

def alchemical_factory_check(reference_system, positions, platform_name=None, precision=None, factory_args=None):
    """
    Compare energies of reference system and fully-interacting alchemically modified system.

    ARGUMENTS

    reference_system : simtk.openmm.System
       The reference System object to compare with
    positions : simtk.unit.Quantity of dimentsion [natoms,3] with units compatible with angstroms
       The positions to assess energetics for
    precision : str, optional, default=None
       Precision model, or default if not specified. ('single', 'double', 'mixed')
    factory_args : dict(), optional, default=None
       Arguments passed to AbsoluteAlchemicalFactory.

    """

    # Create a factory to produce alchemical intermediates.
    logger.info("Creating alchemical factory...")
    initial_time = time.time()
    factory = AbsoluteAlchemicalFactory(reference_system, **factory_args)
    final_time = time.time()
    elapsed_time = final_time - initial_time
    logger.info("AbsoluteAlchemicalFactory initialization took %.3f s" % elapsed_time)
    platform = None
    if platform_name:
        platform = openmm.Platform.getPlatformByName(platform_name)
    alchemical_system = factory.createPerturbedSystem()
    compareSystemEnergies(positions, [reference_system, alchemical_system], ['reference', 'alchemical'], platform=platform, precision=precision)
    return

def compare_platforms(system, positions, factory_args=dict()):
    # Create annihilated version of vacuum system.
    factory = AbsoluteAlchemicalFactory(system, **factory_args)
    alchemical_system = factory.createPerturbedSystem()

    def set_lambda(alchemical_system, lambda_value):
        alchemical_state = AlchemicalState(lambda_electrostatics=lambda_value, lambda_sterics=lambda_value, lambda_torsions=lambda_value)
        AbsoluteAlchemicalFactory.perturbSystem(alchemical_system, alchemical_state)

    # Compare energies
    energies = dict()
    platform_names = list()
    for platform_index in range(openmm.Platform.getNumPlatforms()):
        platform = openmm.Platform.getPlatform(platform_index)
        platform_name = platform.getName()
        if platform_name != 'Reference':
            platform_names.append(platform_name)
        energies[platform_name] = dict()
        energies[platform_name]['full'] = compute_energy(system, positions, platform=platform)
        set_lambda(alchemical_system, 1.0)
        energies[platform_name]['lambda = 1'] = compute_energy(alchemical_system, positions, platform=platform)
        set_lambda(alchemical_system, 0.0)
        energies[platform_name]['lambda = 0'] = compute_energy(alchemical_system, positions, platform=platform)

    # Check deviations.
    for platform_name in platform_names:
        for energy_name in ['full', 'lambda = 1', 'lambda = 0']:
            delta = energies[platform_name][energy_name] - energies['Reference'][energy_name]
            if (abs(delta) > MAX_DELTA):
                raise Exception("Maximum allowable deviation on platform %s exceeded (was %.8f kcal/mol; allowed %.8f kcal/mol); test failed." % (platform_name, delta / unit.kilocalories_per_mole, MAX_DELTA / unit.kilocalories_per_mole))

    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

    def __init__(self, store_directory, verbose=False):
        """
        Initialize YANK object with default parameters.

        Parameters
        ----------
        store_directory : str
           The storage directory in which output NetCDF files are read or written.
        verbose : bool, optional, default=False
           If True, will turn on verbose output.

        """

        # Record that we are not yet initialized.
        self._initialized = False

        # Store output directory.
        self._store_directory = store_directory

        # Public attributes.
        self.verbose = verbose
        self.restraint_type = 'flat-bottom' # default to a flat-bottom restraint between the ligand and receptor
        self.randomize_ligand = True
        self.randomize_ligand_sigma_multiplier = 2.0
        self.randomize_ligand_close_cutoff = 1.5 * unit.angstrom # TODO: Allow this to be specified by user.
        self.mc_displacement_sigma = 10.0 * unit.angstroms

        # Set internal variables.
        self._phases = list()
        self._store_filenames = dict()

        # Default alchemical protocols.
        self.default_protocols = dict()
        self.default_protocols['vacuum'] = AbsoluteAlchemicalFactory.defaultVacuumProtocol()
        self.default_protocols['solvent-implicit'] = AbsoluteAlchemicalFactory.defaultSolventProtocolImplicit()
        self.default_protocols['complex-implicit'] = AbsoluteAlchemicalFactory.defaultComplexProtocolImplicit()
        self.default_protocols['solvent-explicit'] = AbsoluteAlchemicalFactory.defaultSolventProtocolExplicit()
        self.default_protocols['complex-explicit'] = AbsoluteAlchemicalFactory.defaultComplexProtocolExplicit()

        # Default options for repex.
        self.default_options = dict()
        self.default_options['number_of_equilibration_iterations'] = 0
        self.default_options['number_of_iterations'] = 100
        self.default_options['verbose'] = self.verbose
        self.default_options['timestep'] = 2.0 * unit.femtoseconds
        self.default_options['collision_rate'] = 5.0 / unit.picoseconds
        self.default_options['minimize'] = False
        self.default_options['show_mixing_statistics'] = True # this causes slowdown with iteration and should not be used for production
        self.default_options['platform_names'] = None
        self.default_options['displacement_sigma'] = 1.0 * unit.nanometers # attempt to displace ligand by this stddev will be made each iteration

        return

def test_softcore_parameters():
    """
    Testing alchemical slave functions
    """
    alchemical_functions = { 'lambda_sterics' : 'lambda', 'lambda_electrostatics' : 'lambda', 'lambda_bonds' : 'lambda', 'lambda_angles' : 'lambda', 'lambda_torsions' : 'lambda' }
    name = 'Lennard-Jones fluid with dispersion correction'
    test_system = copy.deepcopy(test_systems[name])
    reference_system = test_system['test'].system
    positions = test_system['test'].positions
    factory_args = test_system['factory_args']
    factory_args.update({ 'softcore_alpha' : 1.0, 'softcore_beta' : 1.0, 'softcore_a' : 1.0, 'softcore_b' : 1.0, 'softcore_c' : 1.0, 'softcore_d' : 1.0, 'softcore_e' : 1.0, 'softcore_f' : 1.0 })
    factory = AbsoluteAlchemicalFactory(reference_system, **factory_args)
    alchemical_system = factory.createPerturbedSystem()
    compareSystemEnergies(positions, [reference_system, alchemical_system], ['reference', 'alchemical'])

def alchemical_factory_check(
    reference_system,
    positions,
    receptor_atoms,
    ligand_atoms,
    platform_name=None,
    annihilate_electrostatics=True,
    annihilate_sterics=False,
    precision=None,
):
    """
    Compare energies of reference system and fully-interacting alchemically modified system.

    ARGUMENTS

    reference_system (simtk.openmm.System) - the reference System object to compare with
    positions - the positions to assess energetics for
    receptor_atoms (list of int) - the list of receptor atoms
    ligand_atoms (list of int) - the list of ligand atoms to alchemically modify
    precision : str, optional, default=None
       Precision model, or default if not specified. ('single', 'double', 'mixed')

    """

    # Create a factory to produce alchemical intermediates.
    logger.info("Creating alchemical factory...")
    initial_time = time.time()
    factory = AbsoluteAlchemicalFactory(
        reference_system,
        ligand_atoms=ligand_atoms,
        annihilate_electrostatics=annihilate_electrostatics,
        annihilate_sterics=annihilate_sterics,
    )
    final_time = time.time()
    elapsed_time = final_time - initial_time
    logger.info("AbsoluteAlchemicalFactory initialization took %.3f s" % elapsed_time)
    platform = None
    if platform_name:
        platform = openmm.Platform.getPlatformByName(platform_name)
    alchemical_system = factory.createPerturbedSystem()
    compareSystemEnergies(
        positions,
        [reference_system, alchemical_system],
        ["reference", "alchemical"],
        platform=platform,
        precision=precision,
    )
    return

    def make_alchemical_system(self, topology_proposal, direction='insert'):
        """
        Generate an alchemically-modified system at the correct atoms
        based on the topology proposal

        Arguments
        ---------
        topology_proposal : TopologyProposal namedtuple
            Contains old topology, proposed new topology, and atom mapping
        direction : str, optional, default='insert'
            Direction of topology proposal to use for identifying alchemical atoms (allowed values: ['insert', 'delete'])

        Returns
        -------
        unmodified_system : simtk.openmm.System
            Unmodified real system corresponding to appropriate leg of transformation.
        alchemical_system : simtk.openmm.System
            The system with appropriate atoms alchemically modified

        """
        if direction not in ['insert', 'delete']:
            raise Exception("'direction' must be one of ['insert', 'delete']; was '%s' instead" % direction)

        atom_map = topology_proposal.new_to_old_atom_map

        #take the unique atoms as those not in the {new_atom : old_atom} atom map
        if direction == 'delete':
            unmodified_system = topology_proposal.old_system
            alchemical_atoms = [atom for atom in range(unmodified_system.getNumParticles()) if atom not in atom_map.values()]
        elif direction == 'insert':
            unmodified_system = topology_proposal.new_system
            alchemical_atoms = [atom for atom in range(unmodified_system.getNumParticles()) if atom not in atom_map.keys()]
        else:
            raise Exception("direction must be one of ['delete', 'insert']; found '%s' instead" % direction)

        # DEBUG
        #print('alchemical atoms:')
        #print(alchemical_atoms)

        # Create an alchemical factory.
        from alchemy import AbsoluteAlchemicalFactory
        alchemical_factory = AbsoluteAlchemicalFactory(unmodified_system, ligand_atoms=alchemical_atoms, annihilate_electrostatics=True, annihilate_sterics=True, alchemical_bonds=None, alchemical_angles=None)

        # Return the alchemically-modified system in fully-interacting form.
        alchemical_system = alchemical_factory.createPerturbedSystem()
        return [unmodified_system, alchemical_system]

def test_ncmc_alchemical_integrator_stability_molecules():
    """
    Test NCMCAlchemicalIntegrator

    """
    molecule_names = ['pentane', 'biphenyl', 'imatinib']
    #if os.environ.get("TRAVIS", None) == 'true':
    #    molecule_names = ['pentane']

    for molecule_name in molecule_names:
        from perses.tests.utils import createSystemFromIUPAC
        [molecule, system, positions, topology] = createSystemFromIUPAC(molecule_name)

        # Eliminate half of the molecule
        # TODO: Use a more rigorous scheme to make sure we are really cutting the molecule in half and not just eliminating hydrogens or something.
        alchemical_atoms = [ index for index in range(int(system.getNumParticles()/2)) ]

        # Create an alchemically-modified system.
        from alchemy import AbsoluteAlchemicalFactory
        alchemical_factory = AbsoluteAlchemicalFactory(system, ligand_atoms=alchemical_atoms, annihilate_electrostatics=True, annihilate_sterics=True)

        # Return the alchemically-modified system in fully-interacting form.
        alchemical_system = alchemical_factory.createPerturbedSystem()

        # Create an NCMC switching integrator.
        from perses.annihilation.ncmc_switching import NCMCVVAlchemicalIntegrator
        temperature = 300.0 * unit.kelvin
        nsteps = 10 # number of steps to run integration for
        functions = { 'lambda_sterics' : 'lambda', 'lambda_electrostatics' : 'lambda^0.5', 'lambda_torsions' : 'lambda', 'lambda_angles' : 'lambda^2' }
        ncmc_integrator = NCMCVVAlchemicalIntegrator(temperature, alchemical_system, functions, direction='delete', nsteps=nsteps, timestep=1.0*unit.femtoseconds)

        # Create a Context
        context = openmm.Context(alchemical_system, ncmc_integrator)
        context.setPositions(positions)

        # Run the integrator
        ncmc_integrator.step(nsteps)

        # Check positions are finite
        positions = context.getState(getPositions=True).getPositions(asNumpy=True)
        if np.isnan(np.any(positions / positions.unit)):
            raise Exception('NCMCAlchemicalIntegrator gave NaN positions')
        if np.isnan(ncmc_integrator.getLogAcceptanceProbability(context)):
            raise Exception('NCMCAlchemicalIntegrator gave NaN logAcceptanceProbability')

        del context, ncmc_integrator

def check_waterbox(platform=None, precision=None, nonbondedMethod=openmm.NonbondedForce.CutoffPeriodic):
    """Compare annihilated states in vacuum and a large box.
    """
    platform_name = platform.getName()
    from openmmtools import testsystems
    testsystem = testsystems.WaterBox()
    system = testsystem.system
    positions = testsystem.positions

    # Use reaction field
    for force in system.getForces():
        if force.__class__.__name__ == 'NonbondedForce':
            force.setNonbondedMethod(nonbondedMethod)

    factory_args = {'ligand_atoms' : [], 'receptor_atoms' : [],
        'annihilate_sterics' : False, 'annihilate_electrostatics' : True }

    # Create alchemically-modified system
    factory = AbsoluteAlchemicalFactory(system, **factory_args)
    alchemical_system = factory.createPerturbedSystem()

    # Compare energies
    system_energy = compute_energy(system, positions, platform=platform, precision=precision)
    alchemical_1_energy = compute_energy(alchemical_system, positions, platform=platform, precision=precision)

    # Set lambda = 0
    lambda_value = 0.0
    alchemical_state = AlchemicalState(lambda_electrostatics=lambda_value, lambda_sterics=lambda_value, lambda_torsions=lambda_value)
    AbsoluteAlchemicalFactory.perturbSystem(alchemical_system, alchemical_state)
    alchemical_0_energy = compute_energy(alchemical_system, positions, platform=platform, precision=precision)

    # Check deviation.
    logger.info("========")
    logger.info("Platform %s" % platform_name)
    logger.info("Alchemically-modified WaterBox with no alchemical atoms")
    logger.info('real system : %8.3f kcal/mol' % (system_energy / unit.kilocalories_per_mole))
    logger.info('lambda = 1  : %8.3f kcal/mol' % (alchemical_1_energy / unit.kilocalories_per_mole))
    logger.info('lambda = 0  : %8.3f kcal/mol' % (alchemical_0_energy / unit.kilocalories_per_mole))
    delta = alchemical_1_energy - alchemical_0_energy
    logger.info("ERROR       : %8.3f kcal/mol" % (delta / unit.kilocalories_per_mole))
    if (abs(delta) > MAX_DELTA):
        raise Exception("Maximum allowable deviation on platform %s exceeded (was %.8f kcal/mol; allowed %.8f kcal/mol); test failed." % (platform_name, delta / unit.kilocalories_per_mole, MAX_DELTA / unit.kilocalories_per_mole))

    def make_alchemical_system(self, unmodified_system, topology_proposal, direction='create'):
        """
        Generate an alchemically-modified system at the correct atoms
        based on the topology proposal

        Arguments
        ---------
        unmodified_system : openmm.System object
            The unmodified system to get alchemical modifications
        topology_proposal : TopologyProposal namedtuple
            Contains old topology, proposed new topology, and atom mapping
        direction : str, optional, default='insert'
            Direction of topology proposal to use for identifying alchemical atoms.

        Returns
        -------
        alchemical_system : openmm.System object
            The system with appropriate atoms alchemically modified

        """
        atom_map = topology_proposal.new_to_old_atom_map
        n_atoms = unmodified_system.getNumParticles()

        #take the unique atoms as those not in the {new_atom : old_atom} atom map
        if direction == 'delete':
            alchemical_atoms = [atom for atom in range(n_atoms) if atom not in atom_map.values()]
        elif direction == 'create':
            alchemical_atoms = [atom for atom in range(n_atoms) if atom not in atom_map.keys()]
        else:
            raise Exception("direction must be one of ['delete', 'create']; found '%s' instead" % direction)

        logging.debug(alchemical_atoms)
        # Create an alchemical factory.
        from alchemy import AbsoluteAlchemicalFactory
        alchemical_factory = AbsoluteAlchemicalFactory(unmodified_system, ligand_atoms=alchemical_atoms)

        # Return the alchemically-modified system.
        alchemical_system = alchemical_factory.createPerturbedSystem()
        return alchemical_system

    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

def testAlchemicalFactory(
    reference_system,
    coordinates,
    receptor_atoms,
    ligand_atoms,
    platform_name="CUDA",
    annihilateElectrostatics=True,
    annihilateLennardJones=False,
):
    """
    Compare energies of reference system and fully-interacting alchemically modified system.

    ARGUMENTS
    
    reference_system (simtk.openmm.System) - the reference System object to compare with
    coordinates - the coordinates to assess energetics for
    receptor_atoms (list of int) - the list of receptor atoms 
    ligand_atoms (list of int) - the list of ligand atoms to alchemically modify

    """

    import simtk.unit as units
    import simtk.openmm as openmm
    import time

    # Create a factory to produce alchemical intermediates.
    print "Creating alchemical factory..."
    initial_time = time.time()
    factory = AbsoluteAlchemicalFactory(reference_system, ligand_atoms=ligand_atoms)
    final_time = time.time()
    elapsed_time = final_time - initial_time
    print "AbsoluteAlchemicalFactory initialization took %.3f s" % elapsed_time

    platform = openmm.Platform.getPlatformByName(platform_name)

    delta = 0.001
    delta = 1.0e-6

    compareSystemEnergies(
        coordinates,
        [
            reference_system,
            factory.createPerturbedSystem(
                AlchemicalState(
                    0,
                    1 - delta,
                    1,
                    1,
                    annihilateElectrostatics=annihilateElectrostatics,
                    annihilateLennardJones=annihilateLennardJones,
                )
            ),
        ],
        ["reference", "partially discharged"],
        platform=platform,
    )
    compareSystemEnergies(
        coordinates,
        [
            factory.createPerturbedSystem(
                AlchemicalState(
                    0,
                    delta,
                    1,
                    1,
                    annihilateElectrostatics=annihilateElectrostatics,
                    annihilateLennardJones=annihilateLennardJones,
                )
            ),
            factory.createPerturbedSystem(
                AlchemicalState(
                    0,
                    0.0,
                    1,
                    1,
                    annihilateElectrostatics=annihilateElectrostatics,
                    annihilateLennardJones=annihilateLennardJones,
                )
            ),
        ],
        ["partially charged", "discharged"],
        platform=platform,
    )
    compareSystemEnergies(
        coordinates,
        [
            factory.createPerturbedSystem(
                AlchemicalState(
                    0,
                    0,
                    1,
                    1,
                    annihilateElectrostatics=annihilateElectrostatics,
                    annihilateLennardJones=annihilateLennardJones,
                )
            ),
            factory.createPerturbedSystem(
                AlchemicalState(
                    0,
                    0,
#.........这里部分代码省略.........

def setup_binding_amber(args):
    """
    Set up ligand binding free energy calculation using AMBER prmtop/inpcrd files.

    Parameters
    ----------
    args : dict
       Command-line arguments dict from docopt.

    Returns
    -------
    alchemical_phases : list of AlchemicalPhase
       Phases (thermodynamic legs) of the calculation.

    """
    verbose = args['--verbose']
    setup_directory = args['--setupdir']  # Directory where prmtop/inpcrd files are to be found
    system_parameters = {}  # parameters to pass to prmtop.createSystem

    # Implicit solvent
    if args['--gbsa']:
        system_parameters['implicitSolvent'] = getattr(app, args['--gbsa'])

    # Select nonbonded treatment
    if args['--nbmethod']:
        system_parameters['nonbondedMethod'] = getattr(app, args['--nbmethod'])

    # Constraints
    if args['--constraints']:
        system_parameters['constraints'] = getattr(app, args['--constraints'])

    # Cutoff
    if args['--cutoff']:
        system_parameters['nonbondedCutoff'] = process_unit_bearing_arg(args, '--cutoff', unit.nanometers)

    # Determine if this will be an explicit or implicit solvent simulation
    if ('nonbondedMethod' in system_parameters and
                system_parameters['nonbondedMethod'] != app.NoCutoff):
        phases_names = ['complex-explicit', 'solvent-explicit']
        protocols = [AbsoluteAlchemicalFactory.defaultComplexProtocolExplicit(),
                     AbsoluteAlchemicalFactory.defaultSolventProtocolExplicit()]
    else:
        phases_names = ['complex-implicit', 'solvent-implicit']
        protocols = [AbsoluteAlchemicalFactory.defaultComplexProtocolImplicit(),
                 AbsoluteAlchemicalFactory.defaultSolventProtocolImplicit()]

    # Prepare Yank arguments
    alchemical_phases = [None, None]
    setup_directory = os.path.join(setup_directory, '')  # add final slash character
    system_files_paths = [[setup_directory + 'complex.inpcrd', setup_directory + 'complex.prmtop'],
                          [setup_directory + 'solvent.inpcrd', setup_directory + 'solvent.prmtop']]
    for i, phase_name in enumerate(phases_names):
        positions_file_path = system_files_paths[i][0]
        topology_file_path = system_files_paths[i][1]

        logger.info("Reading phase {}".format(phase_name))
        alchemical_phases[i] = pipeline.prepare_phase(positions_file_path, topology_file_path, args['--ligand'],
                                                      system_parameters, verbose=verbose)
        alchemical_phases[i].name = phase_name
        alchemical_phases[i].protocol = protocols[i]

    return alchemical_phases

 def _AutoAlchemyStates(self, phase, real_R_states=None, real_A_states=None, real_E_states=None, real_C_states=None, alchemy_source=None):
     #Generate the real alchemical states automatically.
     if alchemy_source: #Load alchemy from an external source
         import imp
         if alchemy_source[-3:] != '.py': #Check if the file or the folder was provided
             alchemy_source = os.path.join(alchemy_source, 'alchemy.py')
         alchemy = imp.load_source('alchemy', alchemy_source)
         AAF = alchemy.AbsoluteAlchemicalFactory
     else: #Standard load
         from alchemy import AbsoluteAlchemicalFactory as AAF
     if phase is 'vacuum':
         protocol = AAF.defaultVacuumProtocol()
     elif phase is 'complex':
         protocol = AAF.defaultComplexProtocolExplicit()
     #Determine which phases need crunched
     if real_R_states is None:
         real_R_states = list()
         crunchR = True
     else:
         crunchR = False
     if real_A_states is None:
         real_A_states = list()
         crunchA = True
     else:
         crunchA = False
     if real_E_states is None:
         real_E_states = list()
         real_PMEFull_states = list()
         crunchE = True
     else:
         crunchE = False
     #Detect for the cap basis property
     if numpy.all([hasattr(state, 'ligandCapToFull') for state in protocol]) and real_C_states is None:
         real_C_states = list()
         crunchC = True
     else:
         crunchC = False
     #Import from the alchemy file if need be
     for state in protocol: #Go through each state
         if crunchE:
             real_E_states.append(state.ligandElectrostatics)
             try:
                 real_PMEFull_states.append(state.ligandPMEFull)
             except:
                 real_PMEFull_states.append(None)
         if crunchR:
             real_R_states.append(state.ligandRepulsion)
         if crunchA:
             real_A_states.append(state.ligandAttraction)
         if crunchC:
             real_C_states.append(state.ligandCapToFull)
     if numpy.all([i is None for i in real_PMEFull_states]): #Must put [...] around otherwise it creates the generator object which numpy.all evals to True
         self.PME_isolated = False
     else:
         self.PME_isolated = True
     #Determine cutoffs
     self.real_E_states = numpy.array(real_E_states)
     self.real_PMEFull_states = numpy.array(real_PMEFull_states)
     self.real_R_states = numpy.array(real_R_states)
     self.real_A_states = numpy.array(real_A_states)
     self.real_C_states = numpy.array(real_C_states)
     indicies = numpy.array(range(len(real_E_states)))
     #Determine Inversion
     if numpy.any(self.real_E_states < 0) or numpy.any(numpy.logical_and(self.real_PMEFull_states < 0,numpy.array([i is not None for i in self.real_PMEFull_states]))):
         self.Inversion = True
     else:
         self.Inversion = False
     #Set the indicies, trap TypeError (logical_and false everywhere) as None (i.e. state not found in alchemy)
     if crunchC: #Check for the cap potential
         print "Not Coded Yet!"
         exit(1)
         #Create the Combinations
         basisVars = ["E", "A", "R", "C"]
         mappedStates = [self.real_E_states, self.real_R_states, self.real_C_states, self.real_A_states]
         nBasis = len(basisVars)
         coupled_states = {}
         decoupled_states = {}
         for iBasis in xrange(nBasis):
             coupled_states[basisVars[iBasis]] = numpy.where(mappedStates[iBasis] == 1.00)[0] #need the [0] to extract the array from the basis
             decoupled_states[basisVars[iBasis]] = numpy.where(mappedStates[iBasis] == 0.00)[0]
         self.coupled_states = coupled_states
         self.decoupled_states = decoupled_states
         self.basisVars = basisVars
     else:
         if self.PME_isolated: #Logic to solve for isolated PME case
             try: #Figure out the Fully coupled state
                 self.real_EAR = int(indicies[ numpy.logical_and(numpy.logical_and(self.real_E_states == 1, self.real_PMEFull_states == 1), numpy.logical_and(self.real_R_states == 1, self.real_A_states == 1)) ])
             except TypeError:
                 self.real_EAR = None
             try:
                 self.real_AR = int(indicies[ numpy.logical_and(numpy.logical_and(self.real_E_states == 0, self.real_PMEFull_states == 0), numpy.logical_and(self.real_R_states == 1, self.real_A_states == 1)) ])
             except TypeError:
                 self.real_AR = None
             try:
                 self.real_R = int(indicies[ numpy.logical_and(numpy.logical_and(self.real_E_states == 0, self.real_PMEFull_states == 0), numpy.logical_and(self.real_R_states == 1, self.real_A_states == 0)) ])
             except TypeError:
                 self.real_R = None
             try:
                 self.real_alloff = int(indicies[ numpy.logical_and(numpy.logical_and(self.real_E_states == 0, self.real_PMEFull_states == 0), numpy.logical_and(self.real_R_states == 0, self.real_A_states == 0)) ])
             except:
#.........这里部分代码省略.........

    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.
#.........这里部分代码省略.........