Previous Section - LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Commands - Next Section

### 4. How-to discussions

The following sections describe what commands can be used to perform certain kinds of LAMMPS simulations.

4.1 Restarting a simulation
4.2 2d simulations
4.3 CHARMM and AMBER force fields
4.4 Running multiple simulations from one input script
4.5 Parallel tempering
4.6 Granular models
4.7 TIP3P water model
4.8 TIP4P water model
4.9 SPC water model
4.10 Coupling LAMMPS to other codes
4.11 Visualizing LAMMPS snapshots
4.12 Non-orthogonal simulation boxes
4.13 NEMD simulations
4.14 Aspherical particles
4.15 Output from LAMMPS

The example input scripts included in the LAMMPS distribution and highlighted in this section also show how to setup and run various kinds of problems.

#### 4.1 Restarting a simulation

There are 3 ways to continue a long LAMMPS simulation. Multiple run commands can be used in the same input script. Each run will continue from where the previous run left off. Or binary restart files can be saved to disk using the restart command. At a later time, these binary files can be read via a read_restart command in a new script. Or they can be converted to text data files and read by a read_data command in a new script. This section discusses the restart2data tool that is used to perform the conversion.

Here we give examples of 2 scripts that read either a binary restart file or a converted data file and then issue a new run command to continue where the previous run left off. They illustrate what settings must be made in the new script. Details are discussed in the documentation for the read_restart and read_data commands.

Look at the in.chain input script provided in the bench directory of the LAMMPS distribution to see the original script that these 2 scripts are based on. If that script had the line

    restart	        50 tmp.restart

added to it, it would produce 2 binary restart files (tmp.restart.50 and tmp.restart.100) as it ran.

This script could be used to read the 1st restart file and re-run the last 50 timesteps:

    read_restart	tmp.restart.50

    neighbor	0.4 bin
    neigh_modify	every 1 delay 1

    fix		1 all nve
    fix		2 all langevin 1.0 1.0 10.0 904297

    timestep	0.012

    run		50

Note that the following commands do not need to be repeated because their settings are included in the restart file: units, atom_style, special_bonds, pair_style, bond_style. However these commands do need to be used, since their settings are not in the restart file: neighbor, fix, timestep.

If you actually use this script to perform a restarted run, you will notice that the thermodynamic data match at step 50 (if you also put a "thermo 50" command in the original script), but do not match at step 100. This is because the fix langevin command uses random numbers in a way that does not allow for perfect restarts.

As an alternate approach, the restart file could be converted to a data file using this tool:

    restart2data tmp.restart.50 tmp.restart.data

Then, this script could be used to re-run the last 50 steps:

    units		lj
    atom_style	bond
    pair_style	lj/cut 1.12
    pair_modify	shift yes
    bond_style	fene
    special_bonds   0.0 1.0 1.0

    read_data	tmp.restart.data

    neighbor	0.4 bin
    neigh_modify	every 1 delay 1

    fix		1 all nve
    fix		2 all langevin 1.0 1.0 10.0 904297

    timestep	0.012

    reset_timestep	50
    run		50

Note that nearly all the settings specified in the original in.chain script must be repeated, except the pair_coeff and bond_coeff commands since the new data file lists the force field coefficients. Also, the reset_timestep command is used to tell LAMMPS the current timestep. This value is stored in restart files, but not in data files.

#### 4.2 2d simulations

Use the dimension command to specify a 2d simulation.

Make the simulation box periodic in z via the boundary command. This is the default.

If using the create box command to define a simulation box, set the z dimensions narrow, but finite, so that the create_atoms command will tile the 3d simulation box with a single z plane of atoms - e.g.

    create box 1 -10 10 -10 10 -0.25 0.25

If using the read data command to read in a file of atom coordinates, set the "zlo zhi" values to be finite but narrow, similar to the create_box command settings just described. For each atom in the file, assign a z coordinate so it falls inside the z-boundaries of the box - e.g. 0.0.

Use the fix enforce2d command as the last defined fix to insure that the z-components of velocities and forces are zeroed out every timestep. The reason to make it the last fix is so that any forces induced by other fixes will be zeroed out.

Many of the example input scripts included in the LAMMPS distribution are for 2d models.

#### 4.3 CHARMM and AMBER force fields

