Narupa Engine
 All Classes Namespaces Functions Variables Enumerations Enumerator Properties Events Pages
User Guide

This page is a guide for running simulations with the Narupa engine.

Installing and Running Narupa

Prerequisites

Linux / OS X

Installation of the mono runtime is required, available for download here.

Running A Simulation

The Narupa engine consists of a single executable. Letting PATH_TO_NARUPA_ENGINE be the location of the engine, the executable is PATH_TO_NARUPA_ENGINE/bin/narupa-server.exe.

On Windows, the executable can be launched by double-clicking on it. On Linux or Mac OS X, it can be started by opening a terminal and running the following commands:

1 cd PATH_TO_NARUPA_ENGINE/bin
2 mono narupa-server.exe

The simulation to run can be configured by editing the ServerConfig.xml file located alongside the executable:

1 <?xml version="1.0" encoding="utf-8"?>
2 <ServerConfig>
3  <Server Name="Methane"
4  FilePath="./Assets/Simulations/TestSims/XML/Methane.xml"
5  OnConnectionState="Play"
6  Options="NoForces, NoCollisions, NoVelocities, VRServer"
7  ConnectionType="Tcp"
8  ServiceDiscoveryType="Ahoy"
9  />
10 </ServerConfig>

The FilePath field is used to specify the simulation to run (see Setting Up a Simulation). The ConnectionType field specifies to use TCP/IP to connect, and the ServiceDiscoveryType field specifies how the simulation will be found by clients. The default settings defined above are best used for running simulations on a local network.

Extending Functionality with Plugins.

You can extend Narupa engine functionality with plugins. We currently support an OpenMM plugin, which enables GPU-accelerated simulation of biomolecules.

Windows

Create a new folder called Plugins in Narupa Engine binary folder (typically PATH_TO_NARUPA_ENGINE/bin), and place plugins in separate folders within that directory. Note: On Windows, before extracting you may need to "unblock" the download link. To do this, right click on the zip file, click Properties and check the Unblock box at the bottom of the window.

Linux

Create a new folder called Plugins in Narupa Engine binary folder (typically PATH_TO_NARUPA_ENGINE/bin), and place plugins in separate folders within that directory. If native binaries are also supplied with the binary (e.g. C++ binaries, .so or .a files) place these in /usr/lib, run the command ldconfig or add the paths to your LD_LIBRARY_PATH (see here for more documentation).

Narupa Simulation Overview

The Narupa engine provides a set of APIs and programs for running simulations and visualising them. The engine described here is designed for molecular simulations, with a particular focus on interactive molecular dynamics.

The engine is designed to be run as a server, running as an independent program either on a local machine, a networked cluster or in the cloud. For visualisation, interaction, or analysis, client programs connect to the server and can interact with the simulation while it is running through communication streams and commands.

Narupa Engine

Narupa can run molecular simulations in several ways. It can either run simulations using the built-in forcefields and integrators, or it can use plugins to connect to 3rd-party simulation packages.

A Note On Units

The units in Narupa are designed to be internally consistent so that there do not have to be unit conversions all over the place. The following units should be used everywhere, whenever possible:

  • distance – nm
  • time – ps
  • mass – atomic mass units
  • charge – proton charge
  • temperature – Kelvin
  • angle – radians
  • energy – kJ/mol

Of course, specific applications may require conversion for display to end-users, and occasionally a force field may use different units (for an excruciating example, see the MM3 force field), but the internal system should endeavour to maintain these units.

Simulation Structure and Theory

The core of a simulation is composed of a surprisingly small number of structures. The most important are the AtomicSystem, the Topology, the Integrator and the Force Fields.

Narupa Engine Structure

Atomic System

The Atomic System (see AtomicSystem) is the container that holds the current atom structure and associated properties. The atomic system includes the following:

  • Particle Collection - The data on each atom, such as positions, velocities, forces, mass.
  • Topology - The connectivity and elements that make up the system.
  • Force Fields - The set of force fields acting on this simulation.
  • Integrator - The integrator that moves simulation forwards in time.
  • Interactions - The set of human interactions currently being applied to the system.
  • Simulation Boundary - The bounding box, or periodic boundary of the simulation.
  • Thermostat(s) - The control schemes used for maintaining temperature (or pressure etc.).

The atomic system is where most of the action happens in the simulation. When an integrator moves the simulation forward in time, it acts on the positions and velocities of the particles in the atomic system. Likewise, when a force field calculates the forces between atoms, it reads from the atomic system state and updates the force vector in the particle collection. Most operations within a simulation involve manipulating the atomic system in some way.

Topology

The Topology can be thought of as the overview of the molecule system. It knows how many particles there are, their elements and the bonds that connect them, but it leaves the details to the Atomic System and Particle Collection. It also contains the list of atomic OSC addresses, which allow for atoms to be selected/identified in a semi-human-readable fashion.

Integrator

A numerical integrator takes the Atomic System and uses the current set of positions, velocities and forces to move the simulation forward in time. Typically an integrator moves a simulation forward by a discrete step, such as one femtosecond. The speed at which a simulation can run is then limited by how fast an integration step can be computed and the largest time step that can be maintained without introducing instability.

Force Fields

A force field is a mathematical description of the way in which a set of atoms interact. If a simulation is not running any force fields, the atoms will simply move under Newtonian motion without interacting with one another. A force field takes the current state of the atomic system and calculates a set of forces based on the how it models the interactions between the atoms it governs. The set of forces are then turned into acceleration and ultimately an update in velocities and position by the integrator.

Interactions

When running a molecular dynamics simulation with the Narupa engine, it is easy to manipulate and guide the dynamics through interactions. Single atoms or whole groups of particles can be selected and moved around. The atoms are moved by applying a force, shaped like a normal distribution, which pulls the atoms towards the interaction site.

Simulation Boundary

The simulation boundary is the bounding box, or periodic boundary conditions, of the simulation. Atoms will either elastically bounce off the boundary, or - if periodic boundary conditions are enabled - wrap around to the other opposite side of the box.

A Typical Simulation Step

