ProtoMol Version 2.0.3 - User Guide

Jesús A. Izaguirre  Thierry Matthey
Department of Computer Science Department of Informatics
University of Notre Dame  University of Bergen
USA  Norway
izaguirr@cse.nd.edu  matthey@ii.uib.no

        http://www.nd.edu/~lcls

protomol@cse.nd.edu

Contributors from University of Notre Dame:

Atul Bahel

Trevor Cickovski Scott Hampton Hong Hu Qun Ma Branden Moore Thomas Slabach Jeffrey Stine George Viamontes Jeremiah Willcock

Edited by:

Jim Bilek

July 2, 2007

Abstract

This user guide details the operation and key features of the ProtoMol program. ProtoMol is designed to be an important tool for Molecular Dynamics (MD) simulations. With a simple configuration setup, compatibility with other popular MD file formats, and the ability to run on parallel systems, Protomol combines high performance and ease of use.

The user guide begins with a short introduction to the useful features of ProtoMol and follows with the necessary setup and commands needed to run the software. Also included is information on getting ProtoMol up and running on parallel systems, procedures for integration with VMD, and instructions on how to download and install the latest version.

GNU GENERAL PUBLIC LICENSE, copyrighted by the University of Notre Dame and the University of Bergen, Norway.

Acknowledgements

• Special thanks to the Center For Applied Mathematics at the University of Notre Dame for helping to fund researchers on the project.

• Parts of the development were performed at the Norwegian super-computing facilities in Bergen through a Norges Forskningsråd grant.

Contents

1 Introduction 4

1.1 Useful Features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.1.1 Compatibility with two versions of CHARMM . . . . . . . . . . . . . . . . . . . . 4 1.1.2 Force field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.1.3 Multiple time-stepping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.4 Fast electrostatic computations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.5 Comparison of force algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.6 Hybrid Monte Carlo sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.7 Ability to interact with VMD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.8 Ability to run on parallel machines . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2 Getting Started 6

2.1 Command Line . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.1 Help Option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.2 Man Option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.3 Version Option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.4 Force Option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.5 Constant Option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.6 Integrator Option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.7 Output Option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.8 Topology Option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.9 Keyword Option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.10 Units Option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.2 Configuration File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.1 Format of ProtoMol keywords in the configuration file . . . . . . . . . . . . . . . 9 2.2.2 Format of the integrator and their arguments . . . . . . . . . . . . . . . . . . . . . . 9 2.2.3 Format of forces and their arguments . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.4 Sample Three-Level Multiple Time-Stepping Integration Scheme . . . . . . . . . . 18 2.2.5 Sample of Force Comparison between PME and plain Ewald . . . . . . . . . . . . . 18 2.2.6 Sample of Force Comparison and Benchmarking PME with Different Accuracy … 19 2.2.7 ProtoMol keywords and descriptions . . . . . . . . . . . . . . . . . . . . . . . . 20

2.3 Required Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3 Input and Output File Types and Formats 26

3.1 Supported File Formats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.1.1 PDB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.1.2 XYZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.1.3 PSF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.1.4 CHARMM Parameter Files: New and Old Versions . . . . . . . . . . . . . . . . . . 27 3.1.5 DCD Trajectory Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.1.6 ProtoMol All Energy Output File . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.1.7 ProtoMol Split Energy Output File . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.1.8 BSDL File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.1.9 User Specific Output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

3.2 Input Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3 Restart Output Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.4 User-Specified Output Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.5 ProtoMol Units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.6 Output of Force Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

4 ProtoMol Configuration Files and Test Molecule Examples 31 4.1 Alanin Configuration File with Leapfrog Integrator . . . . . . . . . . . . . . . . . . . . . . 31 4.2 20 Å diameter water droplet With A Two-Level Integration Scheme . . . . . . . . . . . . . 33 4.3 BPTI with Hybrid Monte Carlo Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

5 Parallelization 37

6 Interaction of ProtoMol with VMD 38

7 Availability and Installation of ProtoMol 40 7.1 How to Download ProtoMol . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 7.2 Platforms that ProtoMol is Able to Run On . . . . . . . . . . . . . . . . . . . . . . . . . 40 7.3 Compiling ProtoMol . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 7.4 Documentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

8 The Future of ProtoMol 42

A Background 43

A.1 Integrators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 A.1.1 Standard integrators in ProtoMol . . . . . . . . . . . . . . . . . . . . . . . . . . 43 A.1.2 Distinction between multiple-time-stepping (MTS) and single-time-stepping (STS)

integrators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 A.1.3 The Verlet-I/r-RESPA Impulse Method [15] . . . . . . . . . . . . . . . . . . . . . 44 A.1.4 Mollified Impulse Methods (MOLLY ) [15] . . . . . . . . . . . . . . . . . . . . . 44

A.2 Forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 A.2.1 Bond force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 A.2.2 Angle force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 A.2.3 Dihedral and improper forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 A.2.4 Non-bonded forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 A.3 Switching Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 A.4 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 A.5 Cell Manager . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

B ProtoMol Molecular Dynamcis Software License 52 B.1 Conditions and Regulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 B.2 Contact Information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

List of Figures

4.1 Picture of alanin, generated by the BOOGA application raytrace . . . . . . . . . . . . . . 33 4.2 Picture of 10 Å diameter water droplet, generated by the BOOGA application raytrace . . 35

6.1 PHANTOM Desktop^{T M} Haptic Device. . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

A.1 Partition of cells, where cell no. 5 is the origin and 1-4, 6-9 the neighbor cells. Cell no. 6, 3, 2 and 1 describe the interaction path. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

Chapter 1

Introduction

Molecular Dynamics (MD) describes a molecular system as a function of time based on integration of equations of motion and interacting forces. The most computationally expensive part is the force calculation among the atoms. There have been many implementations to solve one given problem very efficiently, but they usually have a lack of flexibility when trying to implement new methods or approaches to solve the problem. What ProtoMol provides is a generic, object-oriented component framework for MD simulations. To meet these high performance expectations, ProtoMol uses cell algorithms, grid techniques and well-optimized libraries to challenge the most computationally expensive forces. The design of ProtoMol also includes parallelization, based on components to distribute the work and data. The approach follows an incremental and partial parallelization scheme, which allows the developer to start with a sequential implementation and then do step by step parallelization.

The overall framework of ProtoMol is designed for non-bonded, bonded, short-range and long-range forces for systems with tens of thousands of atoms representing water and several large molecules. It is designed for high flexibility, ease in extension and maintainence, and high performance demands.

1.1 Useful Features

1.1.1 Compatibility with two versions of CHARMM