There are many different ways to compute forces in the CHARMM and AMBER molecular dynamics codes, only some of which are available as options in LAMMPS. A force field has 2 parts: the formulas that define it and the coefficients used for a particular system. Here we only discuss formulas implemented in LAMMPS. Setting coefficients is done in the input data file via the read_data command or in the input script with commands like pair_coeff or bond_coeff. See this section for additional tools that can use CHARMM or AMBER to assign force field coefficients and convert their output into LAMMPS input.

See (MacKerell) for a description of the CHARMM force field. See (Cornell) for a description of the AMBER force field.

These style choices compute force field formulas that are consistent with common options in CHARMM or AMBER. See each command's documentation for the formula it computes.

- bond_style harmonic
- angle_style charmm
- dihedral_style charmm
- pair_style lj/charmm/coul/charmm
- pair_style lj/charmm/coul/charmm/implicit
- pair_style lj/charmm/coul/long

- special_bonds charmm
- special_bonds amber

#### 4.4 Running multiple simulations from one input script

This can be done in several ways. See the documentation for individual commands for more details on how these examples work.

If "multiple simulations" means continue a previous simulation for more timesteps, then you simply use the run command multiple times. For example, this script

    units lj
    atom_style atomic
    read_data data.lj
    run 10000
    run 10000
    run 10000
    run 10000
    run 10000

would run 5 successive simulations of the same system for a total of 50,000 timesteps.

If you wish to run totally different simulations, one after the other, the clear command can be used in between them to re-initialize LAMMPS. For example, this script

    units lj
    atom_style atomic
    read_data data.lj
    run 10000
    clear
    units lj
    atom_style atomic
    read_data data.lj.new
    run 10000

would run 2 independent simulations, one after the other.

For large numbers of independent simulations, you can use variables and the next and jump commands to loop over the same input script multiple times with different settings. For example, this script, named in.polymer

    variable d index run1 run2 run3 run4 run5 run6 run7 run8
    cd $d
    read_data data.polymer
    run 10000
    cd ..
    clear
    next d
    jump in.polymer

would run 8 simulations in different directories, using a data.polymer file in each directory. The same concept could be used to run the same system at 8 different temperatures, using a temperature variable and storing the output in different log and dump files, for example

    variable a loop 8
    variable t index 0.8 0.85 0.9 0.95 1.0 1.05 1.1 1.15
    log log.$a
    read data.polymer
    velocity all create $t 352839
    fix 1 all nvt $t $t 100.0
    dump 1 all atom 1000 dump.$a
    run 100000
    next t
    next a
    jump in.polymer

All of the above examples work whether you are running on 1 or multiple processors, but assumed you are running LAMMPS on a single partition of processors. LAMMPS can be run on multiple partitions via the "-partition" command-line switch as described in this section of the manual.

In the last 2 examples, if LAMMPS were run on 3 partitions, the same scripts could be used if the "index" and "loop" variables were replaced with universe-style variables, as described in the variable command. Also, the "next t" and "next a" commands would need to be replaced with a single "next a t" command. With these modifications, the 8 simulations of each script would run on the 3 partitions one after the other until all were finished. Initially, 3 simulations would be started simultaneously, one on each partition. When one finished, that partition would then start the 4th simulation, and so forth, until all 8 were completed.

#### 4.5 Parallel tempering

The temper command can be used to perform a parallel tempering or replica-exchange simulation where multiple copies of a simulation are run at different temperatures on different sets of processors, and Monte Carlo temperature swaps are performed between pairs of copies.

Use the -procs and -in command-line switches to launch LAMMPS on multiple partitions.

In your input script, define a set of temperatures, one for each processor partition, using the variable command:

    variable t proc 300.0 310.0 320.0 330.0

Define a fix of style nvt or langevin to control the temperature of each simulation:

    fix myfix all nvt $t $t 100.0

Use the temper command in place of a run command to perform a simulation where tempering exchanges will take place:

    temper 100000 100 $t myfix 3847 58382

#### 4.6 Granular models

To run a simulation of a granular model, you will want to use the following commands:

- atom_style granular
- fix nve/gran
- fix gravity
- thermo_style gran

Use one of these 3 pair potentials:

- pair_style gran/history
- pair_style gran/no_history
- pair_style gran/hertzian

These commands implement fix options specific to granular systems:

- fix freeze
- fix gran/diag
- fix pour
- fix viscous
- fix wall/gran

The fix style freeze zeroes both the force and torque of frozen atoms, and should be used for granular system instead of the fix style setforce.

