sDMD

KO | EN

gitlite — search

#javascript #python #hacktoberfest #react #ai #typescript #llm #go #golang #android #machine-learning #rust #deep-learning #linux

★ 16 Open GitHub ↗

Discontinuous Molecular Dynamics (DMD) Simulation Package

// repository documentation

Was this content helpful?

★ 0 (0 ratings)

**sDMD** - *Version 0.95*
=====================================

# Discontinuous Molecular Dynamics (DMD) Simulation Package

**Reference & Citing**

If you use this code in your academic work, please cite:
- Zheng, S.; Javidpour, L.; Sahimi, M.; Shing, K. S.; Nakano, A. 
 "sDMD: An open source program for discontinuous molecular dynamics simulation of protein folding and aggregation" 
 *Comput. Phys. Commun.* 2020, 247, 106873-106886.
 
- Zheng, S.; Sajib, M. S. J.; Wei, Y.; Wei, T. 
 "Discontinuous Molecular Dynamics Simulations of Biomolecule Interfacial Behavior: Study of Ovispirin-1 Adsorption on a Graphene Surface" 
 *J. Chem. Theory Comput.* 2021, 17 (3), 1874-1882.

**Table of Contents**
=====================

+ [Introduction](#introduction)

+ [Units](#units)

+ [Models and Interactions](#models-and-interactions)

+ [Models](#models)

+ [Intermediate-resolution Coarse Grained Model](#intermediate-resolution-coarse-grained-model)

+ [High-resolution All-atom Model](#high-resolution-all-atom-model)

+ [Interactions](#interactions)

+ [Algorithms](#algorithms)

+ [Energy Minimization](#energy-minimization)

+ [Confined Walls](#confined-walls)

+ [Replica Exchange Molecular Dynamics Simulation (REMD)](#replica-exchange-molecular-dynamics-simulation-remd)

+ [Program Structure](#program-structure)

+ [Program Usage](#program-usage)

+ [Installation](#installation)

+ [Run Parameters](#run-parameters)

+ [Non-REMD simulation](#non-remd-simulation)

+ [REMD Simulation](#remd-simulation)

+ [Analysis](#analysis)

+ [Format](#format)

+ [Coordinate File](#coordinate-file)

+ [GRO File](#gro-file)

+ [PDB File](#pdb-file)

+ [Force Field File](#force-field-file)

+ [Parameters of amino acid](#parameters-of-amino-acid)

+ [Potential Table](#potential-table)

+ [Functions](#functions)

+ [Simulation Source Code](#simulation-source-code)

+ [DMD.c](#dmd.c)

+ [DataSave.c](#datasave.c)

+ [Event.c](#event.c)

+ [Initialization.c](#initialization.c)

+ [List.c](#list.c)

+ [Models.c](#models.c)

+ [REMD.c](#remd.c)

+ [SGThread.c](#sgthread.c)

+ [ThreadProcess.c](#threadprocess.c)

+ [TimePrediction.c](#timeprediction.c)

+ [ToolFunctions.c](#toolfunctions.c)

+ [sREMD.c](#sremd.c)

+ [sServer.c](#sserver.c)

+ [Analysis Source Code](#analysis-source-code)

+ [Analysis.c](#analysis.c)

+ [FileManage.c](#filemanage.c)

+ [SystemInformation.c](#systeminformation.c)

+ [Cluster.c](#cluster.c)

+ [Energy.c](#energy.c)

+ [HBRamach.c](#hbramach.c)

+ [PBCAdjust.c](#pbcadjust.c)

+ [REMD.c](#remd.c-1)

+ [Tools.c](#tools.c)

**Introduction**
================

Discrete molecular dynamics (DMD) is an extremely fast alternative to the traditional molecular dynamics. It was first introduced in 1959 by Alder and Wainwright[1](#c11) for simulations of hard spheres. Later it was used by Rapaport[2-4](#c12) for simulation of polymer chains. Nowadays it was adopted for simulations of protein-like polymers. It is extremely fast and suitable for the simulations of large systems on long time scales.

During DMD simulation, potentials applied on particles are approximated by discontinuous step-functions of inter-particle distance *r* (see Figure *1-1*). Thus no forces will be exerted on particles until their distance becomes equal to the point of a discontinuity on the potential. When particle encounters a potential discontinuity, we call it have an "event". Between events, particles move at constant velocities. Because the energy change during each event is known, the post-event velocities can be calculated by solving the conservation of momentum and the conservation of energy simultaneously. Thus, the trajectory of a particle can be simulated discontinuously between events. The simulation code for DMD is significantly different from that for the traditional molecular dynamics due to the discontinuous movement of particles from one event to the next at known velocities, and the algorithm could become much more complex due to the event linking, time scheduling and especially the explicit hydrogen bonding mechanism.

*Figure 1-1 DMD potential (a) hard sphere collision (b) attractive square well interaction (c) repulsive soft sphere interaction (d) covalent bond and auxiliary bond interaction*

Here we proudly introduce *sDMD*, a simulation package based on the DMD technique with a high-resolution all-atom molecular model. The package is written by *C* language from scratch and has been highly optimized. Its modularized structure has it highly expandable: new molecular groups can be introduced by simply adding new formatted connection and potential data files, and no modification on the code would be necessary. The package, including the simulation executable and analysis executable, is also very easy to use: the configurations are controlled by an easy-to-read configuration file and several command flags.

Currently the package is capable to run simulations under various situations as follows,

-   NVT and NVE ensembles

-   systems in bulk

-   system in confined environments such as parallel pores, cylindrical tubes and spherical cavities

-   replica exchange molecular dynamics simulation (REMD)

-   dynamic confined systems as the confined walls expand, compress and pulse.

-   molecule flow through a tube between two reservoirs

and the analysis executable can be used to analyze the following,

-   remove periodic boundary conditions of trajectory

-   compute the number of different types of hydrogen bonds

-   compute Ramachandran plot

-   compute the contact maps between pairs of atoms and pairs of amino acids

-   compute system temperature, kinetic energy and potential energy

-   compute number of clusters during aggregation

-   analyze data from T-REMD by weighted histogram analysis method (WHAM)

This manual will describe the theories behind the DMD simulation, as well as the models of molecules and the algorithms of interactions. It will also describe the function of each file in the source code and demonstrate the ways to install and use the package.

*Reference*

[1] Alder, B. J.; Wainwright, T. E., Studies in Molecular Dynamics. I.General Method. *The Journal of Chemical Physics* **1959,** *31* (2), 459-466.[↩](#c11)

[2] Rapaport, D. C., Molecular dynamics simulation of polymer chains with excluded volume. *Journal of Physics A: Mathematical and General* **1978,** *11* (8), L213.[↩](#c12)

[3] Rapaport, D. C., Molecular dynamics study of a polymer chain in solution. *The Journal of Chemical Physics* **1979,** *71* (8), 3299-3303.[↩](#c12)

[4] Rapaport, D., The event scheduling problem in molecular dynamic simulation. *Journal of Computational Physics* **1980,** *34* (2), 184-201.[↩](#c12)

**Units**
=========

*Table 2-1 Basic units*

**Quantity** | **Symbol** | **Unit**
------------ | ---------- | --------
length | L | Å = 10-10 m
mass | m | u (atomic mass unit) = 1.6605402(10) x 10-27 kg (1/12 the mass of a 12C atom)
energy | E | kcal mol-1 = 4.184 kJ mol-1 = 4184 m2 kg s-1 mol-1

*Table 2-2 Derived units*

**Quantity** | **Symbol** | **Unit**
------------ | ---------- | --------
time | T | 50 fs = 50 x 10-15 s
velocity | V | Å (50 fs)-1 = 2 x 103 m s-1
temperature | T | 503 K
force | F | u Å (50 fs)-2 = 66.4 pN

*Table 2-3 Constants*

**Symbol** | **Name** | **Unit**
---------- | -------- | --------
*NA* | Avogadro's number | 6.0221367(36) x 1023 mol-1
*kB* | Boltzmann's constant | 1.987204118 x 10-3 kcal mol-1 K-1

**Models and Interactions**
===========================

This section will describe the molecular models and interaction algorithms that used in sDMD.

**Models**
----------

Different protein models can be used in sDMD ranging from high-resolution all-atom model to low-resolution coarse-grained model. The latest version of library only uses the high-resolution all-atom model. The early version uses coarse grained model. Both models share a lot of characteristics. Users can add customized molecular models into the library as desired. Just need to follow the required format of the model file and there is no need to modify the source code.

Here will give a description of the models starting from the low-resolution and then upgrade to the high-resolution.

### **Intermediate-resolution Coarse Grained Model**

In the early version, sDMD uses an intermediate-resolution coarse grained model, named PRIME, developed by Smith and Hall[1-4](#c31). It correctly reproduces the backbone geometry and has been highly successful in reproducing several important and experimentally determined features of proteins under bulk conditions. In PRIME, every amino acid is represented by four beads: N, Cα, C, and R (see *Figure 3-1*). The N bead represents the amide nitrogen and its hydrogen; the Cα bead represents the α-carbon and its hydrogen; the C bead represents the carbonyl carbon and oxygen; the R bead represent the side chain, all of which are assumed to have the same diameter of CH3 (alanine). The inter-peptide bond is assumed to be in the *trans* configurations, all the backbone bonds' lengths and bond angles are fixed at their experimentally measure values, and the distance between the consecutive Cα beads is also fixed according to the experimental data. *Table 3-1* presents all the relevant parameters of the model. Solvent molecules are NOT explicitly included in the model. The effect of solvent is factored into the energy function as a potential of mean force.

*Figure 3-1 The four-bead model of the protein backbone. Bold lines show covalent bond, dashed lies show auxiliary bonds which helps maintain correct backbone geometry, and thick broken lines show possible hydrogen bonds.[5](#c35)*

*Table 3-1 Parameter values for the proteins used in the DMD simulations*

**Beads** | **Diameter (Å)**
--------- | ----------------
N | 3.300
Cα | 3.700
C | 4.000
R | 4.408
 
 
**Bonds** | **Length (Å)**
--------- | --------------
Ni-Cα,i | 1.460
Cα,i-Ni+1 | 1.510
Ci-Ni+1 | 1.330
Cα,i-Ri | 1.531
 
 
**Auxiliary Bonds** | **Length (Å)**
------------------- | --------------
Ni-Ci | 2.450
Cα,i-Ni+1 | 2.410
Ci-Cα,i+1 | 2.450
Cα,i-Cα,i+1 | 3.800
Ni-Ri | 2.440
Ci-Ri | 2.490
 
 
**Bond Angles** | **Angle (deg)**
--------------- | ---------------
∠Ni-Cα,i-Ci | 111.0
∠Cα,i-Ci-Ni+1 | 116.0
∠Ci-Ni+1-Cα,i+1 | 122.0
∠Ri-Cα,i-Ci | 109.6
∠Ri-Cα,i-Ni | 110.1

Beads in the protein are subject to five different types of forces during events: (1) infinite repulsion due to excluded volume effect (hard-sphere collision), (2) infinite attraction between the beads connected with covalent bonds and auxiliary bonds, (3) attraction between pairs of backbone beads during hydrogen bond formation, (4) attraction between hydrophobic side chains, and (5) repulsion between the neighboring beads in hydrogen bond (soft-sphere collision). All these forces can be represented by the following step-functions,

This is for a hard-sphere potential. *r* is the distance between beads *i* and *j*; *sigmaij* is the sum of bead radii.

This is for a square-well or shoulder potential depending on the sign of *epsilon* (positive for well depth, attractive; negative for shoulder height, repulsive). *lammda_sigmaij* is the well/shoulder width.

This is for a bond interaction. *l* is the bond length and *delta* is the bond vibration tolerance.

Note that, as explained in the work by Takada et al,[6](#c36) interaction involving a pair of beads that are very close on the chain are more faithfully represented by an interaction between the atoms themselves, not the united atoms or the coarse grained beads. In other words, for example, neighboring carbonyl carbon atoms experience interactions based on the diameter of a carbon atom, not on the diameter of the C bead. Consequently, for interaction between pairs of beads that separated by three or fewer bonds, the beads are allowed to overlap in their simulations by up to 25% of their diameters; for those separated by four bonds, they are allowed to overlap up to 15%.

Backbone beads as well as Cα and R beads are connected by covalent bonds. Bonds are allowed to move freely over a small range between *l(1 - delta)* and *l(1 + delta)*. The choice of *delta* defines the acceptable range of bond vibration. In the simulation, *delta* is chosen to be 0.02375. Auxiliary bonds are used between next-neighboring beads along the backbone of the chain to hold backbone angles fixed. For example, an auxiliary is placed between Ni and Ci in *Figure (3-1)* to force Ni-Cα,i-Ci angle to be near its ideal value. Auxiliary bonds are also included between neighboring pairs of Cα beads to maintain their distances near the experimentally determined value. The combination of Cα,i-Ni+1, Ci-Cα,i+1, and Cα,i-Cα,i+1 auxiliary bonds and the Cα,i-Ci, Ci-Ni+1, and Ni+1-Cα,i+1 covalent bonds holds the inter-peptide group (Cα,i-Ci-Ni+1-Cα,i+1) in the *trans* configuration. Auxiliary bonds are also placed between side-chains and N and C beads on backbone. This is to hold the side-chain beads fixed relative to the backbone so that all model residues are *L*-isomers.

Hydrogen bond (HB) can be formed between N and C beads on backbones. Four criteria have to be met (see *Figure 3-2*): (1) the distance between Ni and Cj beads is 4.2 Å; (2) the neighboring beads, Ci-1 and Cα,i, Cα,j and Nj+1, are at the appropriate positions; (3) neither N nor C bead is already involved in a hydrogen bond; (4) N and C beads are separated by at least three intervening residues or they are in the different peptides. To ensure that criteria (2) is satisfied, we require that the four atom pairs Ni-Cα,j, Ni-Nj+1, Cj-Cα,i, and Cj-Ci-1 shown in *Figure (3-2)* be separated by a distance greater than *dij*, which is chosen to maintain the hydrogen bond orientations. Their values are given in *Table 3-2*.

*Figure 3-2 Hydrogen bond structure of coarse grained model. Thick dashed lines indicate hydrogen bond; thin dashed lines indicate auxiliary bonds helping to maintain hydrogen bond orientation, and bold lines show the covalent bonds of the backbone.*

*Table 3-2 Parameter values for hydrogen bond*

**Pairs** | ***dij(Å)***
--------- | -------------
Ni-Cα,j | 5.00
Ni-Nj+1 | 4.74
Cj-Cα,i | 4.86
Cj-Ci-1 | 4.83

When two beads Ni and Cj approach each other at the distance *dij* = 4.2 Å, we will evaluate the total potential energy change by checking the four auxiliary interactions. The potential energy change can be -*epsilonHB*, *0*, *epsilonHB*, 2*epsilonHB* and 3*epsilonHB*, depending on the orientation of the Ni, Cj, and their neighbors. If the kinetic energy is enough to overcome the total potential energy change, the forward reaction will happen. Ni and Cj beads will change their type into Ni' and Cj', respectively, and interact with each other according to the parameters related to their new types (the other interactions with Ni and Cj will not be affected). Otherwise, the two atoms Ni and Cj do not change their types and undergo a regular hard-sphere collision.

The thermal fluctuation distorts the orientation of hydrogen bond and large fluctuations may break the hydrogen bond. Once the two atoms Ni' and Cj' come again to the exact distance *dij* = 4.2 Å, a reverse reaction may occur. We will check the potential energy change by using the similar way we check during HB formation, and decide if the HB will break so that Ni' and Cj' will change back to their original types, Ni and Cj, respectively.

To maintain the temperature around the target value in NVT ensemble, we implement Andersen thermostat.[7](#c37) With this procedure, all beads in the simulation are subjected to random, infrequent collisions with ghost particles having the same mass as themselves. The post-event velocity of a bead colliding with a ghost particle is chosen randomly from a Maxwell-Boltzmann distribution centered at the simulation temperature. The collision frequency is controlled at 1% of the total events number.

The solvent molecules are not explicitly included. They are modeled implicitly by means of a square well attraction (potential of mean force) between hydrophobic side-chain beads. The side-chain beads, R beads, in a peptide are assigned as either "H" or "P": "H" means it is hydrophobic and such R beads will attract with each other when their distance is equal to 1.5*sigmaij* and as long as they are separated by at least three intervening residues in the same peptide or they are in the different peptides; "P" means the side-chain bead is polar. The interactions between two polar beads or between a polar and a hydrophobic side-chain beads will be a normal hard sphere collision.

### **High-resolution All-atom Model**

In order to apply sDMD to more complex systems, we upgraded the molecular and potential models from coarse grained to all-atom. Large amounts of codes were rewritten and optimized to fulfill the more modularized and efficient architecture. Details of the format of the library files will be described in Section 5.

The all-atom models are developed by Ding et al[8](#c38), and thanks for the helps of Dr. Ding, we got the updated force field data (the data on the paper was outdated). The basic concepts and algorithms used in the all-atom models are similar to those in the coarse grained models: it still uses auxiliary bonds to maintain the configuration of molecules, and it still uses the step functions to represent the forces. But now, after the updates, the model represents all the atoms explicitly (except for some light hydrogens), see *Figure 3-3*; the step function is more realistic and closer to the continuous force, see *Figure 3-4(a)*. Thanks to the updated multi-step potential profile, now most interactions combine both attraction and repulsion, unlike the coarse grained model that only has either one of them. It also adds a new type of force, the dihedral constraint. Although it is also maintained by an auxiliary bond, it has more complex profile than a regular bond, see *Figure 3-4(b)*.

*Figure 3-3 Schematic diagram for the all-atom protein model. The solid lines represent the covalent bonds and the dashed lines represent the auxiliary bonds.*

<img src="img/Potential.jpg" alt="Potential" height="300"/>
(a)

<img src="img/Dihedral.jpg" alt="Dihedral" height="300"/>
(b)

*Figure 3-4 (a) Non-bonded interactions in all-atom model. The continuous dashed line corresponds to the VDW and solvation interaction between two carbon atoms. The solid step function is the discretized potential in DMD. (b) Potential profile of dihedral constraint.*

Since the all-atom model represents the oxygen and hydrogen atoms explicitly, the structure of hydrogen bond is also optimized. The new hydrogen bond structure is shown in *Figure 3-5*.

*Figure 3-5 Hydrogen bond structure of all-atom model. Hi is the hydrogen atom; Aj is the acceptor and could be different types of atoms; Di is the donor and Xj is the heavy atom directly bonded to Aj. The thick solid lines indicate the covalent bonds between the atoms; the thin solid lines indicate the hydrogen bond; the dashed lines indicate the auxiliary bonds for helping maintain the hydrogen bond orientation.*

By adjusting the type of Aj atom, the hydrogen bond can be formed between backbone-backbone, backbone-side chain, and side chain-side chain. The donors include backbone amide hydrogen atoms and side chain polar hydrogen atoms of His, Trp, Tyr, Asn, Gln, Arg, and Lys. The acceptors include backbone carbonyl oxygens, side chain oxygen of Asp, Glu, Ser, Thr, and Tyr, and the side chain nitrogen of His.

The mechanism to model the hydrogen bond forming and breaking is similar to the coarse grained model: once the hydrogen, Hi, and the acceptor, Aj, reach the interaction distance, the distances between Hi-Xj and Di-Aj, which define the orientation of the hydrogen bond, will be evaluated; the total potential energy change, ΔE, between Hi/Aj and the other surrounding atoms are evaluated before and after the putative hydrogen bond formation; if these distances satisfy the preset range, and the total kinetic energy is enough to overcome the potential energy change, the hydrogen bond will be allowed to form, and forbad otherwise; after the formation of hydrogen bond, the atoms Hi and Aj will become Hi' and Aj', respectively, and they will interact with the other atoms by using the new pairs of potential.

**Interactions**
----------------

Here will describe the basic algorithms behind the DMD simulation. It will also give a description about the mechanisms in energy minimization, wall interaction, and REMD.

### **Algorithms**

Here we will describe how DMD works.

Assume we have two particles, *i* and *j*, each with diameter *sigma*. The coordinates of the particles are the vector *vector_ri* = [xi,yi,zi] = *vector_rj* = [xj,yj,zj] and they have velocity vector *vector_vi* and *vector_vj*. The collision event and collision time can be determined by using the relative position and relative velocity, *vector_rij* and *vector_vij*.

*Figure 3-6 Sample System of a collision of a pair of atoms (2D)*

A collision event is dependent of whether the distance between particles *i* and *j* will become equal to the average molecular dimeter
*sigmaij = (sigmai + sigmaj) / 2* sometime in the future. Using superscript "*o*" to represent an initial condition, the relative position vector at any time in the future is found by the following equation,

So if *vector_rij* and *vector_vij* have opposite sense, the particles will be approaching and a collision is possible. We may check to see if the dot product *vector_rij ∙ vector_rij* becomes equal to sigmaij2 in the future at some time increment *(t - to)*,

Here we define a variable *bij* as a dot product,

If *bij* > *0*, the relative position |*rij*| will be increasing with time. The two particles will be further away with each other. If *bij* = *0*, the relative velocity and relative position vectors are perpendicular. If *bij* < *0*, the relative position |*rij*| will be decreasing and there may be a collision. So after removing the superscript "*o*" on
*vector_rij* and *vector_vij*, *Eqn 3-7* can be rewritten as follows,

Solving for time increment by using the quadratic formula, we can have,

where *s_time = +/- 1*. The discriminant, *bij2 - vij2(rij2 - sigmaij2)*, in *Eqn 3-11* has to be non-negative for a collision to occur. If it is negative, then the particles will miss each other. The smaller real root *stime = - 1* is used for the solution to the quadratic as shown in *Figure 3-7*. The smaller real root occurs when the particles collide at *vector_rjs-*. The location labeled *vector_rjs+* is the root that will not occur physically.

*Figure 3-7 Only the smaller real root, stime = - 1 is used.*

If the particle has a potential shell (See *Figure 3-8*), the equation for the event time can be modified to,

where *dij* is the actual distance of centers for the event.

*Figure 3-8 (a) Approaching particles that will enter the potential shell; (b) particles in the shell that will have a hard sphere collision; (c) particles in the shell that will miss a hard sphere collision, but experience a shell interaction event after they pass;
(d) particles moving away that will experience a potential change event at the shell boundary.*

Calculation of velocity changes upon an event will be a little more complex. We use prime (') to denote the state after the event. The
change in potential energy for the event is indicated by *Δepsilon = epsilon' - epsilon*, which may be zero, positive or negative depending on whether it is a hard sphere collision, or particles enter or escape from a potential shell.

Total energy is conserved upon an event between particles *i* and *j*,

which can be rearranged,

When the two particles contact, the relative position will be,

The change in momentum for an event is proportional to the collision
vector *vector_rijc*,

*phi* is the parameter we will try to determine. Insert *Eqn 5-14* into
*Eqn 5-11*,

Rearrange *Eqn 5-14* also leads to,

and

Plugging *Eqn 5-16* and *Eqn 5-17* into *Eqn 5-15* to eliminate the
primed variables,

Define
*dij2 = vector_rijc* ∙ vector_rij*c*, *vector_vij = vector_vi - vector_vj* and *bij*c = vector_rij*c ∙ vector_vij*,
and rearrange *Eqn 5-18*,

Applying the quadratic formula,

Define the reduced mass, <img src="img/otherEqn.png" alt="otherEqn" height="50"/>, then,

where the sign of *svel* depends on the context of the event (See *Figure 3-9*). For an event where *Δepsilon* = *0*, *svel* = *- 1* is the only reasonable solution. For a hard sphere collision event, the distance is the inter-particle radius, *dij* = *sigmaij* and the relation becomes,

Based on the context of the event, the interaction calculation will follow the flowsheet below,

*Figure 3-9 Flowsheet for DMD with wells.[9](#c39)*

### **Energy Minimization**

Theoretically, sDMD can import any PDB file of proteins to create the simulation system. In the PDB file, the original coordinates of
atoms may violate the ranges of bond lengths or the minimum distances between pairs of atoms in sDMD's library. To eliminate these conflicts, sDMD uses steepest descent algorithm.

In DMD, the potential function is a step function, then the first-order derivative would not be possible. However, by computing the relative
positions of a pair of atoms, we could obtain the vectors that can guide the atoms to move closer or further; based on their distance, then, we can find at which step the potential currently locates, as well as to which direction the potential will decrease. From these information, we can decide the next move of an atom to approach the local minimum of potential energy.

In practice, to achieve a fast energy minimization, sDMD will only eliminate the bond length violations and atom diameter overlapping. In other words, the program will only vibrate an atom if it has bond that is too short or too long and if it overlaps any surrounding atoms. User can follow a short simulation under a low temperature to further relax the system.

### **Confined Walls**

sDMD is able to simulate a system in various types of confined environments, such as parallel pores, cylindrical tubes, and spherical
cavities. Currently the confined walls can be only smooth surfaces, which means there are no explicit atoms on the walls. The wall surface is like a multi-layer shell, and each layer represents a potential step.

To simulate the parallel pores is very straightforward: once the positions of a wall is set, the event time in the nearest future
when an atom interact with the wall can be easily determined, and the after-interaction velocity can be computed directly: it is just a one-dimension problem and only the velocity perpendicular to the wall's surface will be changed.

However, for a cylindrical tube or a spherical cavity, the situation will become a little bit more complex. Due to the curvature of surface and the implicit atoms, it would be difficult to determine the exact direction and position on the wall at which the atom will interact at the nearest event, both of which are required to compute the interaction potential and the after-interaction velocities.

To deal with this issue, sDMD uses an indirect way, see *Figure 3-10*: it is easy, instead, to compute the instant distance of an atom from the center of either a cylinder tube or a spherical cavity, *L*. Then the distance between the atom and the wall will be *D = R - L*, where *R* is the diameter of the tube or cavity. Thus, the potential region in the potential step can be determined by *D*. Next, instead of doing a lot of trigonometric calculations or trying to find a pseudo-atom on the wall to model the interaction, sDMD considers this as an interaction between the target atom and a huge pseudo-atom at the center of the tube or cavity. In other words, any interaction between an atom and the wall will be like the atom inside another huge atom interacts internally with the shell of the huge atom. The interaction types, however, must be switched. For example, the well-capture event between the atom and the pseudo-atom will become a well-escape event between the atom and the wall; similarly, a hard-sphere collision event will become a well-bounce event, etc. Thus, sDMD could use the same algorithm to compute both the interactions (a) between a pair of atoms and (b) between an atom and the wall.

*Figure 3-10 Schematic diagram of wall interaction in a cylindrical tube or a spherical cavity. The big solid circle is the wall boundary, which is infinitely repulsive; the big dashed circle is the potential shell of wall. The small black dashed circle is the original position of the target atom; the small black solid circle is the future position of the target atom after time t. The small read circle is the pseudo-atom at the center. The black line is the diameter of the wall, R; the red line is the distance between the target atom and the pseudo-atom, L; the blue line is the distance between the target atom and the wall, D = R - L.*

sDMD also allows users to use the obstruction wall. The obstruction wall is an infinite wall parallel to the *xy-*, *xz-* or *yz-* planes. Multiple obstruction wall can be inserted into the system box at any positions. This can be used to construct close-end cylindrical tubes, cuboid confined environments, and tunnel and reservoirs (combined with tunnel function), see Figure 3-11, or divide the system box into different regions.

*Figure 3-11 Combination of obstruction walls and cylindrical tube confinement can create a tube between two reservoirs, which can be used
to model protein translocations. The solid rectangle represents the boundary of the system box; the black dashed rectangles represent the
obstruction walls; the red dashed rectangle represents the cylindrical tube.*

### **Replica Exchange Molecular Dynamics Simulation (REMD)**

sDMD allows users to perform REMD. The results can be analyzed by T-WHAM[10-11](#c310) in the analysis executable.

REMD can be used to efficiently explore the potential energy landscape of molecular systems. The ruggedness and the slope toward the energy
minimum in the landscape govern sampling efficiency at a given temperature. Although escape out of local minima is accelerated at higher temperatures, the free energy landscape is altered due to large entropic contributions. To efficiently overcome energy barriers while maintaining conformational sampling corresponding to a relevant free energy surface, sDMD employs the replica exchange sampling
scheme[8](#c38),[12-13](#c12).

In REMD, multiple simulations/replicas of the same system are performed in parallel at different temperatures. Individual simulations are
coupled through Monte Carlo-based exchanges of simulation temperatures between replicas at periodic time intervals. Temperatures are exchanged between two replicas, *i* and *j*, maintained at temperatures Ti and Tj and with energies Ei and Ej, according to the canonical Metropolis criterion with the exchange probability *P*, where *P* = *1* if <img src="img/otherEqn2.png" alt="otherEqn2" height="50"/>, and *P* = *e-Δ*, if *Δ* > *0*. To achieve this, sDMD uses a Socket server to send and receive the temperature and energy data between pairs of replicas at the preset time intervals. sDMD would suggest the exchange rate to be every 1 x 103 time units.

*Reference*

[1] Voegler Smith, A.; Hall, C. K., alpha-helix formation: discontinuous molecular dynamics on an intermediate-resolution protein model. *Proteins* **2001,** *44* (3), 344-60.[↩](#c31)

[2] Voegler Smith, A.; Hall, C. K., Bridging the gap between homopolymer and protein models: A discontinuous molecular dynamics study. *The Journal of Chemical Physics* **2000,** *113* (20), 9331-9342.[↩](#c31)

[3] Smith, A. V.; Hall, C. K., Assembly of a tetrameric alpha-helical bundle: computer simulations on an intermediate-resolution protein model. *Proteins* **2001,** *44* (3), 376-91.[↩](#c31)

[4] Smith, A. V.; Hall, C. K., Protein refolding versus aggregation: computer simulations on an intermediate-resolution protein model. *J Mol Biol* **2001,** *312* (1), 187-202.[↩](#c31)

[5] Buldyrev, S. V., Application of Discrete Molecular Dynamics to Protein Folding and Aggregation. In *Aspects of Physical Biology: Biological Water, Protein Solutions, Transport and Replication*, Franzese, G.; Rubi, M., Eds. Springer Berlin Heidelberg: Berlin, Heidelberg, 2008; pp 97-131.[↩](#c35)

[6] S. Takada, Z. L.-S., P. G. Wolynes, Folding dynamics with nonadditive forces: A simulation study of a designed helical protein and a random heteropolymer. *J. Chem. Phys.* **1999,** *110* (23), 11616-11629.[↩](#c36)

[7] Andersen, H. C., Molecular dynamics simulations at constant pressure and/or temperature. *The Journal of Chemical Physics* **1980,** *72* (4), 2384-2393.[↩](#c37)

[8] Ding, F.; Tsao, D.; Nie, H.; Dokholyan, N. V., Ab initio folding of proteins with all-atom discrete molecular dynamics.
*Structure* **2008,** *16* (7), 1010-8.[↩](#c38)

[9] Elliott, J. R.; Lira, C. T., Supplement. In *Introductory chemical engineering thermodynamics*, 2nd ed.; Prentice Hall: Upper Saddle River, NJ, 2012.[↩](#c39)

[10] Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A., The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method. *Journal of Computational Chemistry* **1992,** *13* (8), 1011-1021.[↩](#c310)

[11] Gallicchio, E.; Andrec, M.; Felts, A. K.; Levy, R. M., Temperature weighted histogram analysis method, replica
exchange, and transition paths. *J Phys Chem B* **2005,** *109* (14), 6722-31.[↩](#c310)

[12] Zhou, R.; Berne, B. J.; Germain, R., The free energy landscape for beta hairpin folding in explicit water. *Proc Natl Acad Sci U S A* **2001,** *98* (26), 14931-6.[↩](#c312)

[13] Okamoto, Y., Generalized-ensemble algorithms: enhanced sampling techniques for Monte Carlo and molecular dynamics simulations. *J Mol Graph Model* **2004,** *22* (5), 425-39.[↩](#c312)

**Program Structure**
=====================

This section will briefly describe the structure of the program.

struct AtomStr{} is the most important data structure in sDMD. It stores almost all the information of an individual atom.

```C
struct AtomStr {
    struct PropertyStr *property;
    struct DynamicStr *dynamic;
};
```

"property" stores all the static information, like mass and atom number; "dynamic" stores all the dynamic information, like coordinate and
velocity.

sDMD has been optimized by Pthreads Library for the future development of parallel computing. The program will establish single/multiple threads and each of the thread can perform the calculation by its own. Thus, the data in sDMD have two layers, one is the raw data and the other is the data copy within each thread. *Figure 4-1* shows a diagram of a single-core simulation.

*Figure 4-1 Single-core simulation process*

At the beginning of a simulation, the thread will make a copy of the data of the target atom(s) involving in the nearest event and do the calculation only on this copy, i.e. the raw data will not be touched. Since DMD is an event-driven simulation, the results of the previous events would have voided the prediction of the current. Thus, any hazards that would have affected the current event will be evaluated. If there exist any hazards, the prediction of the current event will be voided. In this case, the program will re-predict the event for the target atom(s). Otherwise, the prediction for the current event will be committed and then the thread will predict a new future event for the target atom(s). After that, the raw data will be updated and the time tree (stores the event sequence) will be re-scheduled. The thread will keep looping this process until finishing the assigned simulation time.

**Program Usage**
=================

**Installation**
----------------

The package includes three folders: simulation source code folder, *src*, analysis source code folder, *analysis*, and an *input* folder. Each source code folder has its dedicated *Make* file, *Makefile*. The *input* folder contains the configuration file, *parameter.txt*, and the force filed library folder, *Library*.

To compile the simulation source code, use the following command,

```bash
$ make all
```

It will generate three executables, *sDMD*, *sServer*, and *sREMD*. *sDMD* is the core program to run the simulation; *sREMD* is the program to perform REMD simulation, only; *sServer* is an assistant program to create a socket server for REMD simulation, and will be executed by sREMD automatically.

If users only would like to perform regular DMD simulations, which means only *sDMD* executable will be required, simply use the command,

```bash
$ make sDMD
```

Note that if users would like to enable the visualization feature, which can visualize the simulation on the fly, he/she would need to install OpenGL utility library independently, and then use the following command to make,

```bash
$ make all VIS=yes
```

```bash
$ make sDMD VIS=yes
```

To compile the analysis source code, use the following command in *analysis* folder,

```bash
$ make
```

It will generate only one executable, *analysis*.

**Parameters**
------------------

### **Non-REMD simulation**

sDMD will need a coordinate file, *\*.gro* or *\*.pdb*, the configuration file, *parameter.txt*, and the force field library folder, *Library*, to start a simulation.

Below is an example of the configuration file,

```txt
Continue?            new
Time                 50
T(reduced)           0.65
Dump_Rate            10
CutOff_R(A)          9.0
Force_Field          Ding
Thermostat           Andersen
sDMD_Method          1
ThreadNo             1
WallExist?           no
WallType             sphere
WallDyn              no
Flow                 charge 1, 150.0 50.0 50.0 1.0 1.0
Obstruct?            4, x 50.0, x 150, y 0.0, y 100
Tunnel?              1, 50 150 50 50 30
SphObst?             1, 30 50 50 20
```

**Item** | **Description**
--------------- | ----------------
continue? | If users would like to perform a new simulation, use "new"; if users would like to extend/continue the previous run, which may have not finished yet, or would like to extend the current simulation, use "continue". Both process will back up the old files. The backup files will have a suffix representing the backing up time and date.
Time | Simulation time.
T(reduced) | Reduced system temperature, *T\** = *T* / 502 K. If set to "no", the program will keep using the temperature set at the last time.
Dump\_Rate | How often the data will dump. Here is each 10 time units.
CutOff\_R(A) | The cutoff radius of VDW interaction. Unit is A.
Force\_Field | The force field employed in the simulation. Currently only support "Ding"[1](#c51).
Thermostat | The thermostat algorithm. Currently only support "Andersen"[2](#c52). If set to "no", the thermostat will be turned off, and the simulation will use a microcanonical (NVE) ensemble instead.
sDMD\_Method | Parallel algorithm. Currently only support "1", which means only one CPU core will be used.
ThreadNo | Number of cores employed to perform the parallel computing. Currently only support "1".
WallExist? | If users would like to perform a simulation under bulk conditions, use "no". If users would like to perform a simulation in confined environments, use "smooth".
WallType | Shape of the wall. Support "parallel", "cylinder", and "sphere".
WallDyn | If the confined walls are fixed, use "no". If the confined walls are dynamic, use "compress", "extend", or "repulse".
Flow | Driving force of the flow. Currently support "velocity", "force" and "charge". If use "velocity", every atom in the system will be applied by a constant velocity. if use "force", every atom in the system will be applied by a constant force during thermostat interaction. if use "charge", the charged atom in the system will be attracted or repulsed by defined charges. Follow the keyword "velocity" follows the velocities on *x-*, *y-* and *z-* axis, unit is 0.02 Å \/ ps. Exp: 50.0 0.0 0.0 Follow the keyword "force" follows the forces on *x-*, *y-* and *z-* axis, unit is pN. Exp: 5.0 0.0 0.0. Follow the keyword "charge" follows the total charge number and the charge coordinates on *x-*, *y-* and *z-* axis and the potential gradient (unit kcal \/ mol) and the interaction distance (unit Å). Exp: 1, 150.0 50.0 50.0 1.0 1.0
Obstruct? | If there exists any impenetrable wall in the system. If no, use "no". Otherwise, the first number represents how many walls will be placed in the system, following the axis type and coordinate of each wall, separating by commas.
Tunnel? | If there will be any tunnel in the system. This is only used for translocation simulations. If no, use "no". Otherwise, the first number represents how many tunnels will be placed in the system, following the start *x-* coordinate, the end *x-* coordinate, the *y-* and *z-* coordinates, and the diameter of tunnel, separating each tunnel by commas.
SphObst? | If there will be any spherical obstacles in the system, following the total obstacle number, the coordinates on *x-*, *y-* and *z-* axis, and the diameter of obstacle, separating each obstacle by commas.

To run the executable sDMD, use the following command,

```bash
$ ./sDMD [FLAGS]
```

Use the flag *-h* or *-help* to show the supported *FLAG*s,

**Flag** | **Description**
--------------- | ---------------
-i | directory of the parameter folder
-o | directory of the output folder
-f | (optional) follow the file name of the input coordinate file
-si | (optional) follow the file name of the input saved data for continuity
-so | (optional) follow the file name of the output saved data for continuity
-trj | (optional) follow the file name of the output trajectory file
-cnt | (optional) follow the file name of the output connection map
-log | (optional) follow the file name of the output Log file 
-sys | (optional) follow the file name of the output SysInfo file
-pot | (optional) output potential energy, follow the file name, otherwise will use the default name
-kin | (optional) output kinetic energy, follow the file name, otherwise will use the default name
-tem | (optional) output temperature, follow the file name, otherwise will use the default name
-HBn | (optional) output hydrogen bond number, follow the file name, otherwise will use the default name
-xyz | (optional) output xyz trajectory, follow the file name, otherwise will use the default name
-pdb | (optional) output pdb trajectory, follow the file name, otherwise will use the default name
-box | (required only if the coordinate file is a .PDB file) the simulation box dimensions, x, y, z
-Wsz | (required only if WallDyn is assigned) follow the maximum changing size of the wall, default 5 A
-Wrt | (required only if WallDyn is assigned) follow the total time needed for the size change, default 200
-ter | (optional) interactive termini selection, instead of charged (default)
-arg | (optional) interactive arginine selection, instead of charged (default)
-asp | (optional) interactive aspartic acid selection, instead of charged (default)
-glu | (optional) interactive glutamic acid acid selection, instead of charged (default)
-his | (optional) interactive histidine selection, instead of charged (default)
-lys | (optional) interactive lysine selection, instead of charged (default)
-AlC | (optional) charge all the atoms with 1 proton
-AlP | (optional) charge the atoms in each peptide with [peptide net charge / peptide atom number]
-npC | (optional) charge source will have no PBC
-COM | (optional) keep the velocity of center of mass of the molecules
-REMD | only use during REMD, follow the server name, the port number, the temperature, the exchange rate, and the suffix name of saving files (this flag will be set automatically by REMD executable)
-visual | (optional) enable real-time visualization, default disabled

### **REMD Simulation**

Besides the required files and commands mentioned above, to perform a REMD simulation, sDMD needs another configuration file, *REMDConfig.txt*, and the other two executables, *sServer* and *sREMD*. The configuration file must be in the *input* folder and the two executables must be in the same folder of executable *sDMD*. The *T(reduced)* item in *parameter.txt* should be set to "no".

Below is an example of the configuration file (the first line is just a comment),

```txt
#REPLICA EXCHANGE MD CONFIGURATION FILE
ReplicaNum   = 3
Temperatures = 275 350 450
SocketPort   = 20202
ExchangeRate = 50
```

**Item** | **Description**
---------- | ----------
ReplicaNum | The number of replicas.
Temperatures | The temperature of each replica. Unit is K.
SocketPort | The port number of the socket server.
ExchangeRate | How often the temperature will be exchanged. Unit is the time unit.

To run the executable *sREMD*, use the following command,

```bash
$ ./sREMD [distribute T or not] -f [configuration file] -args [args of executable without flag -REMD]
```

**Flag** | **Description**
------------ | ---------
-nodist | [yes] Default is yes, if absent. Otherwise, the program will NOT distribute the preset temperature to each replica. Use this flag if users would like to extend/restart the simulation.
-f | follow the file name of the REMD configuration file
-args | flags used in each individual replica

The flags of executable *sServer* and the *-REMD* flag of executable *sDMD* will be set by *sREMD* automatically.

### **Analysis**

To perform an analysis, use the following command,

```bash
$ ./analysis [FLAGS]
```

Use the flag *-h* or *-help* to show the supported *FLAGs*,

**Flag** | **Description**
----------- | -----------
-path | exact path of the data folder
-trj | trajectory input file
-cnt | connect map input file
-sys | (optional) system info input file
-log | (optional) log input file
-obs | (optional) exact path of the obstacle info input file
-nPP | (optional) only analyze this specific number of peptide
-sum | (optional) ignore the old time record in .log file
-rPBC | remove PBC, call at the beginning. require -trj input file
-HB | analyze HB info. require -cnt input file
-En | analyze energy info
-Ag | analyze aggregation info
-RG | calculate RG (not support yet)
-MSD | calculate MSD (not support yet)
-Ramach | (require -HB) plot Ramachandran plots for each amino acid
-ConMap | (require -HB) plot contact maps for atoms and amino acids
-REMD | analyze REMD data. follow the flags below:
-reNo | number of replicas
-bin | (optional) number of bins in histogram. default 10
-aHB | (optional) use alpha HB number as the reaction coordinate
-bHB | (optional) use beta HB number as the reaction coordinate
-tHB | (optional) use total HB number as the reaction coordinate
-maxE | (optional) max value of potential energy (most negative). default -150
-minE | (optional) min value of potential energy (least negative). default 0
-rate | (optional) data dumping rate, unit is the time unit. default 10
-temp | (optional) target temperature. default 300
-st | (optional) start frame of analyzing. default 0
-et | (optional) end frame of analyzing. default -1 (the final frame)

The analysis program is also highly modularized and optimized. It will be easy to add more analysis functions into it.

*Reference*

[1] Ding, F.; Tsao, D.; Nie, H.; Dokholyan, N. V., Ab initio folding of proteins with all-atom discrete molecular dynamics. *Structure* **2008,** *16* (7), 1010-8.[↩](#c51)

[2] Andersen, H. C., Molecular dynamics simulations at constant pressure and/or temperature. *The Journal of Chemical Physics* **1980,** *72* (4), 2384-2393.[↩](#c52)

**Format**
==========

**Coordinate File**
-------------------

sDMD supports *.pdb* and *.gro* files as the coordinate input file. But it recommends to use a *.gro* file as the coordinate file, since the *.pdb* file has too many alternative formats.

Users can install Gromacs independently and use the following command to convert *.pdb* file to *.gro* format.

```bash
$ gmx pdb2gmx -f coor.pdb -o coor.gro -ignh
```

During the converting process, Gromacs would ask users to choose a force field type, i.e. GROMOS, CHARMM, and OPLS, it recommands to choose any GROMOS force field since the defined atom names in the library is similar to those applied in sDMD. After converting, users would need to delete all the solvent molecules in the *.gro* file if there exists and remember to renew the total atom number at the 2nd line of the file.

When sDMD imports the coordinate file, it will first read the original *.pdb* or *.gro* file. During this process, the program will remove the atoms those are not in the model library, such as the light hydrogens and one of the carboxylic oxygen on the C terminal of a peptide. At the same time, the program will extract the atom information and assign them onto the internal variables about protein/peptide, amino acid and atom.

After read through the whole coordinate file, some data will be dumped into a *.xyz* file. The program will then read the *.xyz* file in order to finish the final assignments for some other variables. The reason of such process is that some of the variables can be correctly assigned only after the program has read through and modified the original coordinate file.

The output trajectory file will be of *.gro* format. It consists of multiple continuous snapshots of the atoms' coordinates. The snapshot/data dumping rate is defined in *parameter.txt*.

### **GRO File**

Pay attention that the coordinate unit of *.gro* file is nm. The length unit in sDMD, however, is Å. The program will convert the unit
automatically.

The following is an example of the *.gro* file,

```txt
model
  375
    1ASP      N    1   1.110   2.333   3.994
    1ASP     CA    2   1.234   2.331   3.919
    1ASP      C    3   1.207   2.279   3.778
    1ASP      O    4   1.143   2.175   3.761
    1ASP     CB    5   1.337   2.238   3.984
    1ASP     CG    6   1.284   2.100   4.023
    1ASP    OD1    7   1.196   2.095   4.111
    1ASP    OD2    8   1.333   2.002   3.964
    2ALA      N    9   1.256   2.353   3.680
    2ALA     CA   10   1.237   2.315   3.540
    2ALA      C   11   1.370   2.330   3.467
    2ALA      O   12   1.476   2.292   3.519
    2ALA     CB   13   1.128   2.403   3.479
    2ALA      H   14   1.308   2.439   3.696
    3GLU      N   15   1.363   2.383   3.346
    3GLU     CA   16   1.482   2.402   3.264
    …
   7.00000   7.00000   7.00000
```

The first line is a comment, can be a brief description of the system. Only one line of comment will be allowed.

The second line has only one number represents the total atom number in the system. In this example it has 375 atoms. To save pages here does not show them all.

The line after the second contains the information of each atom. Use the line below as an example,

```txt
    2ALA     N    9    1.256 2.353 3.680
```

The first term describes at which amino acid in a protein/peptide the atom locates. Here is the second amino acid whose name is "ALA". Pay attention that the amino acid number here is not continuous between proteins/peptides, which means that if there are multiple molecules in the system, this number will start from 1 for each individual molecule. This will help the program to judge how many molecules there are in the system. The second term is the atom's name. The third term is the sequence number of this atom in the whole system. This number is continuous in the whole system. The following three terms are the atom's coordinates on *x-*, *y-* and *z-* axis, respectively. If this is a output *.gro* file from the simulation, there will be three more terms following the coordinates. They are the atom's velocities on *x-*, *y-* and *z-* axis, respectively.

The three numbers on the last line of the *.gro* file are the dimensions of the simulation box, on *x-*, *y-* and *z-* axis, respectively.

### **PDB File**

The following is an example of the supported *.pdb* file,

```txt
REMARK    THIS IS A SIMULATION BOX
CRYST1  150.000  150.000  150.000  90.00  90.00  90.00 P 1        1
ATOM      1  N   ASP X   1      60.227  61.745  83.465  0.00  0.00            
ATOM      2  CA  ASP X   1      61.467  61.725  82.715  0.00  0.00            
ATOM      3  C   ASP X   1      61.197  61.205  81.305  0.00  0.00            
ATOM      4  O   ASP X   1      60.557  60.165  81.135  0.00  0.00            
ATOM      5  CB  ASP X   1      62.497  60.795  83.365  0.00  0.00            
ATOM      6  CG  ASP X   1      61.967  59.415  83.755  0.00  0.00            
ATOM      7  OD1 ASP X   1      61.087  59.365  84.635  0.00  0.00            
ATOM      8  OD2 ASP X   1      62.457  58.435  83.165  0.00  0.00            
ATOM      9  N   ALA X   2      61.687  61.945  80.325  0.00  0.00            
ATOM     10  CA  ALA X   2      61.497  61.565  78.925  0.00  0.00
…
END
```

The *.pdb* file can have several lines of comments/descriptions at the head. The useful information is on the lines starting with the keyword "ATOM". Use the following line as an example,

```txt
ATOM      9  N   ALA X   2      61.687  61.945  80.325  0.00  0.00            
```

The first term is the keyword. The second term is the sequence number of this atom in the whole system, which is continuous. The third term is the atom's name. The fourth term is the amino acid name at which the atom locates. The fifth term is just an arbitrary symbol for the protein molecule. The atoms in the same molecule will share the same symbol. Different molecules can have the same symbol, but not the adjoining two. The sixth term is the amino acid's sequence number in the protein. This number is not continuous between different molecules. The following three terms are the atom's coordinates on *x-*, *y-* and *z-* axis, respectively. sDMD will ignore all the terms behind.

To import the coordinates from a *.pdb* format file, users need to specify the dimensions of the simulation box by using -box flag while executing sDMD.

**Force Field File**
--------------------

The force field consists of two types of files: one contains the parameters of amino acid and the other contains the potential table.

### **Parameters of amino acid**

There are 22 parameter files for the 21 types of amino acids. Proline has 2. The file's name indicates which amino acid it represents. Here use ALA.txt as an example,

```txt
RESIDUE ALA
   ATOM N      NZB    14.0  N N
   ATOM CA     CA1    12.0  N N
   ATOM C      CRB    12.0  N N
   ATOM O      OZB    16.0  N N
   ATOM CB     CA3    12.0  N N
   ATOM H      HB     12.0  N N
   BOND N      H      1.000 0.020
   BOND N      CA     1.458 0.020
…
  ANGLE CA     H      2.152 0.020
  ANGLE N      CB     2.460 0.012
…
  INTER CA     +N     2.425 0.014
  INTER C      +N     1.329 0.020
…
CONSTR2 C      +C     2.700 2.846 1.500 2.910 0.500 4.000
…
```

The first line is the name of the amino acid. The lines starting with the keyword "ATOM" list all the atoms in this amino acid. Pay attention that all the amino acids, except for proline "PRO", share the same backbone characteristics of alaline "ALA" (glycine "GLY" skips the sidechain "CB" in "ALA"). So for those amino acids the "ATOM" lines only list their own dedicated atoms and the backbone data will be read from *ALA.txt*. Use the following line as an example,

```txt
   ATOM CA     CA1    12.0  N N
```

The second term is the atom's name. The third term is the atom's type. (All the supporting types of atoms can be find in *InteractionPotentialTable.txt* in the library folder.) The fourth term is the atom's mass. Pay attention that the mass of hydrogen is set as 12.0, as the author suggests.[1](#c61) This can avoid over-vibration of hydrogen atoms due to its small mass and accelerate the simulation speed. If users would like to change the mass back to 1.0, he/she can modify the source code of function *ScanAA* in *initialization.c*. The last two terms are left for the future development. They can be assigned by any values which gives the atom extra properties.

The lines starting with the keyword "BOND" list all the bonds in this amino acid. Again, all the amino acids except "ALA", "GLY" and "GRO" only list their own dedicated bonds. Use the following line as an example,

```txt
   BOND N      H      1.000 0.020
```

This line means there is a covalent bond between atom N and atom H. The bond length is 1.000 Å and can vibrate between *1.0 x (1 - 0.02)* and *1.0 x (1 + 0.02)*.

The format of the lines starting with the keywords "ANGLE" and "INTER" are similar to "BOND" lines. They are however auxiliary bonds. The "+" symbol in front of the atom name means this atom locates in the next amino acid.

The line starting with the keyword "CONSTR" list all the dihedral constraints in the amino acid. Use the following line as an example.

```txt
CONSTR2 C      +C     2.700 2.846 1.500 2.910 0.500 4.000
```

The number after the keyword "CONSTR" represents the number of steps there will be in the potential function. The second and the third terms are the head and tail atoms in the dihedral angle. Starting from the fourth term are the values construct the potential function. The first and the last values represent the minimum and the maximum distance between the two atoms; the values in the middle are presented as pairs: the first values are the discontinuity positions in the step function and the second values are the energy change crossing the steps. In this example, the maximum distance is 4.000 Å; when the two atoms move closer than 2.910 Å, the potential energy will decrease 0.500 kcal / mol, or become less negative by 0.500 kcal / mol; when they move even closer than 2.846 Å, the potential energy will further decrease 1.500 kcal / mol. The distance between the two atoms cannot be smaller than 2.700 Å.

sDMD allows users to use any supporting atom to construct the wall surface. Simply copy the atom info from the library into the file *wall.txt* locates in the library folder. The default is "CA". Users can change it at will. Below is an example,

```txt
Wall
 ATOM CA CA 12.0 N N
```

### **Potential Table**

There are two types of potential table: one is for VDW interaction, called "InteractionPotentialTable.txt", and the other is for hydrogen bond interaction, called *HBPotentialTable.txt*. They share the same format and are similar to that used in "CONSTR" line aforementioned. Any lines starting with "\#" symbol are the comments. Use the following line from *InteractionPotentialTable.txt* as an example,

```txt
HB   CA   2.000000 2.50000 0.20000 4.50000 0.02000 6.50000 0.01000
```

The first two terms represent the pair of atoms that holds the following potential function. The only difference from the "CONSTR" line is that there is no maximum distance here, which means the number of values after the atom names is always odd, while in the "CONSTR" line it is even.

*Reference*

[1] Ding, F.; Tsao, D.; Nie, H.; Dokholyan, N. V., Ab initio folding of proteins with all-atom discrete molecular dynamics. *Structure* **2008,** *16* (7), 1010-8.[↩](#c61)

**Functions**
=============

This section will briefly describe the core functions in the main source code files. Most functions can be identified by their names.

**Simulation Source Code**
--------------------------

#### DMD.c

This file contains the main function.

#### CBT.c
This file contains the functions to maintain the complete binary tree (CBT).

#### DataSave.c

The functions in this file manage naming, backing up, creating, initializing and saving. All the saving processes will call the function *SaveData()*.

#### Event.c

This file contains all the functions used to process the different events. The *DoEvent()* manages which event function will be called. The events like collision, bonding, formation and breaking of hydrogen bonds all call the function *InteractionEvent()*. Hydrogen bonding also has its own event function, *HBEvent()*. Other events have their own event functions: *ThermostatEvent()* for thermostat event, *PBCandCrossCellEvent()* for period boundary crossing and cell crossing event, *WallEvent()* for regular wall interaction event, *ObstEvent()* for additional obstacle wall interaction event, *TunnelEvent()* for tunnel interaction event, *ChargeEvent()* for electrostatic interaction event, and *SphObstEvent()* for spherical obstacle interaction event.

The function *LinkList()* manages the linked list between cells. To increase the calculation speed, sDMD divides the simulation box into many cells and keeps tracking the cell number in which an atom locates. Thus, by using the linked list between the atom number and the cell number, it will be easy to find which atoms are in the target cell. Combining the linked list and the cut off radius, it is easy and fast to scan the 27 cells surrounding the target atom by which its next event can be computed.

#### Initialization.c

The functions in this file are used to read and import all the configurations, parameters and coordinates. They also allocate the memory for most of the variables in the simulation.

The functions are designed based on the formats of the force field.

#### List.c

This file contains the functions to manage the dynamic linked list.

#### Models.c

The functions in this file are used to convert the models' names (of atoms and amino acids) into numbers.

#### REMD.c

The functions in this file are only used to perform the calculations for just one replica in a REMD simulation. *REMD()* initializes the variables. *REMDRun()* does the main jobs. *ReplicaExchange()* manages the data exchange, and *client\_open()* initializes a socket client for this replica.

There are two other files to create replicas and server, *REMDExec.c* and *Server.c*, respectively.

#### SGThread.c

This is the core file to perform a single-core non-REMD simulation. *SingleThread()* initializes the variables and *SingleThreadRun()* does the main jobs.

#### ThreadProcess.c

The functions to perform the main procedures of a DMD simulation are in this file. *FirstRun()* is only used when start a fresh new simulation. It will calculate the coming events for all the atoms for the first time and create a complete binary tree to manage the events sequence. *ProcessEvent()* is used by each thread to process the event, by the function *DoEvent()* (located in *Events.c*), and predict the future event, by the function *Predict()*. These two procedures will be done only if there exists no hazard for the current event, checked by the function *HazardCheck()*. Otherwise, the current event is marked as invalid and need to be re-predicted. *CommitEvent()* is used to commit the updated atom data to the raw data.

#### TimePrediction.c

This file contains all the functions used to predict the event time. Different events have their own time prediction functions. The function name indicates the event type.

#### ToolFunctions.c

This file contains all the helper functions. The function FindPair() is used to find the right potential step at which the target atom pair locates.

#### sREMD.c

The functions in this file are only used to perform REMD simulation. It helps to establish multiple replicas and import arguments and flags into the executables *sDMD* and *sServer*.

#### sServer.c

The functions in this file are only used to perform REMD simulation. It helps to create and maintain a socket server for the data exchanging between replicas.

**Analysis Source Code**
------------------------

#### Analysis.c

This file contains the main function.

#### FileManage.c

The functions in this file help to assign names to files.

#### SystemInformation.c

The functions in this file are used to read the system information data from the *.dat* file and *.log* file output from the simulation.

#### Cluster.c

The functions in this file are used to calculate the cluster evolution during aggregation.

#### Energy.c

Th functions in this file are used to calculate the total kinetic energy by *CalKeEnergy()*, and potential energy by *CalPoEnergy()*, of the individual protein/peptide or the whole system, as well as the potential energy between the proteins/peptides and the confinement walls by *CalWlEnergy()*. *CalTemp()* is used to calculate the system temperature, and *FindPair()* and *RightPair()* are used to find the right potential step at which a pair of atoms locates.

#### HBRamach.c

The functions in this file are used to calculate the numbers of different types of hydrogen bonds. *RamachandranPlot()* is used to compute the Ramachandran plot of amino acids. If the flag *-ConMap* is set, *HBInfo()* will also produce the contact maps for both pairs of amino acids and pairs of atoms.

#### PBCAdjust.c

The functions in this file are used to remove the period boundary condition in the trajectory file.

#### REMD.c

The functions in this file are used to perform WHAM[1-2](#c71) analysis for the results from REMD simulations.

#### Tools.c

This file contains all the helper functions.

*Reference*

[1] Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A., THE weighted histogram analysis method for free-energy calculations on biomolecules. I. The method. *Journal of Computational Chemistry* **1992,** *13* (8), 1011-1021.[↩](#c71)

[2] Gallicchio, E.; Andrec, M.; Felts, A. K.; Levy, R. M., Temperature weighted histogram analysis method, replica exchange, and transition paths. *J Phys Chem B* **2005,** *109* (14), 6722-31.[↩](#c71)

sDMD

Today's picks

Was this content helpful?

Recent Feedback

sDMD

Today's picks

Was this content helpful?

Recent Feedback

Download README