CHARMM is a program developed for the purposes of conducting MD simulations (for more information, see http://www.lobos.nih.gov/Charmm ). Files generated by CHARMM can provide initial information that ProtoMol needs to run. However, there have been several new versions of CHARMM released throughout its history, and thus some different file formats. ProtoMol is compatible with both CHARMM version 19 and the most recent version as of this release, CHARMM version 28a2. By covering both of these versions, ProtoMol covers all possible formats of the files that it needs and thus provides more flexibility for its users.

1.1.2 Force field

ProtoMol supports force fields in the same format that CHARMM uses, including bonded interactions between groups of 2 (bonds), 3 (angles), and 4 (dihedrals and impropers) atoms, as well as electrostatic and van der Waals nonbonded interactions.

1.1.3 Multiple time-stepping

Multiple time-stepping is a useful technique for cost reduction of calculating long-range electrostatic forces. The idea behind multiple time-stepping is to compute at every timestep bonded, van der Waals and shortrange electrostatic forces, while computing long-range electrostatic forces less often. Overtime this will result in improved efficiency because the cost for computing electrostatic interactions will be amortized over the number of timesteps to be run in the simulation [4].

1.1.4 Fast electrostatic computations

Efficient electrostatic force evaluation algorithms such as plain Multi Grid (MG) [5, 25, 20], Ewald [10, 7] and Particle Mesh Ewald (PME) [8, 13, 6, 9] are supported by ProtoMol . The plain Ewald, for example, is able to compute electrostatic forces with a complexity of O(N^{3/2}), and even PME algorithm has a complexity of O(NlogN) in computing these forces, where MG has complextiy of O(N). These algorithms are improvements over the usual quadratic complexity [4] for the direct method. For more details can be found in [20].

1.1.5 Comparison of force algorithms

ProtoMol is able to compare forces measuring the error and the time; e.g. an exact algorithm with a fast approximation.

1.1.6 Hybrid Monte Carlo sampling

The Hybrid Monte Carlo algorithm is one of the most popular techniques of molecular sampling used today. ProtoMol has been made to support this algorithm. The benefits of this algorithm over many others, including conventional Monte Carlo, include: conservation of energy is not required, longer timesteps are used, and global update of the configuration file is allowed with a reasonable acceptance ratio.

1.1.7 Ability to interact with VMD

VMD is a program developed by the University of Illinois that is used for displaying large biomolecular systems in three dimensions. ProtoMol is conveniently able to interact with this useful program. For more information, see:

http://www.ks.uiuc.edu/Research/vmd

1.1.8 Ability to run on parallel machines

The most computationally expensive part of an MD simulation is force evaluation between atoms. Because so much computing power is required for the force evaluation of large systems, i.e., 5000 atoms or more, the ability to run MD simulations on parallel machines will hold a large advantage with regard to speed. ProtoMol has been designed to take advantage of parallel computing in running MD simulations.

Getting Started

This chapter will introduce the commands and settings needed to run ProtoMol . Included are the exact formats of the command line on a UNIX machine and the configuration file containing all initial information to run the MD simulation. In Chapter 4 we will show three sample configuration files.

2.1 Command Line

The MD application of the ProtoMol framework has conveniently been named protomol . At a UNIX prompt, a user types protomol followed by an alternating list of keywords and arguments (see section 2.2.7 for a list of keywords). Thus, the general format for the ProtoMol execution command is the following:

protomol [--keyword1 value1 ] [--keyword2 value2 ] [--keyword3 value3 ] .......

Note that keywords must be preceded by two dashes, where a value can be a list of values - e.g. for the keyword CellBasisVector1. Also note that any keyword-value pair specified on the command line overrides any according pair in the configuration file.

The user may also specify any of these keywords and values in the configuration file that he or she is using. However, at least the exact pathname of the configuration file being used for running ProtoMol must be specified on the command line. For example, the configuration file can be specified one of the following ways:

   protomol [--config] <pathname/myconfigfilename> ......

The --config keyword can be omitted if the configuration file is the first thing specified.

      protomol .....  --config <pathname/myconfigfilename> .....

2.1.1 Help Option

A user may type one of the following two possibilities at a UNIX prompt for help with the ProtoMol command line:

1. protomol -h
2. protomol --help

2.1.2 Man Option

ProtoMol supports man alike command to query options and meaning of a keyword:

1. protomol -m <keyword>
2. protomol --man <keyword>

All registered keywords of ProtoMol :

1. protomol -m --keywords
2. protomol --man --keywords

2.1.3 Version Option

The actual version of ProtoMol :

1. protomol -v
2. protomol --version

2.1.4 Force Option

Output of all actually supported forces in ProtoMol :

1. protomol -f
2. protomol --forces

2.1.5 Constant Option

Output of numerical constants :

1. protomol -c
2. protomol --constants

2.1.6 Integrator Option

Output of all actually supported integrators in ProtoMol :

1. protomol -i
2. protomol --integrators

2.1.7 Output Option

Output of all actually supported output possibilities in ProtoMol :

1. protomol -o
2. protomol --outputs

2.1.8 Topology Option

Output of all actually supported output topologies and theis parameters in ProtoMol :

1. protomol -t
2. protomol --topologies

2.1.9 Keyword Option

Output of all actually supported keywords their default values and aliases in the configuration file:

1. protomol -k
2. protomol --keywords

2.1.10 Units Option

Output of ProtoMol units :

1. protomol -c
2. protomol --constants

2.2 Configuration File

The configuration file is a text file containing a collection of keyword-value pairs specifying the simulation configuration, I/O files and formats, and the definition of the integrator scheme. A comment starts with the pound sign (#) and anything on the same line after the # is ignored.

2.2.1 Format of ProtoMol keywords in the configuration file

The configuration file format for ProtoMol is quite simple, making it convenient for creation and modification of the file. The advantage of straightforward modification of the configuration file is that the user can very easily switch between running the same molecule under different initial conditions, or even switch to a different molecule without much trouble. The general format for a ProtoMol configuration file is a list of keywords and values, with whitespace between each keyword and value and a newline between each new keyword-value pair:

keyword1 value1

keyword2 value2

keyword3 value3 . . .

A list of ProtoMol keywords, possible values and defaults can be found in section 2.2.7

2.2.2 Format of the integrator and their arguments

In addition to the keyword-value pairs in the configuration file, one must set up an integrator in the following manner:

Integrator {

level N-1 <integrator type>{ <integrator arguments> (These will differ depending on the integrator type). <integrator forces> (These are all optional).

} . . .

level 0 <integrator type>{ . . .

} }

Note that the order of definition for each level is not strict, but ProtoMol expects one definition for each level.

Integrator Types

• Multiple Timestep Integrators

1. BSplineMOLLY

2. DihedralHMC

3. DihedralLiftMC

4. EquilibriumMOLLY

5. HBondMOLLY

6. HybridMC (Hybrid Monte Carlo Integrator)

7. Impulse

8. ShadowHMC

9. Umbrella

• Single Timestep Integrators

1. BBK

2. DMDLeapfrog

3. LangevinImpulse

4. Leapfrog (Velocity Leap-Frog Integrator)

5. NoseNVTLeapfrog

6. NPTVerlet

7. PaulTrap

8. PLeapfrog (Position Leap-Frog Integrator)

Integrator Argument Types

1. Impulse cyclelength <length of cycle (integer)>

2. HybridMC:

cyclelength <length of cycle (integer)> warmupcycles <number of warm-up cycles (integer)> temperature <Kelvin temperature (float)> randomCycLen <Use a random adjustment to cyclelength (bool)>

3. EquilibriumMOLLY cyclelength <length of cycle (integer)>

4. BSplineMOLLY cyclelength <length of cycle (integer)> BSplineType <superShort|extraShort|short|shortLinear|long|longLinear|longQuadratic> BSplineMollyForceType <Bond|Angle|Bond+Angle>

5. HBondMOLLY cyclelength <length of cycle (integer)> BSplineType <superShort|extraShort|short|shortLinear|long|longLinear|longQuadratic> HBondMollyForceType <Bond+Angle+LJ|Bond+Angle+LJ+Coulomb> mollifyNonbondedRangeFrom <starting point (float)> mollifyNonbondedRangeTo <ending point (float)> SWon <point (float)> cutoff <cutoff point (float)>

6. Leapfrog timestep <length of step (float)>

7. PLeapfrog timestep <length of step (float)>

8. NoseNVTLeapfrog timestep <length of step (float)> temperature <Kelvin temperature (float)> thermal <thermal inertia (float)>

9. BBK timestep <length of step (float)> temperature <Kelvin temperature (float)> gamma <gamma (float)>

seed <random seed (integer)>

10. LangevinImpulse

timestep <length of step (float)>

temperature <Kelvin temperature (float)>

gamma <gamma (float)>

seed <random seed (integer)>

11. DMDLeapfrog

cyclelength <length of cycle (integer)>

temperature <Kelvin temperature (float)>

gamma <gamma (float)>

numIter <number of iterations (integer)>

General Format of Integrator Forces

force <force1 type>

<force1 arguments> force <force2 type>

<force2 arguments> force <force3 type>

<force3 arguments> …

Force Comparison and Timing

At present, pairs of forces can be compared to determine energy and force errors of new force methods. The comparison is performed on-the-fly, such that the reference force does not affect the current simulation..

force compare <force approximation type>

<force approximation arguments> force compare <force exact type>

<force exact arguments>

For the benchmarking of forces a timer, function is provided to measure the total and average time spent in dedicated force methods.

force time <force to benchmark type>

Both comparisons functions can be nested to evaluate accuracy and run-time performance simultaneously. force compare time <force approximation type>

<force approximation arguments> force compare time <force exact type>

<force exact arguments>

2.2.3 Format of forces and their arguments

1. Bonded Forces

(a) Bond - force on each of two atoms caused by a linear bond between them.

Arguments: None

(b) Angle - force on each of three atoms caused by an angular bond between them.

Arguments: None

(c) Dihedral - force on each of four atoms caused by the interaction between them. In an arrangement of four bonded atoms A-B-C-D, rotation occurs around the center covalent bond (B-C) - this is what dihedrals describe. [3] Arguments: None

(d) Improper - force on each of four atoms caused by the interaction between them. Impropers are like dihedrals but the atoms are not necessarily covalently attached in the pattern given. [3] Arguments: None

2. Nonbonded Forces

(a) LennardJones - van der Waals forces Arguments:

1. -algorithm < NonbondedFull| (Direct method over several images (cutoff), only for periodic boundary condtions) NonbondedSimpleFull| (Direct method, no switching function required) NonbondedCutoff>(Uses a cutoff)

2. -switchingfunction < [Complement]Cutoff| (Simple cutoff) [Complement]Shift| (Cutoff and shift) [Complement]C2 (Continuous second derivative)| [Complement]C1 (Continuous first derivative)| Universal>(This is the default)

3. -cutoff <cutoff point (float)>(This is not necessary if using a universal switching function or NonbondedSimpleFull)

4. -switchon <switchon (float)>(Only for C2 switching function)

5. -blocksize <integer (default is 64)>(This is not necessary if using NonbondedCutoff)

(b) Coulomb - interaction of particles due to electric charge Arguments: Same as LennardJones

Note: for optimization purposes, LennardJones and Coulomb forces can be combined if the same algorithm and switching function is desired, as shown in the sample integration scheme on the following page. ProtoMol may not recognize all possible combination since they may already be covered by others, e.g -algorithm NonbondedSimpleFull and an arbitrary switching function (not Universal) is equivalent (and faster) with NonbondedCuttoff. To see what sort of forces are supported

type protomol -f .

3. Fast electrostatic Forces

(a) Coulomb - Coulomb force using a plain Ewald summation Arguments:

1. -algorithm FullEwald

2. -real (Real-space part of the Ewald sum)

3. -reciprocal (Reciprocal-space part of the Ewald sum)

4. -correction (Point self-, intra-molecular self-, charged system and surface dipole (omitted for the moment) part of the Ewald sum)

5. -alpha <α (float)>(Splitting parameter of the Ewald sum (default: optimal splitting))

6. -accuracy <accuracy (float)>(Splitting parameter of the Ewald sum (default 1e − 5))

7. -j <expansion factor (float)>(This is only necessary for vacuum boundary conditions (default 3)

8. -switchingfunction < [Complement]Cutoff| (This is the default, simple cutoff) [Complement]Shift| (Cutoff and shift) [Complement]C2 (Continuous second derivative)| [Complement]C1 (Continuous first derivative)

9. -switchon <switchon (float)>(Only for C2 switching function)

(b) Coulomb - Coulomb force using the Particle Mesh Ewald method Arguments:

1. -algorithm PMEwald

2. -real (Real-space part of the Ewald sum)

3. -reciprocal (Reciprocal-space part of the Ewald sum)

4. -correction (Point self-, intra-molecular self-, charged system and surface dipole (omitted for the moment) part of the Ewald sum)

5. -interpolation <BSpline|Hermite|Lagrange>

6. -gridsize <nx (integer)><ny (integer)><nz (integer)>

7. -cutoff <cutoff (float)>(Cutoff of the Real-space term of the Ewald sum)

8. -order <order (integer)>(Interpolation order (default 4 (=qubic)))

9. -accuracy <accuracy (float)>(Splitting parameter of the Ewald sum (default 1e − 6))

10. -alpha <α (float)>(Splitting parameter of the Ewald sum (default: optimal splitting))

11. -j <expansion factor(float)>(This is only necessary for vacuum boundary conditions (default 3))

12. -switchingfunction < [Complement]Cutoff| (This is the default, simple cutoff) [Complement]Shift| (Cutoff and shift) [Complement]C2 (Continuous second derivative)| [Complement]C1 (Continuous first derivative)

13. -switchon <switchon (float)>(Only for C2 switching function)

(c) Coulomb - Coulomb force using multigrid method Arguments:

1. -algorithm MultiGrid

2. -direct (Direct part)

3. -smooth (Smooth part)

4. -interpolation <BSpline|Hermite|Lagrange>

5. -kernel <C1|C2|C3|C4>

6. -toplevel <nx (integer)><ny (integer)><nz (integer)>(Number of grid points of the coarsest grid, only used for periodic boundary conditions)

7. -levles <number of levels (integer)>(>0)

8. -s <softening distance (float)>

9. -order <order (integer)>(Interpolation order (default 4 (=qubic)))

10. -ratio <grid ratio fine-coarse (float)>

11. -h <x (float)><y (float)><z (float)>(Meshsize of the finest grid, only used for vacuum boundary conditions)

12. -origin <x (float, default 0)><y (float, default 0)><z (float, default 0)>(Origin of the finest grid, only used for vacuum boundary conditions)

4. Other Forces

(a) Haptic - haptic force Arguments:

1. -port <port id (integer)>

2. -trate <transmission rate (integer)>

3. -timeout <timeout value (integer)>

4. -step inc <step increment value (integer)>

5. -wait <wait value (integer)>

(b) Friction - friction force Arguments:

1. -k <friction constant (float)>

2. -rnd <friction random term, default 0 (float)>

(c) PaulTrap - Paul trap force[12, 24] Arguments:

1. -omegaR <ω_{r} (float) [fs^{−1}]>

2. -omegaZ <ω_{z} (float) [fs^{−1}]>

3. -alpha <α (float, default is 0.0) >

4. -r0 <r_{0} (float, default is 0.0) [Å]>

(d) Gravitation - gravitation force Arguments:

1. -algorithm <NonbondedSimpleFull>

2. -G <gravitation constant (float)>

3. -blocksize <integer (default is 64)> (e) MagneticDipole - magnetic dipole force

Arguments:

1. -algorithm <NonbondedSimpleFull| (Direct method, no switching function required) NonbondedCutoff>(Uses a cutoff)

2. -switchingfunction < [Complement]Cutoff| (Simple cutoff) [Complement]Shift| (Cutoff and shift) [Complement]C1 (Continuous first derivative)| Universal>(This is the default)

3. -cutoff <cutoff point (float)>(This is not necessary if using a universal switching function or NonbondedSimpleFull)

4. -chi <χ (float)>(Effective susceptibility of the spheres)

5. -r <radius (float)>(Radius of the spheres)

6. -omega <ω (float)>(Angular frequency of the spheres)

7. -phi <ϕ (float)>(Initial angel of the field)

8. -H <x (float)><y (float)><z (float)>(Magnetic flux-density in x, y and z-direction)

9. -d <(float)>(Separation between boundaries (glass-plates, 0.0 turns mirroreffects off)

10. -blocksize <integer (default is 64)>(This is not necessary if using NonbondedCutoff)

(f) MagneticDipoleMirror - magnetic miror force Arguments:

1. -chi <χ (float)>(Effective susceptibility of the spheres)

2. -r <radius (float)>(Radius of the spheres)

3. -H <x (float)><y (float)><z (float)>(Magnetic flux-density in x, y and z-direction)

4. -d <(float)>(Separation between boundaries (glass-plates, 0.0 turns mirroreffects off)

(g) Spherical - spherical force/ boundary conditions Arguments:

1. -center <x (float)><y (float)><z (float)>(Center of the spherical boundary )

2. -radius <radius (float)>(Radius of the spherical boundary)

3. -k <force (float)>(Force constant)

4. -j <exponent (integer)>(Exponent for the potential function)

2.2.4 Sample Three-Level Multiple Time-Stepping Integration Scheme

Integrator {
level 2 BSplineMOLLY {
cyclelength 2
BSplineType long
BSplineMollyForceType Bond+Angle
force Coulomb
-algorithm NonbondedSimpleFull
-switchingfunction ComplementC1
-cutoff 6.5
-blocksize 64
}
level 1 Impulse {
cyclelength 3
force Bond
force Angle
force LennardJones Coulomb
-algorithm NonbondedCutoff
-switchingfunction C2
-cutoff 7.5
-switchon 1.0
-blocksize 32
}
level 0 LeapFrog {
timestep 1
force Bond
force Angle
force Dihedral
force Improper
}
}

2.2.5 Sample of Force Comparison between PME and plain Ewald

force compare Coulomb
-algorithm PMEwald
-real -reciprocal  -correction
-interpolation BSpline
-gridsize 8 8 8
-order 4
-cutoff 6.5
force compare Coulomb
-algorithm FullEwald
-real -reciprocal -correction

2.2.6 Sample of Force Comparison and Benchmarking PME with Different Accuracy

force compare time force Coulomb -algorithm PMEwald -real -reciprocal
         -correction -interpolation BSpline  -cutoff 6.5 -gridsize 10 10 10
force compare time force Coulomb -algorithm PMEwald -real -reciprocal
         -correction -interpolation BSpline  -cutoff 6.5 -gridsize 20 20 20
force time LennardJones -algorithm NonbondedCutoff -switchingFunction C1
                           -cutoff 8.0

2.2.7 ProtoMol keywords and descriptions

Keyword and Aliases Type

Description

inputfileprefix string Contains the full pathname of the file prefix for all initial data files. If this is specified, there is no need to specify any initial data files as long as each of them can be defined as this file prefix concatenated with the appropriate flag (.psf, .par, .pdb, .xyz).

configfile string Contains the full pathname of the configuration file to be used in the simulation. Note: this keyword would only be used on the command line.

posfile / coordinates string Contains the full pathname of the initial positions file. ProtoMol supports PDB and XYZ formatted position files. Either this or the inputfileprefix must be specified.

velfile string Contains the full pathname of the initial velocities file. Once again, ProtoMol supports PDB and XYZ formats. If no initial velocity file is specified, random velocities are generated based on the initial temperature and a seed that can also be specified. Therefore one of either the velfile, inputfileprefix, or temperature needs to be specified. A seed is optional.

psffile string Contains the full pathname of the initial topology file in PSF format. Either this or the inputfileprefix must be specified.

parfile / parameters string Contains the full pathname of the initial CHARMM parameter file. Either this or the inputfileprefix must be specified.

usecharmm28parfile boolean (default: false) Specifies if the initial parameter file is of new CHARMM format (after CHARMM22) or old format. CHARMM28 is the most recent version of CHARMM.

numsteps integer Specifies the number of steps for the simulation to run.

firststep

integer (default: 0) Specifies the number of the initial timestep.

exclude “none”, “1-2”, “1-3”, “1-4”, “scaled1-4” (default: “1-3”)

Specifies nonbonded exclusions

temperature float

Specifies the initial Kelvin temperature.

one4scaling / 1-4scaling float (default: 1.0) Value for 1-4 scaling, if necessary. seed unsigned integer Random number seed for velocity genera-

tion. This is only used if an initial velocity is not specified, and in that case if the seed is specified as 0, it defaults to the timer.

boundaryconditions string ProtoMol supports vacuum or periodic boundary conditions, specified with “Vacuum” or “Periodic”, respectively.

Keyword and Aliases Type

Description

cellbasisvector1 x y z (x, y, and z are floats) Basis vector 1, for periodic boundaries cellbasisvector2 x y z (x, y, and z are floats) Basis vector 2, for periodic boundaries cellbasisvector3 x y z (x, y, and z are floats) Basis vector 3, for periodic boundaries cellorigin x y z (x, y, and z are floats) Center of periodic cell, for periodic bound-

aries

cellmanager string Cell manager - current ProtoMol only supports a cubic cell manager, specified by “Cubic”. The reason why it needs to be specified even though there is only one option is to allow for future flexibility.

cellsize float Specifies the size of cell used by the cell manager.

hgroupmincutoff  float (default: 2.5)  Specifies minimal distance for a H-group
hgroupmaxcutoff  float (default: 4.0)  Specifies maximal distance for a H-group

useshadow

boolean (default: false) Specifies if shadow Hamiltonian is applied.

extendedcoordinates boolean (default: false) Specifies if the user likes position outputs in extended coordinates, per default the minimal image convention is applied.

commotion boolean Specifies if the user would like center of mass motion removed from the velocities.

restartfile string Contains the file prefix of the restart output files, if the user desires them. This option allows a configuration file (prefix.conf), a PDB positions file (prefix.pos.pdb), a PSF topology file (prefix.psf) and a PDB velocities file (prefix.vel.pdb) to be printed out on several occasions throughout the simulation, the exact frequency specified by restartfreq (below). This way the user can view more closely how the molecule changes throughout the simulation, rather than simply viewing initial and final data.

restartfreq integer (default: 24000) Specifies the frequency in timesteps for

restart files to be written, if desired. This applies only if restart files are desired, but must be specified if they are.

dorestartfiles boolean (default: false) Specifies if the user desires restart files to be written.

Keyword and Aliases Type

Description

xyzposfile string Contains the full pathname of the positions trajectory file in XYZ format, once again if desired.

doxyzposfile boolean (default: false) Specifies if the user would like an XYZ tra-

jectory file for positions to be written.

xyzposoutputfreq integer (default: 1) Specifies the frequency in outputfreq of the XYZ trajectory file for positions to be written, if desired.

xyzvelfile string Contains the full pathname of the XYZ velocities trajectory file, if desired.

doxyzvelfile boolean (default: false) Specifies if the user would like an XYZ ve-

locities trajectory file to be written.

xyzveloutputfreq integer (default: 1) Specifies the frequency in outputfreq of the XYZ trajectory file for velocities to be written, if desired.

dcdfile string Contains the full pathname of the DCD trajectory file to be written, if desired.

dodcdfile boolean (default: false) Specifies if the user would like a DCD tra-

jectory file to be written.

dcdoutputfreq integer (default: 1) Specifies the frequency in outputfreq of the DCD trajectory file to be written, if desired.

xyzforcesfile string Contains the full pathname of the XYZ forces trajectory file to be written, if desired.

doxyzforcesfile boolean (default: false) Specifies if the user would like an XYZ

forces file to be written.

xyzforcesoutputfreq integer (default: 1) Specifies the frequency in outputfreq of the

force trajectory file to be written, if desired.

allenergiesfile string Contains the full pathname of the energies file to be written, again if desired.

doallenergiesfile boolean (default: true) Specifies if the user would like an energies file to be written, with all energies (bond, angle, dihedral... whatever was forced in the integrator section) in one file.

splitenergiesfile string Contains the full pathname of the split energies file prefix, if split energy files are desired. In this case, there will be different files generated for different energies attached with the appropriate flag (bond = .bond.dat, angle = .angle.dat, dihedral = .dihedral.dat, improper = .improper.dat, coulomb = .coulomb.dat, Lennard Jones = .lennardjones.dat, total = .total.dat, kinetic = .kinetic.dat, potential = .potential.dat, temperature = .temperature.dat).

dosplitenergiesfile boolean (default: false) Specifies if the user desires split energy output files.

Keyword and Aliases Type

Description

bsdlfile string Contains the full pathname of the BSDL output files, if desired.

dobsdlfile

boolean (default: false) Specifies if BSDL output files are desired.

bsdloutputfreq integer (default: 1) Specifies the frequency in outputfreq of the

BSDL files to be written, if desired.

dobsdlshowwater boolean (default: true) Specifies if water molecules should be sup-

pressed for the output of the BSDL files.

dobsdlminimalimage boolean (default: false) Specifies if boundary conditions are ap-

plied for the output of the BSDL files.

dobsdlfixedsize boolean (default: true) Specifies if the output of the BSDL files

should have same scene view.

outputfreq integer Specifies the frequency in timesteps for the writing of energy data to the console.

samplefile string Contains the full pathname of the sample file to be written, if desired. The sample file is an example output of a user implemented output routine. The actual implementation computes the diffusion and the volume.

dosamplefile boolean (default: false) Specifies if the user would like a sample

file to be written.

sampleoutputfreq integer (default: 1) Specifies the frequency in outputfreq of the

sample file to be written, if desired.

paulfile string Contains the full pathname of the Paul trap output file to be written, if desired. The output contains the kinetic energy, temperature, energy difference of Coulomb minus twice the Paul trap, and a histogram in function of the distance (position) to the origin.

dopaulfile boolean (default: false) Specifies if the user would like a Paul trap

output file to be written.

pauloutputfreq integer (default: 1) Specifies the frequency in outputfreq of the Paul trap output file to be written, if desired.

paulhistsize integer (default: 15) Specifiesthe size of the histogram. momentumfile string Contains the full pathname of the momentum output file to be written, if desired. The output contains the momentum for each dimension..

domomentumfile boolean (default: false) Specifies if the user would like a momen-

tum output file to be written.

momentumoutputfreq integer (default: 1) Specifies the frequency in outputfreq of the tmomentum file to be written, if desired.

Keyword and Aliases Type

Description

temperaturefile string Contains the full pathname of the temperature output file to be written, if desired. The output contains the kinetic energy and temperature for the total system and for all oxygen.

dotemperaturefile boolean (default: false) Specifies if the user would like a tempera-

ture output file to be written.

temperatureoutputfreq integer (default: 1) Specifies the frequency in outputfreq of the temperature output file to be written, if desired.

finpdbposfile string Contains the full pathname of the final positions file in PDB format, if the user desires.

dofinpdbposfile boolean (default: false) Specifies if the user would like a final PDB

positions file to be written.

finxyzposfile string Contains the full pathname of the final positions file in XYZ format, if the user desires.

dofinxyzposfile boolean (default: false) Specifies if the user would like a final XYZ

positions file to be written.

finpdbvelfile string Contains the full pathname of the final PDB velocities file, if so desired.

dofinpdbvelfile boolean (default: false) Specifies if the user would like a final PDB

velocities file to be written.

finxyzvelfile string Contains the full pathname of the final XYZ velocities file, if desired.

dofinxyzvelfile boolean (default: false) Specifies if the user would like a final XYZ

velocities file to be written.

usebarrier boolean (default: true) Specifies if explicit synchronization (bar-

rier) is desired before global communication.

masterslave integer (default: 3) Specifies if the user would like to use

master-slave or static distribution. The integer defines the minimal number of nodes required for master-slave distribution. A negative number forces the usage of static distribution for any number of nodes.

parallelpipe integer (default: 0) Specifies if the depth of the pipe to assign work to each slave.

• Type boolean can be specified as yes/no or true/false.

• Xfile defines the path and the filename of the output X, where doXfile is the flag to turn on/off the output. The default of doXfile is true, if Xfile is defined, otherwise false.

• Xoutputfreq defines the output frequency in relation to outputfreq .

2.3 Required Parameters

As mentioned earlier, the command line must contain either the full pathname of the configuration file or an input file prefix. These are the only restrictions specifically applied to the command line.

The following parameters MUST be specified on either the command line or in the configuration file specified on the command line (you may view Chapter 3 for a list of supported ProtoMol files and their formats):

• One of either an initial positions file in PDB or XYZ format, or an input file prefix.

• One of an initial velocities file in PDB or XYZ format, an input file prefix, or an initial temperature.

• One of either a PSF topology file or an input file prefix.

• One of either a CHARMM Par parameters file or an input file prefix.

• If the CHARMM Par parameters file follows a version of CHARMM after version 22, USECHARMM28PARFILE must be set to true.

• If the user specifies any output files that they want written, an output frequency must be specified.

• If the user specifies that restart files are desired, both a restart frequency and restart file prefix must be specified.

• If the user specifies any specific output files that should be written, they must specify a filename.

• The number of steps for the simulation.

• A cubic cell manager and either vacuum or periodic boundary conditions.

• An integrator.

Input and Output File Types and Formats

3.1 Supported File Formats

3.1.1 PDB

ProtoMol supports position and velocity files written in Protein Data Bank (PDB) format, both for input and output. This format is one of the most commonly used by other MD programs as well. The site http://www.rcsb.org/pdb/docs/format/pdbguide2.2/part 62.html gives a full description of the PDB file.

3.1.2 XYZ

In addition to PDB files, ProtoMol also supports XYZ position and velocity files for input and output. There is no universal format for XYZ files. ProtoMol will be compatible with XYZ files of the following format: # of atoms

<one comment>

atom1 x y z atom2 x y z atom3 x y z . . . .

Full description of the XYZ format file:

http://hackberry.chem.trinity.edu/IJC/Text/xmolxyz.html

3.1.3 PSF

ProtoMol supports topology files in PSF P roteinStructureF ile format for both input and output. PSF files can be generated by CHARMM, or they can be built by hand. Generally speaking, PSF files contain information about atoms, bonds, angles, dihedrals, impropers, hydrogen donors, hydrogen acceptors, and nonbonded parameters. Information on PSF files, including their structure and how to generate, create and modify them, can be found at the following site:

http://www.sinica.edu.tw/~scimath/msi/insight2K/
charmm principles/Ch02 model build.FM5.html#444511

3.1.4 CHARMM Parameter Files: New and Old Versions

CHARMM parameter files are also supported by ProtoMol for both input and output. There have been several versions of CHARMM released throughout the past two decades. With regard to the parameter file, there have been two main formats that we will refer to as “old” and “new” CHARMM parameter formats. The “old” format was abolished after CHARMM22 and was changed to the “new” format. More information on this can be found at the CHARMM website: http://www.lobos.nih.gov/Charmm. Future CHARMM versions may require updates if the parameter file format changes. However, chances are the parameters will stay the same. Thus it is convenient to know what exactly a parameter file holds. There are seven sections of a CHARMM parameter file, one each for bonds, angles, dihedrals, impropers, nonbondeds, nbfix and hydrogen bonds. Here are what each individual structure should contain as data members:

bonds : bond number, two atom numbers involved, force constant and distance

angles : angle number, three atom numbers involved, force constant, angle value, and the Urey-Bradley constants if they exist.

dihedrals : dihedral number, four atom numbers involved, multiplicity (old), force constant, periodicity, and phase shift.

impropers : improper number, four atom numbers involved, multiplicity (old), force constant, periodicity, and phase shift.

nonbondeds : nonbonded number, atom number, polarizability (if applicable), epsilon (if applicable), number of effective electrons (if applicable), radius, and room for van der Waals parameters if they exist. In the new CHARMM version, whether or not polarizability exists depends on if the value read after it is negative or positive - if the next value is negative, polarizability does not exist and the first value is ignored, and the second value is epsilon. If the next value is positive, the first term is polarizability and the second term is the number of electrons. However, future formats may not be like this, and so it helps to know what needs to be stored and under what conditions.

nbfix : number, two atom numbers involved, epsilon, sigma, epsilon14, sigma14

hydrogen bonds : bond number, two atom numbers involved, emin and rmin

3.1.5 DCD Trajectory Files

ProtoMol is able to write DCD trajectory files throughout a simulation. A DCD file is binary and contains coordinate information for the atoms in the molecule being simulated. The format of a DCD file can be found at the following site:

http://www.arl.hpc.mil/PET/cta/ccm/software/fmd/Docs/html/fmd ug.html#SEC9.

3.1.6 ProtoMol All Energy Output File

A ProtoMol all energy file includes all different types of energy in the same file. Each row in the all energies file represents one timestep. There may be several rows of data depending on the output frequency specified in the configuration file, since that specifies the number of timesteps between data writes.

The format of one row is as follows, with each entry separated by whitespace:

<time><potential energy><kinetic energy><total energy><Kelvin temperature><bond energy><angle energy><dihedral energy><improper energy><van der Waals energy><electrostatic energy><other/user defined energy types>

3.1.7 ProtoMol Split Energy Output File

If ProtoMol split energy files are desired, ten separate files will be generated, one for each different type of energy shown above in the ProtoMol All Energy Output File. Each file will have the same format, composed of rows of the following depending on the output frequency:

<time><value>

where <value>is the value of the appropriate energy parameter at the given time.

3.1.8 BSDL File

A BSDL (BOOGA Scene Description Language) file provides a three-dimensional scene description of the simulation. BOOGA, which stands for “Berne’s Object-Oriented Graphics Architecture” is an object-oriented graphics framework aimed at a wide range of application areas within the domain of computer graphics, providing data abstractions and mechanisms for two and three-dimensional objects as well as including a component layer to model high-level operations that generate or manipulate graphical objects of different types. [1]

BSDL is the language supported by the BOOGA framework and is used to define two or three-dimensional scenes, which is very close the definition for rayshade or similar tray tracer. Most of the BOOGA-applications are capable of reading or writing BSDL files. For more information, see http://www.iam.unibe.ch/~booga .

3.1.9 User Specific Output

To give the user an idea how he/she can implement their own output routines, ProtoMol provides a simple example, which in the actual implementation computes the diffusion.

The user may modify this example or make is his/her own implementation, which involves the source files

framework/frontend/OutputData.h  , framework/frontend/OutputData.C  and
framework/frontend/SimulationConfiguration.C  .
   grep -i sample framework/*/*.[Ch]  will point out the locations.

3.2 Input Files

Input files are specified on the command line and/or the configuration file. The following guidelines need to be followed:

Initial Positions File: PDB or XYZ format
Initial Velocities File: PDB or XYZ format
Topology File:  PSF format
Parameter File:  CHARMM Par format

3.3 Restart Output Files

If the user desires restart files to be printed throughout the simulation, here is what will be printed every certain number of timesteps, given by the restart frequency:

Positions File:  PDB format - filename <restartprefix.pos.pdb>
Velocities File:  PDB format - filename <restartprefix.vel.pdb>
Topology File:  PSF format - filename <restartprefix.psf>
Configuration File: ProtoMol configuration file format - filename <restartprefix.conf>

3.4 User-Specified Output Files

There are several supported types for output files, but the advantage here is that the user can have the same information printed in multiple formats if multiple formats are supported for this information, unlike the input files where there could only be one file for each different initial information type. the following are the supported formats:

Positions File:  PDB and XYZ format
Velocities File:  PDB and XYZ format
Positions Trajectory File:  XYZ and DCD format
Velocities Trajectory File: XYZ format
Forces File:  XYZ format
All/Split Energies File:  ProtoMol energies file format
BSDL File:  BSDL format
Sample File:  Example of user defined output

3.5 ProtoMol Units

Quantity  ProtoMol Unit SI  Comment
Time  fs  10^{−15} s
Length  Å  10^{−10} m
Velocity  10^{3}\sqrt{4.184} m s^{−1} PDB factor^{1}: 20.45482706
Energy  kcal mol^{−1}  \frac{4184}{6.022045e+23} J
Force  kcal mol^{−1} Å^{−1} \frac{4184}{6.022045e+13} N
Mass  amu  1.6605e − 27 Kg
Temperature K  K
Charge  e  1.602e − 19 C

1 Velocities in PDB files are scaled to get a more accurate representation since the PDB format has a limited representation of floating numbers.

3.6 Output of Force Comparison

ProtoMol provides the option to compare forces pairwise measuring the error and the time spent for each force. The output is computed as followed:

Name  Definition  Comment
absF2 max_{i} ||F_{i} − F̃_{i}||
rFavg  \begin{matrix}∑ & _{i} & \sqrt{\frac{||F_{i}−F̃_{i}||^{2}}{m_{i}}} \\ ∑ & _{i} & \sqrt{\frac{||F_{i}||^{2}}{m_{i}}} \end{matrix}
rFmax N \begin{matrix}max & i & \sqrt{\frac{||F_{i}−F̃_{i}||^{2}}{m_{i}}} \\ ∑ & _{i} & \sqrt{\frac{||F_{i}||^{2}}{m_{i}}} \end{matrix} N: number of atoms
rPE |\frac{E_{p}−Ẽ_{p}}{E_{p}} | E_{p}: Potential Energy

where m_{i} is the mass of atom i; F is the exact algorithm/ reference and F̃ is the approximation.

ProtoMol Configuration Files and Test Molecule Examples

The ProtoMol source release includes a folder of example simulation configurations and select test molecules. Examples are organized by test molecule and each folder contains one or more simulation configuration files. To provide the user with a uniform set of example simulation configurations, the solvated BPTI test molecule with 14281 atoms includes the most robust set of configuration files. Each test molecule example include a README file which describes the origination of the molecular parameters and any preprocessing (equilibration, minimization, etc...) which has been performed. The following three subsections provide the reader with a more indepth look at representative simulation configuration files.

4.1 Alanin Configuration File with Leapfrog Integrator

Following is an example of a configuration file for the protein molecule alanin. In this case, ProtoMol will run 10000 steps starting at time 0, and random velocities will be generated based on the initial temperature of 300.0 and the random seed of 1234. A cubic cell manager is used (this must be true) with a cell size of 6.5. Output files, in this case just an all energies file and a DCD trajectory file will be written every 100 timesteps. Also a final PDB position and velocity file will be written. Restart files alanin.out.pos.pdb, alanin.out.vel.pdb, alanin.out.psf and alanin.out.conf will be written every 1000 timesteps. The initial parameter file alanin.par uses the newer CHARMM format. Vacuum boundary conditions are present, and the integrator is a single-timestep leapfrog integrator, with all bonded and two nonbonded forces present. The van der Waals and electrostatic forces both use a nonbonded cutoff algorithm and have cutoffs at 6.5. The switching function for the van der Waals potential function makes the second derivative continuous at the cutoff point, while the electrostatic potential function has a continuous first derivative due to its switching function. The switchon option has been used for van der Waals (recall that it can only be used with the C2 switching function). And the haptic force is operating at port 2001, with a transmission rate of 1 and a timeout value of 1000, a step increment value of 100 and a wait value of 100, which will allow time for an IMD connection, since the haptic force allows for interaction with the program VMD (see chapter 6).

temperature 300.0
firststep 0
numsteps 10000
cellsize 6.5
outputfreq 100
restartfreq 1000
seed 1234
posfile alanin.pdb
psffile alanin.psf
parfile alanin.par
usecharmm28parfile yes
finpdbposfile alanin.out.pos.pdb
finpdbvelfile alanin.out.vel.pdb
dcdfile alanin.out.dcd
allenergiesfile alanin.out.energy
restartfile alanin.out
boundaryConditions vacuum
cellManager Cubic
Integrator {
     level 0 Leapfrog {
               timestep 1
       force Improper
       force Dihedral
       force Bond
       force Angle
       force LennardJones
                 -algorithm NonbondedCutoff
                 -switchingFunction C2
         -switchon 1.0
                 -cutoff 12
       force Coulomb
                 -algorithm NonbondedCutoff
                 -switchingFunction C1
                 -cutoff 12
       force Haptic
         -port 2001
         -trate 1
         -timeout 1000
         -step_inc 100
         -wait 100
     }
}

Figure 4.1: Picture of alanin, generated by the BOOGA application raytrace .

4.2 20 Å diameter water droplet With A Two-Level Integration Scheme

This sample of water is being run using a two-level integration scheme, starting with an MTS Impulse Integrator with a cycle length of 5 and finishing with an STS PLeapFrog Integrator that has a timestep of 0.5. The initial temperature is 300 and the first timestep is 0, and the simulation will be run for 200 steps. Note the cubic cell manager and periodic boundary conditions, and that all cell basis vectors have been provided. Note also that though the user specified a restart file prefix and a restart frequency, restart files will never be written because restartfreq at 1000 is higher than the total number of steps for the simulation. Output files, however, will be written, including a DCD trajectory file and an all energies file, every 10 steps. The final files written will be both a position and velocity file in XYZ format, and also a final BSDL file. Note that all necessary initial data files have been provided.

temperature 300.0
firststep 0
numsteps 200
cellsize 6.5
outputfreq 10
restartfreq 1000
seed 1234
posfile equil298K_01.pos.pdb
velfile equil298K_01.vel.pdb
psffile equil298K_01.psf
parfile equil298K_01.par
finXYZPOSFILE          equil298K_01.out.fin.pos.xyz
finXYZVELFILE          equil298K_01.out.fin.vel.xyz
DCDFILE                 equil298K_01.out.trajectory.pos.dcd
ALLENERGIESFILE        equil298K_01.out.energies
BSDLFILE                equil298K_01.out.fin.bsdl
RESTARTFILE             equil298K_01.out.restart
doRestartFile no
usecharmm28parfile no
cellBasisVector1      28.0 0.0 0.0
cellBasisVector2      0.0 28.0 0.0
cellBasisVector3      0.0 0.0 28.0
cellorigin             0.0 0.0  0.0
boundaryConditions Periodic
cellManager Cubic
Integrator {
     level 1 Impulse {
                cyclelength 5
        force Coulomb
                  -algorithm NonbondedSimpleFull
                  -switchingFunction ComplementC1
                  -cutoff 6.5
      }
      level 0 PLeapfrog {
                timestep .5
        force Bond, Angle
        force Coulomb
                  -algorithm NonbondedCutoff
                  -switchingFunction C1
                  -cutoff 6.5
        force LennardJones
                  -algorithm NonbondedCutoff
                  -switchingFunction C2
                  -cutoff 6.5
                  -switchon 0.1
     }
}

Figure 4.2: Picture of 10 Å diameter water droplet, generated by the BOOGA application raytrace .

4.3 BPTI with Hybrid Monte Carlo Sampling

Here is an example of the Hybrid Monte Carlo algorithm used to sample a BPTI molecule - a 2-level Hybrid Monte Carlo MTS integrator is used. Once again, the initial timestep is at time 0, but the number of timesteps to run this time is 100. 1-2 exclusions are desired, and the cellsize is once again 6.5. There are now however periodic boundary conditions, and as a result three cell basis vectors along with a cell origin are required. Output files will be written to every 10 timesteps - this time they include XYZ position and trajectory files, a DCD trajectory file, an XYZ forces file, and a split energies file (thus there will be 10 different energy files generated - one for each type of energy). The restart files will be written three times, at 0, 50 and 100 timesteps. Note that a final XYZ position file will be written but NOT a final XYZ velocities file even though the filename was specified; setting dofinxyzvelfile to “no” overrides the filename. Also a final BSDL file - bpti.out.bsdl3 will be written. Since commotion is set to “yes”, center of mass motion will be removed when calculating velocities. Once again notice the cubic cell manager, and also notice that the HMCIntegrator does not have forces. The integrator will run 5 Hybrid Monte Carlo cycles, with a cycle length of 10 and 50 warm-up cycles, at a temperature of 300 K. Four bonded forces (bond, angle, improper, dihedral) are present, and one nonbonded (FullEwald, the Coulomb force solved for using Ewald summation).

firststep 0
numsteps 100
exclude 1-2
cellsize 6.5
cellbasisvector1  32.0     0     0
cellbasisvector2      0 32.0     0
cellbasisvector3      0     0 32.0
cellorigin             0     0     0
outputfreq 10
restartfreq 50
posfile bpti.pos.pdb
velfile bpti.vel.xyz
psffile bpti.psf
parfile bpti.par
finxyzposfile bpti.out.pos.xyz
finxyzvelfile bpti.out.vel.xyz
dofinxyzvelfile no
xyzposfile bpti.out.trajectory.pos.xyz
xyzvelfile bpti.out.trajectory.vel.xyz
dcdfile bpti.out.dcd
xyzforcesfile bpti.out.forces.xyz
doallenergiesfile no
splitenergiesfile bpti.out.energy
bsdlfile bpti.out
restartfile bpti.out
boundaryConditions Periodic
commotion yes
cellManager Cubic
Integrator {
       level 1 HybridMC {
                cyclelength 10
                warmupcycles 50
                temperature 300.0
   }
       level 0 PLeapFrog {
                timestep 1
       force Bond, Angle, Dihedral, Improper
       force Coulomb
          -algorithm FullEwald
          -real
          -reciprocal
          -correction
   }
}

Parallelization

Parallel communication for ProtoMol is based on Message Passing Interface (MPI). What ProtoMol uses is an incremental parallelization scheme for clusters with a moderate number of nodes, e.g. up to 64 CPUs. The design is based on a force decomposition with a master-slave paradigm.

The appropriate steps for running ProtoMol on parallel machines depends greatly on the machine being used. On SGI you need to add “mpirun -np <number of nodes> ” at the beginning of the command line, and on the

IBM Regatta you start ProtoMol by “poe ./protomol <protomol arugements> -procs <number of nodes>  ”.

Otherwise, the appropriate flags may be set in the compilation of ProtoMol to generate an MPI compilation.

On some machines a queuing system like LSF or PBS may be necessary. For more information on this ask your system administrator.

Interaction of ProtoMol with VMD

You may start VMD normally before you begin running the desired simulation to be viewed in ProtoMol . Select the “molecule” option from the main window, and in the resulting molecule window select “Load From Files”. Then in the files window that opens, under the menu “Molecule File Types” select psf and pdb, then write the full pathnames of the PSF and PDB files used in the simulation that will be run in the provided spaces. Now VMD has the initial state of the molecule to be simulated.

Next, in the configuration file that will be used, make sure there is a Haptic force in the integrator. Check to see that a positive nonzero value is supplied for the “-wait” value under the haptic force. This will force ProtoMol to wait some time for an IMD connection, 100 would be a decent value so that you have time to set VMD to simulate. Also note the port number, which will be used later.

Now run ProtoMol . It should display “Waiting for IMD connection.....” after a short time. At this point, go to

the VMD console window and type “imd connect <hostname> <port number>  ”, where <hostname> is

the name of the machine that ProtoMol is running on and <port number> is the port number value specified in the configuration file that is being used.

At this point, VMD will view the molecule throughout the simulation as ProtoMol runs. Also from the VMD window IMD commands can be run that will send messages to ProtoMol . For example, typing “imd trate <integer value> ” will change the transmission rate for information from ProtoMol to VMD.

As an example, let us return to the integrator section of the alanin configuration file from section 4.1:

Integrator {
        level 0 Leapfrog {
                  timestep 1
          force Improper
          force Dihedral
          force Bond
          force Angle
          force LennardJones
                     -algorithm NonbondedCutoff
                     -switchingFunction C2
             -switchon 0.1
                     -cutoff 12
          force Coulomb
                     -algorithm NonbondedCutoff

        -switchingFunction C1
        -cutoff 12
force Haptic
-port 2001
-trate 1
-timeout 1000
-step_inc 100
-wait 100
}
}

The only area that counts as far as the VMD simulation is concerned is the “force Haptic ” section at the end. Notice first that a -wait value is given that is greater than zero, thus ProtoMol will wait for an IMD connection, which is what we want. The port ID is 2001. The other three values are IMD variables.

In this case as soon as ProtoMol begins waiting for the IMD connection, a user would type

“imd connect <name of the machine on which ProtoMol is running> 2001  ” and the visual

simulation will begin. Now the user can run IMD commands in the IMD window while the simulation is running.

Figure 6.1: PHANTOM Desktop^{T M} Haptic Device.

Availability and Installation of ProtoMol

7.1 How to Download ProtoMol

The download site for ProtoMol is a link from the page for the Laboratory for Computational Life Sciences at the University of Notre Dame: http://protomol.sourceforge.net/. Executable code is available online, and source code will be available upon request. An agreement to a Non-Exclusive, Non-Commercial Use License is required (this exact license can be found in the appendix of this user guide, appendix B) for download.

7.2 Platforms that ProtoMol is Able to Run On

ProtoMol executables are available for the following three platforms:

1. Sun OS 5.8

2. Linux 2.4

3. IRIX 6.5 (optional with MPI)

4. AIX 5.1 (optional with MPI)

5. HP-UX (optional with MPI)

6. Windows

7.3 Compiling ProtoMol

Because executables are available online, compiling ProtoMol should not be necessary unless the user would like to change features. Here are the steps for generating the appropriate makefiles on different platforms:

1. Change to the main protomol directory.

2. View the README file for specific directions on how to compile ProtoMol correctly on the apporpriate platform.

Once the makefiles have been generated, the user types “make ” from the protomol directory. A “make clean ” will clear all object files and executables recursively throughout the directory tree, if necessary.

After everything has compiled successfully, the user may run ProtoMol from the

protomol/applications/protomol-app  directory as described in Chapter 2.

7.4 Documentation

Available documentation for ProtoMol can be downloaded from:

http://protomol.sourceforge.net/

The Future of ProtoMol

Work is already under way to further exapand ProtoMol . In the next few years, various new algorithms will be placed in the back end, including new MOLLY integrators and semi-implicit integrators using iterative solvers and sparse matrices, rigid body dynamics, and table lookups. The goal is to acheive larger time steps and thus decrease the amount of time necessary to run a simulation.

The ProtoMol design team is also looking into automated empirical optimization of the software (AEOS) for optimal parameters and algorithm selection. With this implementation, the software could automatically generate an extremely efficient simulation, enhancing both usability and speed.

Additionally, the next release of ProtoMol will hopefully include a graphical interface, again increasing ease of use. The program is currently being ported to Windows and will be released in the near future.

Appendix A

Background

A.1 Integrators

MD describes a molecular system as a function of time based on the integration of Newton’s equations of motion and interacting forces. The integrators are the part of the program that that solve the differential equations that describe the system. Specifically, the integrators provide a set of forces, velocities, positions, etcetera, that describe the system at each time step. Thus, it is easy to see that the central part of the entire system is the integrators: the front end reads in data and passes it to the integrators, then the integrators call the back end’s high performance functions to do the massive computations. The way those computations are assembled into useful data describes the integrator.

A.1.1 Standard integrators in ProtoMol

For ProtoMol , a standard integrator is one that updates the velocities by half a step, updates the positions, calculates forces, and finally updates the velocities by another half a step. This method is noted in ProtoMol as:

• halfkick();

• doDriftOrNextIntegrator();

• calculateForces();

• halfkick();

A.1.2 Distinction between multiple-time-stepping (MTS) and single-time-stepping (STS) integrators

Because of the massive amounts of computation for very large systems (as are found in MD), it is essential to have more sophisticated ways of integration than to just use a single time unit per step. With an STS integrator, one typically takes timesteps of one femtosecond. This small of a timestep would however not be feasable for long simulations with many atoms since useful simulations are typically run for somewhere in the range of one microsecond to one second. This can mean literally months of computations on the computer. Thus, MTS integrators are desired since they are capable of taking longer timesteps.

When running ProtoMol , the integrator at level 0 (or terminating integrator) must be STS. Looking at the sequence of standard integrator operations in section A.1.1, when an MTS integrator calls doDriftOrNextIntegrator(), it looks for an MTS or STS integrator that is next in the chain to run. But when an STS integrator calls this function, it actually updates the positions, calculates the forces, updates the velocities and returns. Therefore because an STS integrator does not look for a following integrator, it must not have one. Thus an STS integrator in ProtoMol must be at level 0, while an MTS integrator can never be at that level.

A.1.3 The Verlet-I/r-RESPA Impulse Method [15]

The Impulse method is a MTS method for solving Newton’s equations of motion and is supported by ProtoMol . The first step is to separate the system force into a fast and slow part, each with their own distinct timesteps. Bonded terms have the fastest normal modes in the system and are thus classified as fast. Nonbonded terms, both Lennard Jones and electrostatic, are split into fast and slow depending on range (short-range is fast, long-range is slow). With appropriate bonded conditions, the fast part becomes cheap to integrate compared to the cost of integrating the full system, since only short-range forces are considered. Long-range forces (electrostatic forces beyond 8 Å ) need to be evaulated sparingly for efficiency purposes. The impulse method “approximates the slow force by a sequence of impulses with weights Δt dictated by consistency.” Analytical evaluation of the fast forces can be avoided using a similar method - a sequence of suitably weighted impulses that are more closely spaced in time than the slow force impulses. The method as a whole can be described by the following equation:

M \frac{d^{2}}{dt^{2}} X = ∑_{n′=−∞}^{∞} δtδ(t − n′Δt)F^{fast}(X) + ∑_{n=−∞}^{∞} Δtδ(t − nΔt)F^{slow}(X)  (A.1)

A.1.4 Mollified Impulse Methods (MOLLY ) [15]

MOLLY is a family of integrators [11] that counteracts the instabilities and accuracy reduction present in the multiple time stepping Verlet-I/r-RESPA integrator. This is accomplished by perturbing the potential using time averaged positions

U^{slow}(X) → U^{slow}(A(X)),  (A.2)

with the force defined as a gradient of this averaged potential,

−∇U^{slow}(X) → −A_{X}(X)^{T}∇U^{slow}(A(X)).  (A.3)

This perturbation is supposed to compensate for finite Δt artifacts. Perturbing the potential rather than the force ensures that the numerical integrator remains symplectic [23]. The force used by MOLLY is the gradient of the perturbed potential [19]. MOLLY can be seen as a filter that eliminates components of the slow force impulse in the directions of the fast forces, and thus improves the stability of Verlet-I/r-RESPA. Different averaging functions give rise to MOLLY integrators with different stability and accuracy properties.

• Equilibrium MOLLY - [18] Here we present a time averaging completely eliminate the components of −U_{x}^{slow} in the directions of the fast forces. Our starting point is the backward Euler averaging:

A(x) = x − Δt^{2}M^{−1}U_{x}^{fastest}(A(x)),  (A.4)

which is obtained with one step of backward Euler starting with zero velocities.

Assume that U^{fastest} can be written as U^{fastest} = χ(g(x)), where g(x) is a vector of independent length constraint functions. The elements of g(x) are of the form

g^{k}(x) = x^{T}G_{k}x − l_{k}^{2},  (A.5)

where the G_{k} are symmetric matrices and l_{k} are rest lengths. Both U^{bond} and U^{angle} can be expressed in the form χ(g(x)). Allowing Δt → ∞ gives the equations for A(x) and μ

\{ \begin{matrix}M(A(x) − x) + g_{x}^{T}(A(x))μ = 0, \\ g(A(x)) = 0, \end{matrix} \} .  (A.6)

which are termed Equilibrium*; this method completely eliminates components of −∇U^{slow} in the directions of the constraints g_{x}^{T}. A proof of the stability of Equilibrium* for a general linear model problem under the assumption that all the fast forces are included in the mollification (U^{fastest} = U^{fast}) is given in Ref. [16]. An earlier version is in Ref. [22]. The stability condition for Equilibrium* MOLLY on the longest time step

Δtρ(M^{−1/2}U_{xx}^{slow}M^{−1/2})^{1/2} < 2,  (A.7)

where ρ is the spectral radius of the mass weighted Hessian matrix M^{−1/2}U_{xx}^{slow}M^{−1/2}, is less restrictive than that of the popular leapfrog integrator:

Δtρ(M^{−1/2}U_{xx}M^{−1/2})^{1/2} < 2,  (A.8)

because it is limited by the fastest frequencies in the slow forces (−∇U^{slow}) rather than by the overall fastest frequencies (in −∇U^{fast}).

• BSpline MOLLY - [18] It is possible to use time averagings that consist of numerically integrating an auxiliary, reduced problem:

A(x) = \frac{1}{Δt} ∫_{0}^{∞} ϕ (\frac{t}{Δt}) X̃(t)dt  (A.9)

where ϕ (\frac{t}{Δt}) is a weight function, and X̃(t) solves an auxiliary problem

M \frac{d^{2}}{dt^{2}} X̃ = F^{reduced}(X̃), X̃(0) = x, \frac{d}{dt}X̃(0) = 0.  (A.10)

This approach is computationally feasible if the weight functions ϕ have compact support in time. The paper [11] suggests using B-spline weight functions, which are non-zero over a short interval. The effectiveness of the averagings induced by these weight functions is directly related to the extensiveness of the time averaging. One such B-spline weight function that has been tested is called LongAverage:

ϕ(s) = \begin{cases}0, & s < 0, \\ 1, & 0 ≤ s < 1, \\ \frac{1}{2}, & s = 1 \\ 0, & s > 1. \end{cases}  (A.11)

The coding of A(x) and A_{x}(x) can be done by hand in a systematic manner. First the calculation of A(x) is coded, and then the differentiation, applying the chain rule with respect to each of the components of x to yield code for A_{x}(x). As an example suppose that the leapfrog method with time step δt is coded for the calculation of A(x). This is then differentiated to obtain A_{x}(x). The result is the following code for calculating A(x) and A_{x}(x): Initialization is given by

\begin{matrix}X := x, & X_{x} := I, \\ P := 0, & P_{x} := 0, \\ B := 0, & B_{x} := 0, \\ t := 0, \end{matrix}  (A.12)

and step by step integration by

P := P + \frac{1}{2}δtF^{reduced}(X), P_{x} := P_{x} + \frac{1}{2}δtF_{x}^{reduced}(X)X_{x},
B := B + \frac{1}{2}δtXϕ(t/Δt), B_{x} := B_{x} + \frac{1}{2}δtX_{x}ϕ(t/Δt),
X := X + δtM^{−1}P,  X_{x} := X_{x} + δtM^{−1}P_{x},
t := t + δt,
B := B + \frac{1}{2}δtXϕ(t/Δt), B_{x} := B_{x} + \frac{1}{2}δtX_{x}ϕ(t/Δt),
P := P + \frac{1}{2}δtF^{reduced}(X), P_{x} := P_{x} + \frac{1}{2}δtF_{x}^{reduced}(X)X_{x}.

(A.13)

The value (1/Δt)B is used for A(x) and (1/Δt)B_{x} for A_{x}(x). We continue the above integration until we reach a value of t such that ϕ(t/Δt) is zero at this value and remains zero for larger values of t. In practice, one can choose δt equal to the stepsize of the lowest level integrator. In the above loop, F_{x} = −U_{xx}^{reduced}(x) must be computed efficiently. We have derived the analytical form of the Hessian matrices, U_{xx}^{reduced}(x) [17].

• H-Bond MOLLY - [14] A large part of using MOLLY is to determine which forces need to be included in the time averaging. Systems that are particularly sensitive to instability are those dissolved in water. Water is a hydrogen bonded system, and the presence of hydrogen bonds accounts for many important properties of liquid water, proteins, and their interactions. It has been found that hydrogen bond interactions are very important terms to include in the time averaging. H-Bond MOLLY uses not only the bonds and angles in the averaging, but also the H-Bonds which are modeled by the nonbonded short range Coulomb and LennardJones forces, and thus presents itself as a superior integrator, even to Equilibrium. A straightforward implementation of B-spline MOLLY and H-Bond MOLLY would compute the filter matrix, A(x) directly, which involves sparse matrices and multiplies of those sparse matrices. Another way to implement these two methods is to construct the filter matrix implicitly and completely get rid of sparse matrix multiplies using a linear approximation of A(x). A more detailed description of the efficient implementation of H-Bond MOLLY is available on line

(http://www.cse.nd.edu/~lcls/docs/fastHbond/  ).

A.2 Forces

Forces define the interaction between the atoms. Newton’s equations of motion, which are solved by the integrators for the interacting atoms, evaluate these forces to calculate new positions and velocities. On a given atom at position r⃗_{i},

m_{i} r̈_{i} = F⃗_{i} = ∑_{j=1,j≠i}^{N} −∇U(r⃗_{ij}).  (A.14)

Here m_{i} is the mass of particle i and r_{ij} is the distance between particle i and j. The complexity of the force calculation is governed by the potential function U(r⃗_{ij}).

In ProtoMol , forces are held by the integrators and can be removed or added dynamically. ProtoMol makes distinctions between system forces and extended forces, as well as non-bonded and bonded forces.

A.2.1 Bond force

Bonds [21] describe a linear bond between two atoms. These bonds are described by a simple harmonic spring. The energy of a bond between atoms i and j is given by:

E_{bond} = k (\left|{ \atop r_{ij}}\right| − r_{0})^{2}  (A.15)

where

k  = Spring constant specified in the parameter file for this bond type.
r⃗_{ij} = r⃗_{j} - r⃗_{i}
\left|{ \atop r_{ij}}\right| = \sqrt{(x_{j} − x_{i})^{2} + (y_{j} − y_{i})^{2} + (z_{j} − z_{i})^{2}}
= calculated distance between atoms i and j.
r_{0}  = Rest distance of the bond specified in the parameter file for this bond type.

By differentiating this formula, the force for this bond can be found to be:

F⃗_{bond} = −\frac{dE_{bond}}{dr⃗} = 2k(\left|{ \atop r_{ij}}\right| − r_{0})r̂_{ij}  (A.16)

A.2.2 Angle force

Angles [21] describe angular bonds between three atoms. These bonds are modeled as harmonic angular springs. The energy of such a bond between atoms i, j, and k is given by:

E_{angle} = E_{θ} + E_{ub}  (A.17)
E_{θ} = k_{θ} (θ − θ_{0})^{2}  (A.18)
E_{ub} = k_{ub} (\left|{ \atop r_{ik}}\right| − r_{ub})^{2}  (A.19)

where

k_{θ}  = Force constant specified in the parameter file for this bond type.
θ  = cos^{−1} (\frac{r⃗_{ij} · r⃗_{kj}}{\left|{ \atop r_{ij}}\right|\left|{ \atop r_{kj}}\right|}) , which is equivalent to tan^{−1}(|r⃗_{ij} × r⃗_{kj}|, r⃗_{ij} · r⃗_{kj})
θ_{0}  = Rest angle of this bond specified in the parameter file for this angle type.
k_{ub} = Urey-Bradley constant, which defaults to zero
r⃗_{ik} = r⃗_{k} - r⃗_{i}
\left|{ \atop r_{ik}}\right| = \sqrt{(x_{k} − x_{i})^{2} + (y_{k} − y_{i})^{2} + (z_{k} − z_{i})^{2}}, calculated distance between atoms i and k.
r_{ub} = Rest distance for the Urey-Bradley term.

By differentiating this formula, the force for this bond can be found to be:

F⃗_{angle} = −\frac{dE_{angle}}{dr⃗} = F⃗_{θ} + F⃗_{ub}  (A.20)

where

F⃗_{θ} = −\frac{dE_{θ}}{dθ} \frac{dθ}{dr⃗} = −2k_{θ}(θ − θ_{0}) \frac{dθ}{dr⃗}

F⃗_{ub} = −\frac{dE_{ub}}{dr⃗} = 2k_{ub}(\left|{ \atop r_{ik}}\right| − r_{ub})r̂_{ik}

where r̂_{ik} is the unit vector in the direction of r⃗_{ik}.

A.2.3 Dihedral and improper forces

Dihedral and improper bonds [21] model the interaction between 4 bonded atoms. They are modeled by an angular spring between planes formed by the first three atoms and the second 3 atoms. The energy for a dihedral or improper between atoms i, j, k, and l is given by:

E_{d/i} = \{ \begin{matrix}k (1 + cos(nϕ + δ)) & if n > 0 \\ k (ϕ − δ)^{2} & if n = 0 \end{matrix}  (A.21)

Where

k = Force constant specified in the parameter file for this bond type.
ϕ = Calculated angle between the plane formed by atoms i, j, and k

and the plane formed by atoms j, k, and l.

n = Periodicity of the bond specified in the parameter file for this bond type.
δ = Phase shift specified in the parameter file for this bond type.

The angle ϕ is calculated by first determining the vectors A⃗ and B⃗ where:

A⃗ = r⃗_{ij} × r⃗_{jk}
B⃗ = r⃗_{jk} × r⃗_{kl}

These two vectors can then be used to calculate ϕ using the formula,

ϕ = cos^{−1} (\frac{A⃗ · B⃗}{|A⃗||B⃗|})  (A.22)

To determine the force, the negative gradient of the energy must be found. It can be shown that if n=0 then the force is given by,

F⃗ = −2k (ϕ − δ) (∇ϕ)  (A.23)

and if n > 0 then the force is given by,

F⃗ = n k sin(nϕ + δ) (∇ϕ)  (A.24)

Using the formula for ϕ given, it can be shown that:

∇ϕ = \frac{1}{sin(ϕ)}∇ (\frac{A⃗ · B⃗}{|A⃗||B⃗|})  (A.25)

But this can lead to a singularity in sin(ϕ) goes to 0. When such a situation has the potential to arise, a third vector C⃗ is determined using

C⃗ = r⃗_{jk} × A⃗  (A.26)

If the angle ψ is the angle between C⃗ and B⃗ then it can be shown that

cos(ϕ) = sin(ψ)  (A.27)

and,

ϕ = sin^{−1} (\frac{C⃗ · B⃗}{|C⃗||B⃗|}) + \frac{π}{2}  (A.28)

Therefore

∇ϕ = \frac{1}{cos(ϕ)}∇ (\frac{C⃗ · B⃗}{|C⃗||B⃗|})  (A.29)

A.2.4 Non-bonded forces

Bonded forces define the intramolecular interactions, where non-bonded forces compromise the intermolecular interactions. Non-bonded interactions also play an important role in the behavior of individual molecules. The non-bonded interactions do not depend upon a specific molecular relationship between atoms. They are usually modeled as a function of some inverse power of the distance. Non-bonded terms are usually divided in two groups: one compromising electrostatic and the other van der Waals interactions.

ProtoMol provides a wide range of different non-bonded force evaluations, which can be combined with an arbitrary boundary condition, switching function and potential. Decoupling the algorithm to compute a given force from the boundary condition, switching function and potential, we achieve flexibility and performance. We template the forces for boundary conditions, switching functions and potentials. van der Waals forces [2] van der Waals forces are “relatively weak electric forces that attract neutral molecules to one another in gases, in liquefied and solidified gases, and in almost all organic liquids and solids. The forces are named for the Dutch physicist Johannes van der Waals, who in 1873 first postulated these intermolecular forces in developing a theory to account for the properties of real gases. Solids that are held together by van der Waals forces characteristically have lower melting points and are softer than those held together by the stronger ionic, covalent, and metallic bonds.

Van der Waals forces may arise from three sources. First, the molecules of some materials, although electrically neutral, may be permanent electric dipoles. Because of fixed distortion in the distribution of electric charge in the very structure of some molecules, one side of a molecule is always somewhat positive and the opposite side somewhat negative. The tendency of such permanent dipoles to align with each other results in a net attractive force. Second, the presence of molecules that are permanent dipoles temporarily distorts the electron charge in other nearby polar or nonpolar molecules, thereby inducing further polarization. An additional attractive force results from the interaction of a permanent dipole with a neighboring induced dipole. Third, even though no molecules of a material are permanent dipoles (e.g., in the noble gas argon or the organic liquid benzene), a force of attraction exists between the molecules, accounting for condensing to the liquid state at sufficiently low temperatures.

The nature of the attractive force in molecules, which requires quantum mechanics for its correct description, was first recognized (1930) by the Polish-born physicist Fritz London, who traced it to electron motion within molecules. London pointed out that at any instant the center of negative charge of the electrons and the center of positive charge of the atomic nuclei would not be likely to coincide. Thus, the fluctuation of electrons makes molecules time-varying dipoles, even though the average of this instantaneous polarization over a brief time interval may be zero. Such timevarying dipoles, or instantaneous dipoles, cannot orient themselves into alignment to account for the actual force of attraction, but they do induce properly aligned polarization in adjacent molecules, resulting in attractive forces. These specific interactions, or forces, arising from electron fluctuations in molecules (known as London forces, or dispersion forces) are present even between permanently polar molecules and produce, generally, the largest of the three contributions to intermolecular forces.”

A common and well known approximation for the van der Waal forces is the general Lennard-Jones 6-12 (LJ) potential

u(r_{ij}) = 4ϵ [(\frac{α}{r_{ij}} )^{12} − (\frac{β}{r_{ij}} )^{6}] .  (A.30)

where ϵ, α and β are LJ parameters with α ≈ β, often is α = β, it is called σ. The potential u(r_{ij}) is mostly cutoff at r_{cutoff} without major loss of accuracy, i.e. no interactions are evaluated beyond this distance. Since the cutoff violates assumptions for a simplectic integrator, switching functions are introduced to smooth the potential at the cutoff.

Coulomb force [2] “also called coulomb interaction, attraction or repulsion of particles or objects because of their electric charge. One of the basic physical forces, the electric force is named for a French physicist, Charles-Augustin de Coulomb, who in 1785 published the results of an experimental investigation into the correct quantitative description of this force.

Two like electric charges, both positive or both negative, repel each other along a straight line between their centers. Two unlike charges, one positive, one negative, attract each other along a straight line joining their centers. The electric force is operative between charges down to distances of at least 10^{−16} meter, or approximately one-tenth of the diameter of atomic nuclei. Because of their positive charge, protons within nuclei repel each other, but nuclei hold together because of another basic physical force, the strong interaction, or nuclear force, which is stronger than the electric force. Massive, but electrically neutral, astronomical bodies such as planets and stars are bound together in solar systems and galaxies by still another basic physical force, gravitation, which though much weaker than the electric force, is always attractive and is the dominant force at great distances. At distances between these extremes, including the distances of everyday life, the only significant physical force is the electric force in its many varieties along with the related magnetic force.

The magnitude of the electric force F is directly proportional to the amount of one electric charge, q_{i}, multiplied by the other, q_{j}, and inversely proportional to the square of the distance r between their centers. Expressed in the form of an equation, this relation, called Coulomb’s law, may be written by including the proportionality factor k as F = kq_{i}q_{j}/|r_{i} − r_{j}|. In the centimeter-gram-second system of units, the proportionality factor k in a vacuum is set equal to 1 and unit electric charge is defined by Coulomb’s law.”

A.3 Switching Functions

ProtoMol allows the cutting off of non-bonded forces (electrostatic or van der Waals) to be handled in two ways. One way is to simply truncate the forces at the cutoff value, i.e. no interactions are evaluated beyond this distance. But this method leads to a discontinuity in the force field and it will violate energy conservation for symplectic integrators. The other way of handling cutoffs are switching functions. These functions smoothly bring the forces and energies to 0 at the cutoff distance to avoid any discontinuities in the force field.

A.4 Boundary Conditions

Boundary conditions define the size of the simulation area and how atoms will be treated when reaching or crossing boundaries. The most common area is a cubic box of a finite region.

Vacuum boundary conditions have no restrictions on the dimension of the enclosing lattice. For some applications, the ions and the molecules/atoms can move in every dimension freely. Simulations requiring no periodic boundaries are best suited to simulation in vacuum, such as the conformational study of an isolated polymer molecule. Vacuum boundary conditions are not recommended for studies in a solvent, since evaporation is likely to be a problem. Note that vacuum boundary conditions cannot be used with the Ewald summation method.

Periodic boundary conditions (PBC) replicate the simulation box through space to form an infinite lattice. For some applications, the PBC may be defined only for some dimensions. PBC use the minimal image convention, whenever an atom moves out of the minimal image it moves in on the opposite of the minimal image. PBC in all dimesions can be used with the Ewald summation method.

Spherical boundary conditions (SBC) causes a force that moves atoms back in to (away from) the center of the sphere depending on the force constant. It may be desirable, however, to simulate a surface tension affect, where a small harmonic well around the sphere boundary is desired. To obtain such an effect, two SBC can be superimposed.

A.5 Cell Manager

The cell manager is implemented to achieve better performance. They control which atom belongs to which cell and store spatial information of the atom positions, such that the particle’s neighbors can be quickly located for the force calculation. Thus, the number of interactions evaluated is minimized.

The multi-cell algorithm uses cells together with a cutoff, where the cells are cubic and the dimension of the cells are equal or more than the cutoff.

Figure A.1: Partition of cells, where cell no. 5 is the origin and 1-4, 6-9 the neighbor cells. Cell no. 6, 3, 2 and 1 describe the interaction path.

Appendix B

ProtoMol Molecular Dynamcis Software License

B.1 Conditions and Regulations

ProtoMol GNU GENERAL PUBLIC LICENSE
This Software is copyrighted by the University of Notre Dame and the University
of Bergen, Norway, whom grant you the right to use this software under
the terms of the GNU General Public License. Although you are under no
obligation, the ProtoMol devlopment team would appreciate the
inclusion of the following citation in all research and development
works which utilize this software:
T. Matthey, T. Cickovski, S. Hampton, A. Ko, Q. Ma, T. Slabach
and J. A. Izaguirre, 'ProtoMol, an Object-Oriented Framework for Prototyping
Novel Algorithms for Molecular Dynamics', ACM Transactions on Mathematical
Software, Vol 30, Num 3, pg 237-265, 2004
                 GNU GENERAL PUBLIC LICENSE
   Version 2, June 1991
Copyright c
               1989, 1991 Free Software Foundation, Inc.
59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Everyone is permitted to copy and distribute verbatim copies of this license
document, but changing it is not allowed.
Preamble
The licenses for most software are designed to take away your freedom to
share and change it. By contrast, the GNU General Public License is in-
tended to guarantee your freedom to share and change free software-to
make sure the software is free for all its users. This General Public Li-
cense applies to most of the Free Software Foundation's software and to any
other program whose authors commit to using it. (Some other Free Software
Foundation software is covered by the GNU Library General Public License
instead.) You can apply it to your programs, too.
   When we speak of free software, we are referring to freedom, not price.
Our General Public Licenses are designed to make sure that you have the

freedom to distribute copies of free software (and charge for this service if
you wish), that you receive source code or can get it if you want it, that you
can change the software or use pieces of it in new free programs; and that
you know you can do these things.
   To protect your rights, we need to make restrictions that forbid anyone to
deny you these rights or to ask you to surrender the rights. These restrictions
translate to certain responsibilities for you if you distribute copies of the
software, or if you modify it.
   For example, if you distribute copies of such a program, whether gratis
or for a fee, you must give the recipients all the rights that you have. You
must make sure that they, too, receive or can get the source code. And you
must show them these terms so they know their rights.
   We protect your rights with two steps: (1) copyright the software, and
(2) offer you this license which gives you legal permission to copy, distribute
and/or modify the software.

   Also, for each author's protection and ours, we want to make certain that
everyone understands that there is no warranty for this free software. If the
software is modified by someone else and passed on, we want its recipients
to know that what they have is not the original, so that any problems
introduced by others will not reflect on the original authors' reputations.
   Finally, any free program is threatened constantly by software patents.
We wish to avoid the danger that redistributors of a free program will in-
dividually obtain patent licenses, in effect making the program proprietary.
To prevent this, we have made it clear that any patent must be licensed for
everyone's free use or not licensed at all.
   The precise terms and conditions for copying, distribution and modifi-
cation follow.
Terms and conditions for copying, distribution and
modification
  0. This License applies to any program or other work which contains a
       notice placed by the copyright holder saying it may be distributed
       under the terms of this General Public License. The "Program", be-
       low, refers to any such program or work, and a "work based on the
       Program" means either the Program or any derivative work under
       copyright law: that is to say, a work containing the Program or a
       portion of it, either verbatim or with modifications and/or translated
       into another language. (Hereinafter, translation is included without
       limitation in the term "modification".) Each licensee is addressed as
       "you".
       Activities other than copying, distribution and modification are not
       covered by this License; they are outside its scope. The act of running
       the Program is not restricted, and the output from the Program is
       covered only if its contents constitute a work based on the Program
       (independent of having been made by running the Program). Whether
       that is true depends on what the Program does.
  1. You may copy and distribute verbatim copies of the Program's source
       code as you receive it, in any medium, provided that you conspicuously
       and appropriately publish on each copy an appropriate copyright no-
       tice and disclaimer of warranty; keep intact all the notices that refer
       to this License and to the absence of any warranty; and give any other

recipients of the Program a copy of this License along with the Pro-
gram.
You may charge a fee for the physical act of transferring a copy, and
you may at your option offer warranty protection in exchange for a
fee.
2. You may modify your copy or copies of the Program or any portion
of it, thus forming a work based on the Program, and copy and dis-
tribute such modifications or work under the terms of Section 1 above,
provided that you also meet all of these conditions:
 (a) You must cause the modified files to carry prominent notices stat-
         ing that you changed the files and the date of any change.
 (b) You must cause any work that you distribute or publish, that in
         whole or in part contains or is derived from the Program or any
         part thereof, to be licensed as a whole at no charge to all third
         parties under the terms of this License.
 (c) If the modified program normally reads commands interactively
         when run, you must cause it, when started running for such in-
         teractive use in the most ordinary way, to print or display an
         announcement including an appropriate copyright notice and a
         notice that there is no warranty (or else, saying that you provide
         a warranty) and that users may redistribute the program under
         these conditions, and telling the user how to view a copy of this
         License. (Exception: if the Program itself is interactive but does
         not normally print such an announcement, your work based on
         the Program is not required to print an announcement.)
These requirements apply to the modified work as a whole. If identifi-
able sections of that work are not derived from the Program, and can
be reasonably considered independent and separate works in them-
selves, then this License, and its terms, do not apply to those sections
when you distribute them as separate works. But when you distribute
the same sections as part of a whole which is a work based on the
Program, the distribution of the whole must be on the terms of this
License, whose permissions for other licensees extend to the entire
whole, and thus to each and every part regardless of who wrote it.
Thus, it is not the intent of this section to claim rights or contest your
rights to work written entirely by you; rather, the intent is to exercise
the right to control the distribution of derivative or collective works
based on the Program.
In addition, mere aggregation of another work not based on the Pro-
gram with the Program (or with a work based on the Program) on a
volume of a storage or distribution medium does not bring the other
work under the scope of this License.
3. You may copy and distribute the Program (or a work based on it,
under Section 2) in object code or executable form under the terms of
Sections 1 and 2 above provided that you also do one of the following:
 (a) Accompany it with the complete corresponding machine-readable
      source code, which must be distributed under the terms of Sec-
      tions 1 and 2 above on a medium customarily used for software
      interchange; or,
 (b) Accompany it with a written offer, valid for at least three years,

      to give any third party, for a charge no more than your cost of
      physically performing source distribution, a complete machine-
      readable copy of the corresponding source code, to be distributed
      under the terms of Sections 1 and 2 above on a medium custom-
      arily used for software interchange; or,
 (c) Accompany it with the information you received as to the offer to
      distribute corresponding source code. (This alternative is allowed
      only for noncommercial distribution and only if you received the
      program in object code or executable form with such an offer, in
      accord with Subsection b above.)
The source code for a work means the preferred form of the work for
making modifications to it. For an executable work, complete source
code means all the source code for all modules it contains, plus any
associated interface definition files, plus the scripts used to control
compilation and installation of the executable. However, as a special
exception, the source code distributed need not include anything that
is normally distributed (in either source or binary form) with the major
components (compiler, kernel, and so on) of the operating system on
which the executable runs, unless that component itself accompanies
the executable.
If distribution of executable or object code is made by offering access
to copy from a designated place, then offering equivalent access to
copy the source code from the same place counts as distribution of the
source code, even though third parties are not compelled to copy the
source along with the object code.
4. You may not copy, modify, sublicense, or distribute the Program ex-
cept as expressly provided under this License. Any attempt otherwise
to copy, modify, sublicense or distribute the Program is void, and will
automatically terminate your rights under this License. However, par-
ties who have received copies, or rights, from you under this License
will not have their licenses terminated so long as such parties remain
in full compliance.
5. You are not required to accept this License, since you have not signed
it. However, nothing else grants you permission to modify or distribute
the Program or its derivative works. These actions are prohibited by
law if you do not accept this License. Therefore, by modifying or
distributing the Program (or any work based on the Program), you
indicate your acceptance of this License to do so, and all its terms
and conditions for copying, distributing or modifying the Program or
works based on it.
6. Each time you redistribute the Program (or any work based on the
Program), the recipient automatically receives a license from the orig-
inal licensor to copy, distribute or modify the Program subject to these
terms and conditions. You may not impose any further restrictions on
the recipients' exercise of the rights granted herein. You are not re-
sponsible for enforcing compliance by third parties to this License.
7. If, as a consequence of a court judgment or allegation of patent in-
fringement or for any other reason (not limited to patent issues), con-
ditions are imposed on you (whether by court order, agreement or
otherwise) that contradict the conditions of this License, they do not
excuse you from the conditions of this License. If you cannot distribute

so as to satisfy simultaneously your obligations under this License and
any other pertinent obligations, then as a consequence you may not
distribute the Program at all. For example, if a patent license would
not permit royalty-free redistribution of the Program by all those who
receive copies directly or indirectly through you, then the only way
you could satisfy both it and this License would be to refrain entirely
from distribution of the Program.
If any portion of this section is held invalid or unenforceable under any
particular circumstance, the balance of the section is intended to apply
and the section as a whole is intended to apply in other circumstances.
It is not the purpose of this section to induce you to infringe any
patents or other property right claims or to contest validity of any such
claims; this section has the sole purpose of protecting the integrity of
the free software distribution system, which is implemented by public
license practices. Many people have made generous contributions to
the wide range of software distributed through that system in reliance
on consistent application of that system; it is up to the author/donor
to decide if he or she is willing to distribute software through any other
system and a licensee cannot impose that choice.
This section is intended to make thoroughly clear what is believed to
be a consequence of the rest of this License.
8. If the distribution and/or use of the Program is restricted in certain
countries either by patents or by copyrighted interfaces, the original
copyright holder who places the Program under this License may add
an explicit geographical distribution limitation excluding those coun-
tries, so that distribution is permitted only in or among countries not
thus excluded. In such case, this License incorporates the limitation
as if written in the body of this License.
9. The Free Software Foundation may publish revised and/or new ver-
sions of the General Public License from time to time. Such new
versions will be similar in spirit to the present version, but may differ
in detail to address new problems or concerns.
Each version is given a distinguishing version number. If the Program
specifies a version number of this License which applies to it and "any
later version", you have the option of following the terms and condi-
tions either of that version or of any later version published by the
Free Software Foundation. If the Program does not specify a version
number of this License, you may choose any version ever published by
the Free Software Foundation.
10. If you wish to incorporate parts of the Program into other free pro-
grams whose distribution conditions are different, write to the author
to ask for permission. For software which is copyrighted by the Free
Software Foundation, write to the Free Software Foundation; we some-
times make exceptions for this. Our decision will be guided by the two
goals of preserving the free status of all derivatives of our free software
and of promoting the sharing and reuse of software generally.
                                 NO WARRANTY
11. Because the Program is licensed free of charge, there is no
  warranty for the Program, to the extent permitted by ap-
  plicable law. except when otherwise stated in writing the
  copyright holders and/or other parties provide the program

"as is" without warranty of any kind, either expressed or im-
plied, including, but not limited to, the implied warranties
of merchantability and fitness for a particular purpose. The
entire risk as to the quality and performance of the Program
is with you. Should the Program prove defective, you assume
the cost of all necessary servicing, repair or correction.
12. In no event unless required by applicable law or agreed to
in writing will any copyright holder, or any other party who
may modify and/or redistribute the program as permitted
above, be liable to you for damages, including any general,
special, incidental or consequential damages arising out of
the use or inability to use the program (including but not
limited to loss of data or data being rendered inaccurate or
losses sustained by you or third parties or a failure of the
Program to operate with any other programs), even if such
holder or other party has been advised of the possibility of
such damages.
           END OF TERMS AND CONDITIONS

B.2 Contact Information

The best contact path for licensing issues is by e-mail to protomol@cse.nd.edu or send correspondence to:

ProtoMol Team c/o Prof. Jesús A. Izaguirre Laboratory for Computational Life Sciences Department of Computer Science and Engineering University of Notre Dame 384 Fitzpatrick Hall of Engineering Notre Dame, Indiana 46556 USA

Bibliography

60