For computational efficiency, you can eliminate needless pairwise computations between frozen atoms by using this command:

- neigh_modify exclude

#### 4.7 TIP3P water model

The TIP3P water model as implemented in CHARMM (MacKerell) specifies a 3-site rigid water molecule with charges and Lennard-Jones parameters assigned to each of the 3 atoms. In LAMMPS the fix shake command can be used to hold the two O-H bonds and the H-O-H angle rigid. A bond style of harmonic and an angle style of harmonic or charmm should also be used.

These are the additional parameters (in real units) to set for O and H atoms and the water molecule to run a rigid TIP3P-CHARMM model with a cutoff. The K values can be used if a flexible TIP3P model (without fix shake) is desired. If the LJ epsilon and sigma for HH and OH are set to 0.0, it corresponds to the original 1983 TIP3P model (Jorgensen).

O mass = 15.9994
H mass = 1.008

O charge = -0.834
H charge = 0.417

LJ epsilon of OO = 0.1521
LJ sigma of OO = 3.188
LJ epsilon of HH = 0.0460
LJ sigma of HH = 0.4000
LJ epsilon of OH = 0.0836
LJ sigma of OH = 1.7753

K of OH bond = 450
r0 of OH bond = 0.9572

K of HOH angle = 55
theta of HOH angle = 104.52

These are the parameters to use for TIP3P with a long-range Coulombic solver (Ewald or PPPM in LAMMPS):

O mass = 15.9994
H mass = 1.008

O charge = -0.830
H charge = 0.415

LJ epsilon of OO = 0.102
LJ sigma of OO = 3.1507
LJ epsilon, sigma of OH, HH = 0.0

K of OH bond = 450
r0 of OH bond = 0.9572

K of HOH angle = 55
theta of HOH angle = 104.52

#### 4.8 TIP4P water model

The four-point TIP4P rigid water model extends the traditional three-point TIP3P model by adding an additional site, usually massless, where the charge associated with the oxygen atom is placed. This site M is located at a fixed distance away from the oxygen along the bisector of the HOH bond angle. A bond style of harmonic and an angle style of harmonic or charmm should also be used.

Currently, only a four-point model for long-range Coulombics is implemented via the LAMMPS pair style lj/cut/coul/long/tip4p. We plan to add a cutoff version in the future. For both models, the bond lengths and bond angles should be held fixed using the fix shake command.

These are the additional parameters (in real units) to set for O and H atoms and the water molecule to run a rigid TIP4P model with a cutoff (Jorgensen). Note that the OM distance is specified in the pair_style command, not as part of the pair coefficients.

O mass = 15.9994
H mass = 1.008

O charge = -1.040
H charge = 0.520

r0 of OH bond = 0.9572
theta of HOH angle = 104.52

OM distance = 0.15

LJ epsilon of O-O = 0.1550
LJ sigma of O-O = 3.1536
LJ epsilon, sigma of OH, HH = 0.0

These are the parameters to use for TIP4P with a long-range Coulombic solver (Ewald or PPPM in LAMMPS):

O mass = 15.9994
H mass = 1.008

O charge = -1.0484
H charge = 0.5242

r0 of OH bond = 0.9572
theta of HOH angle = 104.52

OM distance = 0.1250

LJ epsilon of O-O = 0.16275
LJ sigma of O-O = 3.16435
LJ epsilon, sigma of OH, HH = 0.0

#### 4.9 SPC water model

The SPC water model specifies a 3-site rigid water molecule with charges and Lennard-Jones parameters assigned to each of the 3 atoms. In LAMMPS the fix shake command can be used to hold the two O-H bonds and the H-O-H angle rigid. A bond style of harmonic and an angle style of harmonic or charmm should also be used.

These are the additional parameters (in real units) to set for O and H atoms and the water molecule to run a rigid SPC model with long-range Coulombics (Ewald or PPPM in LAMMPS).

O mass = 15.9994
H mass = 1.008

O charge = -0.820
H charge = 0.410

LJ epsilon of OO = 0.1553
LJ sigma of OO = 3.166
LJ epsilon, sigma of OH, HH = 0.0

r0 of OH bond = 1.0
theta of HOH angle = 109.47

#### 4.10 Coupling LAMMPS to other codes

LAMMPS is designed to allow it to be coupled to other codes. For example, a quantum mechanics code might compute forces on a subset of atoms and pass those forces to LAMMPS. Or a continuum finite element (FE) simulation might use atom positions as boundary conditions on FE nodal points, compute a FE solution, and return interpolated forces on MD atoms.

