# Copyright (c) 2020-2021 The Center for Theoretical Biological Physics (CTBP) - Rice University and Northeastern University
# This file is from the OpenSMOG project, released under the MIT License.
R"""
The :class:`~.OpenSMOG` classes perform molecular dynamics using Structure-Based Models (SBM) for biomolecular simulations.
:class:`~.OpenSMOG` uses force fields generated by SMOG 2, and it allows the simulations of a wide variety of potential forms, including commonly employed C-alpha and all-atom variants.
Details about the default models in SMOG 2 can be found in the following resources:
- **SMOG server**: https://smog-server.org/smog2/
- **C-alpha**: Clementi, C., Nymeyer, H. and Onuchic, J.N., 2000. Topological and energetic factors: what determines the structural details of the transition state ensemble and “en-route” intermediates for protein folding? An investigation for small globular proteins. Journal of molecular biology, 298(5), pp.937-953.
- **All-Atom**: Whitford, P.C., Noel, J.K., Gosavi, S., Schug, A., Sanbonmatsu, K.Y. and Onuchic, J.N., 2009. An all‐atom structure‐based potential for proteins: bridging minimal models with all‐atom empirical forcefields. Proteins: Structure, Function, and Bioinformatics, 75(2), pp.430-441.
"""
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
import os
import numpy as np
import xml.etree.ElementTree as ET
import warnings
from lxml import etree
import sys
#from .OpenSMOG_Reporter import forcesReporter
[docs]class SBM:
R"""
The :class:`~.SBM` class performs Molecular dynamics simulations using structure-based (SMOG) models to investigate a broad range of biomolecular dynamics, including domain rearrangements in proteins, folding and ligand binding in RNA, and large-scale rearrangements in ribonucleoprotein assemblies. In its simplest form, a structure-based model defines a particular structure (usually obtained from X-ray, or NMR, methods) as the energetic global minimum. Find more information about SMOG models and OpenSMOG at http://smog-server.org
The :class:`~.SBM` sets the environment to start the molecular dynamics simulations.
Args:
time_step (float, required):
Simulation time step in units of :math:`\tau`.
collision_rate (float, required):
Friction/Damping constant in units of reciprocal time (:math:`1/\tau`).
r_cutoff (float, required):
Cutoff distance to consider non-bonded interactions in units of nanometers.
pbc (boolean, optional):
Turn PBC on/off. Default value: False
temperature (float, required):
Temperature in reduced units.
name (str):
Name used in the output files. (Default value: :code:`OpenSMOG`).
"""
def __init__(self, time_step, collision_rate, r_cutoff, temperature, pbc = False, name = "OpenSMOG", warn = True):
self.printHeader()
print("\nFor more information, including descriptions of units and an example for how to launch a\nsimulation with OpenSMOG, run the help module. For example, if you named your simulation\nobject \"SMOGMODEL\", then issue the command:\n>SMOGMODEL.help()\n")
self.name = name
self.warn = warn
self.dt = time_step * picoseconds
self.gamma = collision_rate / picosecond
self.rcutoff = r_cutoff * nanometers
if self.warn:
if not time_step in [0.0005, 0.002]:
print('Note: The given time_step value is not the one usually employed in the SBM models. Make sure this value is correct. The suggested values are: time_step=0.0005 for C-alpha and time_step = 0.002 for All-Atoms.')
if collision_rate != 1.0:
print('Note: The given collision_rate value is not the one usually employed in the SBM models. Make sure this value is correct. The suggested value is: collision_rate=1.0.')
if not r_cutoff in [1.1 ,0.65]:
print('Note: The given r_cutoff value is not the one usually employed in the SBM models with OpenSMOG. Make sure this value is correct. The suggested values for r_cutoff are: 1.1 for the default C-alpha model and 0.65 for the all-atom model.')
self.temperature = (temperature / 0.00831446261815) * kelvin
self.forceApplied = False
self.loaded = False
self.folder = "."
self.forceCount = 0
self.pbc=pbc
self.nonbonded_present=False
[docs] def help(self):
R"""Prints information about using OpenSMOG.
"""
print("""
#############################################################################
Information about using OpenSMOG
#############################################################################
Units in OpenSMOG
Time: When using OpenSMOG, it is important to consider the units. Since
OpenSMOG is built on OpenMM libraries, some unit labels can be misleading.
Specifically, SMOG models use reduced units. While it is possible to estimate
the effective timescales in these models, OpenMM reporters will label the
timescales as ns. In terms of accounting, 1 \"ns\" is equal to 1000 reduced
time units. For detailed discussions of units, one can consult any number of
recent papers by the SMOG team, such as Yang et al. JCP 2019, 151, 085102.
Temperature: When initializing your simulation with the SBM class, the
temperature is provided in reduced units, where a typical protein will fold
at around 1 reduced unit. However, similar to other MD packages, the
Boltzmann constant is built into the software. Accordingly, the reported
temperature will be equal to (reduced temperature)*120.272. So, if you
initialize your system with a temperature of 0.5, the reported temperature
would read 0.5*120.272~60.136.
Energy: OpenSMOG will report all energies in reduced units. However, the
OpenMM libraries will label the energies as units of kJ/mole. This label
should be ignored, unless you calibrated your SMOG model to use this energy
scale.
Usage Example
While a broad range of simulation protocols may be used, here is a very basic
example for how to launch simulations using the OpenSMOG module with OpenMM.
>from OpenSMOG import SBM
Choose some basic runtime settings. We will call our system 2ci2
>SMOGrun = SBM(name='2ci2', time_step=0.002, collision_rate=1.0, r_cutoff=1.2, temperature=0.5)
Select a platform and GPU IDs (if needed)
>SMOGrun.setup_openmm(platform='cuda',GPUindex='default')
Decide where to save your data (here, output_2ci2
>SMOGrun.saveFolder('output_2ci2')
You may optionally save some input file names to variables
>SMOG_grofile = '2ci2.gro'
>SMOG_topfile = '2ci2.top'
>SMOG_xmlfile = '2ci2.xml'
Load your force field data
>SMOGrun.loadSystem(Grofile=SMOG_grofile, Topfile=SMOG_topfile, Xmlfile=SMOG_xmlfile)
Create the simulation, and prepare to run
>SMOGrun.createSimulation()
Decide how frequently to save data
>SMOGrun.createReporters(trajectory=True, energies=True, energy_components=True, interval=10**3)
Launch the simulation
>SMOGrun.run(nsteps=10**6, report=True, interval=10**3)
For more information and help, see the OpenSMOG and SMOG 2 websites.
OpenSMOG: https://opensmog.readthedocs.io/en/latest/
SMOG 2: https://smog-server.org
If you have questions/suggestions, you can also email us at info@smog-server.org
""")
[docs] def minimize(self,tolerance=1.0,maxIterations=0):
R"""Wrapper for minimization simulation.
Args:
tolerance (float, required):
Stopping criteria value between iteration. When the error between iteration is below this value, the minimization stops. (Default value: :code:`1.0`).
maxIteration (int, required):
Number of maximum steps to be performed in the minimization simulation. (Default value: :code:`1.0`).
"""
self.simulation.minimizeEnergy(tolerance=tolerance,maxIterations=maxIterations)
[docs] def setup_openmm(self, platform='opencl', precision='single', GPUindex='default', integrator="langevin"):
R"""Sets up the parameters of the simulation OpenMM platform.
Args:
platform (str, optional):
Platform to use in the simulations. Opitions are *CUDA*, *OpenCL*, *HIP*, *CPU*, *Reference*. (Default value: :code:`OpenCL`).
precision (str, optional):
Numerical precision type of the simulation. Options are *single*, *mixed*, *double*. (Default value: :code:`single`). For details check the `OpenMM Documentation <http://docs.openmm.org/latest/developerguide/developer.html#numerical-precision>`__.
GPUIndex (str, optional):
Set of Platform device index IDs. Ex: 0,1,2 for the system to use the devices 0, 1 and 2. (Use only when GPU != default).
integrator (str):
Integrator to use in the simulations. Options are *langevin*, *variableLangevin*, *verlet*, *variableVerlet* and, *brownian*. (Default value: :code:`langevin`).
"""
precision = precision.lower()
if precision not in ["mixed", "single", "double"]:
raise ValueError("Precision must be mixed, single or double")
properties = {}
properties["Precision"] = precision
if GPUindex.lower() != "default":
properties["DeviceIndex"] = GPUindex
self.properties = properties
if platform.lower() == "opencl":
platformObject = Platform.getPlatformByName('OpenCL')
elif platform.lower() == "reference":
platformObject = Platform.getPlatformByName('Reference')
self.properties = {}
elif platform.lower() == "cuda":
platformObject = Platform.getPlatformByName('CUDA')
elif platform.lower() == "cpu":
platformObject = Platform.getPlatformByName('CPU')
self.properties = {}
elif platform.lower() == "hip":
platformObject = Platform.getPlatformByName('HIP')
else:
self.exit("\n!!!! Unknown platform !!!!\n")
self.platform = platformObject
if integrator.lower() == "langevin":
self.integrator = LangevinIntegrator(self.temperature,
self.gamma, self.dt)
self.integrator_type = integrator
else:
self.integrator = integrator
self.integrator_type = "UserDefined"
self.forceDict = {}
self.forcesDict = {}
[docs] def saveFolder(self, folder):
R"""Sets the folder path to save data.
Args:
folder (str, optional):
Folder path to save the simulation data. If the folder path does not exist, the function will create the directory.
"""
if os.path.exists(folder) == False:
os.mkdir(folder)
self.folder = folder
[docs] def loadSystem(self, Grofile, Topfile, Xmlfile):
R"""Loads the input files in the OpenMM system platform. The input files are generated using SMOG2 software with the flag :code:`-OpenSMOG`. Details on how to create the files can be found in the `SMOG2 User Manual <https://smog-server.org/smog2/>`__.
A tutorial on how to generate the inputs files for default all-atom and C-alpha models can be found `here <https://opensmog.readthedocs.io>`__.
Args:
Grofile (file, required):
Initial structure for the MD simulations in *.gro* file format generated by SMOG2 software with the flag :code:`-OpenSMOG`. (Default value: :code:`None`).
Topfile (file, required):
Topology *.top* file format generated by SMOG2 software with the flag :code:`-OpenSMOG`. The topology file lists the interactions between the system atoms except for the "Native Contacts" potential that is provided to OpenSMOG in a *.xml* file. (Default value: :code:`None`).
Xmlfile (file, required):
The *.xml* file that contains the all information that defines the "Contact" potential. The *.xml* file is generated by SMOG2 software with the flag :code:`-OpenSMOG`, which support custom potential functions. (Default value: :code:`None`).
"""
def _checknames(f1,f2,f3):
fn1 = os.path.basename(f1).rsplit('.', 1)[0]
fn2 = os.path.basename(f2).rsplit('.', 1)[0]
fn3 = os.path.basename(f3).rsplit('.', 1)[0]
if fn1 == fn2 and fn1 ==fn3:
return False
else:
return True
if _checknames(Grofile, Topfile, Xmlfile) and self.warn :
warnings.warn('The Gro, Top and Xml files have different prefixes. Most people use the same name, so this may be a mistake.')
self.inputNames = [Grofile, Topfile, Xmlfile]
self._check_file(Grofile, '.gro')
self.loadGro(Grofile)
self._check_file(Topfile, '.top')
self.loadTop(Topfile)
self._check_file(Xmlfile, '.xml')
self.loadXml(Xmlfile)
print("Loaded force field and config files.")
def _check_file(self, filename, ext):
if not (filename.lower().endswith(ext)):
raise ValueError('Wrong file extension: {} must to be {} extension'.format(filename, ext))
[docs] def loadGro(self, Grofile):
R"""Loads the *.gro* file format in the OpenMM system platform. The input files are generated using SMOG2 software with the flag :code:`-OpenSMOG`. Details on how to create the files can be found in the `SMOG2 User Manual <https://smog-server.org/smog2/>`__.
A tutorial on how to generate the inputs files for the default all-atom and C-alpha models can be found `here <https://opensmog.readthedocs.io>`__.
Args:
Grofile (file, required):
Initial structure for the MD simulations in *.gro* file format generated by SMOG2 software with the flag :code:`-OpenSMOG`. (Default value: :code:`None`).
"""
self.Gro = GromacsGroFile(Grofile)
[docs] def loadTop(self, Topfile):
R"""Loads the *.top* file format in the OpenMM system platform. The input files are generated using SMOG2 software with the flag :code:`-OpenSMOG`. Details on how to create the files can be found in the `SMOG2 User Manual <https://smog-server.org/smog2/>`__.
A tutorial on how to generate the inputs files for the default all-atom and C-alpha models can be found `here <https://opensmog.readthedocs.io>`__.
Args:
Topfile (file, required):
Topology *.top* file format generated by SMOG2 software with the flag :code:`-OpenSMOG`. The topology file defines the interactions between atoms, except for the "Native Contacts" potential that is provided to OpenSMOG in the form of a *.xml* file. (Default value: :code:`None`).
"""
if self.pbc == True:
print('This simulation will use Periodic boundary conditions')
self.Top = GromacsTopFile(Topfile,unitCellDimensions=self.Gro.getUnitCellDimensions())
self.system = self.Top.createSystem(nonbondedMethod=CutoffPeriodic,nonbondedCutoff=self.rcutoff)
else:
print('This simulation will not use Periodic boundary conditions')
self.Top = GromacsTopFile(Topfile)
self.system = self.Top.createSystem(nonbondedMethod=CutoffNonPeriodic,nonbondedCutoff=self.rcutoff)
nforces = len(self.system.getForces())
for force_id, force in enumerate(self.system.getForces()):
force.setForceGroup(force_id)
self.forcesDict[force.__class__.__name__] = force
self.forceCount +=1
def _splitForces_contacts(self):
#Contacts
cont_data=self.data['contacts']
n_forces = len(cont_data[3])
forces = {}
for n in range(n_forces):
forces[cont_data[3][n]] = [cont_data[0][n], cont_data[1][n], cont_data[2][n]]
self.contacts = forces
def _splitForces_nb(self):
#Contacts
nb_data=self.data['nonbond']
n_forces = len(nb_data[0])
forces = {}
for n in range(n_forces):
forces[nb_data[0][n]] = [nb_data[1][n], nb_data[2][n], nb_data[3][n]]
self.nonbond = forces
def _customSmogForce(self, name, data):
#first set the equation
contacts_ff = CustomBondForce(data[0])
#second set the number of variable
for pars in data[1]:
contacts_ff.addPerBondParameter(pars)
#third, apply the bonds from each pair of atoms and the related variables.
pars = [pars for pars in data[1]]
for iteraction in data[2]:
atom_index_i = int(iteraction['i'])-1
atom_index_j = int(iteraction['j'])-1
parameters = [float(iteraction[k]) for k in pars]
contacts_ff.addBond(atom_index_i, atom_index_j, parameters)
#forth, if the are global variables, add them to the force
if self.constants_present==True:
for const_key in self.data['constants']:
contacts_ff.addGlobalParameter(const_key,self.data['constants'][const_key])
self.forcesDict[name] = contacts_ff
contacts_ff.setForceGroup(self.forceCount)
self.forceCount +=1
def _customSmogForce_nb(self, name, data):
#first set the equation
nonbond_ff = CustomNonbondedForce(data[0])
#Define per particle parameters
nonbond_ff.addPerParticleParameter('q')
nonbond_ff.addPerParticleParameter('type')
#If the are global variables, add them to the force
if self.constants_present==True:
for const_key in self.data['constants']:
nonbond_ff.addGlobalParameter(const_key,self.data['constants'][const_key])
#Add cutoff as global parameter
nonbond_ff.addGlobalParameter('r_c',self.rcutoff.value_in_unit(nanometer))
#Load old nonbonded force for later use
original_nonbonded=self.system.getForce(0)
original_customnonbonded=self.system.getForce(1)
#Get atoms types and number of types
atom_types=[]
for i in range(len(data[2])):
atom_types.append(data[2][i]['type1'])
self.atom_types=np.unique(np.array(atom_types))
natom_types=len(self.atom_types)
#Generate tables for each parameter
tables={}
for par in data[1]:
tables[par]=np.ones([natom_types,natom_types])*np.nan
#Fill tables with nonbond_param
for nonbond_params in data[2]:
#Get atoms id (name to number)
type1=np.where(self.atom_types==nonbond_params['type1'])[0][0]
type2=np.where(self.atom_types==nonbond_params['type2'])[0][0]
for par in tables:
tables[par][type1][type2]=nonbond_params[par]
tables[par][type2][type1]=nonbond_params[par]
#Generate Function from tables
table_fun={}
for par in data[1]:
#Check none have nans
if np.sum(tables[par]==np.nan)!=0:
raise ValueError('Nonbonded force parameter: {:} is not defined for all atom interactions'.format(par))
table_fun[par]=Discrete2DFunction(natom_types,natom_types,list(np.ravel(tables[par])))
nonbond_ff.addTabulatedFunction(par,table_fun[par])
#Get exceptions from topfile import
for i in range(original_customnonbonded.getNumExclusions()):
exclusion_id = original_customnonbonded.getExclusionParticles(i)
nonbond_ff.addExclusion(exclusion_id[0],exclusion_id[1])
## ADD PARTICLES TO THE FORCE BASED ON SPECIFYING TYPE AND CHARGE
## CHARGES
## From nonbondedforce we get the charges for each atom
atom_charges=[]
## Loop over every atom
for i in range(self.system.getNumParticles()):
atom_charge=original_nonbonded.getParticleParameters(i)[0]
atom_charges.append(atom_charge.value_in_unit(constants.elementary_charge))
## ATOM TYPES
## From molecule information we get atom types
atom_types=[]
## Get name and multiplicity of each molecule
molecules_keys=[self.Top._molecules[i][0] for i in range(len(self.Top._molecules)) ]
molecules_mul=[self.Top._molecules[i][1] for i in range(len(self.Top._molecules)) ]
## Loop over molecules
for i in range(len(molecules_keys)):
molecule=molecules_keys[i]
mult=molecules_mul[i]
# Loop over atoms in each molecule
for atom in self.Top._moleculeTypes[molecule].atoms:
# If atoms are repeated, then include the same entrie (CL,CL,CL,... or MG,MG,MG,...)
for _ in range(mult):
atom_types.append(atom[1])
for i in range(self.system.getNumParticles()):
# GET ATOM TYPE
at_type=np.where(atom_types[i]==self.atom_types)[0][0]
# GET ATOM CHARGE
charge=atom_charges[i]
# ADD PARTICLE TO EACH FORCE WITH CORRESPONDING CHARGE AND TYPE
nonbond_ff.addParticle([charge,at_type])
#Set cutoff and nonbonded method
if self.pbc == True:
nonbond_ff.setNonbondedMethod(NonbondedForce.CutoffPeriodic)
else:
nonbond_ff.setNonbondedMethod(NonbondedForce.CutoffNonPeriodic)
nonbond_ff.setCutoffDistance(self.rcutoff.value_in_unit(nanometer))
self.forcesDict['Nonbonded'+str(name)] = nonbond_ff
nonbond_ff.setForceGroup(self.forceCount)
self.forceCount +=1
[docs] def loadXml(self, Xmlfile):
R"""Loads the *.xml* file format in the OpenMM system platform. The input files are generated using SMOG2 software with the flag :code:`-OpenSMOG`. Details on how to create the files can be found in the `SMOG2 User Manual <https://smog-server.org/smog2/>`__.
A tutorial on how to generate the inputs files for default all-atom and C-alpha models can be found `here <https://opensmog.readthedocs.io>`__.
Args:
Xmlfile (file, required):
The *.xml* file that contains all information that defines the "Contact" potentials. The *.xml* file is generated by SMOG2 software with the flag :code:`-OpenSMOG`, which support custom potential energy functions. (Default value: :code:`None`).
"""
def validate(Xmlfile):
path = "share/OpenSMOG.xsd"
pt = os.path.dirname(os.path.realpath(__file__))
filepath = os.path.join(pt,path)
xmlschema_doc = etree.parse(filepath)
xmlschema = etree.XMLSchema(xmlschema_doc)
xml_doc = etree.parse(Xmlfile)
result = xmlschema.validate(xml_doc)
log = xmlschema.error_log
return result,log
def import_xml2OpenSMOG(file_xml):
XML_potential = ET.parse(file_xml)
root = XML_potential.getroot()
xml_data={}
## Constants
self.constants_present=False
if root.find('constants') != None:
self.constants_present=True
constants={}
constants_xml = root.find('constants')
for const in constants_xml.iter('constant'):
constants[const.attrib['name']]=float(const.attrib['value'])
xml_data['constants']=constants
## Contacts
Force_Names=[]
Expression=[]
Parameters=[]
Pairs=[]
if root.find('contacts') == None:
raise ValueError("No contacts were found in the XML file")
contacts_xml=root.find('contacts')
for i in range(len(contacts_xml)):
for name in contacts_xml[i].iter('contacts_type'):
Force_Names.append(name.attrib['name'])
for expr in contacts_xml[i].iter('expression'):
Expression.append(expr.attrib['expr'])
internal_Param=[]
for par in contacts_xml[i].iter('parameter'):
internal_Param.append(par.text)
Parameters.append(internal_Param)
internal_Pairs=[]
for atompairs in contacts_xml[i].iter('interaction'):
internal_Pairs.append(atompairs.attrib)
Pairs.append(internal_Pairs)
xml_data['contacts']=[Expression,Parameters,Pairs,Force_Names]
#Launch contact force function
self.nonbond_present=False
if root.find('nonbond') != None:
self.nonbond_present=True
nonbond_xml=root.find('nonbond')
NonBond_Num=[]
NBExpression=[]
NBExpressionParameters=[]
NBParameters=[]
nonbond_xml=root.find('nonbond')
for i in range(len(nonbond_xml)):
for types in nonbond_xml[i].iter('nonbond_bytype'):
NonBond_Num.append(i)
for expr in nonbond_xml[i].iter('expression'):
NBExpression.append(expr.attrib['expr'])
internal_Param=[]
for par in nonbond_xml[i].iter('parameter'):
internal_Param.append(par.text)
NBExpressionParameters.append(internal_Param)
internal_NBParam=[]
for nbpar in nonbond_xml[i].iter('nonbond_param'):
internal_NBParam.append(nbpar.attrib)
NBParameters.append(internal_NBParam)
xml_data['nonbond']=[NonBond_Num,NBExpression,NBExpressionParameters,NBParameters]
else:
self.nonbond_present=False
return xml_data
if not (self.forceApplied):
result,log=validate(Xmlfile)
if not result:
raise ValueError("The xml file is not adhere to the required schema (same as that used by SMOG 2). Error message:\n\n"+str(log)+"\n\n This typically means your xml file was corrupted, or you are using an incompatible version of SMOG 2")
self.data = import_xml2OpenSMOG(Xmlfile)
self._splitForces_contacts()
for force in self.contacts:
print("Creating Contacts force {:} from xml file".format(force))
self._customSmogForce(force, self.contacts[force])
self.system.addForce(self.forcesDict[force])
if self.nonbond_present==True:
self._splitForces_nb()
for force in self.nonbond:
print("Creating Nonbonded force {:} from xml file".format(force))
self._customSmogForce_nb(force, self.nonbond[force])
self.system.addForce(self.forcesDict['Nonbonded'+str(force)])
## REMOVE OTHER NONBONDED FORCES
self.system.removeForce(0)
self.system.removeForce(0)
self.forceApplied = True
else:
print('\n OpenSMOG forces already applied!!! \n')
[docs] def createSimulation(self):
R"""Creates the simulation context and loads into the OpenMM platform.
"""
if not self.loaded:
self.simulation = Simulation(self.Top.topology, self.system, self.integrator, self.platform, self.properties)
self.simulation.context.setPositions(self.Gro.positions)
self.simulation.context.setVelocitiesToTemperature(self.temperature)
self.loaded = True
else:
print('\n Simulations context already created! \n')
def _checkFile(self,filename):
if os.path.isfile(filename):
i = 1
ck = True
while i <= 10 and ck:
newname = filename + "_" + str(i)
if not os.path.isfile(newname):
print("{:} already exists. Backing up to {:}".format(filename,newname))
os.rename(filename, newname)
ck = False
else:
i += 1
[docs] def createReporters(self, trajectory=True, trajectoryName=None, trajectoryFormat='dcd', energies=True, energiesName=None, energy_components=False, energy_componentsName=None, interval=1000):
R"""Creates the reporters to provide the output data.
Args:
trajectory (bool, optional):
Whether to save the trajectory *.dcd* file containing the position of the atoms as a function of time. (Default value: :code:`True`).
trajectoryName (str, optional):
Name of the trajectory file.
trajectoryFormat (str, optional):
File format of the trajectory file. Options are dcd, pdb, pdbx, hdf5, netcdf, and xtc. Saving the trajectory in the file formats hdf5, netcdf, and xtc require the MDTraj package. (Default value: :code:`dcd`).
energies (bool, optional):
Whether to save the energies in a *.txt* file containing five columns, comma-delimited. The header of the files shows the information of each collum: #"Step","Potential Energy (kJ/mole)","Kinetic Energy (kJ/mole)","Total Energy (kJ/mole)","Temperature (K)". (Default value: :code:`True`).
energiesName (str, optional):
Name of the energy file.
energy_components (bool, optional):
Whether to save the potential energy for each applied force in a *.txt* file containing several columns, comma-delimited. The header of the files shows the information of each column. An example of the header is: #"Step","electrostatic","Non-Contacts","Bonds","Angles","Dihedrals","contact_1-10-12". (Default value: :code:`False`).
energy_componentsName (str, optional):
Name of the energy_components file.
interval (int, required):
Frequency to write the data to the output files. (Default value: :code:`10**3`)
"""
self.outputNames = []
if trajectory:
trajectoryFormat=trajectoryFormat.lower()
if trajectoryName is None:
trajfile = os.path.join(self.folder, self.name + '_trajectory.'+trajectoryFormat)
else:
trajfile = os.path.join(self.folder, trajectoryName + "."+trajectoryFormat)
self._checkFile(trajfile)
self.outputNames.append(trajfile)
for i in ['hdf5', 'xtc', 'netcdf']:
if i == trajectoryFormat:
print("""
The """+str(trajectoryFormat)+""" trajectory format requires mdtraj to be loaded.
Will try to import mdtraj...""")
import mdtraj as md
if trajectoryFormat == 'dcd':
self.simulation.reporters.append(DCDReporter(trajfile, interval))
elif trajectoryFormat == 'pdb':
self.simulation.reporters.append(PDBReporter(trajfile, interval))
elif trajectoryFormat == 'pdbx':
self.simulation.reporters.append(PDBxReporter(trajfile, interval))
elif trajectoryFormat == 'hdf5':
self.simulation.reporters.append(md.reporters.HDF5Reporter(trajfile, interval))
elif trajectoryFormat == 'netcdf':
self.simulation.reporters.append(md.reporters.NetCDFReporter(trajfile, interval))
elif trajectoryFormat == 'xtc':
self.simulation.reporters.append(md.reporters.XTCReporter(trajfile, interval))
else:
raise ValueError("Trajectory format "+str(trajectoryFormat)+" not recognized")
if energies:
if energiesName is None:
energyfile = os.path.join(self.folder, self.name+ '_energies.txt')
else:
energyfile = os.path.join(self.folder, energiesName + ".txt")
self._checkFile(energyfile)
self.outputNames.append(energyfile)
self.simulation.reporters.append(StateDataReporter(energyfile, interval, step=True,
potentialEnergy=True, kineticEnergy=True,
totalEnergy=True,temperature=True, separator=","))
if energy_components:
if energy_componentsName is None:
forcefile = os.path.join(self.folder, self.name + '_forces.txt')
else:
forcefile = os.path.join(self.folder, energy_componentsName + '.txt')
self._checkFile(forcefile)
self.outputNames.append(forcefile)
self.simulation.reporters.append(forcesReporter(forcefile, interval, self.forcesDict, step=True))
[docs] def run(self, nsteps, report=True, interval=10**4):
R"""Run the molecular dynamics simulation.
Args:
nsteps (int, required):
Number of steps to be performed in the simulation. (Default value: :code:`10**7`)
report (bool, optional):
Whether to print the simulation progress. (Default value: :code:`True`).
interval (int, required):
Frequency to print the simulation progress. (Default value: :code:`10**4`)
"""
if report:
self.simulation.reporters.append(StateDataReporter(sys.stdout, interval, step=True, remainingTime=True,
progress=True, speed=True, totalSteps=nsteps, separator="\t"))
self._createLogfile()
self.simulation.step(nsteps)
def _createLogfile(self):
import platform
import datetime
logFilename = os.path.join(self.folder, 'OpenSMOG.log')
self._checkFile(logFilename)
self.outputNames.append(logFilename)
with open(logFilename, 'w') as f:
ori = sys.stdout
sys.stdout = f
self.printHeader()
sys.stdout = ori
#system_information
f.write('\nSystem Information:\n')
f.write('-------------------\n')
f.write('Date and time: {:}\n'.format(datetime.datetime.now()))
f.write('Machine information: {:} : {:}, {:} : {:}\n'.format("System", platform.uname().system, "Version", platform.uname().version))
f.write('Platform: {:}\n'.format(self.platform.getName()))
if (self.platform.getName() in ["CUDA", "OpenCL", "HIP"]):
f.write('Precision: {:}\n'.format(self.properties['Precision']))
f.write('Integrator: {:}\n'.format(self.integrator_type))
f.write('Savefolder: {:}\n'.format(self.folder))
f.write('\nSimulation Information:\n')
f.write('-----------------------\n')
f.write('Name: {:}\n'.format(self.name))
f.write('Time step: {:}\n'.format(self.dt/picoseconds))
f.write('Collision Rate: {:}\n'.format(self.gamma*picosecond))
f.write('r_Cutoff: {:}\n'.format(self.rcutoff/nanometers))
f.write('Temperature: {:}\n'.format(self.temperature * 0.008314/kelvin))
f.write('\nInput Files:\n')
f.write('------------\n')
f.write('GroFile: {:}\n'.format(self.inputNames[0]))
f.write('TopFile: {:}\n'.format(self.inputNames[1]))
f.write('XmlFile: {:}\n'.format(self.inputNames[2]))
f.write('\nOutput Files:\n')
f.write('-------------\n')
for n in self.outputNames:
f.write(n+"\n")
sys.stdout = ori