Python AbsoluteAlchemicalFactory.perturbSystem方法代码示例

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import perturbSystem [as 别名]
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))

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import perturbSystem [as 别名]
"""
Create alchemical intermediates for default alchemical protocol for one water in a water box.

"""
from alchemy import AbsoluteAlchemicalFactory, AlchemicalState
from openmmtools import testsystems

# Create a reference system.
print "Creating a water box..."
waterbox = testsystems.WaterBox()
[reference_system, positions] = [waterbox.system, waterbox.positions]

# Create a factory to produce alchemical intermediates.
print "Creating an alchemical factory..."
factory = AbsoluteAlchemicalFactory(reference_system, ligand_atoms=[0, 1, 2])

# Create a perturbed systems using this protocol.
print "Creating a perturbed system..."
alchemical_state = AlchemicalState()
alchemical_system = factory.createPerturbedSystem(alchemical_state)

# Perturb this system.
print "Perturbing the system..."
alchemical_state = AlchemicalState(lambda_sterics=0.90, lambda_electrostatics=0.90)
factory.perturbSystem(alchemical_system, alchemical_state)

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

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

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import perturbSystem [as 别名]
def test_annihilated_states(platform_name=None, precision=None):
    """Compare annihilated states in vacuum and a large box.
    """
    from openmmtools import testsystems
    testsystem = testsystems.TolueneVacuum()
    vacuum_system = testsystem.system
    positions = testsystem.positions

    factory_args = {'ligand_atoms' : range(0,15), 'receptor_atoms' : [],
        'annihilate_sterics' : False, 'annihilate_electrostatics' : True }

    # Create annihilated version of vacuum system.
    factory = AbsoluteAlchemicalFactory(vacuum_system, **factory_args)
    vacuum_alchemical_system = factory.createPerturbedSystem()

    # Make copy of system that has periodic boundaries and uses reaction field.
    periodic_system = copy.deepcopy(vacuum_system)
    box_edge = 18.5 * unit.angstroms
    from simtk.openmm import Vec3
    periodic_system.setDefaultPeriodicBoxVectors(Vec3(box_edge,0,0), Vec3(0,box_edge,0), Vec3(0,0,box_edge))
    for force in periodic_system.getForces():
        if force.__class__.__name__ == 'NonbondedForce':
            force.setNonbondedMethod(openmm.NonbondedForce.PME)
            force.setCutoffDistance(9.0 * unit.angstroms)
            force.setUseDispersionCorrection(False)
            force.setReactionFieldDielectric(1.0)
    factory = AbsoluteAlchemicalFactory(periodic_system, **factory_args)
    periodic_alchemical_system = factory.createPerturbedSystem()

    # Compare energies
    platform = None
    if platform_name:
        platform = openmm.Platform.getPlatformByName(platform_name)

    vacuum_alchemical_1_energy = compute_energy(vacuum_alchemical_system, positions, platform=platform, precision=precision)
    periodic_alchemical_1_energy = compute_energy(periodic_alchemical_system, positions, platform=platform, precision=precision)

    #compareSystemEnergies(positions, [vacuum_alchemical_system, periodic_alchemical_system], ['vacuum (fully interacting)', 'periodic (fully interacting)'], 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(vacuum_alchemical_system, alchemical_state)
    AbsoluteAlchemicalFactory.perturbSystem(periodic_alchemical_system, alchemical_state)

    #compareSystemEnergies(positions, [vacuum_alchemical_system, periodic_alchemical_system], ['vacuum (noninteracting)', 'periodic (noninteracting)'], platform=platform, precision=precision)

    vacuum_alchemical_0_energy = compute_energy(vacuum_alchemical_system, positions, platform=platform, precision=precision)
    periodic_alchemical_0_energy = compute_energy(periodic_alchemical_system, positions, platform=platform, precision=precision)

    logger.info('vacuum   lambda = 1 : %8.3f kcal/mol' % (vacuum_alchemical_1_energy / unit.kilocalories_per_mole))
    logger.info('vacuum   lambda = 0 : %8.3f kcal/mol' % (vacuum_alchemical_0_energy / unit.kilocalories_per_mole))
    logger.info('difference          : %8.3f kcal/mol' % ((vacuum_alchemical_1_energy - vacuum_alchemical_0_energy) / unit.kilocalories_per_mole))

    logger.info('periodic lambda = 1 : %8.3f kcal/mol' % (periodic_alchemical_1_energy / unit.kilocalories_per_mole))
    logger.info('periodic lambda = 0 : %8.3f kcal/mol' % (periodic_alchemical_0_energy / unit.kilocalories_per_mole))
    logger.info('difference          : %8.3f kcal/mol' % ((periodic_alchemical_1_energy - periodic_alchemical_0_energy) / unit.kilocalories_per_mole))

    delta = (vacuum_alchemical_1_energy - vacuum_alchemical_0_energy) - (periodic_alchemical_1_energy - periodic_alchemical_0_energy)
    if (abs(delta) > MAX_DELTA):
        raise Exception("Maximum allowable difference lambda=1 energy and lambda=0 energy in vacuum and periodic box exceeded (was %.8f kcal/mol; allowed %.8f kcal/mol); test failed." % (delta / unit.kilocalories_per_mole, MAX_DELTA / unit.kilocalories_per_mole))

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import perturbSystem [as 别名]
 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)