LAMMPS can be coupled to other codes in at least 3 ways. Each has advantages and disadvantages, which you'll have to think about in the context of your application.

(1) Define a new fix command that calls the other code. In this scenario, LAMMPS is the driver code. During its timestepping, the fix is invoked, and can make library calls to the other code, which has been linked to LAMMPS as a library. This is the way the POEMS package that performs constrained rigid-body motion on groups of atoms is hooked to LAMMPS. See the fix_poems command for more details. See this section of the documentation for info on how to add a new fix to LAMMPS.

(2) Define a new LAMMPS command that calls the other code. This is conceptually similar to method (1), but in this case LAMMPS and the other code are on a more equal footing. Note that now the other code is not called during the timestepping of a LAMMPS run, but between runs. The LAMMPS input script can be used to alternate LAMMPS runs with calls to the other code, invoked via the new command. The run command facilitates this with its every option, which makes it easy to run a few steps, invoke the command, run a few steps, invoke the command, etc.

In this scenario, the other code can be called as a library, as in (1), or it could be a stand-alone code, invoked by a system() call made by the command (assuming your parallel machine allows one or more processors to start up another program). In the latter case the stand-alone code could communicate with LAMMPS thru files that the command writes and reads.

See this section of the documentation for how to add a new command to LAMMPS.

(3) Use LAMMPS as a library called by another code. In this case the other code is the driver and calls LAMMPS as needed. Or a wrapper code could link and call both LAMMPS and another code as libraries. Again, the run command has options that allow it to be invoked with minimal overhead (no setup or clean-up) if you wish to do multiple short runs, driven by another program.

This section of the documentation describes how to build LAMMPS as a library. Once this is done, you can interface with LAMMPS either via C++, C, or Fortran (or any other language that supports a vanilla C-like interface, e.g. a scripting language). For example, from C++ you could create one (or more) "instances" of LAMMPS, pass it an input script to process, or execute individual commands, all by invoking the correct class methods in LAMMPS. From C or Fortran you can make function calls to do the same things. Library.cpp and library.h contain such a C interface with the functions:

    void lammps_open(int, char **, MPI_Comm, void **);
    void lammps_close(void *);
    void lammps_file(void *, char *);
    char *lammps_command(doivd *, char *);

The functions contain C++ code you could write in a C++ application that was invoking LAMMPS directly. Note that LAMMPS classes are defined wihin a LAMMPS namespace (LAMMPS_NS) if you use them from another C++ application.

Two of the routines in library.cpp are of particular note. The lammps_open() function initiates LAMMPS and takes an MPI communicator as an argument. It returns a pointer to a LAMMPS "object". As with C++, the lammps_open() function can be called mutliple times, to create multiple instances of LAMMPS.

LAMMPS will run on the set of processors in the communicator. This means the calling code can run LAMMPS on all or a subset of processors. For example, a wrapper script might decide to alternate between LAMMPS and another code, allowing them both to run on all the processors. Or it might allocate half the processors to LAMMPS and half to the other code and run both codes simultaneously before syncing them up periodically.

Library.cpp contains a lammps_command() function to which the caller passes a single LAMMPS command (a string). Thus the calling code can read or generate a series of LAMMPS commands (e.g. an input script) one line at a time and pass it thru the library interface to setup a problem and then run it.

A few other sample functions are included in library.cpp, but the key idea is that you can write any functions you wish to define an interface for how your code talks to LAMMPS and add them to library.cpp and library.h. The routines you add can access any LAMMPS data. The examples/couple directory has example C++ and C codes which show how a stand-alone code can link LAMMPS as a library, run LAMMPS on a subset of processors, grab data from LAMMPS, change it, and put it back into LAMMPS.

#### 4.11 Visualizing LAMMPS snapshots

LAMMPS itself does not do visualization, but snapshots from LAMMPS simulations can be visualized (and analyzed) in a variety of ways.

LAMMPS snapshots are created by the dump command which can create files in several formats. The native LAMMPS dump format is a text file (see "dump atom" or "dump custom") which can be visualized by the xmovie program, included with the LAMMPS package. This produces simple, fast 2d projections of 3d systems, and can be useful for rapid debugging of simulation geoemtry and atom trajectories.

