Python Molecule.write方法代码示例

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import write [as 别名]
def _writeInputsFunction(i, f, epoch, inputpath, coorname):
    regex = re.compile('(e\d+s\d+)_')
    frameNum = f.frame
    piece = f.piece
    if f.sim.parent is None:
        currSim = f.sim
    else:
        currSim = f.sim.parent

    traj = currSim.trajectory[piece]
    if currSim.input is None:
        raise NameError('Could not find input folder in simulation lists. Cannot create new simulations.')

    wuName = _simName(traj)
    res = regex.search(wuName)
    if res:  # If we are running on top of adaptive, use the first name part for the next sim name
        wuName = res.group(1)

    # create new job directory
    newName = 'e' + str(epoch) + 's' + str(i + 1) + '_' + wuName + 'p' + str(piece) + 'f' + str(frameNum)
    newDir = path.join(inputpath, newName, '')

    # copy previous input directory including input files
    copytree(currSim.input, newDir, symlinks=False, ignore=ignore_patterns('*.coor', '*.rst', '*.out', *_IGNORE_EXTENSIONS))

    # overwrite input file with new one. frameNum + 1 as catdcd does 1 based indexing

    mol = Molecule(currSim.molfile)  # Always read the mol file, otherwise it does not work if we need to save a PDB as coorname
    mol.read(traj)
    mol.dropFrames(keep=frameNum)  # Making sure only specific frame to write is kept
    mol.write(path.join(newDir, coorname))

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import write [as 别名]
def _filtSim(i, sims, outFolder, filterSel):
    name = _simName(sims[i].trajectory[0])
    directory = path.join(outFolder, name)
    if not path.exists(directory):
        makedirs(directory)

    logger.debug('Processing trajectory ' + name)

    fmolfile = path.join(outFolder, 'filtered.pdb')
    (traj, outtraj) = _renameSims(sims[i].trajectory, name, outFolder)
    if not traj:
        ftrajectory = _listXTCs(path.join(outFolder, name))
        return Sim(simid=sims[i].simid, parent=sims[i], input=None, trajectory=ftrajectory, molfile=fmolfile)

    try:
        mol = Molecule(sims[i].molfile)
    except:
        logger.warning('Error! Skipping simulation ' + name)
        return

    sel = mol.atomselect(filterSel)

    for j in range(0, len(traj)):
        try:
            mol.read(traj[j])
        except IOError as e:
            logger.warning(e.strerror + ', skipping trajectory')
            break

        mol.write(outtraj[j], sel)

    ftrajectory = _listXTCs(path.join(outFolder, name))
    #bar.progress()
    return Sim(simid=sims[i].simid, parent=sims[i], input=None, trajectory=ftrajectory, molfile=fmolfile)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import write [as 别名]
    def write(self, inputdir, outputdir):
        """ Write the equilibration protocol

        Writes the equilibration protocol and files into a folder for execution
        using files inside the inputdir directory

        Parameters
        ----------
        inputdir : str
            Path to a directory containing the files produced by a build process.
        outputdir : str
            Directory where to write the equilibration setup files.

        Examples
        --------
        >>> md = Equilibration()
        >>> md.write('./build','./equil')
        """
        self._findFiles(inputdir)
        self._amberFixes()

        from htmd.units import convert
        numsteps = convert(self.timeunits, 'timesteps', self.runtime, timestep=self.acemd.timestep)

        pdbfile = os.path.join(inputdir, self.acemd.coordinates)
        inmol = Molecule(pdbfile)

        if self.constraintsteps is None:
            constrsteps = int(numsteps / 2)
        else:
            constrsteps = int(self.constraintsteps)

        tcl = list(self.acemd.TCL)
        tcl[0] = tcl[0].format(NUMSTEPS=numsteps, KCONST=self.fb_k,
                               REFINDEX=' '.join(map(str, inmol.get('index', self.fb_reference))),
                               SELINDEX=' '.join(map(str, inmol.get('index', self.fb_selection))),
                               BOX=' '.join(map(str, self.fb_box)),
                               NVTSTEPS=self.nvtsteps, CONSTRAINTSTEPS=constrsteps, TEMPERATURE=self.temperature)
        self.acemd.TCL = tcl[0] + tcl[1]

        if self.acemd.celldimension is None and self.acemd.extendedsystem is None:
            coords = inmol.get('coords', sel='water')
            if coords.size == 0:  # It's a vacuum simulation
                coords = inmol.get('coords', sel='all')
                dim = np.max(coords, axis=0) - np.min(coords, axis=0)
                dim = dim + 12.
            else:
                dim = np.max(coords, axis=0) - np.min(coords, axis=0)
            self.acemd.celldimension = '{} {} {}'.format(dim[0], dim[1], dim[2])
        if self.useconstantratio:
            self.acemd.useconstantratio = 'on'
        self.acemd.setup(inputdir, outputdir, overwrite=True)

        # Adding constraints
        inmol.set('occupancy', 0)
        inmol.set('beta', 0)
        for sel in self.constraints:
            inmol.set('beta', self.constraints[sel], sel)
        outfile = os.path.join(outputdir, self.acemd.coordinates)
        inmol.write(outfile)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import write [as 别名]
