Python AbsoluteAlchemicalFactory.defaultComplexProtocolExplicit方法代码示例

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

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import defaultComplexProtocolExplicit [as 别名]
def setup_binding_gromacs(args):
    """
    Set up ligand binding free energy calculation using gromacs 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']

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

    # Select nonbonded treatment
    # TODO: Carefully check whether input file is periodic or not.
    if args['--nbmethod']:
        nonbondedMethod = getattr(app, args['--nbmethod'])
    else:
        nonbondedMethod = None

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

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

    # COM removal
    removeCMMotion = False

    # Prepare phases of calculation.
    phase_prefixes = ['solvent', 'complex'] # list of calculation phases (thermodynamic legs) to set up
    components = ['ligand', 'receptor', 'solvent'] # components of the binding system
    systems = dict() # systems[phase] is the System object associated with phase 'phase'
    topologies = dict() # topologies[phase] is the Topology object associated with phase 'phase'
    positions = dict() # positions[phase] is a list of coordinates associated with phase 'phase'
    atom_indices = dict() # ligand_atoms[phase] is a list of ligand atom indices associated with phase 'phase'
    setup_directory = args['--setupdir'] # Directory where prmtop/inpcrd files are to be found
    for phase_prefix in phase_prefixes:
        if verbose: logger.info("reading phase %s: " % phase_prefix)
        # Read gromacs input files.
        gro_filename = os.path.join(setup_directory, '%s.gro' % phase_prefix)
        top_filename = os.path.join(setup_directory, '%s.top' % phase_prefix)
        if verbose: logger.info('reading gromacs .gro file: %s' % gro_filename)
        gro = app.GromacsGroFile(gro_filename)
        if verbose: logger.info('reading gromacs .top file "%s" using gromacs include directory "%s"' % (top_filename, args['--gromacsinclude']))
        top = app.GromacsTopFile(top_filename, unitCellDimensions=gro.getUnitCellDimensions(), includeDir=args['--gromacsinclude'])
        # Assume explicit solvent.
        # TODO: Modify this if we can have implicit solvent.
        is_periodic = True
        phase_suffix = 'explicit'
        # Adjust nonbondedMethod.
        # TODO: Ensure that selected method is appropriate.
        if nonbondedMethod == None:
            if is_periodic:
                nonbondedMethod = app.CutoffPeriodic
            else:
                nonbondedMethod = app.NoCutoff
        # TODO: Check to make sure both prmtop and inpcrd agree on explicit/implicit.
        phase = '%s-%s' % (phase_prefix, phase_suffix)
        systems[phase] = top.createSystem(nonbondedMethod=nonbondedMethod, nonbondedCutoff=nonbondedCutoff, constraints=constraints, removeCMMotion=removeCMMotion)
        topologies[phase] = top.topology
        positions[phase] = gro.getPositions(asNumpy=True)
        # Check to make sure number of atoms match between prmtop and inpcrd.
        prmtop_natoms = systems[phase].getNumParticles()
        inpcrd_natoms = positions[phase].shape[0]
        if prmtop_natoms != inpcrd_natoms:
            raise Exception("Atom number mismatch: prmtop %s has %d atoms; inpcrd %s has %d atoms." % (prmtop_filename, prmtop_natoms, inpcrd_filename, inpcrd_natoms))

        # Find ligand atoms and receptor atoms.
        ligand_dsl = args['--ligand'] # MDTraj DSL that specifies ligand atoms
        atom_indices[phase] = find_components(systems[phase], top.topology, ligand_dsl)

    phases = systems.keys()

    alchemical_phases = [None, None]
    protocols = {'complex-explicit': AbsoluteAlchemicalFactory.defaultComplexProtocolExplicit(),
                 'solvent-explicit': AbsoluteAlchemicalFactory.defaultSolventProtocolImplicit()}
    for i, name in enumerate(phases):
        alchemical_phases[i] = AlchemicalPhase(name, systems[name], topologies[name],
                                               positions[name], atom_indices[name],
                                               protocols[name])

    return alchemical_phases

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

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

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import defaultComplexProtocolExplicit [as 别名]
    def __init__(self, store_directory, mpicomm=None, **kwargs):
        """
        Initialize YANK object with default parameters.

        Parameters
        ----------
        store_directory : str
           The storage directory in which output NetCDF files are read or written.
        mpicomm : MPI communicator, optional
           If an MPI communicator is passed, an MPI simulation will be attempted.
        restraint_type : str, optional
           Restraint type to add between protein and ligand. Supported types are
           'flat-bottom' and 'harmonic'. The second one is available only in
           implicit solvent (default: 'flat-bottom').
        randomize_ligand : bool, optional
           Randomize ligand position when True. Not available in explicit solvent
           (default: False).
        randomize_ligand_close_cutoff : simtk.unit.Quantity (units: length), optional
           Cutoff for ligand position randomization (default: 1.5*unit.angstrom).
        randomize_ligand_sigma_multiplier : float, optional
           Multiplier for ligand position randomization displacement (default: 2.0).
        mc_displacement_sigma : simtk.unit.Quantity (units: length), optional
           Maximum displacement for Monte Carlo moves that augment Langevin dynamics
           (default: 10.0*unit.angstrom).

        Other Parameters
        ----------------
        **kwargs
           More options to pass to the ReplicaExchange or AlchemicalFactory classes
           on initialization.

        See Also
        --------
        ReplicaExchange.default_parameters : extra parameters accepted.

        """

        # Copy kwargs to avoid modifications
        parameters = copy.deepcopy(kwargs)

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

        # Store output directory.
        self._store_directory = store_directory

        # Save MPI communicator
        self._mpicomm = mpicomm

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

        # Store Yank parameters
        for option_name, default_value in self.default_parameters.items():
            setattr(self, '_' + option_name, parameters.pop(option_name, default_value))

        # Store repex parameters
        self._repex_parameters = {par: parameters.pop(par) for par in
                                  ModifiedHamiltonianExchange.default_parameters
                                  if par in parameters}

        # Store AlchemicalFactory parameters
        self._alchemy_parameters = {par: parameters.pop(par) for par in
                                    inspect.getargspec(AbsoluteAlchemicalFactory.__init__).args
                                    if par in parameters}

        # Check for unknown parameters
        if len(parameters) > 0:
            raise TypeError('got an unexpected keyword arguments {}'.format(
                ', '.join(parameters.keys())))

# 需要导入模块: from alchemy import AbsoluteAlchemicalFactory [as 别名]
# 或者: from alchemy.AbsoluteAlchemicalFactory import defaultComplexProtocolExplicit [as 别名]
options['collision_rate'] = collision_rate
options['platform'] = platform

options['restraint_type'] = None
if not is_periodic:
    options['restraint_type'] = 'harmonic'

# Turn off MC ligand displacement.
options['mc_displacement_sigma'] = None

# Prepare phases of calculation.
phase_prefixes = ['solvent', 'complex'] # list of calculation phases (thermodynamic legs) to set up
components = ['ligand', 'receptor', 'solvent'] # components of the binding system
phase_prefixes = ['complex'] # DEBUG, since 'solvent' doesn't work yet
if is_periodic:
    protocols = {'complex': AbsoluteAlchemicalFactory.defaultComplexProtocolExplicit(),
                 'solvent': AbsoluteAlchemicalFactory.defaultSolventProtocolExplicit()}
else:
    protocols = {'complex': AbsoluteAlchemicalFactory.defaultComplexProtocolImplicit(),
                 'solvent': AbsoluteAlchemicalFactory.defaultSolventProtocolImplicit()}
alchemical_phases = []  # alchemical phases of the calculations
for phase_prefix in phase_prefixes:
    # Retain the whole system
    if is_periodic:
        phase_suffix = 'explicit'
    else:
        phase_suffix = 'implicit'

    # Form phase name.
    phase = '%s-%s' % (phase_prefix, phase_suffix)
    logger.info("phase %s: " % phase)