With all the components as detailed above, a typical simulation step proceeds as follows:

  • For each force field acting on the system, the integrator calculates updated forces based on current particle positions. These forces include both the internal forces between atoms and any external interactive forces.
  • The integrator updates particle velocities and positions using the newly calculated forces.
  • If thermostats/barostats are enabled, they will act to ensure that temperature/pressure are maintained at the desired values.
  • If the system is subject to constraints, a constraint solver is applied to ensure constraints are satisfied.

Setting Up a Simulation

This section describes the simulation input file and describes in detail the myriad of options available.

Overview of A Simulation File

The snippet below shows an example simulation input file for a toy simulation of two methane molecules and serves as an illustrative example. Paths are relative to the directoy in which the narupa-engine was launched. The ^ symbol has a special meaning, referring to the directory the narupa-engine executable is stored in. The ~ symbol can be used to access the home directory.

1 <?xml version="1.0" encoding="utf-8"?>
2 <Simulation Name="Methane">
3 
4  <SystemProperties>
5  <SimulationBoundary SimulationBox="1,1,1" PeriodicBoundary="false" />
6  <Thermostat Type="BerendsenThermostat" EquilibriumTemperature="390" MaximumTemperature="10000" BerendsenCoupling="0.008" />
7  <Integrator Type="VelocityVerlet" TimeStep="0.001" />
8  <Logging>
9  <Logger LogPath="^/Logs/Methane" Type="XYZ" Positions="true" Velocities="true" Forces="true" WriteFrequency="100" UseDateTimeStamp="true" />
10  <Logger Type="StateLogger" Time="true" PotentialEnergy="true" KineticEnergy="true" Temperature="true" WriteFrequency="10" UseDateTimeStamp="true" />
11  <Logger LogPath="^/Logs/Methane" Type="Interaction" WriteFrequency="100" UseDateTimeStamp="true" />
12  </Logging>
13  </SystemProperties>
14  <Topology>
15  <Templates>
16  <Residue Name="Methane">
17  <Atoms>
18  <Atom Element="Carbon" Position="0,0,0" />
19  <Atom Element="Hydrogen" Position="0.05288, 0.01610, 0.09359" />
20  <Atom Element="Hydrogen" Position="0.02051, 0.08240, -0.06786" />
21  <Atom Element="Hydrogen" Position="0.03345, -0.09314, -0.04496" />
22  <Atom Element="Hydrogen" Position="-0.10685, -0.00537, 0.01921" />>
23  </Atoms>
24  <Bonds>
25  <Bond A="1" B="0" />
26  <Bond A="1" B="2" />
27  <Bond A="1" B="3" />
28  <Bond A="1" B="4" />
29  </Bonds>
30  <ForceFields>
31  <InteractiveGaussianForceField GradientScaleFactor="1000">
32  <ParticleInteraction ScaleFactor="0" AtomPath="0" />
33  <ParticleInteraction ScaleFactor="0" AtomPath="2" />
34  <ParticleInteraction ScaleFactor="0" AtomPath="3" />
35  <ParticleInteraction ScaleFactor="0" AtomPath="4" />
36  </InteractiveGaussianForceField>
37  <MM3ForceField>
38  <MM3AtomMappings>
39  <MM3AtomMapping AtomPath="1" Type="1" />
40  <MM3AtomMapping AtomPath="0" Type="5" />
41  <MM3AtomMapping AtomPath="2" Type="5" />
42  <MM3AtomMapping AtomPath="3" Type="5" />
43  <MM3AtomMapping AtomPath="4" Type="5" />>
44  </MM3AtomMappings>
45  </MM3ForceField>
46  <LennardJonesForceField>
47  <LennardJonesAtomMappings>
48  <LennardJonesAtomMapping AtomPath="1" MM3Type="1" />
49  <LennardJonesAtomMapping AtomPath="0" MM3Type="5" />
50  <LennardJonesAtomMapping AtomPath="2" MM3Type="5" />
51  <LennardJonesAtomMapping AtomPath="3" MM3Type="5" />
52  <LennardJonesAtomMapping AtomPath="4" MM3Type="5" />
53  </LennardJonesAtomMappings>
54  </LennardJonesForceField>
55  </ForceFields>
56  </Residue>
57  </Templates>
58  <Spawners>
59  <Spawner Name="Methane" Template="Methane" Count="2" />
60  </Spawners>
61  </Topology>
62 </Simulation>

System Properties

The system properties are the set of options that will affect the whole molecular system. These include the choice of an integrator, simulation box dimensions and any thermostats.

Simulation Boundary

The simulation boundary defines the edges of the simulation. Currently, it must be included in the simulation properties. The simulation boundary chosen will affect how the molecules are spawned, so make sure the box is big enough to fit all the molecules in.

The following options are available:

  • SimulationBox - A 3d vector that results in a box centred at zero with the lengths specified. Not to be used in conjunction with SimulationBoxMin and SimulationBoxMax.
  • SimulationBoxMin - The position of the lower extremity of the box. To be used in conjunction with SimulationBoxMax.
  • SimulationBoxMax - The position of the upper extremity of the box.
  • PeriodicBoundary - If set to true, the simulation will use periodic boundary conditions, otherwise particles will elastically collide with the boundary.
  • IsVariableBox - If set to true, the bounding box can be resized at runtime through user input, otherwise, the simulation boundary can only be changed by the simulation itself.
  • MinimumBoxVolume - If the IsVariableBox has been set, a minimum box volume should be specified to prevent overlap and instability.

Thermostat

Adding a thermostat to the system properties means that the simulation will try to maintain a constant temperature (otherwise known as the NVT ensemble). For most simulations in Narupa, it makes sense to have a thermostat as the system is subject to external forces - the human interaction. This external force results in the addition of energy to the system, and a thermostat is a way to have that energy leak back to the thermal bath.

Note that you can run a simulation without a thermostat, in which case the system will retain whatever energy it starts with. In this case, user interaction should be minimised to prevent the continuous addition of energy resulting in instability.

All thermostats share the following options:

  • EquilibriumTemperature - The temperature the system will try to maintain (in degrees Kelvin).
  • Type - The name of the specific thermostat to be used. At the time of writing, the only option is the BerendsenThermostat.
  • MaximumTemperature - The maximum temperature a user can change the temperature to in this simulation.
  • MinimumTemperature - The minimum temperature a user can change the temperature to in this simulation.

