Python AbsoluteAlchemicalFactory.createPerturbedSystem方法代码示例

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

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

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

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

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

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

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

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import createPerturbedSystem [as 别名]
    def __init__(self, **kwargs):
        super(LoopSoftening, self).__init__(**kwargs)
        self.description = 'Alchemical Loop Softening script'

        padding = 9.0*unit.angstrom
        explicit_solvent_model = 'tip3p'
        setup_path = 'data/mtor'

        # Create topology, positions, and system.
        from pkg_resources import resource_filename
        gaff_xml_filename = resource_filename('sams', 'data/gaff.xml')
        system_generators = dict()
        ffxmls = [gaff_xml_filename, 'amber99sbildn.xml', 'tip3p.xml']
        forcefield_kwargs={ 'nonbondedMethod' : app.CutoffPeriodic, 'nonbondedCutoff' : 9.0 * unit.angstrom, 'implicitSolvent' : None, 'constraints' : app.HBonds, 'rigidWater' : True }

        # Load topologies and positions for all components
        print('Creating mTOR test system...')
        forcefield = app.ForceField(*ffxmls)
        from simtk.openmm.app import PDBFile, Modeller
        pdb_filename = resource_filename('sams', os.path.join(setup_path, 'mtor_pdbfixer_apo.pdb'))

        pdbfile = PDBFile(pdb_filename)
        modeller = app.Modeller(pdbfile.topology, pdbfile.positions)
        print('Adding solvent...')
        modeller.addSolvent(forcefield, model=explicit_solvent_model, padding=padding)
        self.topology = modeller.getTopology()
        self.positions = modeller.getPositions()
        print('Creating system...')
        self.system = forcefield.createSystem(self.topology, **forcefield_kwargs)

        # DEBUG: Write PDB
        outfile = open('initial.pdb', 'w')
        PDBFile.writeFile(self.topology, self.positions, outfile)
        outfile.close()

        # Atom Selection using MDtraj
        res_pairs = [[403, 483], [1052, 1109]]
        t = md.load(pdb_filename)
        alchemical_atoms = set()
        for x in res_pairs:
            start = min(t.top.select('residue %s' % min(x)))
            end = max(t.top.select('residue %s' % max(x))) + 1
            alchemical_atoms.union(set(range(start, end)))


        # Create thermodynamic states.
        print('Creating alchemically-modified system...')
        temperature = 300 * unit.kelvin
        pressure = 1.0 * unit.atmospheres

        from alchemy import AbsoluteAlchemicalFactory
        factory = AbsoluteAlchemicalFactory(self.system, ligand_atoms=alchemical_atoms, annihilate_electrostatics=True,
                                            alchemical_torsions=True, annihilate_sterics=True,
                                            softcore_beta=0.0)  # turn off softcore electrostatics
        self.system = factory.createPerturbedSystem()
        print('Setting up alchemical intermediates...')
        from sams import ThermodynamicState
        self.thermodynamic_states = list()
        for state in range(26):
            parameters = {'lambda_sterics' : 1.0, 'lambda_electrostatics' : (1.0 - float(state)/25.0) }
            self.thermodynamic_states.append( ThermodynamicState(system=self.system, temperature=temperature,
                                                                 parameters=parameters) )
        for state in range(1,26):
            parameters = {'lambda_sterics' : (1.0 - float(state)/25.0), 'lambda_electrostatics' : 0.0 }
            self.thermodynamic_states.append( ThermodynamicState(system=self.system, temperature=temperature,
                                                                 parameters=parameters) )

        #minimize(self.system, self.positions)
        minimize(self.system)

        # Create SAMS samplers
        print('Setting up samplers...')
        from sams.samplers import SamplerState, MCMCSampler, ExpandedEnsembleSampler, SAMSSampler
        thermodynamic_state_index = 0  # initial thermodynamic state index
        thermodynamic_state = self.thermodynamic_states[thermodynamic_state_index]
        sampler_state = SamplerState(positions=self.system.positions)
        self.mcmc_sampler = MCMCSampler(sampler_state=sampler_state, thermodynamic_state=thermodynamic_state,
                                        ncfile=self.ncfile)
        self.mcmc_sampler.pdbfile = open('output.pdb', 'w')
        self.mcmc_sampler.topology = self.topology
        self.mcmc_sampler.verbose = True
        self.exen_sampler = ExpandedEnsembleSampler(self.mcmc_sampler, self.thermodynamic_states)
        self.exen_sampler.verbose = True
        self.sams_sampler = SAMSSampler(self.exen_sampler)
        self.sams_sampler.verbose = True

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import createPerturbedSystem [as 别名]
def testAlchemicalFactory(reference_system, coordinates, receptor_atoms, ligand_atoms, platform_name='Reference', 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

    # Create an alchemically-perturbed state corresponding to nearly fully-interacting.
    # NOTE: We use a lambda slightly smaller than 1.0 because the AlchemicalFactory does not use Custom*Force softcore versions if lambda = 1.0 identically.
    lambda_value = 1.0 - 1.0e-6
    alchemical_state = AlchemicalState(0.00, lambda_value, lambda_value, lambda_value)
    alchemical_state.annihilateElectrostatics = annihilateElectrostatics
    alchemical_state.annihilateLennardJones = annihilateLennardJones

    #platform_name = 'Reference' # DEBUG
    platform = openmm.Platform.getPlatformByName(platform_name)
    
    # Create the perturbed system.
    print "Creating alchemically-modified state..."
    initial_time = time.time()
    alchemical_system = factory.createPerturbedSystem(alchemical_state)    
    final_time = time.time()
    elapsed_time = final_time - initial_time
    # Compare energies.
    timestep = 1.0 * units.femtosecond
    print "Computing reference energies..."
    reference_integrator = openmm.VerletIntegrator(timestep)
    reference_context = openmm.Context(reference_system, reference_integrator, platform)
    reference_context.setPositions(coordinates)
    reference_state = reference_context.getState(getEnergy=True)
    reference_potential = reference_state.getPotentialEnergy()    
    print "Computing alchemical energies..."
    alchemical_integrator = openmm.VerletIntegrator(timestep)
    alchemical_context = openmm.Context(alchemical_system, alchemical_integrator, platform)
    alchemical_context.setPositions(coordinates)
    alchemical_state = alchemical_context.getState(getEnergy=True)
    alchemical_potential = alchemical_state.getPotentialEnergy()
    delta = alchemical_potential - reference_potential 
    print "reference system       : %24.8f kcal/mol" % (reference_potential / units.kilocalories_per_mole)
    print "alchemically modified  : %24.8f kcal/mol" % (alchemical_potential / units.kilocalories_per_mole)
    print "ERROR                  : %24.8f kcal/mol" % ((alchemical_potential - reference_potential) / units.kilocalories_per_mole)
    print "elapsed alchemical time  %.3f s" % elapsed_time

    return delta

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import createPerturbedSystem [as 别名]
def test_overlap():
    """
    BUGS TO REPORT:
    * Even if epsilon = 0, energy of two overlapping atoms is 'nan'.
    * Periodicity in 'nan' if dr = 0.1 even in nonperiodic system
    """

    # Create a reference system.    
    import testsystems

    print "Creating Lennard-Jones cluster system..."
    #[reference_system, coordinates] = testsystems.LennardJonesFluid()
    #receptor_atoms = [0]
    #ligand_atoms = [1]

    [reference_system, coordinates] = testsystems.LysozymeImplicit()
    receptor_atoms = range(0,2603) # T4 lysozyme L99A
    ligand_atoms = range(2603,2621) # p-xylene

    import simtk.unit as units
    unit = coordinates.unit
    coordinates = units.Quantity(numpy.array(coordinates / unit), unit)

    factory = AbsoluteAlchemicalFactory(reference_system, ligand_atoms=ligand_atoms)
    alchemical_state = AlchemicalState(0.00, 0.00, 0.00, 1.0)

    # Create the perturbed system.
    print "Creating alchemically-modified state..."
    alchemical_system = factory.createPerturbedSystem(alchemical_state)    
    # Compare energies.
    import simtk.unit as units
    import simtk.openmm as openmm
    timestep = 1.0 * units.femtosecond
    print "Computing reference energies..."
    integrator = openmm.VerletIntegrator(timestep)
    context = openmm.Context(reference_system, integrator)
    context.setPositions(coordinates)
    state = context.getState(getEnergy=True)
    reference_potential = state.getPotentialEnergy()    
    del state, context, integrator
    print reference_potential
    print "Computing alchemical energies..."
    integrator = openmm.VerletIntegrator(timestep)
    context = openmm.Context(alchemical_system, integrator)
    dr = 0.1 * units.angstroms # TODO: Why does 0.1 cause periodic 'nan's?
    a = receptor_atoms[-1]
    b = ligand_atoms[-1]
    delta = coordinates[a,:] - coordinates[b,:]
    for k in range(3):
        coordinates[ligand_atoms,k] += delta[k]
    for i in range(30):
        r = dr * i
        coordinates[ligand_atoms,0] += dr
          
        context.setPositions(coordinates)
        state = context.getState(getEnergy=True)
        alchemical_potential = state.getPotentialEnergy()    
        print "%8.3f A : %f " % (r / units.angstroms, alchemical_potential / units.kilocalories_per_mole)
    del state, context, integrator

    return

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import createPerturbedSystem [as 别名]
    reference_system.addForce(barostat)

# Identify ligand indices by residue name
print('Identifying ligand atoms to be alchemically modified...')
reference = md.load(pdb_filename)
alchemical_atoms = reference.topology.select(ligand_dsl_selection) # these atoms will be alchemically softened
alchemical_atoms = [ int(index) for index in alchemical_atoms ] # recode as Python int
print("MDTraj DSL selection '%s' identified %d atoms" % (ligand_dsl_selection, len(alchemical_atoms)))

# Create alchemically-modified system using fused softcore electrostatics and sterics
print('Creating alchemically modified system...')
print('lambda schedule: %s' % str(alchemical_lambdas))
from alchemy import AbsoluteAlchemicalFactory
from sams import ThermodynamicState
factory = AbsoluteAlchemicalFactory(reference_system, ligand_atoms=alchemical_atoms, annihilate_electrostatics=True, annihilate_sterics=False)
system = factory.createPerturbedSystem()

# Add umbrella restraint with global variable to control umbrella position
print('umbrella schedule between atoms %d and %d: %s' % (umbrella_atoms[0], umbrella_atoms[1], str(umbrella_distances)))
energy_function = '(umbrella_K/2.0)*(r-umbrella_r0)^2'
umbrella_force = openmm.CustomBondForce(energy_function)
umbrella_force.addGlobalParameter('umbrella_K', 0.0) # spring constant
umbrella_force.addGlobalParameter('umbrella_r0', 0.0) # umbrella distance
umbrella_force.addBond(umbrella_atoms[0], umbrella_atoms[1], [])
umbrella_K = kT/umbrella_sigma**2
system.addForce(umbrella_force)

# Create thermodynamic states
thermodynamic_states = list()
for alchemical_lambda in alchemical_lambdas:
    # Umbrella off state

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import createPerturbedSystem [as 别名]
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,
#.........这里部分代码省略.........

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import createPerturbedSystem [as 别名]
def benchmark(
    reference_system,
    positions,
    receptor_atoms,
    ligand_atoms,
    platform_name=None,
    annihilate_electrostatics=True,
    annihilate_sterics=False,
    nsteps=500,
    timestep=1.0 * unit.femtoseconds,
):
    """
    Benchmark performance of alchemically modified system relative to original system.

    Parameters
    ----------
    reference_system : simtk.openmm.System
       The reference System object to compare with
    positions : simtk.unit.Quantity with units compatible with nanometers
       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.
    platform_name : str, optional, default=None
       The name of the platform to use for benchmarking.
    annihilate_electrostatics : bool, optional, default=True
       If True, electrostatics will be annihilated; if False, decoupled.
    annihilate_sterics : bool, optional, default=False
       If True, sterics will be annihilated; if False, decoupled.
    nsteps : int, optional, default=500
       Number of molecular dynamics steps to use for benchmarking.
    timestep : simtk.unit.Quantity with units compatible with femtoseconds, optional, default=1*femtoseconds
       Timestep to use for benchmarking.

    """

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

    # Create an alchemically-perturbed state corresponding to nearly fully-interacting.
    # NOTE: We use a lambda slightly smaller than 1.0 because the AlchemicalFactory does not use Custom*Force softcore versions if lambda = 1.0 identically.
    lambda_value = 1.0 - 1.0e-6
    alchemical_state = AlchemicalState(
        lambda_coulomb=lambda_value, lambda_sterics=lambda_value, lambda_torsions=lambda_value
    )

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

    # Create the perturbed system.
    logger.info("Creating alchemically-modified state...")
    initial_time = time.time()
    alchemical_system = factory.createPerturbedSystem(alchemical_state)
    final_time = time.time()
    elapsed_time = final_time - initial_time
    # Compare energies.
    logger.info("Computing reference energies...")
    reference_integrator = openmm.VerletIntegrator(timestep)
    if platform:
        reference_context = openmm.Context(reference_system, reference_integrator, platform)
    else:
        reference_context = openmm.Context(reference_system, reference_integrator)
    reference_context.setPositions(positions)
    reference_state = reference_context.getState(getEnergy=True)
    reference_potential = reference_state.getPotentialEnergy()
    logger.info("Computing alchemical energies...")
    alchemical_integrator = openmm.VerletIntegrator(timestep)
    if platform:
        alchemical_context = openmm.Context(alchemical_system, alchemical_integrator, platform)
    else:
        alchemical_context = openmm.Context(alchemical_system, alchemical_integrator)
    alchemical_context.setPositions(positions)
    alchemical_state = alchemical_context.getState(getEnergy=True)
    alchemical_potential = alchemical_state.getPotentialEnergy()
    delta = alchemical_potential - reference_potential

    # Make sure all kernels are compiled.
    reference_integrator.step(1)
    alchemical_integrator.step(1)

    # Time simulations.
    logger.info("Simulating reference system...")
    initial_time = time.time()
    reference_integrator.step(nsteps)
    reference_state = reference_context.getState(getEnergy=True)
    reference_potential = reference_state.getPotentialEnergy()
    final_time = time.time()
    reference_time = final_time - initial_time
#.........这里部分代码省略.........

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import createPerturbedSystem [as 别名]
def lambda_trace(reference_system, positions, platform_name=None, precision=None, nsteps=100, factory_args=None):
    """
    Compute potential energy as a function of lambda.

    """

    # Create a factory to produce alchemical intermediates.
    factory = AbsoluteAlchemicalFactory(reference_system, **factory_args)

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

    if precision:
        if platform_name == 'CUDA':
            platform.setDefaultPropertyValue('CudaPrecision', precision)
        elif platform_name == 'OpenCL':
            platform.setDefaultPropertyValue('OpenCLPrecision', precision)

    # Take equally-sized steps.
    delta = 1.0 / nsteps

    def compute_potential(system, positions, platform=None):
        timestep = 1.0 * unit.femtoseconds
        integrator = openmm.VerletIntegrator(timestep)
        if platform:
            context = openmm.Context(system, integrator, platform)
        else:
            context = openmm.Context(system, integrator)
        context.setPositions(positions)
        state = context.getState(getEnergy=True)
        potential = state.getPotentialEnergy()
        del integrator, context
        return potential

    # Compute unmodified energy.
    u_original = compute_potential(reference_system, positions, platform)

    # Scan through lambda values.
    lambda_i = np.zeros([nsteps+1], np.float64) # lambda values for u_i
    u_i = unit.Quantity(np.zeros([nsteps+1], np.float64), unit.kilocalories_per_mole) # u_i[i] is the potential energy for lambda_i[i]
    for i in range(nsteps+1):
        lambda_value = 1.0-i*delta # compute lambda value for this step
        alchemical_system = factory.createPerturbedSystem(AlchemicalState(lambda_electrostatics=lambda_value, lambda_sterics=lambda_value, lambda_torsions=lambda_value))
        lambda_i[i] = lambda_value
        u_i[i] = compute_potential(alchemical_system, positions, platform)
        logger.info("%12.9f %24.8f kcal/mol" % (lambda_i[i], u_i[i] / unit.kilocalories_per_mole))

    # Write figure as PDF.
    import pylab
    from matplotlib.backends.backend_pdf import PdfPages
    import matplotlib.pyplot as plt
    with PdfPages('lambda-trace.pdf') as pdf:
        fig = plt.figure(figsize=(10, 5))
        ax = fig.add_subplot(111)
        plt.plot(1, u_original / unit.kilocalories_per_mole, 'ro', label='unmodified')
        plt.plot(lambda_i, u_i / unit.kilocalories_per_mole, 'k.', label='alchemical')
        plt.title('T4 lysozyme L99A + p-xylene : AMBER96 + OBC GBSA')
        plt.ylabel('potential (kcal/mol)')
        plt.xlabel('lambda')
        ax.legend()
        rstyle(ax)
        pdf.savefig()  # saves the current figure into a pdf page
        plt.close()

    return