def _filterPDBPSF(sim, outfolder, filtsel):
    try:
        mol = Molecule(sim.molfile)
    except IOError as e:
        raise NameError('simFilter: ' + e.strerror + ' Cannot create filtered.pdb due to problematic pdb: ' + sim.molfile)

    if not path.isfile(path.join(outfolder, 'filtered.pdb')):
        mol.write(path.join(outfolder, 'filtered.pdb'), filtsel)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import write [as 别名]
def _filterPDBPSF(sim, outfolder, filtsel):
    try:
        mol = Molecule(sim.molfile)
    except IOError as e:
        raise NameError('simFilter: {}. Cannot create filtered.pdb due to problematic pdb: {}'.format(e, sim.molfile))

    if not path.isfile(path.join(outfolder, 'filtered.pdb')):
        if mol.coords.size == 0:  # If we read for example psf or prmtop which have no coords, just add 0s everywhere
            mol.coords = np.zeros((mol.numAtoms, 3, 1), dtype=np.float32)
        mol.write(path.join(outfolder, 'filtered.pdb'), filtsel)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import write [as 别名]
    def write(self, inputdir, outputdir):
        """ Writes the production protocol and files into a folder.

        Parameters
        ----------
        inputdir : str
            Path to a directory containing the files produced by a equilibration process.
        outputdir : str
            Directory where to write the production setup files.
        """
        self._findFiles(inputdir)
        self._amberFixes()

        from htmd.units import convert
        numsteps = convert(self.timeunits, 'timesteps', self.runtime, timestep=self.acemd.timestep)
        self.acemd.temperature = str(self.temperature)
        self.acemd.langevintemp = str(self.temperature)
        if self.fb_k > 0: #use TCL only for flatbottom
            mol = Molecule(os.path.join(inputdir, self.acemd.coordinates))
            self.acemd.tclforces = 'on'
            tcl = list(self.acemd.TCL)
            tcl[0] = tcl[0].format(NUMSTEPS=numsteps, KCONST=self.fb_k,
                                   REFINDEX=' '.join(map(str, mol.get('index', self.fb_reference))),
                                   SELINDEX=' '.join(map(str, mol.get('index', self.fb_selection))),
                                   BOX=' '.join(map(str, self.fb_box)))
            self.acemd.TCL = tcl[0] + tcl[1]
        else:
            self.acemd.TCL = 'set numsteps {}\n'.format(numsteps)
        if self.useconstraints:
            # Turn on constraints
            self.acemd.constraints = 'on'
            self.acemd.constraintscaling = '1.0'
        else:
            if len(self.constraints) != 0:
                logger.warning('You have setup constraints to {} but constraints are turned off. '
                               'If you want to use constraints, define useconstraints=True'.format(self.constraints))
        self.acemd.setup(inputdir, outputdir, overwrite=True)

        # Adding constraints
        if self.useconstraints:
            inmol = Molecule(os.path.join(inputdir, self.acemd.coordinates))
            inmol.set('occupancy', 0)
            inmol.set('beta', 0)
            if len(self.constraints) == 0:
                raise RuntimeError('You have set the production to use constraints (useconstraints=True), but have not '
                               'defined any constraints (constraints={}).')
            else:
                for sel in self.constraints:
                    inmol.set('beta', self.constraints[sel], sel)
                outfile = os.path.join(outputdir, self.acemd.coordinates)
                inmol.write(outfile)
                self.acemd.consref = self.acemd.coordinates

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import write [as 别名]
    def write(self, inputdir, outputdir):
        """ Write the equilibration protocol

        Writes the equilibration protocol and files into a folder for execution
        using files inside the inputdir directory

        Parameters
        ----------
        inputdir : str
            Path to a directory containing the files produced by a build process.
        outputdir : str
            Directory where to write the equilibration setup files.

        Examples
        --------
        >>> md = Equilibration()
        >>> md.write('./build','./equil')
        """
        self._findFiles(inputdir)
        self._amberFixes()

        pdbfile = os.path.join(inputdir, self.acemd.coordinates)
        inmol = Molecule(pdbfile)

        self.acemd.TCL = self.acemd.TCL.replace('NUMSTEPS', str(self.numsteps))
        self.acemd.TCL = self.acemd.TCL.replace('TEMPERATURE', str(self.temperature))
        self.acemd.TCL = self.acemd.TCL.replace('KCONST', str(self.k))
        self.acemd.TCL = self.acemd.TCL.replace('REFINDEX', ' '.join(map(str, inmol.get('index', self.reference))))
        self.acemd.TCL = self.acemd.TCL.replace('SELINDEX', ' '.join(map(str, inmol.get('index', self.selection))))
        self.acemd.TCL = self.acemd.TCL.replace('BOX', ' '.join(map(str, self.box)))
        if self.acemd.celldimension is None and self.acemd.extendedsystem is None:
            coords = inmol.get('coords', sel='water')
            if coords.size == 0:  # It's a vacuum simulation
                coords = inmol.get('coords', sel='all')
                dim = np.max(coords, axis=0) - np.min(coords, axis=0)
                dim = dim + 12.
            else:
                dim = np.max(coords, axis=0) - np.min(coords, axis=0)
            self.acemd.celldimension = '{} {} {}'.format(dim[0], dim[1], dim[2])
        if self.useconstantratio:
            self.acemd.useconstantratio = 'on'
        self.acemd.setup(inputdir, outputdir, overwrite=True)

        # Adding constraints
        inmol.set('occupancy', 0)
        inmol.set('beta', 0)
        for sel in self.constraints:
            inmol.set('beta', self.constraints[sel], sel)
        outfile = os.path.join(outputdir, self.acemd.coordinates)
        inmol.write(outfile)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import write [as 别名]
    def _writeInputs(self, simsframes, epoch=None):
        if epoch is None:
            epoch = self._getEpoch() + 1

        test = glob(path.join(self.inputpath, 'e' + str(epoch) + '*'))
        if len(test) != 0:
            raise NameError('Input dirs of epoch ' + str(epoch) + ' already exists.')

        if path.exists(path.join(self.inputpath, 'e' + str(epoch) + '_writeinputs.log')):
            raise NameError('Epoch logfile already exists. Cant overwrite it.')

        fid = open(path.join(self.inputpath, 'e' + str(epoch) + '_writeinputs.log'), 'w')

        regex = re.compile('(e\d+s\d+)_')
        for i, f in enumerate(simsframes):
            frameNum = f.frame
            piece = f.piece
            #print(frameNum)
            if f.sim.parent is None:
                currSim = f.sim
            else:
                currSim = f.sim.parent

            traj = currSim.trajectory[piece]
            if currSim.input is None:
                raise NameError('Could not find input folder in simulation lists. Cannot create new simulations.')

            wuName = _simName(traj)
            res = regex.search(wuName)
            if res:  # If we are running on top of adaptive, use the first name part for the next sim name
                wuName = res.group(1)

            # create new job directory
            newName = 'e' + str(epoch) + 's' + str(i+1) + '_' + wuName + 'p' + str(piece) + 'f' + str(frameNum)
            newDir = path.join(self.inputpath, newName, '')
            # copy previous input directory including input files
            copytree(currSim.input, newDir, symlinks=False, ignore=ignore_patterns('*.dcd', '*.xtc', '*.coor'))
            # overwrite input file with new one. frameNum + 1 as catdcd does 1 based indexing
            mol = Molecule()
            mol.read(traj)
            mol.frame = frameNum
            mol.write(path.join(newDir, 'input.coor'))

            # write nextInput file
            fid.write('# {0} \n{1} {2}\n'.format(newName, traj, frameNum))

        fid.close()

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import write [as 别名]
    def _analyse(self, mol, pdb, rtf, prm, traj, ftraj):
        t = Molecule(pdb)
        t.read(traj)
        t.filter('not water')
        t.write(ftraj)
        m = FFMolecule(filename=mol, rtf=rtf, prm=prm)
        m.read(ftraj)
        torsions = m.getRotatableDihedrals()
        # For each torsion
        for i in range(len(torsions)):
            # Create title
            title = '{}-{}-{}-{}'.format(m.name[torsions[i][0]], m.name[torsions[i][1]], m.name[torsions[i][2]],
                                         m.name[torsions[i][3]])
            # Measure
            (r, theta) = self._measure_torsion(torsions[i], m.coords)

            self._plot_scatter(r, theta, title)
            self._plot_hist(theta, title)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import write [as 别名]