Once a thermostat type is specified, additional options may be given. The Berendsen Thermostat has the following additional option:

  • BerendsenCouplingFactor - This value affects how strongly the thermostat will attempt to enforce the temperature. A higher value means a longer delay in restoring equilibrium temperature.

The Berendsen thermostat adjusts all velocities in the system equally, which can result in a dampening of the vibrational motion of a system when interacting with it. The AndersenThermostat may be a better choice in situtations where this is not desirable, and has the following additional option:

  • CollisionFrequency - This is the frequency (in ps-1) at which a collision with the bath will occur. A higher value results in more frequent dampening.

Integrator

To perform a simulation, an integrator is required to govern how the atom positions and velocities update in time.

All integrators share the following options which are required:

  • TimeStep - The size of the integration step (in picoseconds).
  • Type - The name of the Integrator to be used. The can be either the Velocity Verlet (VelocityVerlet) or Verlet (Verlet). Both integrators are similar, but the Verlet supports constraints.

You may also specify the following:

  • RemoveCMMotion - If set to true, centre-of-mass motion is removed from the simulation. For many simulations, this is a useful property as it prevents the simulation from drifting.

Logging

For scientific applications, it is usually desirable to store some data from the trajectory. Narupa has several loggers that can be specified to store various kinds of data. All loggers should be placed within the Logging node.

XYZ Logger

The XYZ logger (XYZ) will save the coordinates (and optionally velocities and forces) of the atoms of the trajectory. It has the following options:

  • LogPath - The directory in which to write files.
  • UseDateTimeStamp - Whether to append a date-time stamp to filenames. Default: false.
  • WriteFrequency - How often (in time steps) to write out the XYZ. Default: 1.
  • Positions - Whether to write positions. Default: true.
  • Velocities - Whether to write velocities. Default: false.
  • Forces - Whether to write forces. Default: false.
  • PotentialEnergy - Whether to write the potential energy in the comment line. Default: true.
  • NarupaUnits - Whether to output data in Narupa units (nm etc.), or the standard XYZ units (Angstroms, kcal/mol etc.). Default: false.

State Logger

