Input options
The overall simulations parameters need to be set in an input file. Standard this file is assumed to be named input.pif. However, every other name is possible, although, when starting PumMa, an extra argument has to be supplied. The input.pif is the first file being parsed by PumMa.
Depending on the compiled features of PumMa, some input options might be ignored. For instance, in the default PumMa code, the macroscopic quantitity analysis and the Lennard-Jones energy analysis are disabled. However, special binaries with these options do exist.
The format of the input file is simple. Every option is specified on a separate line, with first the option name and second the option value (separated by at least one space). Blank lines are ignored and lines starting with a # are treated as comments. A possible third argument supplied on one line is ignored. The input file is not case sensitive with respect to the options. If an option is specified multiple times, PumMa terminates.
Options which represent physical relevant values (such as time, temperature or pressure) are to be given in SI-units. The relevant units are discussed together with the description of the option.
First an example file is shown, thereafter all basic input options are discussed and finally some more specific input options are introduced, which require special compiled PumMa code.
Example file
To give a quick impression on the input file, an example is showed below. This file has been used by the PumMa development group for their bilayer simulations, see the literature section. The meaning of each option is discussed after the example.
# Integrator IT 0 max_IT 500000 Dt 0.012 # Temperature & Pressure T 310 lambdac 0.005 P 1.0 muc 0.000001 Prep 25 anisopress on # File dumping Erep 25 Crep 1000 sysfile sys outpfile output.pof
All options
Below all options of the input file and their default values are discussed. For clarity in this manual, the input options are discussed groupwise.
Integrator
IT 0 # Start iteration max_IT 1 # Final iteration Dt 0.001 # Time step (in picoseconds) integration leapfrog # Integrator type bbkseed 0 # Random seed for Langevin dynamics
These five options control the integrator behavior. Both the starting and final iteration number are to be specified. Hence, a total of max_IT-IT steps are performed. By specifying the starting iteration explicitely it is possible to continue a simulation from that time step, without the need to change any file names. However, possible output files will be overwritten, so backing them up is necessary. The time step is to be specified in picoseconds. The integration option allows the user to choose between for different routines. The default routine is the Verlet leapfrog (leapfrog) integration scheme. Furthermore, the Velocity Verlet (velverlet) and Brunger-Brooks-Karplus (bbk) integration schemes can be chosen. The last possibility for the integration method is minimize, which starts the minimize routine rather than a dynamics routine. The bbkseed is only useful in conjunction with the BBK-integration method. It allows the rerun the simulation with the same conditions. However, when running in parallel this might not occur (due to the parallel decomposition). The default value of zero means that the current computer time will be used as a seed.
Temperature and pressure
T 300.0 # Target temperature (in Kelvin) lambdac 0.0 # Berendsen temperature coupling constant P 1.0 # Target pressure (in bar) muc 0.0 # Berendsen pressure coupling constant anisopress off # Anisotropic pressure scaling Prep 1000 # Pressure scaling frequency reportheat off # Report the heat flow to and from the reservoir
Both the temperature and pressure values control the target temperature or pressure in the Berendsen loose coupling technique. With the BBK-integration method only temperature coupling can be acchieved (through the Langevin damping). Besides these targets also the coupling constants lambdac for the temperature coupling and muc for the pressure coupling in the case of the Berendsen loose coupling technique need to be set. If both lambdac and muc are set to zero, the system is a microcanonical (NVE) ensemble. Specifying only lambdac moves the system towards the canonical (NVT) ensemble, while specifying muc as well moves the system into the isothermal-isobaric (NPT) ensemble. The option anisopress specifies whether the orthorombic box is scaled uniform in all three directions (off) or each direction is scaled independently (on). Furthermore, one can choose to prevent pressure scaling in specific directions by choosing the appropriate value for the anisopress option (either fixx, fixy, fixz, fixxy, fixxz, fixyz or fixxyz). Finally, the option Prep specifies the frequency for scaling the box. Decreasing this value might force you to decrease the value of muc as well, otherwise the box can be scaled to fast, resulting in a crash of the simulation. The reportheat option allows, in the case of Berendsen loose coupling for temperature, to report the heat flow to or from the heat bath reservoir to the system. This heat flow is both reported on screen as well as in the output file.
Simulated Annealing
annealing off # Turning annealing on or off sepheatbath off # Turning separate heat baths on or off Thold 0.0 # Final temperature (in Kelvin) Trate 0.0 # Change in temperature per scaling (in Kelvin) Annrep 1000 # Temperature scaling every Annrep'th iteration
Five options control simulated annealing. First of all two switch options determine whether to use annealing (annealing), and if so, to split the annealing across molecule types or scale the system as one (sepheatbath). In the latter case, the annealing constants for each molecule need to be expressed in a separate file: specials.pif, which can be controlled by the input option specfile. The starting temperature of simulated annealing is specified by the normal temperature option T, see above. The final temperature on the contrary is specified by the Thold option. Although it depends on the combination of Trate and Annrep if the simulation reaches Thold before the simulation finishes. The option Trate specifies the amount of temperature change (either positive or negative) to be applied every Annrep'th iteration. A note of caution is appropriate, since the code does not check if the annealing is moving away from the final temperature instead towards it.
Energy report and file dumping
Erep 1 # Energy report frequency
Epp on # Report total or average per particle energies
Crep 10000 # Frequency of writing configuration to file
verbose off # Do or do not report additional information on screen
silent off # Do or do not report energies on screen
Dumpsize 100000 # How many bytes are broadcasted at once while
# dumping in a parallel run
Output of the simulation is necessary for analysis. Six options in total control the basic output of PumMa. The option Erep specifies the frequency of writing the systems energy information to file and screen. With Epp set to off the total energy of the system is reported, while setting it to on averages the energy per particle (or bond, or angle). A complete configuration can be written to file every Crep iterations, allowing trajectory analysis. The option Dumpsize is important when using in parallel, since it specifies the size of the communication buffer when dumping a configuration. With the two switch options verbose and silent one can control the amount of information that is displayed on screen.
File names
sysfile sys # Base name of configuration files outputcomment # # Comment sign to be used in output files parfile parameters.ppf # Name of parameter file specfile specials.pif # Name of files containing special options outpfile output.pof # Name of energy output file nbtablefile nbtables.ppf # Name of file containing a tabulated potential nbtableoutfile nbtables.pof # Name of file containing all nonbonded potentials steerfile steerwork.pof # Name of file containing work done by steering
To be able to read or to write files, their filenames must be known. The option sysfile is the only option not to specify a complete file name, but instead the base name of the configuration files, for instance sysIT1200.cfg. The parameter file is given by the option parfile, whereas the file containing special requests during the simulation is given by the option specfile. All system energies are written to the output file specified by the option outpfile. The nbtablefile specifies a file from which PumMa can read potential information (not that this file has a ppf extension, which means it is a parameter file), whereas the nbtableoutfile is a dump of all tabulated potentials used by PumMa. The steerfile is used by PumMa, when some atoms are being steered, to dump information (per iteration) on the work done on the system by the steering.
Nonbonded parameters
Rc 2.5 # Van der Waals cutoff radius factor Rl 0.9 # Factor specifying size of additional radius for pairlist RcCB 0.0 # Cutoff radius for Coulomb interactions (in nanometer) fixedrc -1 # Cutoff radius for all particle types (in nanometer) fixedrl -1 # Pairlist radius for all particle types (in nanometer) noshifts off # Disable shift function for nonbonded interactions vlrep 0 # Frequency of building pairlist again EpsR 1.0 # Dielectric constant Excl13 0.0 # Exclude Van der Waals 1-3 interactions Excl14 0.0 # Exclude Van der Waals 1-4 interactions NBtables on # Use tabulated potentials NBprec 7 # Precision of tabulated potentials NBtol 0.1 # Tolerance of tabulated potentials TabPotRl 0.4 # Pairlist radius for file specified potentials
Cut offs are very important in the behavior of the nonbonded interactions. Quite a few options exist in PumMa to control this behavior. The option Rc gives a cut off radius factor for Van der Waals interactions. It is important to notice that it is a factor, and not a radius. Therefore, it has no dimensions. By default PumMa determines a specific cut off radius for every nonbonded interaction pair, based on the cut off factor Rc. For instance, if the collision diameter σij for the atoms i and j is exactly 1 nm, their Van der Waals cut off radius is 2.5 nm. Beyond that point the interaction is assumed to be non existing. To avoid discontinuities, the interaction function is shifted towards zero by the amount of potential energy that is found for the cut off radius, in this example at 2.5 nm. The option Rl is a factor too, which indicates the size of the pairlist (Verlet list) additional to the cut off radius. In the example used before, the additional pairlist radius is thus 0.9 nm, since the collision diameter equals one nm. The total size of the pairlist, therefore, is 3.4 nm (2.5+0.9=3.4). For the electrostatic interactions the same cut off radius as for the Van der Waals interactions is used (determined by the factors Rc and Rl), unless the option RcCB is specified otherwise. With RcCB one specifies the exact cut off radius (also in nanometers), and the pairlist size for the electrostatic interaction is set using the Rl factor specified before. The use of fixed cut off radii is also possible. The option fixedrc gives the cut off radius (in nm), and the option fixedrl the pairlist radius (also in nm). If these two options are negative they are ignored and the values specified by Rc, Rl and RcCB are used. If, for any reason, one wants to avoid shifting of the potential energies at their cut off radius, the option noshifts can be set to on.
By default PumMa determines by itself when the pairlists needs to be updated (based on the motion of particles, in such a way that a particle not present in the pairlist cannot enter the region below the cut off radius). However, one can force the pairlist to update at a certain frequency by setting the option vlrep. A zero value means that PumMa uses its default behavior. For the electrostatic interactions a dielectric constant can be specified by the option EpsR. Sometimes it is desirable to exclude nonbonded interactions across angle or dihedral terms. This can be acchieved by setting the options Excl13 (for angles) and Excl14 (for dihedrals), respectively. Nonbonded interactions for two bonded atoms are excluded by default. The default value of Excl13 and Excl14 is 0.0, which means that nonbonded interactions are calculated, whereas setting them to 1.0 would exclude them completely. Any other value within this range is possible, scaling the nonbonded interactions with this specific factor.
To allow different functions for the potential energy, PumMa uses tabulated potentials. The option NBtables allows to switch from tabulated (on) to functions (off). However, in the latter case only the Lennard-Jones 12-6 potential can be used. By default PumMa uses a precision for the tabulated potential given by the option NBprec and a tolerance given by NBtol. These two option together specify the number of points in the tabulated potential and their spacing. The default precision of 7 is the lower limit with respect to accuracy. Since PumMa also allows tabulated potentials to be read from file and stored in the programs memory to be used for a certain pair interactions, the option TabPotRl allows to manually set the pairlist distance (in nanometers).
Miscellaneous
fixed off # Switch fixed atoms on or off steeredMD off # Switch steering atoms on or off lmrep 1000 # Frequency of correcting linear momentum vlsafety 1.0 # Safety value for memory allocation Verlet list vcsafety 1.0 # Safety value for memory allocation Coulomb list
The first two of these five miscellaneous options let you specify whether some atoms are fixed in the system or if some atoms are being steered (they need to specified separately in the specfile). The third option lmrep allows you to give the frequency of correcting the linear momentum, so the total system has no linear momentum. The two safety options are only important if one wants to study a system which has no homogeneous density (such as a protein in vacuum or a liquid droplet embedded in a gas). They control the amount of memory that should be allocated for the Verlet and Coulomb lists, and default to one. Increasing this value (for instance to 5.0) would satisfy a simulation of a protein in vacuum.
Cell and processor partitioning
NCxp -1 # Number of cells in x direction per processor NCyp -1 # Number of cells in y direction per processor NCz -1 # Number of cells in z direction Stretch 3 # Multiply number of cells by this factor
These option rarely need to be changed by the user, since the code determines their values automatically. However, to allow user freedom, they can be set through the input file as well. Each of the options NCxp, NCyp and NCz gives the number of cells in the corresponding orthogonal direction, where it has to be mentioned that the x- and y-direction are per processor and the z-direction is identical for all processors (due to the way parallelism is integrated in PumMa). The cells are used when building the pairlist, so only neigboring cells need to be searched, instead of the entire system. By default the number of cells depends on the maximum size of any of the pairlists, i.e. by fitting as many cells into the box length as possible for the maximum pairlist size. The option Stretch allows the user to force a cell division with a specific factor given by the options value. A value of 1 would decrease the amount of cells 27 times with respect to the default option 3. In large systems with many cells, a Stretch of 3 can induce some speedup, since the shape of the surrounding cells becomes more spherical with higher Stretch values. This is one way of improving the performance of the code, although for biological systems an optimal Stretch value of 3 is most often encountered. For higher values, the looping over all cells becomes more time consuming than the gain in a more spherical cell surrounding.
Special options for PumMa
Lennard-Jones contribution analyis
ljoutpfile ljoutp.pof # Name of file containing splitted energies per
# Lennard-Jones interaction
When the Lennard-Jones interaction energies are collected pair-wise, the information it written to the file specified by the ljoutpfile option. Note that the use of this option requires a different compiled PumMa code.
Macroscopic quantities
MQstart 0 # Start iteration of MQ collection MQstep 1 # Information collection frequency MQstop -1 # Final iteration of MQ collection MQcells 25 # Number of cells MQfreq -1 # Frequency of MQ dumping MQcoor x # Coordinate for MQ collection MQfile mq.pof # Name of file for MQ report
When the special macroscopic quantities (MQ) are compiled into the PumMa code, there are six options which control the behavior of the collection of macroscopic quantities. The option MQstart gives the iteration at which the MQ collection should start. The frequency at which the information of specific iterations is used in the MQ collection is given by the option MQstep. The default value of 1, means that every configuration is used in the MQ collection. The iteration at which the collection is stopped is given by the option MQstop, where a value of -1 means that the collection is continued until the simulation stops. The number of cells, which indicates the resolution, along the coordinate of interest is given by the option MQcells. Normally, the MQ information is only written to file at the end of the simulation. However, by specifying the MQfreq option one can force the program to report the MQ information every so often. The coordinate along which the MQ collection is performed is given by MQcoor, which can either be x, y or z. Finally the file name to which all the MQ information is reported is given by the option MQfile.
<< Reduced units | Manual | Parameters >>