def _filterTopology(sim, outfolder, filtsel):
    from htmd.util import ensurelist
    try:
        from htmd.molecule.molecule import Molecule
        mol = Molecule(sim.molfile)
    except IOError as e:
        raise RuntimeError('simFilter: {}. Cannot read topology file {}'.format(e, sim.molfile))

    if mol.coords.size == 0:  # If we read for example psf or prmtop which have no coords, just add 0s everywhere
        mol.coords = np.zeros((mol.numAtoms, 3, 1), dtype=np.float32)

    extensions = ['pdb',]  # Adding pdb to make sure it's always written
    for m in ensurelist(sim.molfile):
        extensions.append(os.path.splitext(m)[1][1:])

    for ext in list(set(extensions)):
        filttopo = path.join(outfolder, 'filtered.{}'.format(ext))
        if not path.isfile(filttopo):
            try:
                mol.write(filttopo, filtsel)
            except Exception as e:
                logger.warning('Filtering was not able to write {} due to error: {}'.format(filttopo, e))

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import write [as 别名]
def _filtSim(i, sims, outFolder, filterSel):
    name = _simName(sims[i].trajectory[0])
    directory = path.join(outFolder, name)
    if not path.exists(directory):
        makedirs(directory)

    logger.debug('Processing trajectory ' + name)

    fmolfile = path.join(outFolder, 'filtered.pdb')
    (traj, outtraj) = _renameSims(sims[i].trajectory, name, outFolder)
    if not traj:
        ftrajectory = _autoDetectTrajectories(path.join(outFolder, name))
        numframes = _getNumFrames(sims[i], ftrajectory)
        return Sim(simid=sims[i].simid, parent=sims[i], input=None, trajectory=ftrajectory, molfile=fmolfile, numframes=numframes)

    try:
        from htmd.molecule.molecule import Molecule
        mol = Molecule(sims[i].molfile)
    except:
        logger.warning('Error! Skipping simulation ' + name)
        return

    sel = mol.atomselect(filterSel)

    for j in range(0, len(traj)):
        try:
            mol.read(traj[j])
        except IOError as e:
            logger.warning('{}, skipping trajectory'.format(e))
            break

        mol.write(outtraj[j], sel)

    ftrajectory = _autoDetectTrajectories(path.join(outFolder, name))
    numframes = _getNumFrames(sims[i], ftrajectory)
    return Sim(simid=sims[i].simid, parent=sims[i], input=None, trajectory=ftrajectory, molfile=fmolfile, numframes=numframes)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import write [as 别名]