Several programs included with LAMMPS as auxiliary tools can convert native LAMMPS dump files to other formats. See the Section_tools doc page for details. The first is the ch2lmp tool, which contains a lammps2pdb Perl script which converts LAMMPS dump files into PDB files. The second is the lmp2arc tool which converts LAMMPS dump files into Accelrys's Insight MD program files. The third is the lmp2cfg tool which converts LAMMPS dump files into CFG files which can be read into the AtomEye visualizer.

A Python-based toolkit distributed by our group can read native LAMMPS dump files, including custom dump files with additional columns of user-specified atom information, and convert them to various formats or pipe them into visualization software directly. See the Pizza.py WWW site for details. Specifically, Pizza.py can convert LAMMPS dump files into PDB, XYZ, Ensight, and VTK formats. Pizza.py can pipe LAMMPS dump files directly into the Raster3d and RasMol visualization programs. Pizza.py has tools that do interactive 3d OpenGL visualization and one that creates SVG images of dump file snapshots.

LAMMPS can create XYZ files directly (via "dump xyz") which is a simple text-based file format used by many visualization programs including VMD.

LAMMPS can create DCD files directly (via "dump dcd") which can be read by VMD in conjunction with a CHARMM PSF file. Using this form of output avoids the need to convert LAMMPS snapshots to PDB files. See the dump command for more information on DCD files.

LAMMPS can create XTC files directly (via "dump xtc") which is GROMACS file format which can also be read by VMD for visualization. See the dump command for more information on XTC files.

#### 4.12 Non-orthogonal simulation boxes

By default, LAMMPS uses an orthogonal simulation box to encompass the particles. The boundary command sets the boundary conditions of the box (periodic, non-periodic, etc). If the box size is xprd by yprd by zprd then the 3 mutually orthogonal edge vectors of an orthogonal simulation box are a = (xprd,0,0), b = (0,yprd,0), and c = (0,0,zprd).

LAMMPS also allows non-orthogonal simulation boxes (triclinic symmetry) to be defined with 3 additional "tilt" parameters which change the edge vectors of the simulation box to be a = (xprd,0,0), b = (xy,yprd,0), and c = (xz,yz,zprd). The xy, xz, and yz parameters can be positive or negative. The simulation box must be periodic in both dimensions associated with a tilt factor. For example, if xz != 0.0, then the x and z dimensions must be periodic.

To avoid extremely tilted boxes (which would be computationally inefficient), no tilt factor can skew the box more than half the distance of the parallel box length, which is the 1st dimension in the tilt factor (x for xz). For example, if xlo = 2 and xhi = 12, then the x box length is 10 and the xy tilt factor must be between -5 and 5. Similarly, both xz and yz must be between -(xhi-xlo)/2 and +(yhi-ylo)/2. Note that this is not a limitation, since if the maximum tilt factor is 5 (as in this example), then configurations with tilt = ..., -15, -5, 5, 15, 25, ... are all equivalent.

You tell LAMMPS to use a non-orthogonal box when the simulation box is defined. This happens in one of 3 ways. If the create_box command is used with a region of style prism, then a non-orthogonal domain is setup. See the region command for details. If the read_data command is used to define the simulation box, and the header of the data file contains a line with the "xy xz yz" keyword, then a non-orthogonal domain is setup. See the read_data command for details. Finally, if the read_restart command reads a restart file which was written from a simulation using a triclinic box, then a non-orthogonal box will be enabled for the restarted simulation.

Note that you can define a non-orthogonal box with all 3 tilt factors = 0.0, so that it is initially orthogonal. This is necessary if the box will become non-orthogonal. Alternatively, you can use the change_box command to convert a simulation box from orthogonal to non-orthogonal and vice versa.

One use of non-orthogonal boxes is to model solid-state crystals with triclinic symmetry. The lattice command can be used with non-orthogonal basis vectors to define a lattice that will tile a non-orthogonal simulation box via the create_atoms command. Note that while the box edge vectors a,b,c cannot be arbitrary vectors (e.g. a must be aligned with the x axis), it is possible to rotate any crystal's basis vectors so that they meet these restrictions.

A second use of non-orthogonal boxes is to shear a bulk solid to study the response of the material. The fix deform command can be used for this purpose. It allows dynamic control of the xy, xz, and yz tilt factors as a simulation runs.

Another use of non-orthogonal boxes is to perform non-equilibrium MD (NEMD) simulations, as discussed in the next section.

