Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def run_simulation(pdb,settings,pbc=True,Trajectory=True):
""" Run a NPT simulation and gather statistics. """
simulation, system = create_simulation_object(pdb, settings, pbc, "single")
# Set initial positions.
simulation.context.setPositions(pdb.positions)
print "Minimizing the energy... (starting energy % .3f kJ/mol)" % simulation.context.getState(getEnergy=True).getPotentialEnergy().value_in_unit(kilojoule_per_mole),
simulation.minimizeEnergy()
print "Done (final energy % .3f kJ/mol)" % simulation.context.getState(getEnergy=True).getPotentialEnergy().value_in_unit(kilojoule_per_mole)
# Assign velocities.
velocities = generateMaxwellBoltzmannVelocities(system, temperature)
simulation.context.setVelocities(velocities)
if verbose:
# Print out the platform used by the context
print "I'm using the platform", simulation.context.getPlatform().getName()
# Print out the properties of the platform
printcool_dictionary({i:simulation.context.getPlatform().getPropertyValue(simulation.context,i) for i in simulation.context.getPlatform().getPropertyNames()},title="Platform %s has properties:" % simulation.context.getPlatform().getName())
# Serialize the system if we want.
Serialize = 0
if Serialize:
serial = XmlSerializer.serializeSystem(system)
with wopen('system.xml') as f: f.write(serial)
#==========================================#
# Computing a bunch of initial values here #
#==========================================#
if pbc:
# Show initial system volume.
box_vectors = system.getDefaultPeriodicBoxVectors()
volume = compute_volume(box_vectors)
if verbose: print "initial system volume = %.1f nm^3" % (volume / nanometers**3)
# Determine number of degrees of freedom.
kB = BOLTZMANN_CONSTANT_kB * AVOGADRO_CONSTANT_NA
# The center of mass motion remover is also a constraint.
# Then look at the data values to stretch out the other column widths.
for v in data.values():
for n, f in enumerate(v):
cwidths[n+1] = max(cwidths[n+1], len(str(f)))
for i in range(1, len(cwidths)):
cwidths[i] += 2
if cwidths[0] < 25:
cwidths[0] = 25
cblocks = [['' for i in range(max(crows) - len(cname.split('\n')))] + cname.split('\n') for cnum, cname in enumerate(headings)]
# The formatting line consisting of variable column widths
fline = ' '.join("%%%s%is" % (("-" if i==0 else ""), j) for i, j in enumerate(cwidths))
vline = ' '.join(["%%%is" % j for i, j in enumerate(cwidths) if i > 0])
clines = [fline % (tuple(cblocks[j][i] for j in range(nc))) for i in range(max(crows))]
tlines += clines
PrintDict = OrderedDict([(key, vline % (tuple(val))) for key, val in data.items()])
printcool_dictionary(PrintDict, title='\n'.join(tlines), keywidth=cwidths[0], center=False, leftpad=4, color=color)
self.Targets = []
enable_smirnoff_prints = False # extra prints for SMIRNOFF forcefield
for opts in tgt_opts:
if opts['type'] not in Implemented_Targets:
logger.error('The target type \x1b[1;91m%s\x1b[0m is not implemented!\n' % opts['type'])
raise RuntimeError
elif opts['type'].endswith("SMIRNOFF"):
enable_smirnoff_prints = True
# Create a target object. This is done by looking up the
# Target class from the Implemented_Targets dictionary
# using opts['type'] as the key. The object is created by
# passing (options, opts, forcefield) to the constructor.
if opts["remote"] and self.wq_port != 0: Tgt = forcebalance.target.RemoteTarget(options, opts, forcefield)
else: Tgt = Implemented_Targets[opts['type']](options,opts,forcefield)
self.Targets.append(Tgt)
printcool_dictionary(Tgt.PrintOptionDict,"Setup for target %s :" % Tgt.name)
if len(set([Tgt.name for Tgt in self.Targets])) != len([Tgt.name for Tgt in self.Targets]):
logger.error("The list of target names is not unique!\n")
raise RuntimeError
if enable_smirnoff_prints:
smirnoff_analyze_parameter_coverage(forcefield, tgt_opts)
## The force field (it seems to be everywhere)
self.FF = forcefield
## Initialize the penalty function.
self.Penalty = Penalty(self.penalty_type,forcefield,self.penalty_additive,
self.penalty_multiplicative,self.penalty_hyperbolic_b,
self.penalty_alpha,self.penalty_power)
## Obtain the denominator.
if self.normalize_weights:
self.WTot = np.sum([i.weight for i in self.Targets])
else:
self.WTot = 1.0
# Add barostat.
system.addForce(barostat)
# Create integrator.
integrator = openmm.LangevinIntegrator(temperature, collision_frequency, timestep)
# Create simulation class.
simulation = Simulation(pdb.topology, system, integrator, platform)
# Set initial positions.
simulation.context.setPositions(pdb.positions)
# Assign velocities.
velocities = generateMaxwellBoltzmannVelocities(system, temperature)
simulation.context.setVelocities(velocities)
if verbose:
# Print out the platform used by the context
print "I'm using the platform", simulation.context.getPlatform().getName()
# Print out the properties of the platform
printcool_dictionary({i:simulation.context.getPlatform().getPropertyValue(simulation.context,i) for i in simulation.context.getPlatform().getPropertyNames()},title="Platform %s has properties:" % simulation.context.getPlatform().getName())
# Serialize the system if we want.
Serialize = 0
if Serialize:
serial = openmm.XmlSerializer.serializeSystem(system)
with open('system.xml','w') as f: f.write(serial)
#==========================================#
# Computing a bunch of initial values here #
#==========================================#
if pbc:
# Show initial system volume.
box_vectors = system.getDefaultPeriodicBoxVectors()
volume = compute_volume(box_vectors)
if verbose: print "initial system volume = %.1f nm^3" % (volume / units.nanometers**3)
# Determine number of degrees of freedom.
kB = units.BOLTZMANN_CONSTANT_kB * units.AVOGADRO_CONSTANT_NA
ndof = 3*system.getNumParticles() - system.getNumConstraints()
if AGrad:
self.FF.print_map(vals=self.Gp[key])
logger.info(bar)
PrintDict[key] = (("% 10.5f % 8.3f % 14.5e" %
(self.Xp[key], self.Wp[key],
self.Xp[key]*self.Wp[key])))
for i, q in enumerate(self.quantities):
print_item(q, self.points[0].ref["units"][i])
PrintDict['Total'] = "% 10s % 8s % 14.5e" % ("","", self.Objective)
Title = ("%s Thermodynamic Properties:\n %-20s %40s" %
(self.name, "Property", "Residual x Weight = Contribution"))
printcool_dictionary(PrintDict, color=4, title=Title, keywidth=31)
return
def indicate(self):
printcool_dictionary(self.PrintDict,title="Interaction Energies (kcal/mol), Objective = % .5e\n %-20s %9s %9s %9s %11s" %
(self.energy_part, "Interaction", "Calc.", "Ref.", "Delta", "Term"))
if len(self.RMSDDict) > 0:
printcool_dictionary(self.RMSDDict,title="Geometry Optimized Systems (Angstrom), Objective = %.5e\n %-38s %11s %11s" % (self.rmsd_part, "System", "RMSD", "Term"), keywidth=45)
# MD options; time step (fs), production steps, equilibration steps, interval for saving data (ps)
lipid_timestep = TgtOptions['lipid_timestep']
lipid_nsteps = TgtOptions['lipid_md_steps']
lipid_nequil = TgtOptions['lipid_eq_steps']
lipid_intvl = TgtOptions['lipid_interval']
lipid_fnm = TgtOptions['lipid_coords']
# Number of threads, multiple timestep integrator, anisotropic box etc.
threads = TgtOptions.get('md_threads', 1)
mts = TgtOptions.get('mts_integrator', 0)
force_cuda = TgtOptions.get('force_cuda', 0)
anisotropic = TgtOptions.get('anisotropic_box', 0)
minimize = TgtOptions.get('minimize_energy', 1)
# Print all options.
printcool_dictionary(TgtOptions, title="Options from ForceBalance")
lipid_snapshots = int((lipid_nsteps * lipid_timestep / 1000) / lipid_intvl)
lipid_iframes = int(1000 * lipid_intvl / lipid_timestep)
logger.info("For the condensed phase system, I will collect %i snapshots spaced apart by %i x %.3f fs time steps\n" \
% (lipid_snapshots, lipid_iframes, lipid_timestep))
if lipid_snapshots < 2:
raise Exception('Please set the number of lipid time steps so that you collect at least two snapshots (minimum %i)' \
% (2000 * int(lipid_intvl/lipid_timestep)))
#----
# Loading coordinates
#----
ML = Molecule(lipid_fnm, toppbc=True)
# Determine the number of molecules in the condensed phase coordinate file.
NMol = len(ML.molecules)
#----
def indicate(self):
"""
print information into the output file about the last objective function evaluated
This function should be called after get()
"""
for property_name, details in self._last_obj_details.items():
dict_for_print = {
" %sK %satm %s" % detail[:3]: "%9.3f %14.3f +- %-7.3f % 7.3f % 9.5f % 9.5f % 9.5f " % detail[3:] for
detail in details}
title = '%s %s\nTemperature Pressure Substance Reference Calculated +- Stdev Delta Weight ' \
'Denom Term ' % (self.name, property_name)
printcool_dictionary(dict_for_print, title=title, bold=True, color=4, keywidth=15)
for i in ['temperature', 'pressure', 'nequil', 'nsteps', 'timestep', 'sample', 'threads', 'minimize']:
if i in args:
MDOpts[i] = args[i]
MDOpts['nsave'] = int(1000.0*MDOpts['sample']/MDOpts['timestep'])
if 'save_traj' in TgtOpts:
MDOpts['save_traj'] = TgtOpts['save_traj']
#----
# Print some options.
# At this point, engine and MD options should be SET!
#----
printcool("ForceBalance simulation using engine: %s" % engname.upper(),
color=4, bold=True)
printcool_dictionary(args, title="Options from command line")
printcool_dictionary(EngOpts, title="Engine options")
printcool_dictionary(MDOpts, title="Molecular dynamics options")
#----
# For convenience, assign some local variables.
#----
# Finite difference step size
h = TgtOpts['h']
# Active parameters to differentiate
pgrad = TgtOpts['pgrad']
# Create instances of the MD Engine objects.
Engine = EngineClass(**EngOpts)
click() # Start timer.
# This line runs the condensed phase simulation.
#----
# The molecular dynamics simulation returns a dictionary of properties
# In the future, the properties will be stored as data inside the object
Results = Engine.molecular_dynamics(**MDOpts)
# Read the command line options (they may override the options from file.)
AGrad = args['gradient'] or Sim.gradient
for i in ['temperature', 'pressure', 'nequil', 'nsteps', 'timestep', 'sample', 'threads', 'minimize']:
if i in args:
Sim.MDOpts[i] = args[i]
#----
# Print some options.
# At this point, engine and MD options should be SET!
#----
printcool("ForceBalance simulation using engine: %s" % engname.upper(),
color=4, bold=True)
printcool_dictionary(args, title="Options from command line")
printcool_dictionary(Sim.EngOpts, title="Engine options")
printcool_dictionary(Sim.MDOpts, title="Molecular dynamics options")
#----
# For convenience, assign some local variables.
#----
# Finite difference step size
h = Sim.h
# Active parameters to differentiate
pgrad = Sim.pgrad
# Create instances of the MD Engine objects.
Engine = EngineClass(name=Sim.type, **Sim.EngOpts)
click() # Start timer.
# This line runs the condensed phase simulation.
Engine.molecular_dynamics(**Sim.MDOpts)
logger.info("MD simulation took %.3f seconds\n" % click())
# Extract properties.
Results = Engine.md_extract(OrderedDict([(i, {}) for i in Sim.timeseries.keys()]))