def _filterTopology(sim, outfolder, filtsel):
    try:
        from htmd.molecule.molecule import Molecule
        mol = Molecule(sim.molfile)
    except IOError as e:
        raise RuntimeError('simFilter: {}. Cannot read topology file {}'.format(e, sim.molfile))

    ext = os.path.splitext(sim.molfile)[1][1:]
    filttopo = path.join(outfolder, 'filtered.{}'.format(ext))
    filtpdb = path.join(outfolder, 'filtered.pdb')

    if not path.isfile(filttopo):
        try:
            mol.write(filttopo, filtsel)
        except Exception as e:
            logger.warning('Was not able to write {} due to error: {}'.format(filttopo, e))

        if mol.coords.size == 0:  # If we read for example psf or prmtop which have no coords, just add 0s everywhere
            mol.coords = np.zeros((mol.numAtoms, 3, 1), dtype=np.float32)

        try:
            mol.write(filtpdb, filtsel)
        except Exception as e:
            logger.warning('Was not able to write {} due to error: {}'.format(filtpdb, e))

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import write [as 别名]
    def write(self, inputdir, outputdir):
        """ Writes the production protocol and files into a folder.

        Parameters
        ----------
        inputdir : str
            Path to a directory containing the files produced by a equilibration process.
        outputdir : str
            Directory where to write the production setup files.
        """
        from htmd.molecule.molecule import Molecule

        # Do version consistency check
        if (self._version == 2 and not isinstance(self.acemd, Acemd2)) and \
                (self._version == 3 and not isinstance(self.acemd, Acemd)):
            raise RuntimeError('Acemd object version ({}) inconsistent with protocol version at instantiation '
                               '({})'.format(type(self.acemd), self._version))

        self._findFiles(inputdir)
        self._amberFixes()

        if self._version == 2:
            self.acemd.temperature = str(self.temperature)
            self.acemd.langevintemp = str(self.temperature)
        elif self._version == 3:
            self.acemd.temperature = self.temperature
            self.acemd.thermostattemp = self.temperature

        from htmd.units import convert
        numsteps = convert(self.timeunits, 'timesteps', self.runtime, timestep=self.acemd.timestep)
        if self._version == 3:
            self.acemd.run = str(numsteps)

        pdbfile = os.path.join(inputdir, self.acemd.coordinates)
        inmol = Molecule(pdbfile)

        if np.any(inmol.atomselect('lipids')) and not self.useconstantratio:
            logger.warning('Lipids detected in input structure. We highly recommend setting useconstantratio=True '
                           'for membrane simulations.')

        if self._version == 2:
            if self.restraints:
                raise RuntimeWarning('restraints are only available on {}(_version=3)'.format(self.__class__.__name__))
            if self.fb_k > 0:  # use TCL only for flatbottom
                self.acemd.tclforces = 'on'
                if isinstance(self.acemd.TCL, tuple):
                    tcl = list(self.acemd.TCL)
                    tcl[0] = tcl[0].format(NUMSTEPS=numsteps, TEMPERATURE=self.temperature, KCONST=self.fb_k,
                                           REFINDEX=' '.join(map(str, inmol.get('index', self.fb_reference))),
                                           SELINDEX=' '.join(map(str, inmol.get('index', self.fb_selection))),
                                           BOX=' '.join(map(str, self.fb_box)))
                    self.acemd.TCL = tcl[0] + tcl[1]
                else:
                    logger.warning('{} default TCL was already formatted.'.format(self.__class__.__name__))
            else:
                self.acemd.TCL = 'set numsteps {NUMSTEPS}\n'.format(NUMSTEPS=numsteps)
            if self.useconstraints:
                # Turn on constraints
                self.acemd.constraints = 'on'
                self.acemd.constraintscaling = '1.0'
            else:
                if len(self.constraints) != 0:
                    logger.warning('You have setup constraints to {} but constraints are turned off. '
                                   'If you want to use constraints, define '
                                   'useconstraints=True'.format(self.constraints))
        elif self._version == 3:
            if self.restraints is not None:
                logger.info('Using user-provided restraints and ignoring constraints and fb_potential')
                self.acemd.restraints = self.restraints
            else:
                restraints = list()
                if self.fb_k > 0:
                    logger.warning('Converting fb_potential to restraints. This is a convenience '
                                   'functional conversion. We recommend start using restraints with '
                                   '{}(_version=3)'.format(self.__class__.__name__))
                    restraints += self._fb_potential2restraints(inputdir)
                if self.useconstraints:
                    logger.warning('Converting constraints to restraints. This is a convenience '
                                   'functional conversion. We recommend start using restraints with '
                                   '{}(_version=3)'.format(self.__class__.__name__))
                    restraints += self._constraints2restraints()
                else:
                    if len(self.constraints) != 0:
                        logger.warning('You have setup constraints to {} but constraints are turned off. '
                                       'If you want to use constraints, define '
                                       'useconstraints=True'.format(self.constraints))
                if len(restraints) != 0:
                    self.acemd.restraints = restraints

        if self.useconstantratio:
            self.acemd.useconstantratio = 'on'

        if self.adaptive:
            self.acemd.binvelocities = None

        self.acemd.setup(inputdir, outputdir, overwrite=True)

        if self._version == 2:
            # Adding constraints by writing them to the consref file
            if self.useconstraints:
#.........这里部分代码省略.........

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import write [as 别名]
    def write(self, inputdir, outputdir):
        """ Write the equilibration protocol

        Writes the equilibration protocol and files into a folder for execution
        using files inside the inputdir directory

        Parameters
        ----------
        inputdir : str
            Path to a directory containing the files produced by a build process.
        outputdir : str
            Directory where to write the equilibration setup files.

        Examples
        --------
        >>> md = Equilibration()
        >>> md.write('./build','./equil')
        """
        self._findFiles(inputdir)
        self._amberFixes()

        from htmd.units import convert
        numsteps = convert(self.timeunits, 'timesteps', self.runtime, timestep=self.acemd.timestep)

        pdbfile = os.path.join(inputdir, self.acemd.coordinates)
        inmol = Molecule(pdbfile)

        if np.any(inmol.atomselect('lipids')) and not self.useconstantratio:
            logger.warning('Lipids detected in input structure. We highly recommend setting useconstantratio=True '
                           'for membrane simulations.')

        if self.constraintsteps is None:
            constrsteps = int(numsteps / 2)
        else:
            constrsteps = int(self.constraintsteps)

        if isinstance(self.acemd.TCL, tuple):
            tcl = list(self.acemd.TCL)
            tcl[0] = tcl[0].format(NUMSTEPS=numsteps, KCONST=self.fb_k,
                                   REFINDEX=' '.join(map(str, inmol.get('index', self.fb_reference))),
                                   SELINDEX=' '.join(map(str, inmol.get('index', self.fb_selection))),
                                   BOX=' '.join(map(str, self.fb_box)),
                                   NVTSTEPS=self.nvtsteps, CONSTRAINTSTEPS=constrsteps, TEMPERATURE=self.temperature)
            self.acemd.TCL = tcl[0] + tcl[1]
        else:
            logger.warning('{} default TCL was already formatted.'.format(self.__class__.__name__))

        if self.acemd.celldimension is None and self.acemd.extendedsystem is None:
            coords = inmol.get('coords', sel='water')
            if coords.size == 0:  # It's a vacuum simulation
                coords = inmol.get('coords', sel='all')
                dim = np.max(coords, axis=0) - np.min(coords, axis=0)
                dim = dim + 12.
            else:
                dim = np.max(coords, axis=0) - np.min(coords, axis=0)
            self.acemd.celldimension = '{} {} {}'.format(dim[0], dim[1], dim[2])

        if self.useconstantratio:
            self.acemd.useconstantratio = 'on'

        self.acemd.setup(inputdir, outputdir, overwrite=True)

        # Adding constraints by writing them to the consref file
        inconsreffile = os.path.join(inputdir, self.acemd.consref)
        consrefmol = Molecule(inconsreffile)
        consrefmol.set('occupancy', 0)
        consrefmol.set('beta', 0)
        if len(self.constraints) == 0:
            raise RuntimeError('You have not defined any constraints for the Equilibration (constraints={}).')
        else:
            for sel in self.constraints:
                consrefmol.set('beta', self.constraints[sel], sel)
        outconsreffile = os.path.join(outputdir, self.acemd.consref)
        consrefmol.write(outconsreffile)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import write [as 别名]