#### 4.13 NEMD simulations

Non-equilibrium molecular dynamics or NEMD simulations are typically used to measure a fluid's rheological properties such as viscosity. In LAMMPS, such simulations can be performed by first setting up a non-orthogonal simulation box (see the preceeding Howto section).

A shear strain can be applied to the simuaation box at a desired strain rate by using the fix deform command. The fix nvt/sllod command can be used to thermostat the sheared fluid and integrate the SLLOD equations of motion for the system. Fix nvt/sllod uses compute temp/deform to compute a thermal temperature by subtracting out the streaming velocity of the shearing atoms. The velocity profile or other properties of the fluid can be monitored via the fix ave/spatial command.

As discussed in the previous section on non-orthogonal simulation boxes, the amount of tilt or skew that can be applied is limited by LAMMPS for computational efficiency to be 1/2 of the parallel box length. However, fix deform can continuously strain a box by an arbitrary amount. As discussed in the fix deform command, when the tilt value reaches a limit, the box is re-shaped to the opposite limit which is an equivalent tiling of periodic space. The strain rate can then continue to change as before. In a long NEMD simulation these box re-shaping events may occur many times.

In a NEMD simulation, the "remap" option of fix deform should be set to "remap v", since that is what fix nvt/sllod assumes to generate a velocity profile consistent with the applied shear strain rate.

An alternative method for calculating viscosities is provided via the fix viscosity command.

#### 4.14 Aspherical particles

LAMMPS supports ellipsoidal particles via the atom_style ellipsoid and shape commands. The latter command defines the 3 axes (diameters) of a general ellipsoid. The pair_style gayberne command can be used to define a Gay-Berne (GB) potential for how ellipsoidal particles interact with each other and with spherical particles. The GB potential is like a Lennard-Jones (LJ) potential, generalized for orientiation-dependent interactions.

The orientation of ellipsoidal particles is stored as a quaternion. See the set command for a brief explanation of quaternions and how the orientation of such particles can be initialized. The data file read by the read_data command contains quaternions for each atom in the Atoms section if atom_style ellipsoid is being used. The compute temp/asphere command can be used to calculate the temperature of a group of ellipsoidal particles, taking account of rotational degrees of freedom. The motion of the particles can be integrated via the fix nve/asphere, fix nvt/asphere, or fix npt/asphere commands. All of these commands are part of the ASPHERE package in LAMMPS.

Computationally, the cost for two ellipsoidal particles to interact is 30 times (or more) expensive than for 2 spherical LJ particles. Thus if you are modeling a system with many spherical particles (e.g. as the solvent), then you should insure sphere-sphere interactions are computed with a cheaper potential than GB. This can be done by setting the particle's 3 shape parameters to all be equal (a sphere). Additionally, the corresponding GB potential coefficients can be set so the GB potential will treat the pair of particles as LJ spheres. Details are given in the doc page for the pair_style gayberne. Alternatively, the pair_style hybrid potential can be used, with the sphere-sphere interactions computed by another pair potential, such as pair_style lj/cut.

#### 4.15 Output from LAMMPS

Aside from restart files, there are two basic kinds of LAMMPS output. The first is thermodynamic output, which is a list of quantities printed every few timesteps to the screen and logfile. The second is dump files, which contain snapshots of atoms and various per-atom values and are written at a specified frequency. A simulation prints one set of thermodynamic output; it may generate zero, or one, or multiple dump files. LAMMPS gives you a variety of ways to determine what quantities are computed and printed when thermodynamic info or dump files are output. There are also three fixes which can do their own output of user-defined quantities: fix ave/time for time averaging, fix ave/spatial for spatial averaging, and fix print. These are described below.

The frequency and format of thermodynamic output is set by the thermo, thermo_style, and thermo_modify commands. The thermo_style command also specifies what values are calculated and written out. Pre-defined keywords can be specified (e.g. press, etotal, etc) which include time-averaged versions of temperature, pressure, and a few other variables (tave, pave, etc). Three addtional kinds of keywords can also be specified (c_ID, f_ID, v_name), where a compute or fix or variable provides the value(s) to be output. Each of these are described in turn.

In LAMMPS, a compute comes in two flavors: ones that compute one or more global values (e.g. temperature, kinetic energy tensor) and ones that compute one or more per-atom values. There is a compute sum command which sums per-atom quantities into a global scalar or vector.

