Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
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.
ndof = 3*system.getNumParticles() - system.getNumConstraints() - 3
# Compute total mass.
mass = compute_mass(system).in_units_of(gram / mole) / AVOGADRO_CONSTANT_NA # total system mass in g
# Initialize statistics.
data = dict()
GD2 += 2*(flat(np.mat(GDz)*col(Dz))/N - avg(Dz)*(np.mean(GDz,axis=1))) - Beta*(covde(Dz**2) - 2*avg(Dz)*covde(Dz))
GEps0 = prefactor*(GD2/avg(V) - mBeta*covde(V)*D2/avg(V)**2)/T
Sep = printcool("Dielectric constant: % .4e +- %.4e\nAnalytic Derivative:" % (Eps0, Eps0_err))
FF.print_map(vals=GEps0)
if FDCheck:
GEps0_fd = property_derivatives(mvals, h, pdb, FF, Xyzs, Settings, Sim, kT, calc_eps0, {'d_':Dips,'v_':V}, Boxes)
Sep = printcool("Numerical Derivative:")
FF.print_map(vals=GEps0_fd)
Sep = printcool("Difference (Absolute, Fractional):")
absfrac = ["% .4e % .4e" % (i-j, (i-j)/j) for i,j in zip(GEps0,GEps0_fd)]
FF.print_map(vals=absfrac)
## Print the final force field.
pvals = FF.make(mvals)
with wopen(os.path.join('npt_result.p')) as f: lp_dump((Rhos, Volumes, Energies, Dips, G, [GDx, GDy, GDz], mEnergies, mG, Rho_err, Hvap_err, Alpha_err, Kappa_err, Cp_err, Eps0_err, NMol),f)
#throw error if more than one script
#write script into .txt file and parse as text
fflist = list(self.ffdata[ffname].iter())
scriptElements = [elem for elem in fflist if elem.tag=='Script']
if len(scriptElements) > 1:
logger.error('XML file'+ffname+'contains more than one script! Consolidate your scripts into one script!\n')
raise RuntimeError
elif len(scriptElements)==1:
Script = scriptElements[0].text
ffnameList = ffname.split('.')
ffnameScript = ffnameList[0]+'Script.txt'
absScript = os.path.join(self.root, self.ffdir, ffnameScript)
if os.path.exists(absScript):
logger.error('XML file '+absScript+' already exists on disk! Please delete it\n')
raise RuntimeError
wfile = forcebalance.nifty.wopen(absScript)
wfile.write(Script)
wfile.close()
self.addff(ffnameScript, xmlScript=True)
os.unlink(absScript)
for e in self.ffdata[ffname].getroot().xpath('//@parameterize/..'):
parameters_to_optimize = [i.strip() for i in e.get('parameterize').split(',')]
for p in parameters_to_optimize:
if p not in e.attrib:
logger.error("Parameter \'%s\' is not found for \'%s\', please check %s" % (p, e.get('type'), ffname) )
raise RuntimeError
pid = self.Readers[ffname].build_pid(e, p)
self.map[pid] = self.np
# offxml file later than v0.3 may have unit strings in the field
quantity_str = e.get(p)
res = re.search(r'^[-+]?[0-9]*\.?[0-9]*([eEdD][-+]?[0-9]+)?', quantity_str)
line_out.append(line)
# Sanity checks: If frcmod and mol2 files are provided to this function,
# they should be in the leap.cmd file as well. There should be exactly
# one PDB file being loaded.
for i in frcmod:
if i not in have_fmod:
warn_press_key("WARNING: %s is not being loaded in %s" % (i, fnm))
for i in mol2:
if i not in have_mol2:
warn_press_key("WARNING: %s is not being loaded in %s" % (i, fnm))
fout = fnm+'_'
line_out.append('saveamberparm %s %s.prmtop %s.inpcrd\n' % (ambername, prefix, prefix))
line_out.append('quit\n')
if os.path.exists(fout): os.remove(fout)
with wopen(fout) as f: print(''.join(line_out), file=f)
GD2 += 2*(flat(np.mat(GDz)*col(Dz))/N - avg(Dz)*(np.mean(GDz,axis=1))) - Beta*(covde(Dz**2) - 2*avg(Dz)*covde(Dz))
GEps0 = prefactor*(GD2/avg(V) - mBeta*covde(V)*D2/avg(V)**2)/T
Sep = printcool("Dielectric constant: % .4e +- %.4e\nAnalytic Derivative:" % (Eps0, Eps0_err))
FF.print_map(vals=GEps0)
if FDCheck:
GEps0_fd = property_derivatives(mvals, h, FF, args.liquid_xyzfile, args.liquid_keyfile, kT, calc_eps0, {'d_':Dips,'v_':V})
Sep = printcool("Numerical Derivative:")
FF.print_map(vals=GEps0_fd)
Sep = printcool("Difference (Absolute, Fractional):")
absfrac = ["% .4e % .4e" % (i-j, (i-j)/j) for i,j in zip(GEps0,GEps0_fd)]
FF.print_map(vals=absfrac)
## Print the final force field.
pvals = FF.make(mvals)
with wopen(os.path.join('npt_result.p')) as f: lp_dump((Rhos, Volumes, Potentials, Energies, Dips, G, [GDx, GDy, GDz], mPotentials, mEnergies, mG, Rho_err, Hvap_err, Alpha_err, Kappa_err, Cp_err, Eps0_err, NMol),f)
# Average over all points
Objective += np.sum(Objs, axis=0)
Gradient += np.sum(Grads, axis=0)
Hessian += np.sum(Hess, axis=0)
# Store gradients and setup print map
GradMapPrint = [["#Point"] + self.FF.plist]
for pt in self.points:
temp = pt.temperature
press = pt.pressure
GradMapPrint.append([' %8.2f %8.1f' % (temp, press)] +
["% 9.3e" % i for i in grads[pt.idnr-1]])
o = wopen('gradient_%s.dat' % quantity)
for line in GradMapPrint:
print >> o, ' '.join(line)
o.close()
printer = OrderedDict([(" %-5d %-12.2f %-8.1f"
% (pt.idnr, pt.temperature, pt.pressure),
("% -10.3f % -10.3f +- %-8.3f % -8.3f % -9.5f % -9.5f"
% (Exp[pt.idnr-1], values[pt.idnr-1],
errors[pt.idnr-1], Deltas[pt.idnr-1],
Weights[pt.idnr-1], Objs[pt.idnr-1])))
for pt in self.points])
return { "X": Objective, "G": Gradient, "H": Hessian, "info": printer }
# 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()]))
# Calculate energy and dipole derivatives if needed.
if AGrad:
Results['derivatives'] = energy_derivatives(Engine, FF, mvals, h, pgrad, dipole='dipole' in Sim.timeseries.keys())
# Dump results to file
logger.info("Writing final force field.\n")
pvals = FF.make(mvals)
logger.info("Writing all simulation data to disk.\n")
with wopen('md_result.p') as f:
lp_dump(Results, f)
def writechk(self):
""" Write the checkpoint file for the main optimizer. """
if self.wchk_fnm is not None:
logger.info("Writing the checkpoint file %s\n" % self.wchk_fnm)
with wopen(os.path.join(self.root,self.wchk_fnm), binary=True) as f: pickle.dump(self.chk, f)