def build(mol, ff=None, topo=None, param=None, prefix='structure', outdir='./build', caps=None, ionize=True, saltconc=0,
          saltanion=None, saltcation=None, disulfide=None, tleap='tleap', execute=True):
    """ Builds a system for AMBER

    Uses tleap to build a system for AMBER. Additionally it allows the user to ionize and add disulfide bridges.

    Parameters
    ----------
    mol : :class:`Molecule <htmd.molecule.molecule.Molecule>` object
        The Molecule object containing the system
    ff : list of str
        A list of leaprc forcefield files. Default: ['leaprc.lipid14', 'leaprc.ff14SB', 'leaprc.gaff']
    topo : list of str
        A list of topology `prepi` files.
    param : list of str
        A list of parameter `frcmod` files.
    prefix : str
        The prefix for the generated pdb and psf files
    outdir : str
        The path to the output directory
        Default: './build'
    caps : dict
        A dictionary with keys segids and values lists of strings describing the caps for a particular protein segment.
        e.g. caps['P'] = ['ACE', 'NME'] or caps['P'] = ['none', 'none']. Default: will apply ACE and NME caps to every
        protein segment.
    ionize : bool
        Enable or disable ionization
    saltconc : float
        Salt concentration to add to the system after neutralization.
    saltanion : {'Cl-'}
        The anion type. Please use only AMBER ion atom names.
    saltcation : {'Na+', 'K+', 'Cs+'}
        The cation type. Please use only AMBER ion atom names.
    disulfide : np.ndarray
        If None it will guess disulfide bonds. Otherwise provide a 2D array where each row is a pair of atom indexes that makes a disulfide bond
    tleap : str
        Path to tleap executable used to build the system for AMBER
    execute : bool
        Disable building. Will only write out the input script needed by tleap. Does not include ionization.

    Returns
    -------
    molbuilt : :class:`Molecule <htmd.molecule.molecule.Molecule>` object
        The built system in a Molecule object

    Example
    -------
    >>> ffs = ['leaprc.lipid14', 'leaprc.ff14SB', 'leaprc.gaff']
    >>> molbuilt = amber.build(mol, ff=ffs, outdir='/tmp/build', saltconc=0.15)
    """
    # Remove pdb bonds!
    mol = mol.copy()
    mol.bonds = np.empty((0, 2), dtype=np.uint32)
    if shutil.which(tleap) is None:
        raise NameError('Could not find executable: `' + tleap + '` in the PATH. Cannot build for AMBER.')
    if not os.path.isdir(outdir):
        os.makedirs(outdir)
    _cleanOutDir(outdir)
    if ff is None:
        ff = ['leaprc.lipid14', 'leaprc.ff14SB', 'leaprc.gaff']
    if topo is None:
        topo = []
    if param is None:
        param = []
    if caps is None:
        caps = _defaultProteinCaps(mol)

    _missingSegID(mol)
    _checkMixedSegment(mol)

    logger.info('Converting CHARMM membranes to AMBER.')
    mol = _charmmLipid2Amber(mol)

    #_checkProteinGaps(mol)
    _applyProteinCaps(mol, caps)

    f = open(path.join(outdir, 'tleap.in'), 'w')
    f.write('# tleap file generated by amber.build\n')

    # Printing out the forcefields
    if isinstance(ff, str):
        ff = [ff]
    for force in ff:
        f.write('source ' + force + '\n')
    f.write('\n')

    # Loading TIP3P water parameters
    f.write('# Loading ions and TIP3P water parameters\n')
    f.write('loadamberparams frcmod.ionsjc_tip3p\n\n')

    # Loading user parameters
    f.write('# Loading parameter files\n')
    for p in param:
        try:
            shutil.copy(p, outdir)
            f.write('loadamberparams ' + path.basename(p) + '\n')
        except:
            f.write('loadamberparams ' + p + '\n')
            logger.info("Path {:s} not found, assuming a standard AmberTools file.".
                        format(p))
#.........这里部分代码省略.........