Only global quantities from a compute can be used for thermodynamic output. The user-defined ID of the compute is used along with an optional subscript as part of the thermo_style command. E.g. c_myTemp outputs the single scalar value generated by the compute; c_myTemp[2] would output the 2nd vector value.

Fixes can also generate global scalar or vector values which can be output with thermodynamic output, e.g. the energy of an indenter's interaction with the simulation atoms. These values are accessed via the same format as a compute's values, as f_ID or f_ID[N]. See the doc pages for individual fix commands to see which ones generate global values that can be output with thermodynamic info. The fix ave/time command generates time-averaged global quantities which can be accessed for thermodynamic output.

Input script variables of various kinds are defined by the variable command. All kinds except the atom-style variable can be used for thermodynamic output. A variable with name "abc" is referenced in a thermo_style command as v_abc.

The variable formula defined in the input script can contain math functions (add, exp, etc), atom values (x[N], fx[N]), groups quantities (mass(), vcm(), etc), references to thermodynamic quantities (e.g. temp, volume, etc), or references to other variables or computes or fixes. Thus a variable is the most general way to define some quantity you want calculated and output with thermodynamic info.

Dump file output is specified by the dump and dump_modify commands. There are several pre-defined formats (dump atom, dump xtc, etc). There is also a dump custom format where you specify what values are output with each atom. Pre-defined keywords can be specified (e.g. tag, type, x, etc). Two additional kinds of keywords can also be specified (c_ID, f_ID), where a compute or fix provides the values to be output.

Computes that generate per-atom values can be accessed by the dump custom command. These are computes that have the word "atom" in their style name, e.g. ke/atom, stress/atom, etc. The values are accessed as c_myKE for a scalar per-atom quantity or as c_myStress[2] for a component of a vector per-atom quantity. The compute variable/atom command takes a user-defined atom-style variable as input and calculates its value for each atom. Since this compute can be accessed by the dump custom command, this is a general way to define some quantity you want calculated and output in a dump file.

Fixes can also generate per-atom values to output to dump files. For example, the fix ave/atom command calculates time-averages of compute quantities. As indicated in the preceeding paragraph, a compute quantity can be a calculated value such as energy or stress or it can be a value calculated by an atom-style variable, or it can be an atom attribute such as velocity or force. These per-atom fix values are accessed by the dump custom command as f_myKE for a scalar per-atom quantity or as f_myStress[2] for a component of a vector per-atom quantity.

Three other fixes are of particular note for output: fix ave/time, fix ave/spatial, and fix print.

The fix ave/time command enables time-averaging of global quantities like temperature or pressure. The global quantities are calculated by a compute or a fix. The compute or fix must generate global scalar or vector quantities. Note that this includes the "compute sum" command which computes global values by summing per-atom quantities. The time-averaged values generated by fix ave/time can be written directly to a file and/or accessed by any output command that uses fixes as a source of input, e.g. the thermo_style custom command. Fix ave/time options allow for running cummulative averages or moving time-windowed averages to be output.

The fix ave/spatial command enables spatial-averaging of per-atom quantities like per-atom energy or stress. The per-atom quantities can be atom density (mass or number) or be calculated by a compute or a fix. The compute or fix must generate per-atom scalar or vector quantities. Note that if you use the fix ave/atom command with fix ave/spatial, it means you are effectively calculating a time average of a spatial average of a time-averaged per-atom quantity. The time-averaged values generated by fix ave/spatial can be written directly to a file and/or accessed by any output command that uses fixes as a source of input, e.g. the thermo_style custom command. Fix ave/spatial options allow for running cummulative averages or moving time-windowed averages to be output.

The fix print command can generate a line of output written to the screen and log file periodically during a running simulation. Since the line can contain one or more variable quantities, this command is a means to output desired calculated quantities that are not part of thermodynamic or dump file output.

(Cornell) Cornell, Cieplak, Bayly, Gould, Merz, Ferguson, Spellmeyer, Fox, Caldwell, Kollman, JACS 117, 5179-5197 (1995).

(Horn) Horn, Swope, Pitera, Madura, Dick, Hura, and Head-Gordon, J Chem Phys, 120, 9665 (2004).

(MacKerell) MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field, Fischer, Gao, Guo, Ha, et al, J Phys Chem, 102, 3586 (1998).

(Jorgensen) Jorgensen, Chandrasekhar, Madura, Impey, Klein, J Chem Phys, 79, 926 (1983).