The state logger (##StateLogger##) will store information such as temperature and energy in CSV format, and is useful for keeping an eye on your system behaviour. It has the following options:

  • LogPath - The directory to write this log in.
  • UseDateTimeStamp - Whether to append a date-time stamp to filenames. Default: false.
  • WriteFrequency - How often (in time steps) to write out the log. Default: 1.
  • Time - The simulation time (ps). Default: true.
  • PotentialEnergy - The potential energy (kJ/mol). Default: true.
  • KineticEnergy - The kinect energy (kJ/mol). Default: false.
  • Temperature - The instantaneous value of the temperature (K). Default: true.

Interaction Logger

The interaction logger (Interaction) will store data relating to interactive biasing applied by the user. This logger is more complex as information about the interaction such as force, time applied etc needs to be store, but also the set of atoms it was applied to. The logger achieves this by storing all the time dependent data about a reaction such as force etc in a CSV file, but stores the set of atoms interacted with in a JSON file, and creates a unique ID for mapping between the two files. The ID is of the form time_atomCount_hash, where time is the simulation time, atomCount is the number of atoms in the set, and hash is a hash function to ensure uniqueness. An example of the CSV file:

1 timeStep, time (ps), player ID, input ID, atom, energy (kJ/mol), net force X (kJ/(mol*nm)), , net force Y (kJ/(mol*nm)), , net force Z (kJ/(mol*nm)), interaction position X (nm), interaction position Y (nm), interaction position Z (nm) , center of mass of selected atoms X (nm), center of mass of selected atoms Y (nm), center of mass of selected atoms Z (nm)
2 151, 0.151, 56, 1, 0.15_1_3316, -19169.19, 582.3007, 237.8996, -1200.039, -1.507557, -0.8658723, -2.071692, -2.34, 1.47, -0.249
3 ...
4 1101, 1.101, 56, 1, 1.1_1_582, -19190.84, 349.2877, -940.3788, -47.67678, -1.948206, 0.3591044, -0.1219435, -2.34, 1.47, -0.249
5 ...
6 4361, 4.361, 56, 2, 4.36_30_19425, -12943.02, 2013.339, 1402.158, -3486.888, -1.57274, 1.31896, -1.207936, -1.728294, 1.210627, -0.9385331
7 4362, 4.362, 56, 2, 4.36_30_19425, -12943.14, 2015.886, 1399.105, -3486.244, -1.57274, 1.31896, -1.207936, -1.728489, 1.210864, -0.9385852
8 ...

and its corresponding json file:

1 {
2  "0.15_1_3316": [
3  3316
4  ],
5  "1.1_1_582": [
6  582
7  ],
8  "4.36_30_19425": [
9  639,
10  633,
11  634,
12  635,
13  636
14  ]
15 }

The interaction logger has the following options:

  • LogPath - The directory to write this log in.
  • UseDateTimeStamp - Whether to append a date-time stamp to filenames. Default: false.
  • WriteFrequency - How often (in time steps) to write out the XYZ. Default: 1.

The following options define what will be written out on each line:

  • Time - Whether to store the simulation time (ps). Default: true.
  • PlayerID - Default: true.
  • SelectedAtoms - Whether to store the selected atoms hash. Default: true.
  • InputID - Whether to store input IDs, e.g. controllers. Default: true.
  • Energy - Whether to store the energy of the interaction (kJ/mol). Default: true.
  • NetForce - Whether to store the net force applied to atoms (kJ/(mol*nm)). Default: true.
  • ForceMagnitude - Whether to store the length of the net force vector. Default: false.
  • Position - Whether to store the position of the interaction. Default: true.
  • SelectedAtomsCOM - Whether to store the centre of mass of the selected atoms in the interaction. Default: true.

Topology

The Topology section of the input file defines the set of molecules that are going to be simulated. The section subdivided into two further sections, Templates and Spawners.

Templates

A template is a description of how to simulate a particular set of atoms (typically a molecule, or macro-molecule). You may use existing templates, or create your own, and mix and match them to create a simulation using Spawners. Within the simulation file, a template can be explicitly described inline as in the example above, or as a reference to another file. For example, to load an existing nitrogen template, one might use:

1 <ResidueFile Name="N2" FilePath="^/Assets/Templates/N2.xml" />

The nitrogen molecule template file then looks like this:

1 <?xml version="1.0" encoding="utf-8"?>
2 <Template>
3  <Residue Name="N2">
4  <Atoms>
5  <Atom Element="Nitrogen" Position="0.06433,0,0" />
6  <Atom Element="Nitrogen" Position="-0.06433,0,0" />
7  </Atoms>
8  <Bonds>
9  <Bond A="0" B="1" />
10  </Bonds>
11  <ForceFields>
12  <InteractiveGaussianForceField GradientScaleFactor="1000" />
13  <MM3ForceField>
14  <MM3AtomMappings>
15  <MM3AtomMapping AtomPath="0" Type="107" />
16  <MM3AtomMapping AtomPath="1" Type="107" />
17  </MM3AtomMappings>
18  </MM3ForceField>
19  <LennardJonesForceField>
20  <LennardJonesAtomMappings>
21  <LennardJonesAtomMapping AtomPath="0" MM3Type="107" />
22  <LennardJonesAtomMapping AtomPath="1" MM3Type="107" />
23  </LennardJonesAtomMappings>
24  </LennardJonesForceField>
25  </ForceFields>
26  </Residue>
27 </Template>

Note The residue has to be given a name, which is then referenced by the spawner. It is advised to use the same name for a template everywhere.

Alternatively, an example of placing the definition of the Residues in line (useful for scientific simulations) is:

1 <IMDSimulation Name="UCH-L1">
2  <IMDSettings Port="54321" />
3  <Topology>
4  <Residue Name="2ETL">
5  <File Path="^/Assets/Simulations/TestSims/PDB/2ETL.pdb" />
6  </Residue>
7  <ExternalForceFields>
8  <InteractiveGaussianForceField GradientScaleFactor="1000" />
9  </ExternalForceFields>
10  </Topology>
11 </IMDSimulation>

A template consists of a series of Residues (aka molecules) which may combine together to form one large molecule or simply be a set of molecules that you want to be associated with each other. For each residue, the set of atoms, the bonds between them, and the forcefields that apply to them must be defined. The atoms and bonds may come from a chemical file format or may be written explicitly.

To import the atom positions, elements, names and bonds from a file, use a line such as:

1 <File Path="^/Assets/Simulations/TestSims/PDB/2ETL.pdb" />

Currently, PDB and XYZ files are supported. A topology will be created for the residue using the information available in the file. Note: The force fields still need to be defined!

If loading in a PDB file, it is possible to centre the atom coordinates to zero (this stops the protein from spawning outside of the simulation box). In order to do this, simply add a "Centre" attribute to the File tag:

1 <File Path="^/Assets/Simulations/TestSims/PDB/2ETL.pdb" Centre="true"/>

Note: Centering atom coordinates is only implemented for .pdb files and will not work for .xyz.

If your chemical file of choice isn't currently supported, then the Atoms and Bonds needs to be added manually (or with a script!).

In the Atoms section, each atom in the molecule is defined, with an Element string and its position within the molecule. Additionally, an atom may be given a name which can then be used to reference it. In the Bonds section, any topological bonds between the atoms defined in the Atoms section can be specified. A bond is defined by detailing the indices of the two atoms that are bonded, or by referencing the path to the atoms, as shown below. Note that this definition of the bond is purely for visualisation and topological information, the forces between the atoms need to be described.

1 <Residue Name="Heroin">
2  <Residue Name="HeroinAcetyl">
3  <Atoms>
4  ...
5  <Atom Name="h40" Element="Carbon" Position="-0.180434,-0.3335661,-0.247448" />
6  <Atom Name="h41" Element="Carbon" Position="-0.080434,-0.2335661,-0.247448" />
7  ...
8  </Atoms>
9  </Residue>
10  <Residue Name="HeroinMain">
11  <Atoms>
12  ...
13  <Atom Name="h0" Element="Oxygen" Position="-0.233834,-0.253666,-0.152648" />
14  ...
15  </Atoms>
16  </Residue>
17  <Bonds>
18  <Bond A="/HeroinMain/h0" B="/HeroinAcetyl/h40" />
19  <Bond A="/HeroinMain/h0" B="/HeroinAcetyl/h41" />
20  ...
21  </Bonds>
22  ...
23 </Residue>

In the ForceFields section, the force fields that are compatible with this molecule and any additional data that is required for that force field to simulate the molecule is given. In the nitrogen template above the Interactive Gaussian force field, the MM3 force field and the Lennard Jones force fields are defined. For simple molecules, these are likely to be the force fields you will want to use. For larger biomolecules, consider using the OpenMM force fields. More details on force fields and how to set them up can be found in Narupa Force Fields

Spawners

Once you have a set of templates in place, you need to decide how the molecules should be spawned. A Spawner takes a template and copies of it in the simulation. All spawners require a name and at least one template to be specified. There are currently the following spawners available:

  • Random Spawner (RandomSpawner or Spawner) - Places a template at random within the bounds of the box, such that it does not overlap with any other molecule. Requires the option ##Template## to be defined, pointing to the template to spawn, and Count, how many of that template to attempt to spawn.
  • Instance Spawner - Places the template in the simulation directly, using the coordinates given in the template file. Requires the option ##Template## to be defined.
  • Lattice Spawner - Spawns a lattice of the given template. Requires the option Template to be defined, as well as one of the following lattice types using the keyword Type: ** FFC - Spawns on a Front Facing Cube, requires the LatticeConstant option. ** PrimitiveCubic - Spawns on a Primitive Cubic Lattice, requires the LatticeConstant option. ** HCP - Spawns on a Hexagonal Close Pack Lattice, requires lattice constants A and C to be specified.

Note: If no spawners are specified, each Residue is spawned as in InstanceSpawner.

Loading in selection data from an external file

Selection data can be loaded in externally from a .json file. This has many advantages - namely, if sending a simulation to multiple front ends, this selection will be automatically loaded in for all players. To do this, add a SelectionData node to the Topology node like so:

1 <Simulation Name="Simulation">
2  <Topology>
3  <Templates>
4  <Residue Name="Residue">
5  <File Path="^/Assets/Simulations/TestSims/PDB/residue.pdb" />
6  <ForceFields>
7  ...
8  </Forcefields>
9  </Residue>
10  </Templates>
11  <Spawners>
12  <InstanceSpawner Name="Residue" Template="Residue" Count="1" />
13  </Spawners>
14  <SelectionData FilePath="Simulation.json" />
15  </Topology>
16 </Simulation>

The FilePath attribute should point to where the .json file is stored.

Troubleshooting

Below is a list of frequently encountered problems encountered while making a simulation, and some details about why they happen.

1 "Could not find a location to spawn the residue [x] in the box without it overlapping with other residues. Try making the simulation box bigger, and check that the residue initial positions are correct."

When this occurs, the simulation is trying to spawn a residue/molecule in the box specified (SimulationBox), but it can't fit the molecule in the box without the risk of being too close too another molecule. You should make the box bigger, or make your molecule smaller (reduce empty space in it).

1 Could not load file or assembly 'OpenMMNET, Version=1.0.0.0, Culture=neutral, PublicKeyToken=null' or one of its dependencies. An attempt was made to load a program with an incorrect format.

This error occurs when a simulation is using the OpenMM force field but OpenMM failed to load. Ensure that the program is compiled in 64-bit mode.

1 Failed to spawn topology into the simulation.
2 Spawner references unknown residue template [x].

This error occurs when the spawner does not recognise the name of the template you have spawned. Make sure the name given to a template is consistent.

Assuming your simulation gets past the loading phase, the main error one will encounter is a forced reset, with one of the following error messages:

1 "Propagation of positions resulted in 'not a number' exception in Positions for atom: "
1 "Calculation of forces resulted in 'not a number' exception in Forces for atom: "
1 "Velocity update resulted in 'not a number' exception in Velocities for atom:"

All these error messages are indicative of the same underlying problem: there is numerical instability in your simulation. There are several reasons this could occur:

  • Time step is too large. For systems with lots of vibrations, a timestep of less than 1fs (0.001ps) is recommended.
  • Particles got too close to one another during the simulation. This could be as a result of your simulation setup (check your units!) or due to user interaction. Try lowering the GradientScaleFactor in the Interactive Gaussian force field, or putting a cap on the maximum temperature.
  • A force field is misconfigured.

Sometimes it is possible to introduce numerical instability without triggering the exception. The symptoms of this are unrealistically long bonds, particles being 'stuck' in place, or wildly flickering in and out of existence. The triggers for this behaviour are usually the same as those listed above.

Narupa Force Fields

This page details the current set of force fields available in Narupa, and how to use them.

Interactive Gaussian

The interactive gaussian force field is the primary method in which a user can interact with a molecule in the system. It places a spherical Gaussian at the user interaction point, but only affects the nearest atom. It results in a smooth attraction of the atom to the interaction point.

Input File

Adding the interactive gaussian force field to a simulation is simple. Below is a template file for the simulation of Methane.

1 <Template>
2  <Residue Name="CH4">
3  <Atoms>
4  <Atom Element="Carbon" Position="0,0,0" />
5  <Atom Element="Hydrogen" Position="0.05288, 0.01610, 0.09359" />
6  <Atom Element="Hydrogen" Position="0.02051, 0.08240, -0.06786" />
7  <Atom Element="Hydrogen" Position="0.03345, -0.09314, -0.04496" />
8  <Atom Element="Hydrogen" Position="-0.10685, -0.00537, 0.01921" />>
9  </Atoms>
10  <Bonds>
11  <Bond A="0" B="1" />
12  <Bond A="0" B="2" />
13  <Bond A="0" B="3" />
14  <Bond A="0" B="4" />
15  </Bonds>
16 
17  <ForceFields>
18  <InteractiveGaussianForceField GradientScaleFactor="2000" >
19  <ParticleInteraction AtomPath="1" GradientScaleFactor="0"/>
20  <ParticleInteraction AtomPath="2" GradientScaleFactor="0"/>
21  <ParticleInteraction AtomPath="3" GradientScaleFactor="0"/>
22  <ParticleInteraction AtomPath="4" GradientScaleFactor="0"/>
23  </InteractiveGaussianForceField>
24 
25  <MM3ForceField>
26  <MM3AtomMappings>
27  <MM3AtomMapping AtomPath="0" Type="1" />
28  <MM3AtomMapping AtomPath="1" Type="5" />
29  <MM3AtomMapping AtomPath="2" Type="5" />
30  <MM3AtomMapping AtomPath="3" Type="5" />
31  <MM3AtomMapping AtomPath="4" Type="5" />>
32  </MM3AtomMappings>
33  </MM3ForceField>
34  <LennardJonesForceField>
35  <LennardJonesAtomMappings>
36  <LennardJonesAtomMapping AtomPath="0" MM3Type="1" />
37  <LennardJonesAtomMapping AtomPath="1" MM3Type="5" />
38  <LennardJonesAtomMapping AtomPath="2" MM3Type="5" />
39  <LennardJonesAtomMapping AtomPath="3" MM3Type="5" />
40  <LennardJonesAtomMapping AtomPath="4" MM3Type="5" />
41  </LennardJonesAtomMappings>
42  </LennardJonesForceField>
43  </ForceFields>
44  </Residue>
45 </Template>

The InteractiveGaussian element is added to the ForceFields node, and the GradientScaleFactor attributed is specified, which adjusts how strongly the molecule will be attracted/repelled. A value of 500 typically results in a pleasing, snappy interaction.

You can modify the interaction of specific atoms within the molecule by adding a ParticleInteraction element to the force field term. In the example above, these were added to so that the hydrogens do not interact with the user. The AtomPath field is the path to the atom within the template, and the GradientScaleFactor term is the scale factor that will be applied to the atom. The total scale factor is calculated by multiplying the atom's scale factor with the global scale factor for the molecule.

Lennard Jones

The Lennard Jones force field is the reference non-bonded interaction implemented in Narupa. It performs a simple evaluation of the Lennard Jones potential and its gradient. The implementation is not in any way optimal and has quadratic computational complexity, so should only be used for small systems. The algorithm requires the sigma and epsilon value of the potential for each atom in the system.

Input File

The Lennard Jones force field requires the sigma and epsilon values for each atom in a template to be specified. See the example file above for syntax details. The values for each atom are represented with the LennardJonesMapping element, which requires an AtomPath to the atom, then either the specific Sigma (nm) and Epsilon (kJ/mol) values, or an MM3Type. If an MM3 type is provided, the sigma and epsilon values from the MM3 database will be used.

MM3

The MM3 force field the reference force field for the simulation of hydrocarbons in Narupa. The reference implementation of the force field includes bonds, angles, torsions and out-of-plane bends.

It should be noted that MM3 was implemented as a starting point for Narupa, but will gradually be replaced by more sophisticated and general force fields. It is primarily parameterised for the simulation of hydrocarbons, so trying to simulate more complex biomolecules will be difficult, and another force field should be used (probably OpenMM).

Input File

The MM3 force field is added to the template like other force fields. For most molecules, the only information required is the MM3 type for each atom. An MM3 type is a number that corresponds the topological condition of the atom. For example, an sp4 carbon has type 1, sp3 carbon has type 2, hydrogen in an alcohol group has type 21 and so on. The full data set of atom types is stored in the mm3.xml file, stored in Assets/Data/mm3.xml. OpenBabel can generate MM3 atom types, but be warned that the atom types generated are typically naive and will not necessarily be correct. Given the types for each atom and the bonds specified in the topology section of the template file, the MM3 force field will attempt to automatically determine the bonds, angles and torsion terms required to simulate the molecule.

An illustrative example is the simulation of H2O + CH3, the product of Methane + OH. Below is the MM3 input for that molecule, taken from an EVB simulation file:

1 <MM3ForceField>
2 <DataFile File= "^/Assets/Data/mm3.xml" />
3 <MM3AtomMappings>
4  <MM3AtomMapping AtomPath="/H2O/0" Type="21" />
5  <MM3AtomMapping AtomPath="/H2O/1" Type="6" />
6  <MM3AtomMapping AtomPath="/CH3/0" Type="2" />
7  <MM3AtomMapping AtomPath="/H2O/2" Type="21" />
8  <MM3AtomMapping AtomPath="/CH3/1" Type="5" />
9  <MM3AtomMapping AtomPath="/CH3/2" Type="5" />
10  <MM3AtomMapping AtomPath="/CH3/3" Type="5" />
11 </MM3AtomMappings>
12 
13 <AdditionalTerms>
14  <OPBend A="/CH3/1" B="/CH3/0" C="/CH3/2" D="/CH3/3" FC="0.150" />
15 </AdditionalTerms>
16 </MM3ForceField>

The DataFile element specifies the location of the MM3 XML data file. The MM3AtomMapping terms define the MM3 types for each atom. The simulation of CH3 requires some additional instruction to the MM3 force field, as it needs an out-of-plane bend term to enforce planarity. Out-of-plane bends are currently the only additional terms that can be added.

Troubleshooting

The most common error one will run into with the MM3 force field is lack of parameterisation for the molecule being simulated. You may have been unable to find a correct MM3 atom type for one of the atoms in the molecule you are simulating, or you did select what seemed like appropriate atom types, and one of the following errors was thrown:

1 "Do not know how to handle bond between atoms {0} and {1}, with MM3 atom types {2} and {3}"
2 "Do not know how to handle angle between atoms {0}, {1}, {2}, with MM3 atom types {3}, {4} and {5}."
3 "No torsion data found for torsion with indices {0}, {1}, {2} and {3} with atom types ..."

Depending on the required accuracy of the simulation, you may simply find the atom type that is 'close enough' to give a qualitatively accurate depiction of the dynamics. If that will not suffice, you should either calculate parameters for that atom type for MM3 (a complicated process) or switch to a different force field.

OpenMM

The OpenMM integration with Narupa gives us the ability to perform research-grade simulations, with hardware-adapted performance. It is the recommended engine for any simulation requiring high accuracy or performance. OpenMM is a GPU-accelerated simulation suite primarily designed for the simulation of biomolecules but is highly flexible and customizable. We currently have support for the following features of OpenMM:

  • Harmonic bond forces
  • Harmonic angle forces
  • Nonbonded forces (all of the flavours including PME).
  • Torsion forces
  • Constraints
  • Centre of Mass motion remover.
  • GBSOAC Implicit Solvent

The details on the various forces can be found in the OpenMM user guide. The OpenMM user guide is a great resource for both the theory and implementation of molecular dynamics.

We have a sophisticated integration with OpenMM due to the C# wrapper we have created. For details on the generation of the wrapper, see the Narupa.OpenMM repository.

Input File

The input for the OpenMM force field has been designed to be plug-and-play with the "Forces" section of the XML output produced by OpenMM itself. Below is an abbreviated example of the amino acid alanine:

1 <OpenMMForceField TargetPlatform="CUDA">
2  <Forces>
3  <Force forceGroup="0" type="HarmonicBondForce" version="1">
4  <Bonds>
5  <Bond d=".1522" k="265265.6" p1="10" p2="4" />
6  ....
7  </Bonds>
8  </Force>
9  <Force forceGroup="0" type="HarmonicAngleForce" version="1">
10  <Angles>
11  <Angle a="1.91113553093" k="418.4" p1="0" p2="4" p3="5" />
12  ....
13  </Angles>
14  </Force>
15  <Force forceGroup="0" type="PeriodicTorsionForce" version="1">
16  <Torsions>
17  <Torsion k=".650844444444" p1="0" p2="4" p3="6" p4="7" periodicity="3" phase="0" />
18  ....
19  </Torsions>
20  </Force>
21  <Force alpha="0" cutoff="1" dispersionCorrection="1" ewaldTolerance=".0005" forceGroup="0" method="0" nx="0" ny="0" nz="0" recipForceGroup="-1" rfDielectric="78.3" switchingDistance="-1" type="NonbondedForce" useSwitchingFunction="0" version="1">
22  <Particles>
23  <Particle eps=".71128" q=".1414" sig=".324999852378" />
24  ....
25  </Particles>
26  <Exceptions>
27  <Exception eps="0" p1="0" p2="1" q="0" sig="1" />
28  ....
29  </Exceptions>
30  </Force>
31  </Forces>
32 </OpenMMForceField>

Hand-writing these input files is not recommended, although the ability to tweak them is useful. For complex simulations, the best way to generate OpenMM files is to use the OpenMM python wrapper and generate the forcefield yourself. Below is a simple script that takes a PDB file for the Villin headpiece, creates an OpenMM simulation of it using the AMBER force field, and prints out the OpenMM XML file.

1 from simtk.openmm.app import *
2 from simtk.openmm import *
3 from simtk.unit import *
4 from sys import stdout
5 
6 inputFile = 'villin_gas.pdb'
7 outputFile = 'villin.xml'
8 
9 pdb = PDBFile(inputFile)
10 forcefield = ForceField('amber99sb.xml')
11 system = forcefield.createSystem(pdb.topology, nonbondedMethod=NoCutoff, constraints=HBonds)
12 integrator = LangevinIntegrator(300 * kelvin, 1 / picosecond, 0.002 * picoseconds)
13 simulation = Simulation(pdb.topology, system, integrator)
14 simulation.context.setPositions(pdb.positions)
15 
16 system_xml = XmlSerializer.serialize(system)
17 
18 with open(outputFile, 'w') as f:
19  f.write(system_xml)

From the resulting file, the Forces section can then be cut-and-pasted into a Narupa template file, or the file can be referenced directly, as shown below. OpenMM supports creating simulations from PDB files, Gromacs, Amber and CHARMM input files, with a good range of force fields.

1 ...
2 <OpenMMForceField SystemFile="villin.xml"/>
3 ...

Troubleshooting

If the potential energy is zero, that is indicative that the OpenMM force field did not load correctly. Ensure that Narupa has been compiled to target 64-bit architectures.

If OpenMM throws an exception, then either your system is in a bizarre configuration or your OpenMM force field was generated using an unsupported version of OpenMM. We currently support OpenMM 7.1.1

EVB

The Empirical Valence Bond (EVB) method, is a way of simulating chemical reactions while maintaining the speed of a classical simulation. It works by simulating all possible configurations of the simulation (e.g. before a bond breaks, and after a bond breaks), then constructing a Hamiltonian that smoothly transitions from one state to another.

EVB Input File

To set up an EVB simulation, you first have to determine how each state will be simulated. Each state can use any of the existing force fields in Narupa. The input file for EVB consists of specifying how each state differs from the reactant (first) state, along with the force fields required to run the state. Below is an abbreviated example of a CH4 + OH reaction:

1 <EVBForceField>
2 
3  <CouplingElements>
4  <EVBCouplingElement State1="0" State2="1" Value="10" />
5  ...
6  </CouplingElements>
7 
8  <States>
9 
10  <EVBState EnergyShift="0" Name="CH4 + OH">
11  <ForceFields>
12  <MM3ForceField>
13  ...
14  </ForceFields>
15  </EVBState>
16 
17 
18  <EVBState EnergyShift="-50" Name="*CH3 + H2O">
19  <Bonds>
20  <RemoveBond A="/CH4/0" B="/CH4/1" />
21  <AddBond A="/OH/1" B="/CH4/1" />
22  </Bonds>
23  <ForceFields>
24  <MM3ForceField>
25  ...
26  </ForceFields>
27  </EVBState>
28  ....
29  </States>
30  </EVBForceField>

The CouplingElements section details how the various states are connected. A value of zero means the simulation cannot transfer between the two states, while a non-zero value affects the shape of the potential energy surface between the two states.

The next section is the list of states. Think of these as mini-templates. Each state takes an optional EnergyShift term that shifts the energy calculation by the specified amount, which can be used to alter the relative energies of the various states. The first one details the reactant state, and simply requires the force field definitions. The following states require any changes in the topology to be detailed via the AddBond and RemoveBond elements. These result in the addition or removal of a bond from the **//first//** state so can be used to simulate the breaking or formation of a bond. Each state may have different force field terms, for example, MM3 atom types may change as a result of a bond breaking.

Particle Restraint Force

It is sometimes desirable to apply a constraint to some atoms, fixing its position in absolute space. This can be achieved in Narupa by adding a ParticleRestraintForce to the ForceFields section of the input file. An example is given below:

1 <ParticleRestraintForce>
2  <ParticleRestraints>
3  <ParticleRestraint AtomIndex="8" K="2000" />
4  <ParticleRestraint AtomIndex="6" RestraintDimensions="XY" />
5  </ParticleRestraints>
6 </ParticleRestraintForce>

A restraint is defined for each atom, with AtomIndex referring to the index of the atom within the template, K referring to a force constant in kJ/mol (with default value 800 kJ/mol) and RestraintDimensions can be used to specify which dimensions to restrain on (by default XYZ).

Connecting Narupa to an Existing Molecular Dynamics Simulation

The simulations supported by the Narupa engine are by no means comprehensive, and you may have existing simulations set up in one of the popular molecular dynamics programs that you would like to visualize and interact with using the Narupa platform.

This is catered for via a plugin for Narupa that implements the VMD Interactive Molecular Dynamics (IMD) API. There are several packages that implement the VMD IMD API directly, including GROMACS and NAMD, but currently this plugin is guaranteed to be compatible with the IMD module in PLUMED, which has been re-enabled in a fork available here. Please check the list of MD packages supported by PLUMED to see if your MD package is supported.

The plugin provides a connection over a network to a simulation that is running the IMD API, in which the positions of the atoms will be sent to Narupa, and any forces applied by the user will be sent to the simulation. This way of running Narupa means that you can run your molecular dynamics simulations as normal in your favourite package, and just make a small change to your input files to launch the IMD.

This plugin is already deployed with the Narupa engine, so to use it only the MD program that you want to run interactive molecular dynamics in needs to be configured.

Using the plugin with PLUMED

Whichever MD code you choose to run needs to be compiled with a special version of plumed2, which is a fork containing the IMD plugin. This fork is available here. You will need to compile plumed and link/patch your MD code, by following PLUMED's documentation.

Note: When configuring this version of plumed, initialise the vmd-imd-api submodule, and be sure to add the IMD flag as a module:

1 cd /path/to/plumed2
2 git submodule update --init --recursive
3 ./configure --enable-modules=imd
4 make
5 make install

If OpenMM is your MD program of choice, then you need to add the OpenMM PLUMED Plugin to your OpenMM installation (see here).

Note for OpenMM users:

If using ubuntu 16.04 you may need to build the OpenMM PLUMED Plugin (and OpenMM itself) with the following compiler flags

1 //Flags used by the compiler during all build types.
2 CMAKE_CXX_FLAGS:STRING=-D_GLIBCXX_USE_CXX11_ABI=0

Once set up, your MD run needs to be configured to run using IMD. Here is a complete example script for OpenMM:

1 from __future__ import print_function
2 from simtk.openmm import app
3 import simtk.openmm as mm
4 from simtk import unit
5 from openmmplumed import PlumedForce
6 import time
7 
8 temp = 300
9 t_equil = 1
10 t_sim = 1
11 fric = 1
12 pdb_str = "file.pdb"
13 imd_port = 54321
14 imd_rate = 10
15 solvent_padding = 1
16 # run IMD if the imd rate is greater than zero.
17 runIMD = imd_rate > 0
18 
19 plumed_script = """
20  IMD HOST=127.0.0.1 PORT={} STRIDE=1 TRATE={} WAIT
21  """.format(imd_port, imd_rate)
22 
23 pdb = app.PDBFile(pdb_str)
24 print("Creating model for pdb")
25 modeller = app.Modeller(pdb.topology, pdb.positions)
26 
27 forcefield = app.ForceField('amber99sbildn.xml', 'tip3p.xml')
28 
29 if solvent_padding > 0:
30  modeller.addSolvent(forcefield, padding=solvent_padding * unit.nanometers)
31 
32 system = forcefield.createSystem(modeller.topology, nonbondedMethod=app.CutoffNonPeriodic,
33  nonbondedCutoff=2.0 * unit.nanometers, constraints=app.HBonds)
34 
35 # Add the plumed IMD force
36 if runIMD:
37  system.addForce(PlumedForce(plumed_script))
38 integrator = mm.LangevinIntegrator(temp * unit.kelvin, fric / unit.picoseconds,
39  0.5 * unit.femtoseconds)
40 integrator.setConstraintTolerance(0.00001)
41 
42 platform = mm.Platform.getPlatformByName('CUDA')
43 platform.loadPluginsFromDirectory("/usr/local/lib/openmm/plugins")
44 properties = {}
45 simulation = app.Simulation(modeller.topology, system, integrator, platform, properties)
46 simulation.context.setPositions(modeller.positions)
47 
48 print("minimizing...")
49 simulation.minimizeEnergy()
50 
51 simulation.context.setVelocitiesToTemperature(temp * unit.kelvin)
52 
53 
54 
55 # equilibrate
56 print("Equilibrating ...")
57 simulation.step(int(t_equil * unit.nanoseconds / (2 * unit.femtoseconds)))
58 
59 # now add the trajectory reporter.
60 print("Running dynamics...")
61 nsteps = int((t_sim * unit.nanoseconds) / (2 * unit.femtoseconds))
62 simulation.step(nsteps)

The key addition is the PLUMED script: IMD HOST=127.0.0.1 PORT=54321 STRIDE=1 TRATE=10 WAIT. Here the host is the IP address you want to connect to (127.0.0.1 means the localhost), and TRATE is how often, in md time steps, the positions of the atoms in the simulation will be transmitted to Narupa. The keyword WAIT results in the MD waiting for a user to connect before running.

With that in place, you just need to set up a small configuration file for Narupa. Here is a complete example:

1 <?xml version="1.0" encoding="utf-8" ?>
2 <IMDSimulation Name="2ETL">
3  <IMDSettings Port="54321">
4  <Command FileName="bash" Arguments="runOpenMMIMDSim.sh" WorkingDir="/path/to/sim/"/>
5  </IMDSettings>
6  <Topology>
7  <Residue Name="2etl">
8  <File Path="^/Assets/Simulations/PLUMED/2etl.pdb" />
9  </Residue>
10  <ExternalForceFields>
11  <InteractiveGaussianForceField GradientScaleFactor="200" HarmonicForce="true"/>
12  </ExternalForceFields>
13  </Topology>
14 </IMDSimulation>

Create this file, and add it to your ServerConfig.xml. Note that you need to specify the Port and provide a pdb file with the topology of the simulation that is running in your MD package. Also observe that the simulation is using the HarmonicForce flag in the interactive force, which replaces the vanilla Gaussian forcefield with a spring-like force which can be easier to control. Narupa can also handle launching of your MD package for you. In the example above, the Command node is used to run the following bash script:

1 /usr/bin/python3 simulation_imd.py 2etl.pdb -imdr 10

With that setup, you can just launch Narupa, which will launch your MD and connect to it automatically. Any script/program can be launched in this manner (be sure to pass any environment variables required for your MD program!).

Advanced Functionality: Resetting Simulations From Narupa

A useful feature of vanilla Narupa is the capability to reset simulations back to their initial conditions. This is not so straightforward to achieve with the VMD IMD plugin, as we do not have direct control over the MD package, and cannot reset positions/velocities. However, we can reset the simulation back to it's initial state by restarting the script used to launch the MD.

Due to a fun feature of sockets, for responsive resetting we need to launch the IMD on a different port each time we reset (unless we improve the protocol). Below is an example of setting this up for OpenMM (using the same python script as above):

Our bash script is modified to accept a port as the first argument:

1 /usr/bin/python3 simulation_imd.py 2etl.pdb -imdr 10 -imdp $1

And our NSB simulation XML file is modified to indicate that the first argument of the script can be interpreted as a port:

1 <?xml version="1.0" encoding="utf-8" ?>
2 <IMDSimulation Name="2ETL">
3  <IMDSettings Port="54321">
4  <Command FileName="bash" SupportsPort="true" Arguments="runOpenMMIMDSim.sh ${PORT}" WorkingDir="/path/to/sim"/>
5  </IMDSettings>
6  <Topology>
7  <Residue Name="2etl">
8  <File Path="^/Assets/Simulations/PLUMED/2etl.pdb" />
9  </Residue>
10  <ExternalForceFields>
11  <InteractiveGaussianForceField GradientScaleFactor="200" />
12  </ExternalForceFields>
13  </Topology>
14 </IMDSimulation>

Whenenever the simulation is reset, NSB will automatically increment the port used, substitute it for ${PORT} in the command string, and relaunch the script. This could be used in interesting ways, for example to return to the last checkpoint of a run rather than the start.