Python PDBFile.writeModel方法代码示例

# 需要导入模块: from simtk.openmm.app import PDBFile [as 别名]
# 或者: from simtk.openmm.app.PDBFile import writeModel [as 别名]
    def update(self):
        """
        Update the sampler with one step of sampling.
        """
        if not self._initialized:
            self._initialize()

        if self.verbose:
            print("." * 80)
            print("MCMC sampler iteration %d" % self.iteration)

        initial_time = time.time()

        # Reset statistics
        if hasattr(self.integrator, 'setGlobalVariableByName'):
            self.integrator.setGlobalVariableByName('naccept', 0)

        # Take some steps
        self.integrator.step(self.nsteps)

        # Get new sampler state.
        self.sampler_state = SamplerState.createFromContext(self.context)

        # Report statistics
        if hasattr(self.integrator, 'getGlobalVariableByName'):
            naccept = self.integrator.getGlobalVariableByName('naccept')
            fraction_accepted = float(naccept) / float(self.nsteps)
            if self.verbose: print("Accepted %d / %d GHMC steps (%.2f%%)." % (naccept, self.nsteps, fraction_accepted * 100))

        final_time = time.time()
        elapsed_time = final_time - initial_time
        self._timing['sample positions'] = elapsed_time

        if self.verbose:
            final_energy = self.context.getState(getEnergy=True).getPotentialEnergy() * self.thermodynamic_state.beta
            print('Final energy is %12.3f kT' % (final_energy))
            print('elapsed time %8.3f s' % elapsed_time)

        if self.ncfile:
            self.ncfile.variables['positions'][self.iteration,:,:] = self.sampler_state.positions[:,:] / unit.nanometers
            for k in range(3):
                self.ncfile.variables['box_vectors'][self.iteration,k,:] = self.sampler_state.box_vectors[k,:] / unit.nanometers
            self.ncfile.variables['potential'][self.iteration] = self.thermodynamic_state.beta * self.context.getState(getEnergy=True).getPotentialEnergy()
            self.ncfile.variables['sample_positions_time'][self.iteration] = elapsed_time

        # Increment iteration count
        self.iteration += 1

        if self.verbose:
            print("." * 80)
        if self.pdbfile is not None:
            print("Writing frame...")
            from simtk.openmm.app import PDBFile
            PDBFile.writeModel(self.topology, self.sampler_state.positions, self.pdbfile, self.iteration)
            self.pdbfile.flush()

# 需要导入模块: from simtk.openmm.app import PDBFile [as 别名]
# 或者: from simtk.openmm.app.PDBFile import writeModel [as 别名]
    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
        """
        if self._nextModel == 0:
            PDBFile.writeHeader(simulation.topology, self._out)
            self._topology = simulation.topology
            self._nextModel += 1
        PDBFile.writeModel(simulation.topology, state.getPositions(), self._out, self._nextModel)
        self._nextModel += 1

# 需要导入模块: from simtk.openmm.app import PDBFile [as 别名]
# 或者: from simtk.openmm.app.PDBFile import writeModel [as 别名]
    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
        """
        if self._nextModel == 0:
            PDBFile.writeHeader(simulation.topology, self._out)
            self._topology = simulation.topology
            self._nextModel += 1
        PDBFile.writeModel(simulation.topology, state.getPositions(), self._out, self._nextModel)
        self._nextModel += 1
        if hasattr(self._out, 'flush') and callable(self._out.flush):
            self._out.flush()

# 需要导入模块: from simtk.openmm.app import PDBFile [as 别名]
# 或者: from simtk.openmm.app.PDBFile import writeModel [as 别名]
    def update(self):
        """
        Update the sampler with one step of sampling.
        """
        if self.verbose:
            print("-" * 80)
            print("Expanded Ensemble sampler iteration %8d" % self.iteration)
        self.update_positions()
        self.update_state()
        self.iteration += 1
        if self.verbose:
            print("-" * 80)

        if self.pdbfile is not None:
            print("Writing frame...")
            from simtk.openmm.app import PDBFile
            PDBFile.writeModel(self.topology, self.sampler.sampler_state.positions, self.pdbfile, self.iteration)
            self.pdbfile.flush()

# 需要导入模块: from simtk.openmm.app import PDBFile [as 别名]
# 或者: from simtk.openmm.app.PDBFile import writeModel [as 别名]
    def report(self, simulation, _):
        """Generate a report.

        Parameters
        ----------
        simulation : Simulation
            The Simulation to generate a report for
        _ : State
            The current state of the simulation
        """
        state = simulation.context.getState(getPositions=True, enforcePeriodicBox=self._enforcePeriodicBox)
        if self._nextModel == 0:
            PDBFile.writeHeader(simulation.topology, self._out)
            self._topology = simulation.topology
            self._nextModel += 1
        PDBFile.writeModel(simulation.topology, state.getPositions(), self._out, self._nextModel)
        self._nextModel += 1
        if hasattr(self._out, 'flush') and callable(self._out.flush):
            self._out.flush()

