Python Molecule.write方法代码示例

# 需要导入模块: from molecule import Molecule [as 别名]
# 或者: from molecule.Molecule import write [as 别名]
_exec("trjconv_d -s eq.tpr -f eq.trr -o eq.gro -ndec 9 -pbc mol -dump 1", stdin="0\n")

# Confirm that constraints are satisfied
OMM_eqgro = app.GromacsGroFile('eq.gro')
OMM_prmtop = app.AmberPrmtopFile('prmtop')
system = OMM_prmtop.createSystem(nonbondedMethod=app.NoCutoff)
integ = mm.VerletIntegrator(1.0*u.femtosecond)
plat = mm.Platform.getPlatformByName('Reference')
simul = app.Simulation(OMM_prmtop.topology, system, integ)
simul.context.setPositions(OMM_eqgro.positions)
simul.context.applyConstraints(1e-12)
state = simul.context.getState(getPositions=True)
pos = np.array(state.getPositions().value_in_unit(u.angstrom)).reshape(-1,3)
M = Molecule('eq.gro')
M.xyzs[0] = pos
M.write('constrained.gro')

# Gromacs calculation
GMX_Energy, GMX_Force, Ecomps_GMX = Calculate_GMX('constrained.gro', 'topol.top', 'shot.mdp')
GMX_Force = GMX_Force.reshape(-1,3)

# Print Gromacs energy components
printcool_dictionary(Ecomps_GMX, title="GROMACS energy components")

# Parse the .mdp file to inform ParmEd
defines, sysargs, mdp_opts = interpret_mdp('shot.mdp')

parm = parmed.amber.AmberParm('prmtop', 'inpcrd')
GmxGro = parmed.gromacs.GromacsGroFile.parse('constrained.gro')
parm.box = GmxGro.box
parm.positions = GmxGro.positions

# 需要导入模块: from molecule import Molecule [as 别名]
# 或者: from molecule.Molecule import write [as 别名]
_exec("trjconv_d -s eq.tpr -f eq.trr -o eq.gro -ndec 9 -pbc mol -dump 1", stdin="0\n")

# Confirm that constraints are satisfied
OMM_eqgro = app.GromacsGroFile("eq.gro")
OMM_prmtop = app.AmberPrmtopFile("prmtop")
system = OMM_prmtop.createSystem(nonbondedMethod=app.NoCutoff)
integ = mm.VerletIntegrator(1.0 * u.femtosecond)
plat = mm.Platform.getPlatformByName("Reference")
simul = app.Simulation(OMM_prmtop.topology, system, integ)
simul.context.setPositions(OMM_eqgro.positions)
simul.context.applyConstraints(1e-12)
state = simul.context.getState(getPositions=True)
pos = np.array(state.getPositions().value_in_unit(u.angstrom)).reshape(-1, 3)
M = Molecule("eq.gro")
M.xyzs[0] = pos
M.write("constrained.gro")

# Gromacs calculation
GMX_Energy, GMX_Force, Ecomps_GMX = Calculate_GMX("constrained.gro", "topol.top", "shot.mdp")
GMX_Force = GMX_Force.reshape(-1, 3)

# Print Gromacs energy components
printcool_dictionary(Ecomps_GMX, title="GROMACS energy components")

# Parse the .mdp file to inform ParmEd
defines, sysargs, mdp_opts = interpret_mdp("shot.mdp")

parm = parmed.amber.AmberParm("prmtop", "inpcrd")
GmxGro = parmed.gromacs.GromacsGroFile.parse("constrained.gro")
parm.box = GmxGro.box
parm.positions = GmxGro.positions

# 需要导入模块: from molecule import Molecule [as 别名]
# 或者: from molecule.Molecule import write [as 别名]
class TinkerReporter(object):
    """TinkerReporter

    To use it, create a TinkerReporter, then add it to the Simulation's list of reporters.
    """

    def __init__(self, file, reportInterval, simulation, xyzfile, pbc=True, pmegrid=None, tinkerpath = ''):
        """Create a TinkerReporter.

        Parameters:
         - tinkerpath (string) The
         - file (string) The file to write to
         - reportInterval (int) The interval (in time steps) at which to write frames
        """
        print "Initializing Tinker Reporter"
        self._reportInterval = reportInterval
        self._openedFile = isinstance(file, str)
        if self._openedFile:
            self._out = open(file, 'w')
        else:
            self._out = file
        self._pbc = pbc
        self._pmegrid = pmegrid
        self._tinkerpath = tinkerpath
        self._simulation = simulation

        self.xyzfile = xyzfile
        self.pdbfile = os.path.splitext(xyzfile)[0] + '.pdb'
        self.keyfile = os.path.splitext(xyzfile)[0] + '.key'
        #self.prmfile = os.path.splitext(xyzfile)[0] + '.prm'
        self.dynfile = os.path.splitext(xyzfile)[0] + '.dyn'

        print "Loading Tinker xyz file"
        if not os.path.exists(self.xyzfile):
            raise IOError('You need a Tinker .xyz')
        if not os.path.exists(self.keyfile):
            raise IOError('You need a Tinker .key file with the same base name as the .xyz file')
        #if not os.path.exists(self.prmfile):
        #    raise IOError('You need a Tinker .prm file with the same base name as the .xyz file')
        if not os.path.exists(self.pdbfile):
            raise IOError('You need a .pdb file with the same base name as the .xyz file (the same one you used to start the simulation)')
        self.M = Molecule(self.xyzfile,ftype='tinker')
        self.comm = open(self.xyzfile).readlines()[0]
        print "Done"

    def describeNextReport(self, simulation):
        """Get information about the next report this object will generate.

        Parameters:
         - simulation (Simulation) The Simulation to generate a report for
        Returns: A five element tuple.  The first element is the number of steps until the
        next report.  The remaining elements specify whether that report will require
        positions, velocities, forces, and energies respectively.
        """
        steps = self._reportInterval - simulation.currentStep%self._reportInterval
        return (steps, True, True, True, True)

    def report(self, simulation, state):
        """Generate a report.

        Parameters:
         - simulation (Simulation) The Simulation to generate a report for
         - state (State) The current state of the simulation
        """

        RunDynamic = 0
        pos = state.getPositions() / angstrom
        if RunDynamic:
            # Write a dyn file.
            dynout = open(self.dynfile,'w')
            print >> dynout, ' Number of Atoms and Title :'
            print >> dynout, self.comm,
            print >> dynout, ' Periodic Box Dimensions :'
            print >> dynout, pvec([state.getPeriodicBoxVectors()[i][i] / angstrom for i in range(3)])
            print >> dynout, pvec([90,90,90])
            print >> dynout, ' Current Atomic Positions :'
            for i in range(len(pos)):
                print >> dynout, pvec(pos[i])
            print >> dynout, ' Current Atomic Velocities :'
            vel = state.getVelocities() / angstrom * picosecond
            for i in range(len(vel)):
                print >> dynout, pvec(vel[i])
            print >> dynout, ' Current Atomic Accelerations :'
            for i in range(len(vel)):
                print >> dynout, pvec([0,0,0])
            print >> dynout, ' Alternate Atomic Accelerations :'
            for i in range(len(vel)):
                print >> dynout, pvec([0,0,0])
            dynout.close()
            o = _exec('%s %s 1 1e-5 1 1' % (os.path.join(self._tinkerpath,'dynamic'), self.xyzfile), print_command = False)
            for line in o.split('\n'):
                if 'Total Energy' in line:
                    total = float(line.split()[2]) * 4.184
                elif 'Potential Energy' in line:
                    pot = float(line.split()[2]) * 4.184
                elif 'Kinetic Energy' in line:
                    kin = float(line.split()[2]) * 4.184
                elif 'Temperature' in line:
                    temp = float(line.split()[1])
                elif 'Pressure' in line:
#.........这里部分代码省略.........

# 需要导入模块: from molecule import Molecule [as 别名]
# 或者: from molecule.Molecule import write [as 别名]

#.........这里部分代码省略.........
        # and "scf_guess_mix 5" which allows us to find broken
        # symmetry states.
        self.coreguess = True
        # Error message if the calculation failed for a known reason
        self.errmsg = ''
        # qcdsav is "known-good qcdir for this object",
        # used to restore from failed calcs (e.g. SCF failure)
        self.qcdsav = self.qcdir+'sav'
        #--------
        # The clean option makes sure nothing on the disk influences this calculation.
        # This can be a bit confusing.  There are two modes of usage:
        # 1) Clean OFF.  Calculation uses whatever is in qcdir and backs it up to qcdsav on successful calcs.
        # 2) Clean ON.  qcdir is always cleared, and copied over from qcdsav (if exist) prior to calling Q-Chem.
        # This allows us to save the state of a good calculation without worrying about outside interference.
        # - Use case 1: AnalyzeReaction.py does not like to read SCF guesses from previous calculations so we use clean = True.
        # - Use case 2: Growing string does like to read SCF guesses so we use clean = False.
        # - Use case 3: IRC calculation requires Hessian from a previous calculation, so again we use clean = False.
        self.clean = clean
        # If readsave is set, then copy it to self.qcdsav and it will be used
        # to initialize this calculation.  Otherwise self.qcdsav will be removed.
        self.readsave = readsave
        if isinstance(self.readsave, str):
            if not os.path.isdir(self.readsave):
                raise RuntimeError('Tried to initialize Q-Chem reading from a save folder but does not exist')
	    if self.readsave == self.qcdsav: pass
            elif os.path.exists(self.qcdsav):
                shutil.rmtree(self.qcdsav)
                shutil.copytree(self.readsave, self.qcdsav)
        elif isinstance(self.readsave, int) and self.readsave: pass
        elif os.path.exists(self.qcdsav): shutil.rmtree(self.qcdsav)
        # Remove self.qcdir; it will be restored from self.qcdsav right before calling Q-Chem.
        if os.path.exists(self.qcdir): shutil.rmtree(self.qcdir)

    def write(self, *args, **kwargs):
        """ Write the Molecule object to a file. """
        self.M.write(*args, **kwargs)

    def write_qcin(self):
        """ Write Q-Chem input file. """
        rems = OrderedDict([('jobtype', self.jobtype)])
        rems['scf_convergence'] = 8
        # If not the first calculation, read SCF guess from the first calculation.
        if self.readguess and os.path.exists(self.qcdsav):
            rems['scf_guess'] = 'read'
            rems['scf_guess_mix'] = None
        elif self.coreguess:
            rems['scf_guess'] = 'core'
            rems['scf_guess_mix'] = 5
        # Add SCF convergence rem variables.
        rems.update(self.remscf)
        # Add job-related rem variables.
        rems.update(self.remextra)
        # If doing stability analysis, loosen SCF convergence tolerance by 1.
        # This is a bootleg solution to our workflow hanging indefinitely
        # when Q-Chem crashes.
        if 'stability_analysis' in rems.keys():
            rems['scf_convergence'] -= 2
        # Create copy of stored Molecule object, update
        # Q-Chem rem variables and write Q-Chem input file.
        M1 = deepcopy(self.M)
        M1.edit_qcrems(rems)
        M1.write(self.qcin, ftype="qcin")

    def DIE(self, errmsg):
        """ Does what it says. """
        raise RuntimeError("Error: Q-Chem calculation failed! (%s)" % errmsg)