1% 2% $Id$ 3% 4\label{sec:nwmd} 5%\newcommand{\mc}[3]{\multicolumn{#1}{#2}{#3}} 6\newcommand{\mc}[3]{\mbox{\bf #3}} 7\newcommand{\vb}[1]{\mbox{\verb.#1.}} 8\newcommand{\none}{\multicolumn{2}{|c|}{ }} 9%%%%%%%\renewcommand{\thetable}{\Roman{table}} 10\newcommand{\mcc}[1]{\multicolumn{2}{c}{#1}} 11\def\bmu{\mbox{\boldmath $\mu$}} 12\def\bE{\mbox{\bf E}} 13\def\br{\mbox{\bf r}} 14\def\tT{\tilde{T}} 15\def\t{\tilde{1}} 16\def\ip{i\prime} 17\def\jp{j\prime} 18\def\ipp{i\prime\prime} 19\def\jpp{j\prime\prime} 20\def\etal{{\sl et al.}} 21\def\nwchem{{\bf NWChem}} 22\def\prepare{{\bf prepare}} 23\def\nwtop{{\bf nwtop}} 24\def\nwrst{{\bf nwrst}} 25\def\nwsgm{{\bf nwsgm}} 26 27\section{Introduction} 28 29\subsection{Spacial decomposition} 30The molecular dynamics module of \nwchem\ uses a distribution of data 31based on a spacial decomposition of the molecular system, offering 32an efficient parallel implementation in terms of both memory 33requirements and communication costs, especially for simulations of 34large molecular systems. 35 36Inter-processor communication using the global array tools and the 37design of a data structure allowing distribution based on spacial 38decomposition are the key elements in taking advantage of 39the distribution of memory requirements and computational work with 40minimal communication. 41 42In the spacial decomposition approach, the physical simulation 43volume is divided into rectangular cells, each of which is 44assigned to a processor. Depending on the conditions of the 45calculation and the number of available processors, each processor 46contains one or more of these spacially grouped cells. 47The most important aspects of this decomposition are the dependence 48of the cell sizes and communication cost on the number of processors 49and the shape of the cells, the frequent reassignment of atoms to 50cells leading to a fluctuating number of atoms per cell, and the 51locality of communication which is the main reason for the efficiency 52of this approach for very large molecular systems. 53 54To improve efficiency, molecular systems are broken up into separately 55treated solvent and solute parts. Solvent molecules are assigned to 56the domains according to their center of geometry and are always owned 57by a one node. This avoids solvent--solvent bonded interactions 58crossing node boundaries. Solute molecules are broken up into 59segments, with each segment assigned to a processor based on its 60center of geometry. This limits the number of solute bonded 61interactions that cross node boundaries. The processor to which a 62particular cell is assigned is responsible for the calculation of all 63interactions between atoms within that cell. For the calculation of 64forces and energies in which atoms in cells assigned to different 65processors are involved, data are exchanged between processors. The 66number of neighboring cells is determined by the size and shape of the 67cells and the range of interaction. The data exchange that takes place 68every simulation time step represents the main communication 69requirements. Consequently, one of the main efforts is to design 70algorithms and data structures to minimize the cost of this 71communication. However, for very large molecular systems, memory 72requirements also need to be taken into account. 73 74To compromise between these requirements exchange of data is performed 75in successive point to point communications rather than using the 76shift algorithm which reduces the number of communication calls 77for the same amount of communicated data. 78 79For inhomogeneous systems, the computational load of evaluating 80atomic interactions will generally differ between cell pairs. 81This will lead to load imbalance between processors. 82Two algorithms have been implemented that allow for dynamically 83balancing the workload of each processor. 84One method is the dynamic resizing of cells such that cells gradually 85become smaller on the busiest node, thereby reducing the computational 86load of that node. Disadvantages of this method are that the 87efficiency depends on the solute distribution in the simulation volume 88and the redistribution of work depends on the number of nodes which 89could lead to results that depend on the number of nodes used. 90The second method is based on the dynamic redistribution of intra-node 91cell-cell interactions. This method represents a more coarse load 92balancing scheme, but does not have the disadvantages of the cell 93resizing algorithm. For most molecular systems the cell pair 94redistribution is the more efficient and preferred method. 95 96The description of a molecular system consists of static and dynamic 97information. The static information does not change during a 98simulation and includes items such as connectivity, excluded and third 99neighbor lists, equilibrium values and force constants for all 100bonded and non-bonded interactions. The static information is called 101the topology of the molecular system, and is kept on a separate 102topology file. The dynamic information includes coordinates and 103velocities for all atoms in the molecular system, and is kept in a 104so-called restart file. 105 106\subsection{Topology} 107\label{sec:nwatopology} 108The static information about a molecular system that is needed for 109a molecular simulation is provided to the simulation module in a 110topology file. 111Items in this file include, among many other things, 112a list of atoms, their non-bonded parameters for van der Waals and 113electrostatic interactions, and the complete connectivity in terms 114of bonds, angles and dihedrals. 115 116In molecular systems, a distinction is made between 117{\it solvent} and {\it solute}, which are treated separately. 118A solvent molecule is defined only once in the topology file, 119even though many solvent molecules usually are included in the 120actual molecular system. In the current implementation only one 121solvent can be defined. Everything that is not solvent in the 122molecular system is solute. Each solute atom in the system must 123be explicitly defined in the topology. 124 125Molecules are defined in terms of one or more {\it segment}s. 126Typically, repetitive parts of a molecule are each defined as a single 127segment, such as the amino acid residues in a protein. 128Segments can be quite complicated to define and are, therefore, 129collected in a set of database files. 130The definition of a molecular system in terms of segments is a 131{\it sequence}. 132 133Topology files are created using the \prepare\ module. 134 135\subsection{Files} 136\label{sec:nwafilenames} 137 138File names used have the form \verb+$system$_$calc$.$ext$+, with 139exception of the topology file (Section \ref{sec:nwatopology}), which is named 140\verb+$system$.top+. 141Anything that refers to the definition of the chemical system can be used 142for \verb+$system$+, as long as no periods or underlines are used. 143The identifier \verb+$calc$+ can be anything that refers to the type of 144calculation to be performed for the system with the topology defined. 145This file naming convention allows for the creation of a single 146topology file \verb+$system$.top+ that can be used for a number of 147different calculations, each identified with a different \verb+$calc$+. 148For example, if {\tt crown.top} is the name of the topology file for 149a crown ether, {\tt crown\_em}, {\tt crown\_md}, {\tt crown\_ti} could 150be used with appropriate extensions for the filenames for energy 151minimization, molecular dynamics simulation and multi-configuration 152thermodynamic integration, respectively. All of these calculations 153would use the same topology file {\tt crown.top}. 154 155\label{sec:nwaextensions} 156 157The extensions \verb+<ext>+ identify the kind of information on a file, 158and are pre-determined. 159\begin{table}[htbp] 160\begin{center} 161\begin{tabular}{ll} 162%{\bf acf} & free energy correlation data file\\ 163%{\bf cnv} & free energy convergence file\\ 164{\bf dbg} & debug file\\ 165%{\bf dre} & distance restraints file\\ 166%{\bf fet} & free energy step contribution file\\ 167{\bf frg} & fragment file\\ 168{\bf gib} & free energy data file\\ 169{\bf mri} & free energy multiple run input file\\ 170{\bf mro} & free energy multiple run output file\\ 171{\bf nw} & \nwchem\ input file\\ 172{\bf nwout} & \nwchem\ output file\\ 173{\bf out} & molecular dynamics output file\\ 174{\bf pdb} & PDB formatted coordinate file\\ 175{\bf prp} & property file\\ 176{\bf qrs} & quenched restart file, resulting from an energy minimization\\ 177%{\bf rdf} & radial distribution function output file\\ 178%{\bf rdi} & radial distribution function input file\\ 179{\bf rst} & restart file, used to start or restart a simulation \\ 180{\bf seq} & sequence file, describing the system in segments\\ 181{\bf sgm} & segment file, describing segments\\ 182{\bf syn} & synchronization time file\\ 183{\bf tst} & test file\\ 184{\bf tim} & timing analysis file\\ 185{\bf top} & topology file, contains the static description of a system\\ 186{\bf trj} & trajectory file\\ 187\end{tabular} 188\end{center} 189\caption{List of file extensions for nwchem chemical system files.} 190\end{table} 191 192 193\subsection{Databases} 194Database file supplied with NWChem and used by the \prepare\ module are found in 195directories with name \verb+$ffield$_$level$+, \\ 196where \verb+$ffield$+ is any of the 197supported force fields (Section \ref{sec:nwaforcefields}). 198The source of the data is identified by \verb+$level$+, and can be 199\begin{center} 200\begin{tabular}{ll} 201\hline 202level & Description \\ 203{\bf s} & original published data \\ 204{\bf x} & additional published data \\ 205{\bf q} & contributed data \\ 206{\bf u} & user preferred data \\ 207{\bf t} & user defined temporary data\\ 208\hline 209\end{tabular} 210\end{center} 211 212The user is can replace these directories or add additional database files 213by specifying them in the $.nwchemrc$ file. or in the \verb+prepare+ input file. 214 215The extension \verb+1-9+ defines the priority of database file. 216 217\begin{table}[htbp] 218\begin{center} 219\begin{tabular}{ll} 220{\bf frg} & fragments\\ 221{\bf par} & parameters\\ 222{\bf seq} & sequences\\ 223{\bf sgm} & segments\\ 224\end{tabular} 225\end{center} 226\caption{List of database file extensions.} 227\end{table} 228 229 230The paths of the different database directories should be defined in a file 231{\tt .nwchemrc} in a user's home directory, and provides the user the 232option to select which database files are scanned. 233 234\subsection{Force fields} 235\label{sec:nwaforcefields} 236Force fields recognized are 237\begin{center} 238\begin{tabular}{lll} 239\hline 240Keyword & Force field & Status \\ 241{\tt amber} & AMBER99 & AMBER95,GLYCAM also available\\ 242{\tt charmm} & CHARMM & \\ 243%{\tt cvff} & CVFF & \\ 244%{\tt gromos} & GROMOS87 & \\ 245%{\tt oplsa} & OPLS/AMBER3.0 & \\ 246%{\tt oplsg} & OPLS/GROMOS87 & \\ 247\hline 248\end{tabular} 249\end{center} 250The units for lengths, angles, and energies are correspondingly 251nanometers, radians, and kJ/mol. 252 253\section{Format of fragment files} 254Fragment files contain the basic information needed to specify all 255interactions that need to be considered in a molecular simulation. 256The format of the fragment files is described in Table \ref{tbl:nwmdfrg}. 257Normally these files are created by the \prepare\ module. Manual 258editing is needed when, for example, the \prepare\ module could not 259complete atom typing, or when modified charges are required. 260 261\section{Creating segment files} 262\label{sec:nwanwsgm} 263The \prepare\ module is used to generate segment files 264from corresponding fragment files. A segment file contains all 265information for the calculation of bonded and non-bonded interactions 266for a given chemical system using a specific force field. 267 268Which atoms form a fragment is specified in the coordinate file, 269currently only in PDB format. 270the restriction is that bonded interactions may only involve atoms on at 271most two segments does no longer exist as of \nwchem\ release 3.2.1. 272The segment entries define three sets of parameters 273for each interaction. 274 275Free energy perturbations can be performed using set 1 for the 276generation of the ensemble while using sets 2 and/or 3 277as perturbations. Free energy multiconfiguration thermodynamic 278integration and multistep thermodynamic perturbation calculations are 279performed by gradually changing the interactions in the system from 280parameter set 2 to parameter set 3. These modifications can be 281edited into the segment files manually, or introduced directly into 282the topology file using the \verb+modify+ commands in the input for 283the \prepare\ module. 284 285The format of a segment is 286described in Tables \ref{tbl:nwmdseg1}--\ref{tbl:nwmdseg6}. 287 288\section{Creating sequence files} 289A sequence file describes a molecular system in terms of segments. This 290file is generated by the \prepare\ module for the molecular system 291provided on a PDB-formatted coordinate file. 292The file format is given in Table \ref{tbl:nwmdseq} 293 294\section{Creating topology files} 295\label{sec:nwanwtop} 296 297The topology (Section \ref{sec:nwatopology}) describes all static information 298that describes a molecular system. This includes the connectivity in 299terms of bond-stretching, angle-bending and torsional interactions, as well as 300the non-bonded van der Waals and Coulombic interactions. 301 302The topology of a molecular system is generated by the \prepare\ module 303from the sequence in terms of segments as specified on the PDB file. 304For each unique segment specified in this file the 305segment database directories are searched for the segment definition. 306For segments not found in one of the database directories a segment definition 307is generated in the temporary directory if a fragment file was found. 308If a fragment file could not be found, it is generated by the \prepare\ module 309base on what is found on the PDB file. 310 311When all segments are found or created, the parameter substitutions are 312performed, using force field parameters taken from the parameter 313databases. After all lists have been generated the 314topology is written to a local topology file \verb+$system$.top+. 315 316\section{Creating restart files} 317\label{sec:nwanwrst} 318 319Restart files contain all dynamical information about a molecular 320system and are created by the \prepare\ module if a topology file 321is available. The \prepare\ module will automatically generate 322coordinates for hydrogen atoms and monatomic counter ions 323not found on the PDB formatted coordinate file, if no fragment or 324segment files were generated using that PDB file. 325 326The \prepare\ module has a number of other optional input command, 327including solvation. 328 329\section{Molecular simulations} 330The type of molecular dynamics simulation is specified by the 331\nwchem\ task directive. 332\begin{verbatim} 333task md [ energy | optimize | dynamics | thermodynamics ] 334\end{verbatim} 335where the theory keyword {\tt md} specifies use of the molecular 336dynamics module, and the operation keyword is one of 337\begin{description} 338\item 339{\tt energy} for single configuration energy evaluation 340\item 341{\tt optimize} for energy minimization 342\item 343{\tt dynamics} for molecular dynamics simulations and single step 344thermodynamic perturbation free energy molecular dynamics simulations 345\item 346{\tt thermodynamics} for combined multi-configuration thermodynamic 347integration and multiple step thermodynamic perturbation free 348energy molecular dynamics simulations. 349\end{description} 350 351\section{System specification} 352 353The chemical system for a calculation is specified in the topology 354and restart files. These files should be created using the utilities 355\nwtop\ and \nwrst\ before a simulation can be performed. 356The names of these files are determined from the required \verb+system+ 357directive. 358 359\begin{verbatim} 360system <string systemid>_<string calcid> 361\end{verbatim} 362 363where the strings \verb+systemid+ and \verb+calcid+ are user defined names 364for the chemical system and the type of calculation to ber performed, 365respectively. These names are used to derive the filenames used for the 366calculation. The topoly file used will be \verb+systemid.top+, while all 367other files are named \verb+systemid_calcid.ext+. 368 369\section{Restarting and continuing simulations} 370\begin{description} 371\item 372\begin{verbatim} 373finish 374\end{verbatim} 375specifies that the current job will finish a previous, incomplete 376simulation, using the input data that have been recorded by that 377previous run in the restart file. Most of the input in the current 378md input block will be ignored. 379 380\item 381\begin{verbatim} 382resume 383\end{verbatim} 384specifies that the current job will be an extension of a previous 385simulation, using most of the input data that have been recorded by that 386previous run in the restart file. Typically the input in the current 387md input block defines a larger number of steps than the previous job. 388 389\section{Parameter set} 390 391\item 392\begin{verbatim} 393set <integer iset> 394\end{verbatim} 395specifies the use of parameter set \verb+<iset>+ for the 396molecular dynamics simulation. 397The topology file contains three separate parameters sets that can 398be used. The default for \verb+<iset>+ is 1. 399 400\item 401\begin{verbatim} 402lambda <integer ilambda> <integer ilambda> 403\end{verbatim} 404specifies the use of parameter set for the \verb+ilambda+-th 405of \verb+mlambda+ steps. 406 407\item 408\begin{verbatim} 409pset <integer isetp1> [<integer isetp2>] 410\end{verbatim} 411specifies the parameter sets to be used as perturbation potentials 412in single step thermodynamic perturbation free energy evaluations, 413where \verb+<isetp1>+ specifies the first perturbation parameter set and 414\verb+<isetp2>+ specifies the second perturbation parameter set. Legal 415values for \verb+<isetp1>+ are 2 and 3. Legal value for \verb+<isetp2>+ is 4163, in which case \verb+<isetp1>+ can only be 2. If specified, \verb+<iset>+ 417is automatically set to 1. 418 419\item 420\begin{verbatim} 421pmf [ equilharm <integer npmfc> | scale <real facpmf>] 422\end{verbatim} 423specifies that any potential of mean force functions defined in the 424topology files are to be used. If \verb+equilharm+ is specified, the 425first \verb+npmfc+ dynamics steps will use a harmonic potential 426in stead of any pmf constraint. If \verb+scale+ is specified, all 427pmf force constants are scaled by a factor \verb+facpmf+. 428 429\item 430\begin{verbatim} 431distar [draver [<integer ndaver default 1>]] 432 [scale <real drsscl>] 433 [after <integer nfdrss>] 434\end{verbatim} 435specifies that any distance restraint functions defined in the 436topology files are to be used. 437 438%Optional keywords are \verb+draver+ 439%to specify that average distances over \verb+ndaver+ steps should be 440%used, \verb+scale+ to scale the specified force constants, and 441%\verb+after+ 442 443\item 444\begin{verbatim} 445qhop [<integer nfhop default 10>] 446 [<real rhop default 0.35>] 447 [<real thop default 0.02>] 448\end{verbatim} 449specifies that a Q-HOP simulation is to be carried out with attempted 450proton hops every $nfhop$ steps, a cutoff for the donor-acceptor pair 451distance of $rhop$ nm, and a minimum time before back hopping can occur 452of $thop$ ps. 453 454 455\end{description} 456 457\section{Energy minimization algorithms} 458The energy minimization of the system as found in the restart file 459is performed with the following directives. If both are specified, 460steepest descent energy minimization precedes conjugate gradient 461minimization. 462 463%\begin{description} 464\begin{description} 465\item 466\begin{verbatim} 467sd <integer msdit> [init <real dx0sd>] [min <real dxsdmx>] \ 468 [max <real dxmsd>] 469\end{verbatim} 470specifies the variables for steepest descent energy minimizations, 471where \verb+<msdit>+ is the maximum number of steepest descent steps taken, 472for which the default is 100, \verb+<dx0sd>+ is the initial step size in nm 473for which the default is 0.001, \verb+<dxsdmx>+ is the threshold for the 474step size in nm for which the default is 0.0001, and \verb+<dxmsd>+ is the 475maximum allowed step size in nm for which the default is 0.05. 476\item 477\begin{verbatim} 478cg <integer mcgit> [init <real dx0cg>] [min <real dxcgmx>] \ 479 [cy <integer ncgcy>] 480\end{verbatim} 481specifies the variables for conjugate gradient energy minimizations, 482where \verb+<mcgit>+ is the maximum number of conjugate gradient steps 483taken, for which the default is 100, \verb+<dx0cg>+ is the initial search 484interval size in nm for which the default is 0.001, \verb+<dxcgmx>+ is the 485threshold for the step size in nm for which the default is 0.0001, and 486\verb+<ncgcy>+ is the number of conjugate gradient steps after which the 487gradient history is discarded for which the default is 10. If conjugate 488gradient energy minimization is preceded by steepest descent energy 489minimization, the search interval is set to twice the final step of the 490steepest descent energy minimization. 491\end{description} 492%\end{description} 493 494\section{Multi-configuration thermodynamic integration} 495The following keywords control free energy difference simulations. 496Multi-configuration thermodynamic integrations are always combined 497with multiple step thermodynamic perturbations. 498\begin{description} 499 500\item 501\begin{verbatim} 502(forward | reverse) [[<integer mrun> of] <integer maxlam>] 503\end{verbatim} 504specifies the direction and number of integration steps in free 505energy evaluations, with {\tt forward} being the default direction. 506\verb+<mrun>+ is the number of ensembles that will be generated in 507this calculation, and \verb+<maxlam>+ is the total number of ensembles 508to complete the thermodynamic integration. The default value for 509\verb+<maxlam>+ is 21. The default value of \verb+<mrun>+ is the 510value of \verb+<maxlam>+. 511 512\item 513\begin{verbatim} 514error <real edacq> 515\end{verbatim} 516specifies the maximum allowed statistical error in each generated 517ensemble, where \verb+<edacq>+ is the maximum error allowed in the 518ensemble average derivative of the Hamiltonian with respect to 519$\lambda$ with a default of 5.0 kJ~mol$^{-1}$. 520 521\item 522\begin{verbatim} 523drift <real ddacq> 524\end{verbatim} 525specifies the maximum allowed drift in the free energy result, 526where \verb+<ddacq>+ is the maximum drift allowed in the 527ensemble average derivative of the Hamiltonian with respect to 528$\lambda$with a default of 5.0 kJ~mol$^{-1}$ps$^{-1}$. 529 530\item 531\begin{verbatim} 532factor <real fdacq> 533\end{verbatim} 534specifies the maximum allowed change in ensemble size 535where \verb+<fdacq>+ is the minimum size of an ensemble relative to the 536previous ensemble in the calculation with a default value of 0.75. 537 538\item 539\begin{verbatim} 540decomp 541\end{verbatim} 542specifies that a free energy decomposition is to be carried out. 543Since free energy contributions are path dependent, results from a 544decomposition analysis can no be unambiguously interpreted, and 545the default is not to perform this decomposition. 546 547\item 548\begin{verbatim} 549sss [delta <real delta>] 550\end{verbatim} 551specifies that atomic non-bonded interactions describe a dummy atom 552in either the initial or final state of the thermodynamic calculation 553will be calculated using separation-shifted scaling, where \verb+<delta>+ 554is the separation-shifted scaling factor with a default of 0.075 nm$^2$. 555This scaling method prevents problems associated with singularities in 556the interaction potentials. 557 558\item 559\begin{verbatim} 560new | renew | extend 561\end{verbatim} 562specifies the initial conditions for thermodynamic calculations. 563{\tt new} indicates that this is an initial mcti calculation, which 564is the default. {\tt renew} instructs to obtain the initial 565conditions for each $\lambda$ from the {\bf mro}-file from a previous 566mcti calculation, which has to be renamed to an {\bf mri}-file. The 567keyword {\tt extend} will extend a previous mcti calculation from the 568data read from an {\bf mri}-file. 569\end{description} 570 571\section{Time and integration algorithm directives} 572Following directives control the integration of the equations of motion. 573\begin{description} 574 575\item 576\begin{verbatim} 577leapfrog | leapfrog_bc 578\end{verbatim} 579specifies the integration algorithm, 580where {\tt leapfrog} specifies the default leap frog integration, and 581{\tt leapfrog\_bc} specifies the Brown-Clarke leap frog integrator. 582 583\item 584\begin{verbatim} 585guided [<real fguide default 0.2> [<real tguide default 0.2>]] 586\end{verbatim} 587specifies the use of the guided molecular dynamics simulation 588technique. Variable $fguide$ defines the fraction of the averaged 589forces $g$ to be added to the forces $f^{f}$ evaluated using the force 590field functions to obtain the forces $f$ used to advance the coordinates. 591\begin{equation} 592f_i=f^{f}_i+fguide * g_{i-1} 593\end{equation} 594Variable $tguide$ defines the length of the averaging relative to the 595timestep $\Delta t$. 596\begin{equation} 597g_i = {\Delta t\over tguide} f_i + \left(1- {\Delta t\over tguide}\right) 598g_{i-1} 599\end{equation} 600The current implementation is still under 601development. 602 603\item 604\begin{verbatim} 605equil <integer mequi> 606\end{verbatim} 607specifies the number of equilibration steps \verb+<mequi>+, with a default 608of 100. 609 610\item 611\begin{verbatim} 612data <integer mdacq> [over <integer ldacq>]] 613\end{verbatim} 614specifies the number of data gathering steps \verb+<mdacq>+ with a 615default of 500. In multi-configuration thermodynamic integrations 616\verb+<mequi>+ and \verb+<mdacq>+ are for each of the ensembles, and 617variable \verb+<ldacq>+ specifies the minimum number of data gathering steps 618in each ensemble. In regular molecular dynamics simulations \verb+<ldacq>+ 619is not used. The default value for \verb+<ldacq>+ is the value of \verb+<mdacq>+. 620 621\item 622\begin{verbatim} 623time <real stime> 624\end{verbatim} 625specifies the initial time \verb+<stime>+ of a molecular simulation in ps, 626with a default of 0.0. 627 628\item 629\begin{verbatim} 630step <real tstep> 631\end{verbatim} 632specifies the time step \verb+<tstep>+ in ps, with 0.001 as the default value. 633\end{description} 634 635\section{Ensemble selection} 636Following directives control the ensemble type. 637 638\begin{description} 639 640\item 641\begin{verbatim} 642isotherm [<real tmpext> [<real tmpext2>]] [trelax <real tmprlx> [<real tmsrlx>]] \ 643 [anneal [<real tann1>] <real tann2>] 644\end{verbatim} 645specifies a constant temperature ensemble using Berendsen's thermostat, 646where \verb+<tmpext>+ is the external temperature with a default of 298.15~K, 647and \verb+<tmprlx>+ and \verb+<tmsrlx>+ are temperature relaxation times in ps 648with a default of 0.1. If only \verb+<tmprlx>+ is given the complete system 649is coupled to the heat bath with relaxation time \verb+<tmprlx>+. If both 650relaxation times are supplied, solvent and solute are independently coupled 651to the heat bath with relaxation times \verb+<tmprlx>+ and \verb+<tmsrlx>+, 652respectively. If keyword \verb+anneal+ is specified, the external temperature 653will change from \verb+tmpext+ to \verb+tempext2+ between simulation time 654\verb+tann1+ and \verb+tann2+ 655 656\item 657\begin{verbatim} 658isobar [<real prsext>] [trelax <real prsrlx> ] \ 659 [compress <real compr>] [anisotropic] [xy | z | xy-z] 660\end{verbatim} 661specifies a constant pressure ensemble using Berendsen's piston, 662where \verb+<prsext>+ is the external pressure with a default of 1.025~10$^5$ Pa, 663\verb+<prsrlx>+ is the pressure relaxation time in ps with a default of 0.5, and 664\verb+<compr>+ is the system compressibility in m$^2$N$^{-1}$ with a 665default of 4.53E-10. Optional keywords \verb+xy+, \verb+z+ and \verb+xy-z+ 666may be used to specify that pressure scaling is to be applied in 667the \verb+x+ and \verb+y+ dimension only, the \verb+z+ dimension only, or, 668in all three dimensions with identical scaling in the \verb+x+ and \verb+y+ 669dimension. The last option requires that \verb+anisotropic+ is also specified. 670 671\end{description} 672 673\section{Velocity reassignments} 674Velocities can be periodically reassigned to reflect a certain temperature. 675\begin{description} 676\item 677\begin{verbatim} 678vreass <integer nfgaus> <real tgauss> 679 [fraction [<real frgaus default 0.5]] 680 [once] 681 [(first | initial)] [(last | final)] 682\end{verbatim} 683specifies that velocities will be reassigned every \verb+<nfgaus>+ molecular 684dynamics steps, reflecting a temperature of \verb+<tgauss>+~K. The default 685is not to reassign velocities, i.e.\ \verb+<nfgaus>+ is 0. Keyword 686\verb+fraction+ allows the specification of the fraction of the new 687velocities are random. Keyword \verb+once+ specifies that velocity 688reassignment only should be done in the first step. Keywords \verb+first+ 689or \verb+initial+ and \verb+last+ or \verb+final+ specify that 690velocity reassigment should only be applied in the first and last 691window of multiple run simulations. 692 693\end{description} 694 695\section{Cutoff radii} 696Cutoff radii can be specified for short range and long range interactions. 697\begin{description} 698\item 699\begin{verbatim} 700cutoff [short] <real rshort> [long <real rlong>] \ 701 [qmmm <real rqmmm>] 702\end{verbatim} 703specifies the short range cutoff radius \verb+<rshort>+, and the long range 704cutoff radius \verb+<rlong>+ in nm. If the long range cutoff radius 705is larger than the short range cutoff radius the twin range method will 706be used, in which short range forces and energies are evaluated every 707molecular dynamics step, and long range forces and energies with a 708frequency of \verb+<nflong>+ molecular dynamics steps. Keyword 709\verb+qmmm+ specifies the radius of the zone around quantum atoms 710defining the QM/MM bare charges. 711The default value for \verb+<rshort>+, \verb+<rlong>+ and \verb+<rqmmm>+ 712is 0.9~nm. 713\end{description} 714 715\section{Polarization} 716First order and self consistent electronic polarization models have 717been implemented. 718\begin{description} 719\item 720\begin{verbatim} 721polar (first | scf [[<integer mpolit>] <real ptol>]) 722\end{verbatim} 723specifies the use of polarization potentials, 724where the keyword {\tt first} specifies the first order polarization 725model, and {\tt scf} specifies the self consistent polarization field 726model, iteratively determined with a maximum of \verb+<mpolit>+ 727iterations to within a tolerance of \verb+<ptol>+ D in the generated 728induced dipoles. The default is not to use polarization models. 729\end{description} 730 731\section{External electrostatic field} 732\begin{description} 733 734\item 735\begin{verbatim} 736field <real xfield> [freq <real xffreq>] [vector <real xfvect(1:3)>] 737\end{verbatim} 738specifies an external electrostatic field, 739where \verb+<xfield>+ is the field strength, \verb+<xffreq>+ is the 740frequency in MHz and \verb+<xfvect>+ is the external field vector. 741\end{description} 742 743\section{Constraints} 744Constraints are satisfied using the SHAKE 745coordinate resetting procedure. 746\begin{description} 747 748\item 749\begin{verbatim} 750shake [<integer mshitw> [<integer mshits>]] \ 751 [<real tlwsha> [<real tlssha>]] 752\end{verbatim} 753specifies the use of SHAKE constraints, 754where \verb+<mshitw>+ is the maximum number of solvent SHAKE iterations, 755and \verb+<mshits>+ is the maximum number of solute SHAKE iterations. If 756only \verb+<mshitw>+ is specified, the value will also be used for \verb+<mshits>+. 757The default maximum number of iterations is 100 for both. 758\verb+<tlwsha>+ is the solvent SHAKE tolerance in nm, and \verb+<tlssha>+ is 759the solute SHAKE tolerance in nm. If only \verb+<tlwsha>+ is specified, the 760value given will also be used for \verb+<tlssha>+. The default tolerance 761is 0.001~nm for both. 762 763\item 764\begin{verbatim} 765noshake (solvent | solute) 766\end{verbatim} 767disables SHAKE and treats the bonded interaction according to the force 768field. 769 770\end{description} 771 772\section{Long range interaction corrections} 773Long range electrostatic interactions are implemented using the 774smooth particle mesh Ewald technique, for neutral periodic cubic systems in 775the constant volume ensemble, using pair interaction potentials. Particle-mesh 776Ewald long range interactions can only be used in molecular dynamics simulations 777using effective pair potentials, and not in free energy simulations, QMD or 778QM/MM simulations. 779 780\begin{description} 781\item 782\begin{verbatim} 783pme [grid <integer ng>] [alpha <real ealpha>] \ 784 [order <integer morder>] [fft <integer imfft>]\ 785 [procs <integer nprocs>] [solvent] 786\end{verbatim} 787specifies the use of smooth particle-mesh Ewald long range 788interaction treatment, 789where \verb+ng+ is the number of grid points per dimension, 790\verb+ealpha+ is the Ewald coefficient in nm$^{-1}$, with a default 791that leads to a tolerance of $10^{-4}$ at the short range cutoff radius, 792and \verb+morder+ is order of the Cardinal B-spline 793interpolation which must be an even number and at least 4 (default 794value). A platform specific 3D fast Fourier transform is used, if 795available, when \verb+imfft+ is set to 2. \verb+nprocs+ can be used to 796define a subset of processors to be used to do the FFT calculations. 797If \verb+solvent+ is specified, the charge grid will be calculated from 798the solvent charges only. 799 800\item 801\begin{verbatim} 802react [<real dielec default 80.0>] 803\end{verbatim} 804specifies that a simple reaction field correction is used with a 805dielectric constant \verb+dielec+. This is an 806experimental option that has not been well tested. 807 808\end{description} 809 810\section{Fixing coordinates} 811Solvent or solute may be fixed using the following keywords. 812\begin{description} 813\item 814\begin{verbatim} 815( fix | free ) 816 solvent ( [<integer idfirst> [<integer idlast>]] | 817 ( within | beyond) <real rfix> <string atomname> ) | \ 818 solute ( [<integer idfirst> [<integer idlast>]] [ heavy | {<string atomname>}] | 819 ( within | beyond) <real rfix> <string atomname> ) 820 [permanent] 821\end{verbatim} 822For solvent the molecule numbers \verb+idfirst+ and \verb+idlast+may be 823specified to be the first and last molecule to which the directive 824applies. If omitted, the directive applies to all molecules. For solute, 825the segment numbers \verb+idfirst+ and \verb+idlast+may be 826specified to be the first and last segment to which the directive 827applies. If omitted, the directive applies to all segments. In addition, 828the keyword \verb+heavy+ may be specified to apply to all non hydrogen 829atoms in the solute, or a set of atom names may be specified in which 830a wildcard character \verb+?+ may be used. Keyword \verb+permanent+ 831is used to keep the specification on the restart file for subsequent 832simulations. 833\end{description} 834 835\section{Special options} 836\begin{description} 837 838\item 839\begin{verbatim} 840import [<integer impfr default 1> [<integer impto default impfr> \ 841 [<integer nftri default 1>]]] 842\end{verbatim} 843specifies the import of frames \verb+impfr+ to \verb+impto+ with 844frequency \verb+nftri+ from a trajectory file with extension 845\verb+tri+ for which energies and forces are to be recalculated. 846This option only applied to \verb+task md energy+. 847 848\item 849\begin{verbatim} 850detail 851\end{verbatim} 852specifies that moments of inertia and radii of gyration will be part of the 853recorded properties. 854 855\item 856\begin{verbatim} 857profile 858\end{verbatim} 859specifies that execution time profiling data will be part of the 860recorded properties. 861 862\item 863\begin{verbatim} 864scale <real scaleq> 865\end{verbatim} 866specifies that all charges will be scaled by the factro \verb+scaleq+. 867 868\item 869\begin{verbatim} 870collapse [<real fcoll default 10.0> [ segment | z | xy ] 871\end{verbatim} 872specifies that additional forces directed to the origin of the 873simulation cell with strength \verb+fcoll+ will be 874applied to all solute molecules. If \verb+z+ or \verb+xy+ is 875specified, these forces will only apply in the specified dimension(s). 876 877 878\item 879\begin{verbatim} 880include fixed 881\end{verbatim} 882specifies that energies will be evaluated between fixed atoms. 883Normally these interactions are excluded from the pairlists. 884 885 886\item 887\begin{verbatim} 888eqm <real eqm> 889\end{verbatim} 890specifies the zero point of energy in QMD simulations. 891 892\item 893\begin{verbatim} 894atomlist 895\end{verbatim} 896specifies that pairlists will be atom based. Normally pairlist 897are charge group based. 898 899\end{description} 900 901\section{Autocorrelation function} 902For the evaluation of the statistical error of multi-configuration 903thermodynamic integration free energy results a correlated data 904analysis is carried out, involving the calculation of the 905autocorrelation function of the derivative of the Hamiltonian with 906respect to the control variable $\lambda$. 907\begin{description} 908 909\item 910\begin{verbatim} 911auto <integer lacf> [fit <integer nfit>] [weight <real weight>] 912\end{verbatim} 913controls the calculation of the autocorrelation, 914where \verb+<lacf>+ is the length of the autocorrelation function, with 915a default of 1000, \verb+<nfit>+ is the number of functions used in the 916fit of the autocorrelation function, with a default of 15, and 917\verb+<weight>+ is the weight factor for the autocorrelation function, 918with a default value of 0.0. 919\end{description} 920 921\section{Print options} 922Keywords that control print to the output file, with extension {\bf out}. 923Print directives may be combined to a single directive. 924 925\begin{description} 926 927\item 928\begin{verbatim} 929print [topol [nonbond] [solvent] [solute]] \ 930 [step <integer nfoutp> [extra] [energy]] \ 931 [stat <integer nfstat>] \ 932 [energies [<integer nfener>]] \ 933 [forces [<integer nfforce>]] \ 934 [matrix] \ 935 [expect <integer npxpct>] \ 936 [timing] \ 937 [pmf [<integer iprpmf>]] \ 938 [out6] \ 939 [dayout] 940\end{verbatim} 941 942Keyword \verb+topol+ specifies printing the topology information, 943where {\tt nonbond} refers to the non-bonded interaction parameters, 944{\tt solvent} to the solvent bonded parameters, and {\tt solute} to the 945solute bonded parameters. If only {\tt topol} is specified, all 946topology information will be printed to the output file. 947 948Keyword \verb+step+ 949specifies the frequency \verb+nfoutp+ of printing molecular dynamics step 950information to the output file. If the keyword {\tt extra} is specified 951additional energetic data are printed for solvent and solute separately. 952If the keyword {\tt energy} is specified, information is printed for 953all bonded solute interactions. 954The default for \verb+nfoutp+ is 0. For molecular dynamics simulations 955this frequency is in time steps, and for multi-configuration thermodynamic 956integration in $\lambda$-steps. 957 958Keyword \verb+stat+ 959specifies the frequency \verb+<nfstat>+ of printing statistical information 960of properties that are calculated during the simulation. 961For molecular dynamics simulation 962this frequency is in time steps, for multi-configuration thermodynamic 963integration in $\lambda$-steps. 964 965Keyword \verb+energies+ 966specifies the frequency \verb+nfener+ of printing solute bonded energies 967the output file for energy/import calculations. 968The default for \verb+nfener+ is 0. 969 970Keyword \verb+forces+ 971specifies the frequency \verb+nfforc+ of printing solute forces 972the output file for energy/import calculations. 973The default for \verb+nfforc+ is 0. 974 975Keyword \verb+matrix+ specifies that a solute distance matrix is to 976be printed. 977 978Keyword \verb+expect+ is obsolete. 979 980Keyword \verb+timing+ specifies that timing data is printed. 981 982Keyword \verb+pmf+ specifies that pmf data is printed every \verb+iprpmf+ 983steps. 984Keyword \verb+out6+ specifies that output is written to standard out in stead 985of the output file with extension \verb+out+. 986 987Keyword \verb+dayout+ is obsolete. 988 989 990\end{description} 991 992\section{Periodic updates} 993Following keywords control periodic events during a molecular 994dynamics or thermodynamic integration simulation. 995Update directives may be combined to a single directive. 996 997\begin{description} 998\item 999\begin{verbatim} 1000update [pairs <integer nfpair default 1>] \ 1001 [long <integer nflong default 1>] \ 1002 [center <integer nfcntr default 0> [zonly | xyonly] \ 1003 [fraction <integer idscb(1:5)>] \ 1004 [motion <integer nfslow default 0>] \ 1005 [analysis <integer nfanal default 0>] \ 1006 [rdf <integer nfrdf default 0> \ 1007 [range <real rrdf>] [bins <integer ngl>] \ 1008\end{verbatim} 1009 1010Keyword \verb+pairs+ 1011specifies the frequency \verb+<nfpair>+ in molecular dynamics steps of 1012updating the pair lists. The default for the frequency is 1. 1013In addition, pair lists are also updated after each step in which 1014recording of the restart or trajectory files is performed. Updating 1015the pair lists includes the redistribution of atoms that changed 1016domain and load balancing, if specified. 1017 1018Keyword \verb+long+ 1019specifies the frequency \verb+<nflong>+ in molecular dynamics steps 1020of updating the long range forces. The default frequency is 1. 1021The distinction of short range and long range forces is only 1022made if the long range cutoff radius was specified to be larger 1023than the short range cutoff radius. Updating the long range forces 1024is also done in every molecular dynamics step in which the 1025pair lists are regenerated. 1026 1027Keywrod \verb+center+ 1028specifies the frequency \verb+<nfcntr>+ in molecular dynamics steps in 1029which the center of geometry of the solute(s) is translated to the 1030center of the simulation volume. Optional keyword \verb+zonly+ or 1031\verb+xyonly+ can be used to specify that centering will take place in 1032the z-direction or in the xy-plane only. 1033The solute fractions determining the 1034solutes that will be centered are specified by the keyword 1035{\tt fraction} and the vector \verb+<idscb>+, with a maximum of 5 entries. 1036This translation is implemented such that it has no effect on any 1037aspect of the simulation. The default is not to center, i.e. nfcntr is 10380. The default fraction used to center solute is 1. 1039 1040Keyword \verb+motion+ 1041specifies the frequency \verb+<nfslow>+ in molecular dynamics steps of 1042removing the overall rotational and center of mass translational motion. 1043 1044Keyword \verb+analysis+ 1045specifies the frequency \verb+<nfanal>+ in molecular dynamics steps of 1046invoking the analysis module. This option is obsolete. 1047 1048Keyword \verb+rdf+ 1049specifies the frequency \verb+<nfrdf>+ in molecular dynamics steps of 1050calculating contributions to the radial distribution functions. 1051The default is 0. The range of the radial distribution 1052functions is given by \verb+<rrdf>+ in nm, with a default of the short 1053range cutoff radius. Note that radial distribution functions are not 1054evaluated beyond the short range cutoff radius. The number of 1055bins in each radial distribution function is given by \verb+<ngl>+, with 1056a default of 1000. This option is no longer supported. 1057If radial distribution function are to be 1058calculated, a {\bf rdi} files needs to be available in which the 1059contributions are specified as follows. 1060\begin{center} 1061\begin{tabular}{lll} 1062\hline\hline 1063Card & Format & Description \\ \hline 1064I-1 & i & Type, 1=solvent-solvent, 2=solvent-solute, 10653-solute-solute\\ 1066I-2 & i & Number of the rdf for this contribution\\ 1067I-3 & i & First atom number \\ 1068I-4 & i & Second atom number \\ 1069\hline 1070\end{tabular} 1071\end{center} 1072\end{description} 1073 1074\section{Recording} 1075The following keywords control recording data to file. 1076Record directives may be combined to a single directive. 1077 1078\begin{description} 1079 1080\item 1081\begin{verbatim} 1082record [rest <integer nfrest> [keep]] \ 1083 [coord <integer nfcoor default 0>] \ 1084 [wcoor <integer nfwcoo default 0>] \ 1085 [scoor <integer nfscoo default 0>] \ 1086 [veloc <integer nfvelo default 0>] \ 1087 [wvelo <integer nfwvel default 0>] \ 1088 [svelo <integer nfsvel default 0>] \ 1089 [force <integer nfvelo default 0>] \ 1090 [wforc <integer nfwvel default 0>] \ 1091 [sforc <integer nfsvel default 0>] \ 1092 [(prop | prop_average) <integer nfprop default 0>] \ 1093 [free <integer nffree default 1>] \ 1094 [sync <integer nfsync default 0>] \ 1095 [times <integer nftime default 0>] \ 1096 [acf] [cnv] [fet] 1097 [binary] [ascii] [ecce] [argos] 1098\end{verbatim} 1099 1100Keyword \verb+rest+ 1101specifies the frequency \verb+<nfrest>+ in molecular dynamics steps 1102of rewriting the restart file, with extension \verb+rst+. 1103For multi-configuration 1104thermodynamic integration simulations the frequency is in 1105steps in $\lambda$. The default is not to record. The restart 1106file is used to start or restart simulations. The keyword {\tt keep} 1107causes all restart files written to be kept on disk, rather than 1108to be overwritten. 1109 1110Keyword \verb+coord+ 1111specifies the frequency \verb+<nfcoor>+ in molecular dynamics steps 1112of writing coordinates to the trajectory file. This directive redefines 1113previous \verb+coord+, \verb+wcoor+ and \verb+scoor+ 1114directives. The default is not to record. 1115 1116Keyword \verb+wcoor+ 1117specifies the frequency \verb+<nfcoor>+ in molecular dynamics steps 1118of writing solvent coordinates to the trajectory file. This keyword 1119takes precedent over \verb+coord+. This directive redefines 1120previous \verb+coord+, \verb+wcoor+ and \verb+scoor+ 1121directives. 1122The default is not to record. 1123 1124Keyword \verb+scoor+ 1125specifies the frequency \verb+<nfscoo>+ in molecular dynamics steps 1126of writing solute coordinates to the trajectory file. This keyword 1127takes precedent over \verb+coord+. This directive redefines 1128previous \verb+coord+, \verb+wcoor+ and \verb+scoor+ 1129directives. 1130The default is not to record. 1131 1132Keyword \verb+veloc+ 1133specifies the frequency \verb+<nfvelo>+ in molecular dynamics steps 1134of writing velocities to the trajectory file. This directive redefines 1135previous \verb+veloc+, \verb+wvelo+ and \verb+svelo+ 1136directives. 1137The default is not to record. 1138 1139Keyword \verb+wvelo+ 1140specifies the frequency \verb+<nfvelo>+ in molecular dynamics steps 1141of writing solvent velocitiesto the trajectory file. This keyword 1142takes precedent over \verb+veloc+. This directive redefines 1143previous \verb+veloc+, \verb+wvelo+ and \verb+svelo+ 1144directives. 1145The default is not to record. 1146 1147Keyword \verb+svelo+ 1148specifies the frequency \verb+<nfsvel>+ in molecular dynamics steps 1149of writing solute velocities to the trajectory file. This keyword 1150takes precedent over \verb+veloc+. This directive redefines 1151previous \verb+veloc+, \verb+wvelo+ and \verb+svelo+ 1152directives. 1153The default is not to record. 1154 1155Keyword \verb+force+ 1156specifies the frequency \verb+<nfvelo>+ in molecular dynamics steps 1157of writing forces to the trajectory file. This directive redefines 1158previous \verb+vforce+, \verb+wforc+ and \verb+sforc+ 1159directives. 1160The default is not to record. 1161 1162Keyword \verb+wforc+ 1163specifies the frequency \verb+<nfvelo>+ in molecular dynamics steps 1164of writing solvent forcesto the trajectory file. This keyword 1165takes precedent over \verb+force+. This directive redefines 1166previous \verb+vforce+, \verb+wforc+ and \verb+sforc+ 1167directives. 1168The default is not to record. 1169 1170Keyword \verb+sforc+ 1171specifies the frequency \verb+<nfsvel>+ in molecular dynamics steps 1172of writing solute forces to the trajectory file. This keyword 1173takes precedent over \verb+force+. This directive redefines 1174previous \verb+vforce+, \verb+wforc+ and \verb+sforc+ 1175directives. 1176The default is not to record. 1177 1178Keyword \verb+prop+ 1179specifies the frequency \verb+<nfprop>+ in molecular dynamics steps 1180of writing information to the property file, with extension 1181\verb+prp+. The default is not to record. 1182 1183Keyword \verb+prop_average+ 1184specifies the frequency \verb+<nfprop>+ in molecular dynamics steps 1185of writing average information to the property file, with extension 1186\verb+prp+. 1187The default is not to record. 1188 1189Keyword \verb+free+ 1190specifies the frequency \verb+<nffree>+ in multi-configuration 1191thermodynamic integration steps to record data to the 1192free energy data file, with extension \verb+gib+. 1193The default is 1, i.e.\ to record at every $\lambda$. 1194This option is obsolete. All data are required to do the 1195final analysis. 1196 1197%\item 1198%\begin{verbatim} 1199%record cnv 1200%\end{verbatim} 1201%specifies that free energy convergence data will be written to the 1202%free energy convergence file, with extension \verb+cnv+. 1203% 1204%\item 1205%\begin{verbatim} 1206%record acf 1207%\end{verbatim} 1208%specifies that free energy derivative autocorrelation data will be 1209%written to the free energy autocorrelation file, with extension 1210%\verb+acf+. 1211% 1212%\item 1213%\begin{verbatim} 1214%record fet 1215%\end{verbatim} 1216%that free energy vs.\ time data will be recorded to the free energy 1217%data file, with extension \verb+fet+. 1218% 1219 1220Keyword \verb+sync+ 1221specifies the frequency \verb+<nfsync>+ in molecular dynamics steps 1222of writing information to the synchronization file, with extension 1223\verb+syn+. 1224The default is not to record. 1225The information written is the simulation time, the wall clock time 1226of the previous MD step, the wall clock time of the previous force 1227evaluation, the total synchronization time, the largest 1228synchronization time and the node on which the largest synchronization 1229time was found. The recording of synchronization times is part of the 1230load balancing algorithm. Since load balancing is only performed when 1231pair-lists are updated, the frequency \verb+<nfsync>+ is correlated 1232with the frequency of pair-list updates \verb+<nfpair>+. This directive 1233is only needed for analysis of the load balancing performance. For 1234normal use this directive is not used. 1235 1236Keyword \verb+times+ 1237specifies the frequency \verb+<nfsync>+ in molecular dynamics steps 1238of writing information to the timings file, with extension 1239\verb+tim+. 1240The default is not to record. 1241The information written is wall clock time used by each of the 1242processors for the different components in the force evaluation. 1243This directive is only needed for analysis of the wall clock time 1244distribution. For normal use this directive is not used. 1245 1246Keywords \verb+acf+, \verb+cnv+ and \verb+fet+ are obsolete. 1247 1248Keywords \verb+binary+, \verb+ascii+, \verb+ecce+ and \verb+argos+ are obsolete. 1249 1250\end{description} 1251 1252\section{Program control options} 1253\begin{description} 1254\item 1255\begin{verbatim} 1256load [reset] 1257 ( none | 1258 size [<real factld>] | 1259 sizez [<real factld>] | pairs | 1260 (pairs [<integer ldpair>] size [<real factld>]) ) 1261 [last] 1262 [minimum] 1263 [average] 1264 [combination] 1265 [iotime] 1266 [experimental] 1267 1268\end{verbatim} 1269determines the type of dynamic load balancing performed, 1270where the default is {\tt none}. Load balancing option {\tt size} 1271is resizing cells on a node, and {\tt pairs} redistributes the 1272cell-cell interactions over nodes. Keyword \verb+reset+ will reset the 1273load balancing read from the restart file. The level of cell resizing 1274can be influenced with $factld$. The cells on the busiest node are 1275resized with a factor 1276\begin{equation} 1277\left( 1 - factld * { {T_{sync} \over n_p} - t^{min}_{sync} \over t_{wall}} 1278\right)^{1\over 3} 1279\end{equation} 1280where $T_{sync}$ is the accumulated synchronization time of all nodes, 1281$n_p$ is the total number of nodes, $t^{min}_{sync}$ is the synchronization 1282time of the busiest node, and $t_{wall}$ is the wall clock time of the 1283molecular dynamics step. 1284For the combined load balancing, \verb+ldpair+ is the number of successive pair 1285redistribution load balancing steps in which the accumulated synchronization 1286time increases, before a resizing load balancing step will be attempted. 1287Load balancing is only performed in molecular dynamics steps in which the 1288pair-list is updated. The default load balancing is equivalent to specifying\\ 1289\begin{verbatim} 1290load pairs 10 size 0.75 1291\end{verbatim} 1292Keyword \verb+last+ specifies that the load balancing is based on the 1293synchronization times of the last step. This is the default. 1294Keyword \verb+average+ specifies that the load balancing is based on the 1295average synchronization times since the last load balancing step. 1296Keyword \verb+minimum+ specifies that the load balancing is based on the 1297minimum synchronization times since the last load balancing step. 1298Keywords \verb+combination+, \verb+iotime+ and \verb+experimental+ are 1299experimental load balancing options that should not be used in 1300production runs. 1301 1302\item 1303\begin{verbatim} 1304(pack | nopack) 1305\end{verbatim} 1306specifies if data are communicated in packed or unpacked form. The 1307default is \verb+pack+. 1308 1309\item 1310\begin{verbatim} 1311procs <integer npx> <integer npy> <integer npz> 1312\end{verbatim} 1313specifies the distribution of the available processors over the three 1314Cartesian dimensions. The default distribution is chosen such that, 1315\verb+<npx>+$*$\verb+<npy>+$*$\verb+<npz>+=\verb+<np>+ 1316and \verb+<npx>+ $<=$ \verb+<npy>+ $<=$ \verb+<npz>+, 1317where \verb+<npx>+, \verb+<npy>+ and \verb+<npz>+ are the processors in the 1318\verb+x+, \verb+y+ and \verb+z+ dimension respectively, and \verb+<np>+ is the number of processors 1319allocated for the calculation. Where more than one combination 1320of \verb+<npx>+, \verb+<npy>+ and \verb+<npz>+ are possible, the 1321combination is chosen with the minimum value of 1322\verb+<npx>+$+$\verb+<npy>+$+$\verb+<npz>+. To change the default setting 1323the following optional input option is provided. 1324 1325\item 1326\begin{verbatim} 1327cells <integer nbx> <integer nby> <integer nbz> 1328\end{verbatim} 1329specifies the distribution of cells, 1330where \verb+<nbx>+, \verb+<nby>+ and \verb+<nbz>+ are the number of 1331cells in \verb+x+, \verb+y+ and \verb+z+ direction, respectively. 1332The molecular system is decomposed into cells that form the smallest 1333unit for communication of atomic data between nodes. The size of the 1334cells is per default set to the short-range cutoff radius. If 1335long-range cutoff radii are used the cell size is set to half the 1336long-range cutoff radius if it is larger than the short-range cutoff. 1337If the number of cells in a dimension is less than the number of 1338processors in that dimension, the number of cells is set to the number 1339of processors. 1340 1341\item 1342\begin{verbatim} 1343extra <integer madbox> 1344\end{verbatim} 1345sets the number of additional cells for which memory is allocated. 1346In rare events the amount of memory set aside per node is insufficient 1347to hold all atomic coordinates assigned to that node. This leads to 1348execution which aborts with the message that {\tt mwm} or {\tt msa} is too 1349small. Jobs may be restarted with additional space allocated by 1350where \verb+<madbox>+ is the number of additional cells that are allocated 1351on each node. The default for \verb+<madbox>+ is 6. 1352In some cases \verb+<madbox>+ can be reduced to 4 if memory usage is a 1353concern. Values of 2 or less will almost certainly result in memory 1354shortage. 1355 1356\item 1357\begin{verbatim} 1358mwm <integer mwmreq> 1359\end{verbatim} 1360sets the maximum number of solvent molecules \verb+<mwmreq>+ per node, 1361allowing increased memory to be allocated for solvent molecules. This 1362option can be used if execution aborted because \verb+mwm+ was too 1363small. 1364 1365\item 1366\begin{verbatim} 1367msa <integer msareq> 1368\end{verbatim} 1369sets the maximum number of solute atoms \verb+<msareq>+ per node, 1370allowing increased memory to be allocated for solute atoms. This 1371option can be used if execution aborted because \verb+msa+ was too 1372small. 1373 1374\item 1375\begin{verbatim} 1376mcells <integer mbbreq> 1377\end{verbatim} 1378sets the maximum number of cell pairs \verb+<mbbreq>+ per node, 1379allowing increased memory to be allocated for the cell pair lists. 1380This option can be used if execution aborted because \verb+mbbl+ was too 1381small. 1382 1383\item 1384\begin{verbatim} 1385boxmin <real rbox> 1386\end{verbatim} 1387sets the minimum size of a cell. This directive is obsolete. The 1388use of mcells is preferred. 1389 1390\item 1391\begin{verbatim} 1392segmentsize <real rsgm> 1393\end{verbatim} 1394sets the maximum size of a segment. This value is used to determine 1395which segments at the boundary of the cutoff radius should be considered 1396in the generation of the pairlists. This value is also determined by the 1397prepare module and written to the restart file. Use of this directive 1398is not needed for simulations that use the current prepare module to 1399generate the restart file. 1400 1401\item 1402\begin{verbatim} 1403memory <integer memlim> 1404\end{verbatim} 1405sets a limit \verb+<memlim>+ in kB on the allocated amount of memory used by 1406the molecular dynamics module. 1407Per default all available memory is allocated. Use of this command 1408is required for QM/MM simulations only. 1409 1410\item 1411\begin{verbatim} 1412expert 1413\end{verbatim} 1414enables the use of certain combinations of features that are considered 1415unsafe. This directive should not be used for production runs. 1416 1417\item 1418\begin{verbatim} 1419develop <integer idevel> 1420\end{verbatim} 1421enables the use of certain development options specified by the 1422integer \verb+idevel+. This option is for development purposes only, 1423and should not be used for production runs. 1424 1425\item 1426\begin{verbatim} 1427control <integer icntrl> 1428\end{verbatim} 1429enables the use of certain development options specified by the 1430integer \verb+icntrl+. This option is for development purposes only, 1431and should not be used for production runs. 1432 1433\item 1434\begin{verbatim} 1435numerical 1436\end{verbatim} 1437writes out analytical and finite difference forces for test purposes. 1438 1439\item 1440\begin{verbatim} 1441server <string servername> <integer serverport> 1442\end{verbatim} 1443allows monitoring over a socket connection to the specified port on the 1444named server of basic data as a simulation is running. 1445 1446\item 1447For development purposes debug information can be written to the debug 1448file with extension {\bf dbg} with 1449 1450\begin{verbatim} 1451debug <integer idebug> 1452\end{verbatim} 1453where $idebug$ specifies the type of debug information being written. 1454 1455\item 1456For testing purposes test information can be written to the test 1457file with extension {\bf tst} with 1458 1459\begin{verbatim} 1460test <integer itest> 1461\end{verbatim} 1462where $itest$ specifies the number of steps test information is written. 1463 1464\item 1465On some platforms prefetching of data can improve the efficiency. This 1466feature can be turned on using 1467 1468\begin{verbatim} 1469prefetch [<integer nbget>] 1470\end{verbatim} 1471where $nbget$ is the number of outstanding communication operations. 1472 1473\item 1474Application of periodic boundary conditions for the evaluation of 1475forces can be controlled with 1476 1477\begin{verbatim} 1478pbc ( atom | residue | molecule ) 1479\end{verbatim} 1480 1481This option rarely needs to be used. 1482 1483\item 1484Autocorrelation functions for error analysis are controlled using 1485 1486\begin{verbatim} 1487auto [ fit <integer iapprx> | weight <real weight> ] 1488\end{verbatim} 1489 1490This option is disabled in the current release. 1491 1492\item 1493Membrane system equilibration can be made more efficient using 1494 1495\begin{verbatim} 1496membrane [ rotations ] 1497\end{verbatim} 1498 1499\item 1500Constraining the center of mass of solute molecules in the xy plane is 1501accomplished using 1502 1503\begin{verbatim} 1504scmxy [<integer icmopt default 1>] 1505\end{verbatim} 1506 1507where $icmopt$ determines if the constraint is mass weighted (2). 1508 1509\item 1510Radius of gyration calculations are enabled using 1511 1512\begin{verbatim} 1513radius_gyration 1514\end{verbatim} 1515 1516\item 1517Calculations of diffusion coefficients is enabled using 1518 1519\begin{verbatim} 1520diffusion 1521\end{verbatim} 1522This option is disabled in the current release. 1523 1524\item 1525 1526\begin{verbatim} 1527comlim ( on | off ) 1528\end{verbatim} 1529is disabled 1530 1531\item 1532 1533To limit the size of recoding files, new files are opened every $nfnewf$ md steps using 1534 1535\begin{verbatim} 1536batch <integer nfnewf> 1537\end{verbatim} 1538 1539 1540\end{description} 1541