# 需要导入模块: from simtk.openmm.app import PDBFile [as 别名]
# 或者: from simtk.openmm.app.PDBFile import writeModel [as 别名]

#.........这里部分代码省略.........
                context.setParameter(parameter_name, 1)

        # Compute initial potential of alchemical state.
        initial_potential = self.beta * context.getState(getEnergy=True).getPotentialEnergy()
        if np.isnan(initial_potential):
            raise NaNException("Initial potential of 'insert' operation is NaN (unmodified potential was %.3f kT, alchemical potential was %.3f kT before changing lambda)" % (unmodified_potential, alchemical_potential))
        from perses.tests.utils import compute_potential_components
        #print("initial potential before '%s' : %f kT" % (direction, initial_potential))
        #print("initial potential components:   %s" % str(compute_potential_components(context))) # DEBUG

        # Take a single integrator step since all switching steps are unrolled in NCMCAlchemicalIntegrator.
        try:
            # Write PDB file if requested.
            if self.write_pdb_interval:
                if direction == 'insert':
                    topology = topology_proposal.new_topology
                    indices = topology_proposal.unique_new_atoms
                else:
                    topology = topology_proposal.old_topology
                    indices = topology_proposal.unique_old_atoms

                # Write atom indices that are changing
                import pickle
                filename = 'ncmc-%s-%d-atomindices.pkl' % (direction, self.nattempted)
                outfile = open(filename, 'wb')
                pickle.dump(indices, outfile)
                outfile.close()

                from simtk.openmm.app import PDBFile
                filename = 'ncmc-%s-%d.pdb' % (direction, self.nattempted)
                outfile = open(filename, 'w')
                PDBFile.writeHeader(topology, file=outfile)
                modelIndex = 0
                PDBFile.writeModel(topology, context.getState(getPositions=True).getPositions(asNumpy=True), file=outfile, modelIndex=modelIndex)
                try:
                    for step in range(self.nsteps):
                        integrator.step(1)
                        if (step+1)%self.write_pdb_interval == 0:
                            modelIndex += 1
                            PDBFile.writeModel(topology, context.getState(getPositions=True).getPositions(asNumpy=True), file=outfile, modelIndex=modelIndex)
                except ValueError as e:
                    # System is exploding and coordinates won't fit in PDB ATOM fields
                    print(e)

                PDBFile.writeFooter(topology, file=outfile)
                outfile.close()
            else:
                integrator.step(self.nsteps)

        except Exception as e:
            # Trap NaNs as a special exception (allowing us to reject later, if desired)
            if str(e) == "Particle coordinate is nan":
                raise NaNException(str(e))
            else:
                raise e

        # Set final context parameters.
        if direction == 'insert':
            for parameter_name in available_parameters:
                context.setParameter(parameter_name, 1)
        elif direction == 'delete':
            for parameter_name in available_parameters:
                context.setParameter(parameter_name, 0)

        # Compute final potential of alchemical state.
        final_potential = self.beta * context.getState(getEnergy=True).getPotentialEnergy()

# 需要导入模块: from simtk.openmm.app import PDBFile [as 别名]
# 或者: from simtk.openmm.app.PDBFile import writeModel [as 别名]
niterations = 500 # number of iterations
mctrials = 10 # number of Monte Carlo trials per iteration
nsalt = 0 # current number of salt pairs
tol = 1e-6 # constraint tolerance

# Determine number of molecules
nmolecules = 0
for residue in topology.residues():
    nmolecules += 1
print('system originally has %d water molecules' % nmolecules)

# Open PDB file for writing.
from simtk.openmm.app import PDBFile
pdbfile = open('output.pdb', 'w')
PDBFile.writeHeader(topology, file=pdbfile)
PDBFile.writeModel(topology, positions, file=pdbfile, modelIndex=0)

# Simulate
for iteration in range(niterations):
    print('iteration %5d / %5d' % (iteration, niterations))

    # Create a simulation
    from openmmtools.integrators import VelocityVerletIntegrator
    integrator = VelocityVerletIntegrator(timestep)
    integrator.setConstraintTolerance(tol)
    context = openmm.Context(system, integrator)
    context.setPositions(positions)

    # Propagate dynamics at constant counterion number.
    print('propagating dynamics for %d steps...' % nsteps)
    context.setVelocitiesToTemperature(temperature)