diff --git a/doc/Section_howto.html b/doc/Section_howto.html
index bad450038..d7bfb8ced 100644
--- a/doc/Section_howto.html
+++ b/doc/Section_howto.html
@@ -1,1963 +1,1967 @@
<HTML>
<CENTER><A HREF = "Section_accelerate.html">Previous Section</A> - <A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A> - <A HREF = "Section_example.html">Next Section</A>
</CENTER>
<HR>
<H3>6. How-to discussions
</H3>
<P>This section describes how to perform common tasks using LAMMPS.
</P>
6.1 <A HREF = "#howto_1">Restarting a simulation</A><BR>
6.2 <A HREF = "#howto_2">2d simulations</A><BR>
6.3 <A HREF = "#howto_3">CHARMM, AMBER, and DREIDING force fields</A><BR>
6.4 <A HREF = "#howto_4">Running multiple simulations from one input script</A><BR>
6.5 <A HREF = "#howto_5">Multi-replica simulations</A><BR>
6.6 <A HREF = "#howto_6">Granular models</A><BR>
6.7 <A HREF = "#howto_7">TIP3P water model</A><BR>
6.8 <A HREF = "#howto_8">TIP4P water model</A><BR>
6.9 <A HREF = "#howto_9">SPC water model</A><BR>
6.10 <A HREF = "#howto_10">Coupling LAMMPS to other codes</A><BR>
6.11 <A HREF = "#howto_11">Visualizing LAMMPS snapshots</A><BR>
6.12 <A HREF = "#howto_12">Triclinic (non-orthogonal) simulation boxes</A><BR>
6.13 <A HREF = "#howto_13">NEMD simulations</A><BR>
6.14 <A HREF = "#howto_14">Extended spherical and aspherical particles</A><BR>
6.15 <A HREF = "#howto_15">Output from LAMMPS (thermo, dumps, computes, fixes, variables)</A><BR>
6.16 <A HREF = "#howto_16">Thermostatting, barostatting and computing temperature</A><BR>
6.17 <A HREF = "#howto_17">Walls</A><BR>
6.18 <A HREF = "#howto_18">Elastic constants</A><BR>
6.19 <A HREF = "#howto_19">Library interface to LAMMPS</A><BR>
6.20 <A HREF = "#howto_20">Calculating thermal conductivity</A><BR>
6.21 <A HREF = "#howto_21">Calculating viscosity</A> <BR>
<P>The example input scripts included in the LAMMPS distribution and
highlighted in <A HREF = "Section_example.html">Section_example</A> also show how to
setup and run various kinds of simulations.
</P>
<HR>
<HR>
<A NAME = "howto_1"></A><H4>6.1 Restarting a simulation
</H4>
<P>There are 3 ways to continue a long LAMMPS simulation. Multiple
<A HREF = "run.html">run</A> 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 <A HREF = "restart.html">restart</A>
command. At a later time, these binary files can be read via a
<A HREF = "read_restart.html">read_restart</A> command in a new script. Or they can
be converted to text data files and read by a
<A HREF = "read_data.html">read_data</A> command in a new script. <A HREF = "Section_tools.html">This
section</A> discusses the <I>restart2data</I> tool that is
used to perform the conversion.
</P>
<P>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 <A HREF = "read_restart.html">read_restart</A> and
<A HREF = "read_data.html">read_data</A> commands.
</P>
<P>Look at the <I>in.chain</I> input script provided in the <I>bench</I> directory
of the LAMMPS distribution to see the original script that these 2
scripts are based on. If that script had the line
</P>
<PRE>restart 50 tmp.restart
</PRE>
<P>added to it, it would produce 2 binary restart files (tmp.restart.50
and tmp.restart.100) as it ran.
</P>
<P>This script could be used to read the 1st restart file and re-run the
last 50 timesteps:
</P>
<PRE>read_restart tmp.restart.50
</PRE>
<PRE>neighbor 0.4 bin
neigh_modify every 1 delay 1
</PRE>
<PRE>fix 1 all nve
fix 2 all langevin 1.0 1.0 10.0 904297
</PRE>
<PRE>timestep 0.012
</PRE>
<PRE>run 50
</PRE>
<P>Note that the following commands do not need to be repeated because
their settings are included in the restart file: <I>units, atom_style,
special_bonds, pair_style, bond_style</I>. However these commands do
need to be used, since their settings are not in the restart file:
<I>neighbor, fix, timestep</I>.
</P>
<P>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 <A HREF = "fix_langevin.html">fix langevin</A> command
uses random numbers in a way that does not allow for perfect restarts.
</P>
<P>As an alternate approach, the restart file could be converted to a data
file using this tool:
</P>
<PRE>restart2data tmp.restart.50 tmp.restart.data
</PRE>
<P>Then, this script could be used to re-run the last 50 steps:
</P>
<PRE>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
</PRE>
<PRE>read_data tmp.restart.data
</PRE>
<PRE>neighbor 0.4 bin
neigh_modify every 1 delay 1
</PRE>
<PRE>fix 1 all nve
fix 2 all langevin 1.0 1.0 10.0 904297
</PRE>
<PRE>timestep 0.012
</PRE>
<PRE>reset_timestep 50
run 50
</PRE>
<P>Note that nearly all the settings specified in the original <I>in.chain</I>
script must be repeated, except the <I>pair_coeff</I> and <I>bond_coeff</I>
commands since the new data file lists the force field coefficients.
Also, the <A HREF = "reset_timestep.html">reset_timestep</A> command is used to tell
LAMMPS the current timestep. This value is stored in restart files,
but not in data files.
</P>
<HR>
<A NAME = "howto_2"></A><H4>6.2 2d simulations
</H4>
<P>Use the <A HREF = "dimension.html">dimension</A> command to specify a 2d simulation.
</P>
<P>Make the simulation box periodic in z via the <A HREF = "boundary.html">boundary</A>
command. This is the default.
</P>
<P>If using the <A HREF = "create_box.html">create box</A> 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.
</P>
<PRE><A HREF = "create_box.html">create box</A> 1 -10 10 -10 10 -0.25 0.25
</PRE>
<P>If using the <A HREF = "read_data.html">read data</A> 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.
</P>
<P>Use the <A HREF = "fix_enforce2d.html">fix enforce2d</A> 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.
</P>
<P>Many of the example input scripts included in the LAMMPS distribution
are for 2d models.
</P>
<P>IMPORTANT NOTE: Some models in LAMMPS treat particles as extended
spheres, as opposed to point particles. In 2d, the particles will
still be spheres, not disks, meaning their moment of inertia will be
the same as in 3d.
</P>
<HR>
<A NAME = "howto_3"></A><H4>6.3 CHARMM, AMBER, and DREIDING force fields
</H4>
<P>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 that correspond to formulas commonly
used in the CHARMM, AMBER, and DREIDING force fields. Setting
coefficients is done in the input data file via the
<A HREF = "read_data.html">read_data</A> command or in the input script with
commands like <A HREF = "pair_coeff.html">pair_coeff</A> or
<A HREF = "bond_coeff.html">bond_coeff</A>. See <A HREF = "Section_tools.html">Section_tools</A>
for additional tools that can use CHARMM or AMBER to assign force
field coefficients and convert their output into LAMMPS input.
</P>
<P>See <A HREF = "#MacKerell">(MacKerell)</A> for a description of the CHARMM force
field. See <A HREF = "#Cornell">(Cornell)</A> for a description of the AMBER force
field.
</P>
<P>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.
</P>
<UL><LI><A HREF = "bond_harmonic.html">bond_style</A> harmonic
<LI><A HREF = "angle_charmm.html">angle_style</A> charmm
<LI><A HREF = "dihedral_charmm.html">dihedral_style</A> charmm
<LI><A HREF = "pair_charmm.html">pair_style</A> lj/charmm/coul/charmm
<LI><A HREF = "pair_charmm.html">pair_style</A> lj/charmm/coul/charmm/implicit
<LI><A HREF = "pair_charmm.html">pair_style</A> lj/charmm/coul/long
</UL>
<UL><LI><A HREF = "special_bonds.html">special_bonds</A> charmm
<LI><A HREF = "special_bonds.html">special_bonds</A> amber
</UL>
<P>DREIDING is a generic force field developed by the <A HREF = "http://www.wag.caltech.edu">Goddard
group</A> at Caltech and is useful for
predicting structures and dynamics of organic, biological and
main-group inorganic molecules. The philosophy in DREIDING is to use
general force constants and geometry parameters based on simple
hybridization considerations, rather than individual force constants
and geometric parameters that depend on the particular combinations of
atoms involved in the bond, angle, or torsion terms. DREIDING has an
<A HREF = "pair_hbond_dreiding.html">explicit hydrogen bond term</A> to describe
interactions involving a hydrogen atom on very electronegative atoms
(N, O, F).
</P>
<P>See <A HREF = "#Mayo">(Mayo)</A> for a description of the DREIDING force field
</P>
<P>These style choices compute force field formulas that are consistent
with the DREIDING force field. See each command's
documentation for the formula it computes.
</P>
<UL><LI><A HREF = "bond_harmonic.html">bond_style</A> harmonic
<LI><A HREF = "bond_morse.html">bond_style</A> morse
</UL>
<UL><LI><A HREF = "angle_harmonic.html">angle_style</A> harmonic
<LI><A HREF = "angle_cosine.html">angle_style</A> cosine
<LI><A HREF = "angle_cosine_periodic.html">angle_style</A> cosine/periodic
</UL>
<UL><LI><A HREF = "dihedral_charmm.html">dihedral_style</A> charmm
<LI><A HREF = "improper_umbrella.html">improper_style</A> umbrella
</UL>
<UL><LI><A HREF = "pair_buck.html">pair_style</A> buck
<LI><A HREF = "pair_buck.html">pair_style</A> buck/coul/cut
<LI><A HREF = "pair_buck.html">pair_style</A> buck/coul/long
<LI><A HREF = "pair_lj.html">pair_style</A> lj/cut
<LI><A HREF = "pair_lj.html">pair_style</A> lj/cut/coul/cut
<LI><A HREF = "pair_lj.html">pair_style</A> lj/cut/coul/long
</UL>
<UL><LI><A HREF = "pair_hbond_dreiding.html">pair_style</A> hbond/dreiding/lj
<LI><A HREF = "pair_hbond_dreiding.html">pair_style</A> hbond/dreiding/morse
</UL>
<UL><LI><A HREF = "special_bonds.html">special_bonds</A> dreiding
</UL>
<HR>
<A NAME = "howto_4"></A><H4>6.4 Running multiple simulations from one input script
</H4>
<P>This can be done in several ways. See the documentation for
individual commands for more details on how these examples work.
</P>
<P>If "multiple simulations" means continue a previous simulation for
more timesteps, then you simply use the <A HREF = "run.html">run</A> command
multiple times. For example, this script
</P>
<PRE>units lj
atom_style atomic
read_data data.lj
run 10000
run 10000
run 10000
run 10000
run 10000
</PRE>
<P>would run 5 successive simulations of the same system for a total of
50,000 timesteps.
</P>
<P>If you wish to run totally different simulations, one after the other,
the <A HREF = "clear.html">clear</A> command can be used in between them to
re-initialize LAMMPS. For example, this script
</P>
<PRE>units lj
atom_style atomic
read_data data.lj
run 10000
clear
units lj
atom_style atomic
read_data data.lj.new
run 10000
</PRE>
<P>would run 2 independent simulations, one after the other.
</P>
<P>For large numbers of independent simulations, you can use
<A HREF = "variable.html">variables</A> and the <A HREF = "next.html">next</A> and
<A HREF = "jump.html">jump</A> commands to loop over the same input script
multiple times with different settings. For example, this
script, named in.polymer
</P>
<PRE>variable d index run1 run2 run3 run4 run5 run6 run7 run8
shell cd $d
read_data data.polymer
run 10000
shell cd ..
clear
next d
jump in.polymer
</PRE>
<P>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
</P>
<PRE>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
</PRE>
<P>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 <A HREF = "Section_start.html#start_6">this
section</A> of the manual.
</P>
<P>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 <I>universe</I>-style variables, as described in the
<A HREF = "variable.html">variable</A> 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.
</P>
<HR>
<A NAME = "howto_5"></A><H4>6.5 Multi-replica simulations
</H4>
<P>Several commands in LAMMPS run mutli-replica simulations, meaning
that multiple instances (replicas) of your simulation are run
simultaneously, with small amounts of data exchanged between replicas
periodically.
</P>
<P>These are the relevant commands:
</P>
<UL><LI><A HREF = "neb.html">neb</A> for nudged elastic band calculations
<LI><A HREF = "prd.html">prd</A> for parallel replica dynamics
<LI><A HREF = "tad.html">tad</A> for temperature accelerated dynamics
<LI><A HREF = "temper.html">temper</A> for parallel tempering
</UL>
<P>NEB is a method for finding transition states and barrier energies.
PRD and TAD are methods for performing accelerated dynamics to find
and perform infrequent events. Parallel tempering or replica exchange
runs different replicas at a series of temperature to facilitate
rare-event sampling.
</P>
<P>These command can only be used if LAMMPS was built with the "replica"
package. See the <A HREF = "Section_start.html#start_3">Making LAMMPS</A> section
for more info on packages.
</P>
<P>In all these cases, you must run with one or more processors per
replica. The processors assigned to each replica are determined at
run-time by using the <A HREF = "Section_start.html#start_6">-partition command-line
switch</A> to launch LAMMPS on multiple
partitions, which in this context are the same as replicas. E.g.
these commands:
</P>
<PRE>mpirun -np 16 lmp_linux -partition 8x2 -in in.temper
mpirun -np 8 lmp_linux -partition 8x1 -in in.neb
</PRE>
<P>would each run 8 replicas, on either 16 or 8 processors. Note the use
of the <A HREF = "Section_start.html#start_6">-in command-line switch</A> to specify
the input script which is required when running in multi-replica mode.
</P>
<P>Also note that with MPI installed on a machine (e.g. your desktop),
you can run on more (virtual) processors than you have physical
processors. Thus the above commands could be run on a
single-processor (or few-processor) desktop so that you can run
a multi-replica simulation on more replicas than you have
physical processors.
</P>
<HR>
<A NAME = "howto_6"></A><H4>6.6 Granular models
</H4>
<P>Granular system are composed of spherical particles with a diameter,
as opposed to point particles. This means they have an angular
velocity and torque can be imparted to them to cause them to rotate.
</P>
<P>To run a simulation of a granular model, you will want to use
the following commands:
</P>
<UL><LI><A HREF = "atom_style.html">atom_style sphere</A>
<LI><A HREF = "fix_nve_sphere.html">fix nve/sphere</A>
<LI><A HREF = "fix_gravity.html">fix gravity</A>
</UL>
<P>This compute
</P>
<UL><LI><A HREF = "compute_erotate_sphere.html">compute erotate/sphere</A>
</UL>
<P>calculates rotational kinetic energy which can be <A HREF = "Section_howto.html#howto_15">output with
thermodynamic info</A>.
</P>
<P>Use one of these 3 pair potentials, which compute forces and torques
between interacting pairs of particles:
</P>
<UL><LI><A HREF = "pair_style.html">pair_style</A> gran/history
<LI><A HREF = "pair_style.html">pair_style</A> gran/no_history
<LI><A HREF = "pair_style.html">pair_style</A> gran/hertzian
</UL>
<P>These commands implement fix options specific to granular systems:
</P>
<UL><LI><A HREF = "fix_freeze.html">fix freeze</A>
<LI><A HREF = "fix_pour.html">fix pour</A>
<LI><A HREF = "fix_viscous.html">fix viscous</A>
<LI><A HREF = "fix_wall_gran.html">fix wall/gran</A>
</UL>
<P>The fix style <I>freeze</I> zeroes both the force and torque of frozen
atoms, and should be used for granular system instead of the fix style
<I>setforce</I>.
</P>
<P>For computational efficiency, you can eliminate needless pairwise
computations between frozen atoms by using this command:
</P>
<UL><LI><A HREF = "neigh_modify.html">neigh_modify</A> exclude
</UL>
<HR>
<A NAME = "howto_7"></A><H4>6.7 TIP3P water model
</H4>
<P>The TIP3P water model as implemented in CHARMM
<A HREF = "#MacKerell">(MacKerell)</A> specifies a 3-site rigid water molecule with
charges and Lennard-Jones parameters assigned to each of the 3 atoms.
In LAMMPS the <A HREF = "fix_shake.html">fix shake</A> command can be used to hold
the two O-H bonds and the H-O-H angle rigid. A bond style of
<I>harmonic</I> and an angle style of <I>harmonic</I> or <I>charmm</I> should also be
used.
</P>
<P>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
<A HREF = "#Jorgensen">(Jorgensen)</A>.
</P>
<P>O mass = 15.9994<BR>
H mass = 1.008 <BR>
</P>
<P>O charge = -0.834<BR>
H charge = 0.417 <BR>
</P>
<P>LJ epsilon of OO = 0.1521<BR>
LJ sigma of OO = 3.1507<BR>
LJ epsilon of HH = 0.0460<BR>
LJ sigma of HH = 0.4000<BR>
LJ epsilon of OH = 0.0836<BR>
LJ sigma of OH = 1.7753 <BR>
</P>
<P>K of OH bond = 450<BR>
r0 of OH bond = 0.9572 <BR>
</P>
<P>K of HOH angle = 55<BR>
theta of HOH angle = 104.52 <BR>
</P>
<P>These are the parameters to use for TIP3P with a long-range Coulombic
solver (Ewald or PPPM in LAMMPS), see <A HREF = "#Price">(Price)</A> for details:
</P>
<P>O mass = 15.9994<BR>
H mass = 1.008 <BR>
</P>
<P>O charge = -0.830<BR>
H charge = 0.415 <BR>
</P>
<P>LJ epsilon of OO = 0.102<BR>
LJ sigma of OO = 3.188<BR>
LJ epsilon, sigma of OH, HH = 0.0 <BR>
</P>
<P>K of OH bond = 450<BR>
r0 of OH bond = 0.9572 <BR>
</P>
<P>K of HOH angle = 55<BR>
theta of HOH angle = 104.52 <BR>
</P>
<P>Wikipedia also has a nice article on <A HREF = "http://en.wikipedia.org/wiki/Water_model">water
models</A>.
</P>
<HR>
<A NAME = "howto_8"></A><H4>6.8 TIP4P water model
</H4>
<P>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 <I>harmonic</I> and an
angle style of <I>harmonic</I> or <I>charmm</I> should also be used.
</P>
-<P>Currently, only a four-point model for long-range Coulombics is
-implemented via the LAMMPS <A HREF = "pair_lj.html">pair style
-lj/cut/coul/long/tip4p</A>. A cutoff version may be added
-the future. For both models, the bond lengths and bond angles should
-be held fixed using the <A HREF = "fix_shake.html">fix shake</A> command.
+<P>A TIP4P model is run with LAMMPS using two commands:
+</P>
+<UL><LI><A HREF = "pair_lj.html">pair_style lj/cut/coul/long/tip4p</A>
+<LI><A HREF = "kspace_style.html">kspace_style pppm/tip4p</A>
+</UL>
+<P>Note that only a TIP4P model with long-range Coulombics is currently
+implemented. A cutoff version may be added in the future. for both
+models, the bond lengths and bond angles should be held fixed using
+the <A HREF = "fix_shake.html">fix shake</A> command.
</P>
<P>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
<A HREF = "#Jorgensen">(Jorgensen)</A>. Note that the OM distance is specified in
the <A HREF = "pair_style.html">pair_style</A> command, not as part of the pair
coefficients.
</P>
<P>O mass = 15.9994<BR>
H mass = 1.008 <BR>
</P>
<P>O charge = -1.040<BR>
H charge = 0.520 <BR>
</P>
<P>r0 of OH bond = 0.9572<BR>
theta of HOH angle = 104.52 <BR>
</P>
<P>OM distance = 0.15 <BR>
</P>
<P>LJ epsilon of O-O = 0.1550<BR>
LJ sigma of O-O = 3.1536<BR>
LJ epsilon, sigma of OH, HH = 0.0 <BR>
</P>
<P>These are the parameters to use for TIP4P with a long-range Coulombic
solver (Ewald or PPPM in LAMMPS):
</P>
<P>O mass = 15.9994<BR>
H mass = 1.008 <BR>
</P>
<P>O charge = -1.0484<BR>
H charge = 0.5242 <BR>
</P>
<P>r0 of OH bond = 0.9572<BR>
theta of HOH angle = 104.52 <BR>
</P>
<P>OM distance = 0.1250 <BR>
</P>
<P>LJ epsilon of O-O = 0.16275<BR>
LJ sigma of O-O = 3.16435<BR>
LJ epsilon, sigma of OH, HH = 0.0 <BR>
</P>
<P>Wikipedia also has a nice article on <A HREF = "http://en.wikipedia.org/wiki/Water_model">water
models</A>.
</P>
<HR>
<A NAME = "howto_9"></A><H4>6.9 SPC water model
</H4>
<P>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 <A HREF = "fix_shake.html">fix shake</A> command can be used to hold
the two O-H bonds and the H-O-H angle rigid. A bond style of
<I>harmonic</I> and an angle style of <I>harmonic</I> or <I>charmm</I> should also be
used.
</P>
<P>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.
</P>
<P>O mass = 15.9994<BR>
H mass = 1.008 <BR>
</P>
<P>O charge = -0.820<BR>
H charge = 0.410 <BR>
</P>
<P>LJ epsilon of OO = 0.1553<BR>
LJ sigma of OO = 3.166<BR>
LJ epsilon, sigma of OH, HH = 0.0 <BR>
</P>
<P>r0 of OH bond = 1.0<BR>
theta of HOH angle = 109.47 <BR>
</P>
<P>Note that as originally proposed, the SPC model was run with a 9
Angstrom cutoff for both LJ and Coulommbic terms. It can also be used
with long-range Coulombics (Ewald or PPPM in LAMMPS), without changing
any of the parameters above, though it becomes a different model in
that mode of usage.
</P>
<P>The SPC/E (extended) water model is the same, except
the partial charge assignemnts change:
</P>
<P>O charge = -0.8476<BR>
H charge = 0.4238 <BR>
</P>
<P>See the <A HREF = "#Berendsen">(Berendsen)</A> reference for more details on both
the SPC and SPC/E models.
</P>
<P>Wikipedia also has a nice article on <A HREF = "http://en.wikipedia.org/wiki/Water_model">water
models</A>.
</P>
<HR>
<A NAME = "howto_10"></A><H4>6.10 Coupling LAMMPS to other codes
</H4>
<P>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.
</P>
<P>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.
</P>
<P>(1) Define a new <A HREF = "fix.html">fix</A> 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
<A HREF = "http://www.rpi.edu/~anderk5/lab">POEMS</A> package that performs constrained rigid-body motion on
groups of atoms is hooked to LAMMPS. See the
<A HREF = "fix_poems.html">fix_poems</A> command for more details. See <A HREF = "Section_modify.html">this
section</A> of the documentation for info on how to add
a new fix to LAMMPS.
</P>
<P>(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
<A HREF = "run.html">run</A> command facilitates this with its <I>every</I> option, which
makes it easy to run a few steps, invoke the command, run a few steps,
invoke the command, etc.
</P>
<P>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.
</P>
<P>See <A HREF = "Section_modify.html">Section_modify</A> of the documentation for how
to add a new command to LAMMPS.
</P>
<P>(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 <A HREF = "run.html">run</A> 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.
</P>
<P>Examples of driver codes that call LAMMPS as a library are included in
the "couple" directory of the LAMMPS distribution; see couple/README
for more details:
</P>
<UL><LI>simple: simple driver programs in C++ and C which invoke LAMMPS as a
library
<LI>lammps_quest: coupling of LAMMPS and <A HREF = "http://dft.sandia.gov/Quest">Quest</A>, to run classical
MD with quantum forces calculated by a density functional code
<LI>lammps_spparks: coupling of LAMMPS and <A HREF = "http://www.sandia.gov/~sjplimp/spparks.html">SPPARKS</A>, to couple
a kinetic Monte Carlo model for grain growth using MD to calculate
strain induced across grain boundaries
</UL>
<P><A HREF = "Section_start.html#start_4">This section</A> of the documentation
describes how to build LAMMPS as a library. Once this is done, you
can interface with LAMMPS either via C++, C, Fortran, or Python (or
any other language that supports a vanilla C-like interface). 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. See
<A HREF = "Section_python.html">Section_python</A> of the manual for a description
of the Python wrapper provided with LAMMPS that operates through the
LAMMPS library interface.
</P>
<P>The files src/library.cpp and library.h contain the C-style interface
to LAMMPS. See <A HREF = "Section_howto.html#howto_19">Section_howto 19</A> of the
manual for a description of the interface and how to extend it for
your needs.
</P>
<P>Note that the lammps_open() function that creates an instance of
LAMMPS takes an MPI communicator as an argument. This means that
instance of LAMMPS will run on the set of processors in the
communicator. Thus 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. Or it might instantiate multiple
instances of LAMMPS to perform different calculations.
</P>
<HR>
<A NAME = "howto_11"></A><H4>6.11 Visualizing LAMMPS snapshots
</H4>
<P>LAMMPS itself does not do visualization, but snapshots from LAMMPS
simulations can be visualized (and analyzed) in a variety of ways.
</P>
<P>LAMMPS snapshots are created by the <A HREF = "dump.html">dump</A> 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 <A HREF = "Section_tools.html#xmovie">xmovie</A> program, included with the
LAMMPS package. This produces simple, fast 2d projections of 3d
systems, and can be useful for rapid debugging of simulation geometry
and atom trajectories.
</P>
<P>Several programs included with LAMMPS as auxiliary tools can convert
native LAMMPS dump files to other formats. See the
<A HREF = "Section_tools.html">Section_tools</A> doc page for details. The first is
the <A HREF = "Section_tools.html#charmm">ch2lmp tool</A>, which contains a
lammps2pdb Perl script which converts LAMMPS dump files into PDB
files. The second is the <A HREF = "Section_tools.html#arc">lmp2arc tool</A> which
converts LAMMPS dump files into Accelrys' Insight MD program files.
The third is the <A HREF = "Section_tools.html#cfg">lmp2cfg tool</A> which converts
LAMMPS dump files into CFG files which can be read into the
<A HREF = "http://mt.seas.upenn.edu/Archive/Graphics/A">AtomEye</A> visualizer.
</P>
<P>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 <A HREF = "http://www.sandia.gov/~sjplimp/pizza.html">Pizza.py
WWW site</A> for details. Specifically, Pizza.py can convert
LAMMPS dump files into PDB, XYZ, <A HREF = "http://www.ensight.com">Ensight</A>, 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.
</P>
<P>LAMMPS can create XYZ files directly (via "dump xyz") which is a
simple text-based file format used by many visualization programs
including <A HREF = "http://www.ks.uiuc.edu/Research/vmd">VMD</A>.
</P>
<P>LAMMPS can create DCD files directly (via "dump dcd") which can be
read by <A HREF = "http://www.ks.uiuc.edu/Research/vmd">VMD</A> in conjunction with a CHARMM PSF file. Using this
form of output avoids the need to convert LAMMPS snapshots to PDB
files. See the <A HREF = "dump.html">dump</A> command for more information on DCD
files.
</P>
<P>LAMMPS can create XTC files directly (via "dump xtc") which is GROMACS
file format which can also be read by <A HREF = "http://www.ks.uiuc.edu/Research/vmd">VMD</A> for visualization.
See the <A HREF = "dump.html">dump</A> command for more information on XTC files.
</P>
<HR>
<A NAME = "howto_12"></A><H4>6.12 Triclinic (non-orthogonal) simulation boxes
</H4>
<P>By default, LAMMPS uses an orthogonal simulation box to encompass the
particles. The <A HREF = "boundary.html">boundary</A> command sets the boundary
conditions of the box (periodic, non-periodic, etc). The orthogonal
box has its "origin" at (xlo,ylo,zlo) and is defined by 3 edge vectors
starting from the origin given by <B>a</B> = (xhi-xlo,0,0); <B>b</B> =
(0,yhi-ylo,0); <B>c</B> = (0,0,zhi-zlo). The 6 parameters
(xlo,xhi,ylo,yhi,zlo,zhi) are defined at the time the simluation box
is created, e.g. by the <A HREF = "create_box.html">create_box</A> or
<A HREF = "read_data.html">read_data</A> or <A HREF = "read_restart.html">read_restart</A>
commands. Additionally, LAMMPS defines box size parameters lx,ly,lz
where lx = xhi-xlo, and similarly in the y and z dimensions. The 6
parameters, as well as lx,ly,lz, can be output via the <A HREF = "thermo_style.html">thermo_style
custom</A> command.
</P>
<P>LAMMPS also allows simulations to be perfored in non-orthogonal
simulation boxes shaped as a parallelepiped with triclinic symmetry.
The parallelepiped has its "origin" at (xlo,ylo,zlo) and is defined by
3 edge vectors starting from the origin given by <B>a</B> = (xhi-xlo,0,0); <B>b</B>
= (xy,yhi-ylo,0); <B>c</B> = (xz,yz,zhi-zlo). <I>Xy,xz,yz</I> can be 0.0 or
positive or negative values and are called "tilt factors" because they
are the amount of displacement applied to faces of an originally
orthogonal box to transform it into the parallelepiped. Note that in
LAMMPS the triclinic simulation box edge vectors <B>a</B>, <B>b</B>, and <B>c</B> cannot be
arbitrary vectors. As indicated, <B>a</B> must be aligned with the x axis, <B>b</B>
must be in the xy plane, and <B>c</B> is arbitrary. However, this is not a
restriction since it is possible to rotate any set of 3 crystal basis
vectors so that they meet this restriction.
</P>
<P>The 9 parameters (xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz) are defined at the
time the simluation box is created. This happens in one of 3 ways.
If the <A HREF = "create_box.html">create_box</A> command is used with a region of
style <I>prism</I>, then a triclinic box is setup. See the
<A HREF = "region.html">region</A> command for details. If the
<A HREF = "read_data.html">read_data</A> 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 triclinic box is setup. See the
<A HREF = "read_data.html">read_data</A> command for details. Finally, if the
<A HREF = "read_restart.html">read_restart</A> command reads a restart file which
was written from a simulation using a triclinic box, then a triclinic
box will be setup for the restarted simulation.
</P>
<P>Note that you can define a triclinic box with all 3 tilt factors =
0.0, so that it is initially orthogonal. This is necessary if the box
will become non-orthogonal, e.g. due to the <A HREF = "fix_nh.html">fix npt</A> or
<A HREF = "fix_deform.html">fix deform</A> commands. Alternatively, you can use the
<A HREF = "change_box.html">change_box</A> command to convert a simulation box from
orthogonal to triclinic and vice versa.
</P>
<P>As with orthogonal boxes, LAMMPS defines triclinic box size parameters
lx,ly,lz where lx = xhi-xlo, and similarly in the y and z dimensions.
The 9 parameters, as well as lx,ly,lz, can be output via the
<A HREF = "thermo_style.html">thermo_style custom</A> command.
</P>
<P>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 geometrically all
equivalent.
</P>
<P>Triclinic crystal structures are often defined using three lattice
constants <I>a</I>, <I>b</I>, and <I>c</I>, and three angles <I>alpha</I>, <I>beta</I> and
<I>gamma</I>. Note that in this nomenclature, the a, b, and c lattice constants
are the scalar lengths of the edge vectors <B>a</B>, <B>b</B>, and <B>c</B> defined
above. The
relationship between these 6 quantities (a,b,c,alpha,beta,gamma) and
the LAMMPS box sizes (lx,ly,lz) = (xhi-xlo,yhi-ylo,zhi-zlo) and tilt
factors (xy,xz,yz) is as follows:
</P>
<CENTER><IMG SRC = "Eqs/box.jpg">
</CENTER>
<P>The inverse relationship can be written as follows:
</P>
<CENTER><IMG SRC = "Eqs/box_inverse.jpg">
</CENTER>
<P>The values of <I>a</I>, <I>b</I>, <I>c</I> , <I>alpha</I>, <I>beta</I> , and <I>gamma</I> can be printed
out or accessed by computes using the
<A HREF = "thermo_style.html">thermo_style custom</A> keywords
<I>cella</I>, <I>cellb</I>, <I>cellc</I>, <I>cellalpha</I>, <I>cellbeta</I>, <I>cellgamma</I>,
respectively.
</P>
<P>As discussed on the <A HREF = "dump.html">dump</A> command doc page, when the BOX
BOUNDS for a snapshot is written to a dump file for a triclinic box,
an orthogonal bounding box which encloses the triclinic simulation box
is output, along with the 3 tilt factors (xy, xz, yz) of the triclinic
box, formatted as follows:
</P>
<PRE>ITEM: BOX BOUNDS xy xz yz
xlo_bound xhi_bound xy
ylo_bound yhi_bound xz
zlo_bound zhi_bound yz
</PRE>
<P>This bounding box is convenient for many visualization programs and is
calculated from the 9 triclinic box parameters
(xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz) as follows:
</P>
<PRE>xlo_bound = xlo + MIN(0.0,xy,xz,xy+xz)
xhi_bound = xhi + MAX(0.0,xy,xz,xy+xz)
ylo_bound = ylo + MIN(0.0,yz)
yhi_bound = yhi + MAX(0.0,yz)
zlo_bound = zlo
zhi_bound = zhi
</PRE>
<P>These formulas can be inverted if you need to convert the bounding box
back into the triclinic box parameters, e.g. xlo = xlo_bound -
MIN(0.0,xy,xz,xy+xz).
</P>
<P>One use of triclinic simulation boxes is to model solid-state crystals
with triclinic symmetry. The <A HREF = "lattice.html">lattice</A> command can be
used with non-orthogonal basis vectors to define a lattice that will
tile a triclinic simulation box via the
<A HREF = "create_atoms.html">create_atoms</A> command.
</P>
<P>A second use is to run Parinello-Rahman dyanamics via the <A HREF = "fix_nh.html">fix
npt</A> command, which will adjust the xy, xz, yz tilt
factors to compensate for off-diagonal components of the pressure
tensor. The analalog for an <A HREF = "minimize.html">energy minimization</A> is
the <A HREF = "fix_box_relax.html">fix box/relax</A> command.
</P>
<P>A third use is to shear a bulk solid to study the response of the
material. The <A HREF = "fix_deform.html">fix deform</A> command can be used for
this purpose. It allows dynamic control of the xy, xz, yz tilt
factors as a simulation runs. This is discussed in the next section
on non-equilibrium MD (NEMD) simulations.
</P>
<HR>
<A NAME = "howto_13"></A><H4>6.13 NEMD simulations
</H4>
<P>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 preceding Howto section).
</P>
<P>A shear strain can be applied to the simulation box at a desired
strain rate by using the <A HREF = "fix_deform.html">fix deform</A> command. The
<A HREF = "fix_nvt_sllod.html">fix nvt/sllod</A> command can be used to thermostat
the sheared fluid and integrate the SLLOD equations of motion for the
system. Fix nvt/sllod uses <A HREF = "compute_temp_deform.html">compute
temp/deform</A> 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 <A HREF = "fix_ave_spatial.html">fix ave/spatial</A> command.
</P>
<P>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, <A HREF = "fix_deform.html">fix deform</A> can continuously strain
a box by an arbitrary amount. As discussed in the <A HREF = "fix_deform.html">fix
deform</A> 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.
</P>
<P>In a NEMD simulation, the "remap" option of <A HREF = "fix_deform.html">fix
deform</A> should be set to "remap v", since that is what
<A HREF = "fix_nvt_sllod.html">fix nvt/sllod</A> assumes to generate a velocity
profile consistent with the applied shear strain rate.
</P>
<P>An alternative method for calculating viscosities is provided via the
<A HREF = "fix_viscosity.html">fix viscosity</A> command.
</P>
<HR>
<A NAME = "howto_14"></A><H4>6.14 Extended spherical and aspherical particles
</H4>
<P>Typical MD models treat atoms or particles as point masses.
Sometimes, however, it is desirable to have a model with finite-size
particles such as spheres or aspherical ellipsoids. The difference is
that such particles have a moment of inertia, rotational energy, and
angular momentum. Rotation is induced by torque from interactions
with other particles.
</P>
<P>LAMMPS has several options for running simulations with these kinds of
particles. The following aspects are discussed in turn:
</P>
<UL><LI>atom styles
<LI>pair potentials
<LI>time integration
<LI>computes, thermodynamics, and dump output
<LI>rigid bodies composed of extended particles
</UL>
<H5>Atom styles
</H5>
<P>There are 2 <A HREF = "atom_style.html">atom styles</A> that allow for definition of
finite-size particles: sphere and ellipsoid. The peri atom style also
treats particles as having a volume, but that is internal to the
<A HREF = "pair_peri.html">pair_style peri</A> potentials. The dipole atom style is
most often used in conjunction with finite-size particles.
</P>
<P>The sphere style defines particles that are spheriods and each
particle can have a unique diameter and mass (or density). These
particles store an angular velocity (omega) and can be acted upon by
torque. The "set" command can be used to modify the diameter and mass
of individual particles, after then are created.
</P>
<P>The ellipsoid style defines particles that are ellipsoids and thus can
be aspherical. Each particle has a shape, specified by 3 diameters,
and mass (or density). These particles store an angular momentum and
their orientation (quaternion), and can be acted upon by torque. They
do not store an angular velocity (omega), which can be in a different
direction than angular momentum, rather they compute it as needed.
The "set" command can be used to modify the diameter, orientation, and
mass of individual particles, after then are created. It also has a
brief explanation of what quaternions are.
</P>
<P>The dipole style does not define extended particles, but is often
used in conjunction with spherical particles, via a command like
</P>
<PRE>atom_style hybrid sphere dipole
</PRE>
<P>This is because when dipoles interact with each other, they induce
torques, and a particle must be extended (i.e. have a moment of
inertia) in order to respond and rotate. See the <A HREF = "atom_style.html">atom_style
dipole</A> command for details. The "set" command can be
used to modify the orientation and length of the dipole moment of
individual particles, after then are created.
</P>
<P>Note that if one of these atom styles is used (or multiple styles via
the <A HREF = "atom_style.html">atom_style hybrid</A> command), not all particles in
the system are required to be finite-size or aspherical. For example,
if the 3 shape parameters are set to the same value, the particle will
be a sphere rather than an ellipsoid. If the 3 shape parameters are
all set to 0.0 or if the diameter is set to 0.0, it will be a point
particle. If the length of the dipole moment is set to zero, the
particle will not have a point dipole associated with it. The pair
styles used to compute pairwise interactions will typically compute
the correct interaction in these simplified (cheaper) cases.
<A HREF = "pair_hybrid.html">Pair_style hybrid</A> can be used to insure the correct
interactions are computed for the appropriate style of interactions.
Likewise, using groups to partition particles (ellipsoids versus
spheres versus point particles) will allow you to use the appropriate
time integrators and temperature computations for each class of
particles. See the doc pages for various commands for details.
</P>
<P>Also note that for <A HREF = "dimension.html">2d simulations</A>, finite-size
spheres and ellipsoids are still treated as 3d particles, rather than
as circular disks or ellipses. This means they have the same moment
of inertia for a 3d extended object. When their temperature is
coomputed, the correct degrees of freedom are used for rotation in a
2d versus 3d system.
</P>
<H5>Pair potentials
</H5>
<P>When a system with extended particles is defined, the particles will
only rotate and experience torque if the force field computes such
interactions. These are the various <A HREF = "pair_style.html">pair styles</A>
that generate torque:
</P>
<UL><LI><A HREF = "pair_gran.html">pair_style gran/history</A>
<LI><A HREF = "pair_gran.html">pair_style gran/hertzian</A>
<LI><A HREF = "pair_gran.html">pair_style gran/no_history</A>
<LI><A HREF = "pair_dipole.html">pair_style dipole/cut</A>
<LI><A HREF = "pair_gayberne.html">pair_style gayberne</A>
<LI><A HREF = "pair_resquared.html">pair_style resquared</A>
<LI><A HREF = "pair_lubricate.html">pair_style lubricate</A>
</UL>
<P>The <A HREF = "pair_gran.html">granular pair styles</A> are used with spherical
particles. The <A HREF = "pair_dipole.html">dipole pair style</A> is used with
<A HREF = "atom_style.html">atom_style dipole</A>, which could be applied to
spherical or ellipsoidal particles. The <A HREF = "pair_gayberne.html">GayBerne</A>
and <A HREF = "pair_resquared.html">REsquared</A> potentials require ellipsoidal
particles, though they will also work if the 3 shape parameters are
the same (a sphere). The <A HREF = "pair_lubricate.html">lubrication potential</A>
works with spherical particles.
</P>
<H5>Time integration
</H5>
<P>There are 3 fixes that perform time integration on extended spherical
particles, meaning the integrators update the rotational orientation
and angular velocity or angular momentum of the particles:
</P>
<UL><LI><A HREF = "fix_nve_sphere.html">fix nve/sphere</A>
<LI><A HREF = "fix_nvt_sphere.html">fix nvt/sphere</A>
<LI><A HREF = "fix_npt_sphere.html">fix npt/sphere</A>
</UL>
<P>Likewise, there are 3 fixes that perform time integration on
ellipsoids as extended aspherical particles:
</P>
<UL><LI><A HREF = "fix_nve_asphere.html">fix nve/asphere</A>
<LI><A HREF = "fix_nvt_asphere.html">fix nvt/asphere</A>
<LI><A HREF = "fix_npt_asphere.html">fix npt/asphere</A>
</UL>
<P>The advantage of these fixes is that those which thermostat the
particles include the rotational degrees of freedom in the temperature
calculation and thermostatting. Other thermostats can be used with
fix nve/sphere or fix nve/asphere, such as fix langevin or fix
temp/berendsen, but those thermostats only operate on the
translational kinetic energy of the extended particles.
</P>
<P>Note that for mixtures of point and extended particles, you should
only use these integration fixes on <A HREF = "group.html">groups</A> which contain
extended particles.
</P>
<H5>Computes, thermodynamics, and dump output
</H5>
<P>There are 4 computes that calculate the temperature or rotational energy
of extended spherical or aspherical particles (ellipsoids):
</P>
<UL><LI><A HREF = "compute_temp_sphere.html">compute temp/sphere</A>
<LI><A HREF = "compute_temp_asphere.html">compute temp/asphere</A>
<LI><A HREF = "compute_erotate_sphere.html">compute erotate/sphere</A>
<LI><A HREF = "compute_erotate_asphere.html">compute erotate/asphere</A>
</UL>
<P>These include rotational degrees of freedom in their computation. If
you wish the thermodynamic output of temperature or pressure to use
one of these computes (e.g. for a system entirely composed of extended
particles), then the compute can be defined and the
<A HREF = "thermo_modify.html">thermo_modify</A> command used. Note that by
default thermodynamic quantities will be calculated with a temperature
that only includes translational degrees of freedom. See the
<A HREF = "thermo_style.html">thermo_style</A> command for details.
</P>
<P>The <A HREF = "dump.html">dump custom</A> command can output various attributes of
extended particles, including the dipole moment (mu), the angular
velocity (omega), the angular momentum (angmom), the quaternion
(quat), and the torque (tq) on the particle.
</P>
<H5>Rigid bodies composed of extended particles
</H5>
<P>The <A HREF = "fix_rigid.html">fix rigid</A> command treats a collection of
particles as a rigid body, computes its inertia tensor, sums the total
force and torque on the rigid body each timestep due to forces on its
constituent particles, and integrates the motion of the rigid body.
</P>
<P>If any of the constituent particles of a rigid body are extended
particles (spheres or ellipsoids), then their contribution to the
inertia tensor of the body is different than if they were point
particles. This means the rotational dynamics of the rigid body will
be different. Thus a model of a dimer is different if the dimer
consists of two point masses versus two extended sphereoids, even if
the two particles have the same mass. Extended particles that
experience torque due to their interaction with other particles will
also impart that torque to a rigid body they are part of.
</P>
<P>See the "fix rigid" command for example of complex rigid-body models
it is possible to define in LAMMPS.
</P>
<P>Note that the <A HREF = "fix_shake.html">fix shake</A> command can also be used to
treat 2, 3, or 4 particles as a rigid body, but it always assumes the
particles are point masses.
</P>
<HR>
<A NAME = "howto_15"></A><H4>6.15 Output from LAMMPS (thermo, dumps, computes, fixes, variables)
</H4>
<P>There are four basic kinds of LAMMPS output:
</P>
<UL><LI><A HREF = "thermo_style.html">Thermodynamic output</A>, which is a list
of quantities printed every few timesteps to the screen and logfile.
<LI><A HREF = "dump.html">Dump files</A>, which contain snapshots of atoms and various
per-atom values and are written at a specified frequency.
<LI>Certain fixes can output user-specified quantities to files: <A HREF = "fix_ave_time.html">fix
ave/time</A> for time averaging, <A HREF = "fix_ave_spatial.html">fix
ave/spatial</A> for spatial averaging, and <A HREF = "fix_print.html">fix
print</A> for single-line output of
<A HREF = "variable.html">variables</A>. Fix print can also output to the
screen.
<LI><A HREF = "restart.html">Restart files</A>.
</UL>
<P>A simulation prints one set of thermodynamic output and (optionally)
restart files. It can generate any number of dump files and fix
output files, depending on what <A HREF = "dump.html">dump</A> and <A HREF = "fix.html">fix</A>
commands you specify.
</P>
<P>As discussed below, LAMMPS gives you a variety of ways to determine
what quantities are computed and printed when the thermodynamics,
dump, or fix commands listed above perform output. Throughout this
discussion, note that users can also <A HREF = "Section_modify.html">add their own computes and fixes
to LAMMPS</A> which can then generate values that can
then be output with these commands.
</P>
<P>The following sub-sections discuss different LAMMPS command related
to output and the kind of data they operate on and produce:
</P>
<UL><LI><A HREF = "#global">Global/per-atom/local data</A>
<LI><A HREF = "#scalar">Scalar/vector/array data</A>
<LI><A HREF = "#thermo">Thermodynamic output</A>
<LI><A HREF = "#dump">Dump file output</A>
<LI><A HREF = "#fixoutput">Fixes that write output files</A>
<LI><A HREF = "#computeoutput">Computes that process output quantities</A>
<LI><A HREF = "#fixoutput">Fixes that process output quantities</A>
<LI><A HREF = "#compute">Computes that generate values to output</A>
<LI><A HREF = "#fix">Fixes that generate values to output</A>
<LI><A HREF = "#variable">Variables that generate values to output</A>
<LI><A HREF = "#table">Summary table of output options and data flow between commands</A>
</UL>
<H5><A NAME = "global"></A>Global/per-atom/local data
</H5>
<P>Various output-related commands work with three different styles of
data: global, per-atom, or local. A global datum is one or more
system-wide values, e.g. the temperature of the system. A per-atom
datum is one or more values per atom, e.g. the kinetic energy of each
atom. Local datums are calculated by each processor based on the
atoms it owns, but there may be zero or more per atom, e.g. a list of
bond distances.
</P>
<H5><A NAME = "scalar"></A>Scalar/vector/array data
</H5>
<P>Global, per-atom, and local datums can each come in three kinds: a
single scalar value, a vector of values, or a 2d array of values. The
doc page for a "compute" or "fix" or "variable" that generates data
will specify both the style and kind of data it produces, e.g. a
per-atom vector.
</P>
<P>When a quantity is accessed, as in many of the output commands
discussed below, it can be referenced via the following bracket
notation, where ID in this case is the ID of a compute. The leading
"c_" would be replaced by "f_" for a fix, or "v_" for a variable:
</P>
<DIV ALIGN=center><TABLE BORDER=1 >
<TR><TD >c_ID </TD><TD > entire scalar, vector, or array</TD></TR>
<TR><TD >c_ID[I] </TD><TD > one element of vector, one column of array</TD></TR>
<TR><TD >c_ID[I][J] </TD><TD > one element of array
</TD></TR></TABLE></DIV>
<P>In other words, using one bracket reduces the dimension of the data
once (vector -> scalar, array -> vector). Using two brackets reduces
the dimension twice (array -> scalar). Thus a command that uses
scalar values as input can typically also process elements of a vector
or array.
</P>
<H5><A NAME = "thermo"></A>Thermodynamic output
</H5>
<P>The frequency and format of thermodynamic output is set by the
<A HREF = "thermo.html">thermo</A>, <A HREF = "thermo_style.html">thermo_style</A>, and
<A HREF = "thermo_modify.html">thermo_modify</A> commands. The
<A HREF = "thermo_style.html">thermo_style</A> command also specifies what values
are calculated and written out. Pre-defined keywords can be specified
(e.g. press, etotal, etc). Three additional kinds of keywords can
also be specified (c_ID, f_ID, v_name), where a <A HREF = "compute.html">compute</A>
or <A HREF = "fix.html">fix</A> or <A HREF = "variable.html">variable</A> provides the value to be
output. In each case, the compute, fix, or variable must generate
global values for input to the <A HREF = "dump.html">thermo_style custom</A>
command.
</P>
<H5><A NAME = "dump"></A>Dump file output
</H5>
<P>Dump file output is specified by the <A HREF = "dump.html">dump</A> and
<A HREF = "dump_modify.html">dump_modify</A> commands. There are several
pre-defined formats (dump atom, dump xtc, etc).
</P>
<P>There is also a <A HREF = "dump.html">dump custom</A> format where the user
specifies what values are output with each atom. Pre-defined atom
attributes can be specified (id, x, fx, etc). Three additional kinds
of keywords can also be specified (c_ID, f_ID, v_name), where a
<A HREF = "compute.html">compute</A> or <A HREF = "fix.html">fix</A> or <A HREF = "variable.html">variable</A>
provides the values to be output. In each case, the compute, fix, or
variable must generate per-atom values for input to the <A HREF = "dump.html">dump
custom</A> command.
</P>
<P>There is also a <A HREF = "dump.html">dump local</A> format where the user specifies
what local values to output. A pre-defined index keyword can be
specified to enumuerate the local values. Two additional kinds of
keywords can also be specified (c_ID, f_ID), where a
<A HREF = "compute.html">compute</A> or <A HREF = "fix.html">fix</A> or <A HREF = "variable.html">variable</A>
provides the values to be output. In each case, the compute or fix
must generate local values for input to the <A HREF = "dump.html">dump local</A>
command.
</P>
<H5><A NAME = "fixoutput"></A>Fixes that write output files
</H5>
<P>Sevarl fixes take various quantities as input and can write output
files: <A HREF = "fix_ave_time.html">fix ave/time</A>, <A HREF = "fix_ave_spatial.html">fix
ave/spatial</A>, <A HREF = "fix_ave_histo.html">fix ave/histo</A>,
<A HREF = "fix_ave_correlate.html">fix ave/correlate</A>, and <A HREF = "fix_print.html">fix
print</A>.
</P>
<P>The <A HREF = "fix_ave_time.html">fix ave/time</A> command enables direct output to
a file and/or time-averaging of global scalars or vectors. The user
specifies one or more quantities as input. These can be global
<A HREF = "compute.html">compute</A> values, global <A HREF = "fix.html">fix</A> values, or
<A HREF = "variable.html">variables</A> of any style except the atom style which
produces per-atom values. Since a variable can refer to keywords used
by the <A HREF = "thermo_style.html">thermo_style custom</A> command (like temp or
press) and individual per-atom values, a wide variety of quantities
can be time averaged and/or output in this way. If the inputs are one
or more scalar values, then the fix generate a global scalar or vector
of output. If the inputs are one or more vector values, then the fix
generates a global vector or array of output. The time-averaged
output of this fix can also be used as input to other output commands.
</P>
<P>The <A HREF = "fix_ave_spatial.html">fix ave/spatial</A> command enables direct
output to a file of spatial-averaged per-atom quantities like those
output in dump files, within 1d layers of the simulation box. The
per-atom quantities can be atom density (mass or number) or atom
attributes such as position, velocity, force. They can also be
per-atom quantities calculated by a <A HREF = "compute.html">compute</A>, by a
<A HREF = "fix.html">fix</A>, or by an atom-style <A HREF = "variable.html">variable</A>. The
spatial-averaged output of this fix can also be used as input to other
output commands.
</P>
<P>The <A HREF = "fix_ave_histo.html">fix ave/histo</A> command enables direct output
to a file of histogrammed quantities, which can be global or per-atom
or local quantities. The histogram output of this fix can also be
used as input to other output commands.
</P>
<P>The <A HREF = "fix_ave_histo.html">fix ave/correlate</A> command enables direct
output to a file of time-correlated quantities, which can be global
scalars. The correlation matrix output of this fix can also be used
as input to other output commands.
</P>
<P>The <A HREF = "fix_print.html">fix print</A> command can generate a line of output
written to the screen and log file or to a separate file, periodically
during a running simulation. The line can contain one or more
<A HREF = "variable.html">variable</A> values for any style variable except the atom
style). As explained above, variables themselves can contain
references to global values generated by <A HREF = "thermo_style.html">thermodynamic
keywords</A>, <A HREF = "compute.html">computes</A>,
<A HREF = "fix.html">fixes</A>, or other <A HREF = "variable.html">variables</A>, or to per-atom
values for a specific atom. Thus the <A HREF = "fix_print.html">fix print</A>
command is a means to output a wide variety of quantities separate
from normal thermodynamic or dump file output.
</P>
<H5><A NAME = "computeoutput"></A>Computes that process output quantities
</H5>
<P>The <A HREF = "compute_reduce.html">compute reduce</A> and <A HREF = "compute_reduce.html">compute
reduce/region</A> commands take one or more per-atom
or local vector quantities as inputs and "reduce" them (sum, min, max,
ave) to scalar quantities. These are produced as output values which
can be used as input to other output commands.
</P>
<P>The <A HREF = "compute_slice.html">compute slice</A> command take one or more global
vector or array quantities as inputs and extracts a subset of their
values to create a new vector or array. These are produced as output
values which can be used as input to other output commands.
</P>
<P>The <A HREF = "compute_property_atom.html">compute property/atom</A> command takes a
list of one or more pre-defined atom attributes (id, x, fx, etc) and
stores the values in a per-atom vector or array. These are produced
as output values which can be used as input to other output commands.
The list of atom attributes is the same as for the <A HREF = "dump.html">dump
custom</A> command.
</P>
<P>The <A HREF = "compute_property_local.html">compute property/local</A> command takes
a list of one or more pre-defined local attributes (bond info, angle
info, etc) and stores the values in a local vector or array. These
are produced as output values which can be used as input to other
output commands.
</P>
<P>The <A HREF = "compute_atom_molecule.html">compute atom/molecule</A> command takes a
list of one or more per-atom quantities (from a compute, fix, per-atom
variable) and sums the quantities on a per-molecule basis. It
produces a global vector or array as output values which can be used
as input to other output commands.
</P>
<H5><A NAME = "fixoutput"></A>Fixes that process output quantities
</H5>
<P>The <A HREF = "fix_ave_atom.html">fix ave/atom</A> command performs time-averaging
of per-atom vectors. The per-atom quantities can be atom attributes
such as position, velocity, force. They can also be per-atom
quantities calculated by a <A HREF = "compute.html">compute</A>, by a
<A HREF = "fix.html">fix</A>, or by an atom-style <A HREF = "variable.html">variable</A>. The
time-averaged per-atom output of this fix can be used as input to
other output commands.
</P>
<P>The <A HREF = "fix_store_state.html">fix store/state</A> command can archive one or
more per-atom attributes at a particular time, so that the old values
can be used in a future calculation or output. The list of atom
attributes is the same as for the <A HREF = "dump.html">dump custom</A> command,
including per-atom quantities calculated by a <A HREF = "compute.html">compute</A>,
by a <A HREF = "fix.html">fix</A>, or by an atom-style <A HREF = "variable.html">variable</A>.
The output of this fix can be used as input to other output commands.
</P>
<H5><A NAME = "compute"></A>Computes that generate values to output
</H5>
<P>Every <A HREF = "compute.html">compute</A> in LAMMPS produces either global or
per-atom or local values. The values can be scalars or vectors or
arrays of data. These values can be output using the other commands
described in this section. The doc page for each compute command
describes what it produces. Computes that produce per-atom or local
values have the word "atom" or "local" in their style name. Computes
without the word "atom" or "local" produce global values.
</P>
<H5><A NAME = "fix"></A>Fixes that generate values to output
</H5>
<P>Some <A HREF = "fix.html">fixes</A> in LAMMPS produces either global or per-atom or
local values which can be accessed by other commands. The values can
be scalars or vectors or arrays of data. These values can be output
using the other commands described in this section. The doc page for
each fix command tells whether it produces any output quantities and
describes them.
</P>
<H5><A NAME = "variable"></A>Variables that generate values to output
</H5>
<P>Every <A HREF = "variable.html">variables</A> defined in an input script generates
either a global scalar value or a per-atom vector (only atom-style
variables) when it is accessed. The formulas used to define equal-
and atom-style variables can contain references to the thermodynamic
keywords and to global and per-atom data generated by computes, fixes,
and other variables. The values generated by variables can be output
using the other commands described in this section.
</P>
<H5><A NAME = "table"></A>Summary table of output options and data flow between commands
</H5>
<P>This table summarizes the various commands that can be used for
generating output from LAMMPS. Each command produces output data of
some kind and/or writes data to a file. Most of the commands can take
data from other commands as input. Thus you can link many of these
commands together in pipeline form, where data produced by one command
is used as input to another command and eventually written to the
screen or to a file. Note that to hook two commands together the
output and input data types must match, e.g. global/per-atom/local
data and scalar/vector/array data.
</P>
<P>Also note that, as described above, when a command takes a scalar as
input, that could be an element of a vector or array. Likewise a
vector input could be a column of an array.
</P>
<DIV ALIGN=center><TABLE BORDER=1 >
<TR><TD >Command</TD><TD > Input</TD><TD > Output</TD><TD ></TD></TR>
<TR><TD ><A HREF = "thermo_style.html">thermo_style custom</A></TD><TD > global scalars</TD><TD > screen, log file</TD><TD ></TD></TR>
<TR><TD ><A HREF = "dump.html">dump custom</A></TD><TD > per-atom vectors</TD><TD > dump file</TD><TD ></TD></TR>
<TR><TD ><A HREF = "dump.html">dump local</A></TD><TD > local vectors</TD><TD > dump file</TD><TD ></TD></TR>
<TR><TD ><A HREF = "fix_print.html">fix print</A></TD><TD > global scalar from variable</TD><TD > screen, file</TD><TD ></TD></TR>
<TR><TD ><A HREF = "print.html">print</A></TD><TD > global scalar from variable</TD><TD > screen</TD><TD ></TD></TR>
<TR><TD ><A HREF = "compute.html">computes</A></TD><TD > N/A</TD><TD > global/per-atom/local scalar/vector/array</TD><TD ></TD></TR>
<TR><TD ><A HREF = "fix.html">fixes</A></TD><TD > N/A</TD><TD > global/per-atom/local scalar/vector/array</TD><TD ></TD></TR>
<TR><TD ><A HREF = "variable.html">variables</A></TD><TD > global scalars, per-atom vectors</TD><TD > global scalar, per-atom vector</TD><TD ></TD></TR>
<TR><TD ><A HREF = "compute_reduce.html">compute reduce</A></TD><TD > per-atom/local vectors</TD><TD > global scalar/vector</TD><TD ></TD></TR>
<TR><TD ><A HREF = "compute_slice.html">compute slice</A></TD><TD > global vectors/arrays</TD><TD > global vector/array</TD><TD ></TD></TR>
<TR><TD ><A HREF = "compute_property_atom.html">compute property/atom</A></TD><TD > per-atom vectors</TD><TD > per-atom vector/array</TD><TD ></TD></TR>
<TR><TD ><A HREF = "compute_property_local.html">compute property/local</A></TD><TD > local vectors</TD><TD > local vector/array</TD><TD ></TD></TR>
<TR><TD ><A HREF = "compute_atom_molecule.html">compute atom/molecule</A></TD><TD > per-atom vectors</TD><TD > global vector/array</TD><TD ></TD></TR>
<TR><TD ><A HREF = "fix_ave_atom.html">fix ave/atom</A></TD><TD > per-atom vectors</TD><TD > per-atom vector/array</TD><TD ></TD></TR>
<TR><TD ><A HREF = "fix_ave_time.html">fix ave/time</A></TD><TD > global scalars/vectors</TD><TD > global scalar/vector/array, file</TD><TD ></TD></TR>
<TR><TD ><A HREF = "fix_ave_spatial.html">fix ave/spatial</A></TD><TD > per-atom vectors</TD><TD > global array, file</TD><TD ></TD></TR>
<TR><TD ><A HREF = "fix_ave_histo.html">fix ave/histo</A></TD><TD > global/per-atom/local scalars and vectors</TD><TD > global array, file</TD><TD ></TD></TR>
<TR><TD ><A HREF = "fix_ave_correlate.html">fix ave/correlate</A></TD><TD > global scalars</TD><TD > global array, file</TD><TD ></TD></TR>
<TR><TD ><A HREF = "fix_store_state.html">fix store/state</A></TD><TD > per-atom vectors</TD><TD > per-atom vector/array</TD><TD ></TD></TR>
<TR><TD >
</TD></TR></TABLE></DIV>
<HR>
<A NAME = "howto_16"></A><H4>6.16 Thermostatting, barostatting, and computing temperature
</H4>
<P>Thermostatting means controlling the temperature of particles in an MD
simulation. Barostatting means controlling the pressure. Since the
pressure includes a kinetic component due to particle velocities, both
these operations require calculation of the temperature. Typically a
target temperature (T) and/or pressure (P) is specified by the user,
and the thermostat or barostat attempts to equilibrate the system to
the requested T and/or P.
</P>
<P>Temperature is computed as kinetic energy divided by some number of
degrees of freedom (and the Boltzmann constant). Since kinetic energy
is a function of particle velocity, there is often a need to
distinguish between a particle's advection velocity (due to some
aggregate motiion of particles) and its thermal velocity. The sum of
the two is the particle's total velocity, but the latter is often what
is wanted to compute a temperature.
</P>
<P>LAMMPS has several options for computing temperatures, any of which
can be used in thermostatting and barostatting. These <A HREF = "compute.html">compute
commands</A> calculate temperature, and the <A HREF = "compute_pressure.html">compute
pressure</A> command calculates pressure.
</P>
<UL><LI><A HREF = "compute_temp.html">compute temp</A>
<LI><A HREF = "compute_temp_sphere.html">compute temp/sphere</A>
<LI><A HREF = "compute_temp_asphere.html">compute temp/asphere</A>
<LI><A HREF = "compute_temp_com.html">compute temp/com</A>
<LI><A HREF = "compute_temp_deform.html">compute temp/deform</A>
<LI><A HREF = "compute_temp_partial.html">compute temp/partial</A>
<LI><A HREF = "compute_temp_profile.html">compute temp/profile</A>
<LI><A HREF = "compute_temp_ramp.html">compute temp/ramp</A>
<LI><A HREF = "compute_temp_region.html">compute temp/region</A>
</UL>
<P>All but the first 3 calculate velocity biases (i.e. advection
velocities) that are removed when computing the thermal temperature.
<A HREF = "compute_temp_sphere.html">Compute temp/sphere</A> and <A HREF = "compute_temp_asphere.html">compute
temp/asphere</A> compute kinetic energy for
extended particles that includes rotational degrees of freedom. They
both allow, as an extra argument, which is another temperature compute
that subtracts a velocity bias. This allows the translational
velocity of extended spherical or aspherical particles to be adjusted
in prescribed ways.
</P>
<P>Thermostatting in LAMMPS is performed by <A HREF = "fix.html">fixes</A>, or in one
case by a pair style. Four thermostatting fixes are currently
available: Nose-Hoover (nvt), Berendsen, Langevin, and direct
rescaling (temp/rescale). Dissipative particle dynamics (DPD)
thermostatting can be invoked via the <I>dpd/tstat</I> pair style:
</P>
<UL><LI><A HREF = "fix_nh.html">fix nvt</A>
<LI><A HREF = "fix_nvt_sphere.html">fix nvt/sphere</A>
<LI><A HREF = "fix_nvt_asphere.html">fix nvt/asphere</A>
<LI><A HREF = "fix_nvt_sllod.html">fix nvt/sllod</A>
<LI><A HREF = "fix_temp_berendsen.html">fix temp/berendsen</A>
<LI><A HREF = "fix_langevin.html">fix langevin</A>
<LI><A HREF = "fix_temp_rescale.html">fix temp/rescale</A>
<LI><A HREF = "pair_dpd.html">pair_style dpd/tstat</A>
</UL>
<P><A HREF = "fix_nh.html">Fix nvt</A> only thermostats the translational velocity of
particles. <A HREF = "fix_nvt_sllod.html">Fix nvt/sllod</A> also does this, except
that it subtracts out a velocity bias due to a deforming box and
integrates the SLLOD equations of motion. See the <A HREF = "#howto_13">NEMD
simulations</A> section of this page for further details. <A HREF = "fix_nvt_sphere.html">Fix
nvt/sphere</A> and <A HREF = "fix_nvt_asphere.html">fix
nvt/asphere</A> thermostat not only translation
velocities but also rotational velocities for spherical and aspherical
particles.
</P>
<P>DPD thermostatting alters pairwise interactions in a manner analagous
to the per-particle thermostatting of <A HREF = "fix_langevin.html">fix
langevin</A>.
</P>
<P>Any of the thermostatting fixes can use temperature computes that
remove bias for two purposes: (a) computing the current temperature to
compare to the requested target temperature, and (b) adjusting only
the thermal temperature component of the particle's velocities. See
the doc pages for the individual fixes and for the
<A HREF = "fix_modify.html">fix_modify</A> command for instructions on how to assign
a temperature compute to a thermostatting fix. For example, you can
apply a thermostat to only the x and z components of velocity by using
it in conjunction with <A HREF = "compute_temp_partial.html">compute
temp/partial</A>.
</P>
<P>IMPORTANT NOTE: Only the nvt fixes perform time integration, meaning
they update the velocities and positions of particles due to forces
and velocities respectively. The other thermostat fixes only adjust
velocities; they do NOT perform time integration updates. Thus they
should be used in conjunction with a constant NVE integration fix such
as these:
</P>
<UL><LI><A HREF = "fix_nve.html">fix nve</A>
<LI><A HREF = "fix_nve_sphere.html">fix nve/sphere</A>
<LI><A HREF = "fix_nve_asphere.html">fix nve/asphere</A>
</UL>
<P>Barostatting in LAMMPS is also performed by <A HREF = "fix.html">fixes</A>. Two
barosttating methods are currently available: Nose-Hoover (npt and
nph) and Berendsen:
</P>
<UL><LI><A HREF = "fix_nh.html">fix npt</A>
<LI><A HREF = "fix_npt_sphere.html">fix npt/sphere</A>
<LI><A HREF = "fix_npt_asphere.html">fix npt/asphere</A>
<LI><A HREF = "fix_nh.html">fix nph</A>
<LI><A HREF = "fix_press_berendsen.html">fix press/berendsen</A>
</UL>
<P>The <A HREF = "fix_nh.html">fix npt</A> commands include a Nose-Hoover thermostat
and barostat. <A HREF = "fix_nh.html">Fix nph</A> is just a Nose/Hoover barostat;
it does no thermostatting. Both <A HREF = "fix_nh.html">fix nph</A> and <A HREF = "fix_press_berendsen.html">fix
press/bernendsen</A> can be used in conjunction
with any of the thermostatting fixes.
</P>
<P>As with the thermostats, <A HREF = "fix_nh.html">fix npt</A> and <A HREF = "fix_nh.html">fix
nph</A> only use translational motion of the particles in
computing T and P and performing thermo/barostatting. <A HREF = "fix_npt_sphere.html">Fix
npt/sphere</A> and <A HREF = "fix_npt_asphere.html">fix
npt/asphere</A> thermo/barostat using not only
translation velocities but also rotational velocities for spherical
and aspherical particles.
</P>
<P>All of the barostatting fixes use the <A HREF = "compute_pressure.html">compute
pressure</A> compute to calculate a current
pressure. By default, this compute is created with a simple <A HREF = "compute_temp.html">compute
temp</A> (see the last argument of the <A HREF = "compute_pressure.html">compute
pressure</A> command), which is used to calculated
the kinetic componenet of the pressure. The barostatting fixes can
also use temperature computes that remove bias for the purpose of
computing the kinetic componenet which contributes to the current
pressure. See the doc pages for the individual fixes and for the
<A HREF = "fix_modify.html">fix_modify</A> command for instructions on how to assign
a temperature or pressure compute to a barostatting fix.
</P>
<P>IMPORTANT NOTE: As with the thermostats, the Nose/Hoover methods (<A HREF = "fix_nh.html">fix
npt</A> and <A HREF = "fix_nh.html">fix nph</A>) perform time
integration. <A HREF = "fix_press_berendsen.html">Fix press/berendsen</A> does NOT,
so it should be used with one of the constant NVE fixes or with one of
the NVT fixes.
</P>
<P>Finally, thermodynamic output, which can be setup via the
<A HREF = "thermo_style.html">thermo_style</A> command, often includes temperature
and pressure values. As explained on the doc page for the
<A HREF = "thermo_style.html">thermo_style</A> command, the default T and P are
setup by the thermo command itself. They are NOT the ones associated
with any thermostatting or barostatting fix you have defined or with
any compute that calculates a temperature or pressure. Thus if you
want to view these values of T and P, you need to specify them
explicitly via a <A HREF = "thermo_style.html">thermo_style custom</A> command. Or
you can use the <A HREF = "thermo_modify.html">thermo_modify</A> command to
re-define what temperature or pressure compute is used for default
thermodynamic output.
</P>
<HR>
<A NAME = "howto_17"></A><H4>6.17 Walls
</H4>
<P>Walls in an MD simulation are typically used to bound particle motion,
i.e. to serve as a boundary condition.
</P>
<P>Walls in LAMMPS can be of rough (made of particles) or idealized
surfaces. Ideal walls can be smooth, generating forces only in the
normal direction, or frictional, generating forces also in the
tangential direction.
</P>
<P>Rough walls, built of particles, can be created in various ways. The
particles themselves can be generated like any other particle, via the
<A HREF = "lattice.html">lattice</A> and <A HREF = "create_atoms.html">create_atoms</A> commands,
or read in via the <A HREF = "read_data.html">read_data</A> command.
</P>
<P>Their motion can be constrained by many different commands, so that
they do not move at all, move together as a group at constant velocity
or in response to a net force acting on them, move in a prescribed
fashion (e.g. rotate around a point), etc. Note that if a time
integration fix like <A HREF = "fix_nve.html">fix nve</A> or <A HREF = "fix_nh.html">fix nvt</A>
is not used with the group that contains wall particles, their
positions and velocities will not be updated.
</P>
<UL><LI><A HREF = "fix_aveforce.html">fix aveforce</A> - set force on particles to average value, so they move together
<LI><A HREF = "fix_setforce.html">fix setforce</A> - set force on particles to a value, e.g. 0.0
<LI><A HREF = "fix_freeze.html">fix freeze</A> - freeze particles for use as granular walls
<LI><A HREF = "fix_nve_noforce.html">fix nve/noforce</A> - advect particles by their velocity, but without force
<LI><A HREF = "fix_move.html">fix move</A> - prescribe motion of particles by a linear velocity, oscillation, rotation, variable
</UL>
<P>The <A HREF = "fix_move.html">fix move</A> command offers the most generality, since
the motion of individual particles can be specified with
<A HREF = "variable.html">variable</A> formula which depends on time and/or the
particle position.
</P>
<P>For rough walls, it may be useful to turn off pairwise interactions
between wall particles via the <A HREF = "neigh_modify.html">neigh_modify
exclude</A> command.
</P>
<P>Rough walls can also be created by specifying frozen particles that do
not move and do not interact with mobile particles, and then tethering
other particles to the fixed particles, via a <A HREF = "bond_style.html">bond</A>.
The bonded particles do interact with other mobile particles.
</P>
<P>Idealized walls can be specified via several fix commands. <A HREF = "fix_wall_gran.html">Fix
wall/gran</A> creates frictional walls for use with
granular particles; all the other commands create smooth walls.
</P>
<UL><LI><A HREF = "fix_wall_reflect.html">fix wall/reflect</A> - reflective flat walls
<LI><A HREF = "fix_wall.html">fix wall/lj93</A> - flat walls, with Lennard-Jones 9/3 potential
<LI><A HREF = "fix_wall.html">fix wall/lj126</A> - flat walls, with Lennard-Jones 12/6 potential
<LI><A HREF = "fix_wall.html">fix wall/colloid</A> - flat walls, with <A HREF = "pair_colloid.html">pair_style colloid</A> potential
<LI><A HREF = "fix_wall.html">fix wall/harmonic</A> - flat walls, with repulsive harmonic spring potential
<LI><A HREF = "fix_wall_region.html">fix wall/region</A> - use region surface as wall
<LI><A HREF = "fix_wall_gran.html">fix wall/gran</A> - flat or curved walls with <A HREF = "pair_gran.html">pair_style granular</A> potential
</UL>
<P>The <I>lj93</I>, <I>lj126</I>, <I>colloid</I>, and <I>harmonic</I> styles all allow the
flat walls to move with a constant velocity, or oscillate in time.
The <A HREF = "fix_wall_region.html">fix wall/region</A> command offers the most
generality, since the region surface is treated as a wall, and the
geometry of the region can be a simple primitive volume (e.g. a
sphere, or cube, or plane), or a complex volume made from the union
and intersection of primitive volumes. <A HREF = "region.html">Regions</A> can also
specify a volume "interior" or "exterior" to the specified primitive
shape or <I>union</I> or <I>intersection</I>. <A HREF = "region.html">Regions</A> can also be
"dynamic" meaning they move with constant velocity, oscillate, or
rotate.
</P>
<P>The only frictional idealized walls currently in LAMMPS are flat or
curved surfaces specified by the <A HREF = "fix_wall_gran.html">fix wall/gran</A>
command. At some point we plan to allow regoin surfaces to be used as
frictional walls, as well as triangulated surfaces.
</P>
<HR>
<A NAME = "howto_18"></A><H4>6.18 Elastic constants
</H4>
<P>Elastic constants characterize the stiffness of a material. The formal
definition is provided by the linear relation that holds between the
stress and strain tensors in the limit of infinitesimal deformation.
In tensor notation, this is expressed as s_ij = C_ijkl * e_kl, where
the repeated indices imply summation. s_ij are the elements of the
symmetric stress tensor. e_kl are the elements of the symmetric strain
tensor. C_ijkl are the elements of the fourth rank tensor of elastic
constants. In three dimensions, this tensor has 3^4=81 elements. Using
Voigt notation, the tensor can be written as a 6x6 matrix, where C_ij
is now the derivative of s_i w.r.t. e_j. Because s_i is itself a
derivative w.r.t. e_i, it follows that C_ij is also symmetric, with at
most 7*6/2 = 21 distinct elements.
</P>
<P>At zero temperature, it is easy to estimate these derivatives by
deforming the cell in one of the six directions using the command
<A HREF = "displace_box.html">displace_box</A> and measuring the change in the
stress tensor. A general-purpose script that does this is given in the
examples/elastic directory described in <A HREF = "Section_example.html">this
section</A>.
</P>
<P>Calculating elastic constants at finite temperature is more
challenging, because it is necessary to run a simulation that perfoms
time averages of differential properties. One way to do this is to
measure the change in average stress tensor in an NVT simulations when
the cell volume undergoes a finite deformation. In order to balance
the systematic and statistical errors in this method, the magnitude of
the deformation must be chosen judiciously, and care must be taken to
fully equilibrate the deformed cell before sampling the stress
tensor. Another approach is to sample the triclinic cell fluctuations
that occur in an NPT simulation. This method can also be slow to
converge and requires careful post-processing <A HREF = "#Shinoda">(Shinoda)</A>
</P>
<HR>
<A NAME = "howto_19"></A><H4>6.19 Library interface to LAMMPS
</H4>
<P>As described in <A HREF = "Section_start.html#start_4">Section_start 4</A>, LAMMPS
can be built as a library, so that it can be called by another code,
used in a <A HREF = "Section_howto.html#howto_10">coupled manner</A> with other
codes, or driven through a <A HREF = "Section_python.html">Python interface</A>.
</P>
<P>All of these methodologies use a C-style interface to LAMMPS that is
provided in the files src/library.cpp and src/library.h. The
functions therein have a C-style argument list, but contain C++ code
you could write yourself in a C++ application that was invoking LAMMPS
directly. The C++ code in the functions illustrates how to invoke
internal LAMMPS operations. Note that LAMMPS classes are defined
within a LAMMPS namespace (LAMMPS_NS) if you use them from another C++
application.
</P>
<P>Library.cpp contains these 4 functions:
</P>
<PRE>void lammps_open(int, char **, MPI_Comm, void **);
void lammps_close(void *);
void lammps_file(void *, char *);
char *lammps_command(void *, char *);
</PRE>
<P>The lammps_open() function is used to initialize LAMMPS, passing in a
list of strings as if they were <A HREF = "Section_start.html#start_6">command-line
arguments</A> when LAMMPS is run in
stand-alone mode from the command line, and a MPI communicator for
LAMMPS to run under. It returns a ptr to the LAMMPS object that is
created, and which is used in subsequent library calls. The
lammps_open() function can be called multiple times, to create
multiple instances of LAMMPS.
</P>
<P>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. Or it might instantiate multiple
instances of LAMMPS to perform different calculations.
</P>
<P>The lammps_close() function is used to shut down an instance of LAMMPS
and free all its memory.
</P>
<P>The lammps_file() and lammps_command() functions are used to pass a
file or string to LAMMPS as if it were an input script or single
command in an input script. Thus the calling code can read or
generate a series of LAMMPS commands one line at a time and pass it
thru the library interface to setup a problem and then run it,
interleaving the lammps_command() calls with other calls to extract
information from LAMMPS, perform its own operations, or call another
code's library.
</P>
<P>Other useful functions are also included in library.cpp. For example:
</P>
<PRE>void *lammps_extract_global(void *, char *)
void *lammps_extract_atom(void *, char *)
void *lammps_extract_compute(void *, char *, int, int)
void *lammps_extract_fix(void *, char *, int, int, int, int)
void *lammps_extract_variable(void *, char *, char *)
int lammps_get_natoms(void *)
void lammps_get_coords(void *, double *)
void lammps_put_coords(void *, double *)
</PRE>
<P>These can extract various global or per-atom quantities from LAMMPS as
well as values calculated by a compute, fix, or variable. The "get"
and "put" operations can retrieve and reset atom coordinates.
See the library.cpp file and its associated header file library.h for
details.
</P>
<P>The key idea of the library interface is that you can write any
functions you wish to define how your code talks to LAMMPS and add
them to src/library.cpp and src/library.h, as well as to the <A HREF = "Section_python.html">Python
interface</A>. The routines you add can access
or change any LAMMPS data you wish. The couple and python directories
have example C++ and C and Python codes which show how a driver code
can link to LAMMPS as a library, run LAMMPS on a subset of processors,
grab data from LAMMPS, change it, and put it back into LAMMPS.
</P>
<HR>
<A NAME = "howto_20"></A><H4>6.20 Calculating thermal conductivity
</H4>
<P>The thermal conductivity kappa of a material can be measured in at
least 3 ways using various options in LAMMPS. (See <A HREF = "Section_howto.html#howto_21">this
section</A> of the manual for an analogous
discussion for viscosity). The thermal conducitivity tensor kappa is
a measure of the propensity of a material to transmit heat energy in a
diffusive manner as given by Fourier's law
</P>
<P>J = -kappa grad(T)
</P>
<P>where J is the heat flux in units of energy per area per time and
grad(T) is the spatial gradient of temperature. The thermal
conductivity thus has units of energy per distance per time per degree
K and is often approximated as an isotropic quantity, i.e. as a
scalar.
</P>
<P>The first method is to setup two thermostatted regions at opposite
ends of a simulation box, or one in the middle and one at the end of a
periodic box. By holding the two regions at different temperatures
with a <A HREF = "Section_howto.html#howto_13">thermostatting fix</A>, the energy added
to the hot region should equal the energy subtracted from the cold
region and be proportional to the heat flux moving between the
regions. See the paper by <A HREF = "#Ikeshoji">Ikeshoji and Hafskjold</A> for
details of this idea. Note that thermostatting fixes such as <A HREF = "fix_nh.html">fix
nvt</A>, <A HREF = "fix_langevin.html">fix langevin</A>, and <A HREF = "fix_temp_rescale.html">fix
temp/rescale</A> store the cumulative energy they
add/subtract. Alternatively, the <A HREF = "fix_heat.html">fix heat</A> command can
used in place of thermostats on each of two regions, and the resulting
temperatures of the two regions monitored with the "compute
temp/region" command or the temperature profile of the intermediate
region monitored with the <A HREF = "fix_ave_spatial.html">fix ave/spatial</A> and
<A HREF = "compute_ke_atom.html">compute ke/atom</A> commands.
</P>
<P>The second method is to perform a reverse non-equilibrium MD
simulation using the <A HREF = "fix_thermal_conductivity.html">fix
thermal/conductivity</A> command which
implements the rNEMD algorithm of Muller-Plathe. Kinetic energy is
swapped between atoms in two different layers of the simulation box.
This induces a temperature gradient between the two layers which can
be monitored with the <A HREF = "fix_ave_spatial.html">fix ave/spatial</A> and
<A HREF = "compute_ke_atom.html">compute ke/atom</A> commands. The fix tallies the
cumulative energy transfer that it performs. See the <A HREF = "fix_thermal_conductivity.html">fix
thermal/conductivity</A> command for
details.
</P>
<P>The third method is based on the Green-Kubo (GK) formula which relates
the ensemble average of the auto-correlation of the heat flux to
kappa. The heat flux can be calculated from the fluctuations of
per-atom potential and kinetic energies and per-atom stress tensor in
a steady-state equilibrated simulation. This is in contrast to the
two preceding non-equilibrium methods, where energy flows continuously
between hot and cold regions of the simulation box.
</P>
<P>The <A HREF = "compute_heat_flux.html">compute heat/flux</A> command can calculate
the needed heat flux and describes how to implement the Green_Kubo
formalism using additional LAMMPS commands, such as the <A HREF = "fix_ave_correlate.html">fix
ave/correlate</A> command to calculate the needed
auto-correlation. See the doc page for the <A HREF = "compute_heat_flux.html">compute
heat/flux</A> command for an example input script
that calculates the thermal conductivity of solid Ar via the GK
formalism.
</P>
<HR>
<A NAME = "howto_21"></A><H4>6.21 Calculating viscosity
</H4>
<P>The shear viscosity eta of a fluid can be measured in at least 3 ways
using various options in LAMMPS. (See <A HREF = "Section_howto.html#howto_20">this
section</A> of the manual for an analogous
discussion for thermal conductivity). Eta is a measure of the
propensity of a fluid to transmit momentum in a direction
perpendicular to the direction of velocity or momentum flow.
Alternatively it is the resistance the fluid has to being sheared. It
is given by
</P>
<P>J = -eta grad(Vstream)
</P>
<P>where J is the momentum flux in units of momentum per area per time.
and grad(Vstream) is the spatial gradient of the velocity of the fluid
moving in another direction, normal to the area through which the
momentum flows. Viscosity thus has units of pressure-time.
</P>
<P>The first method is to perform a non-equlibrium MD (NEMD) simulation
by shearing the simulation box via the <A HREF = "fix_deform.html">fix deform</A>
command, and using the <A HREF = "fix_nvt_sllod.html">fix nvt/sllod</A> command to
thermostat the fluid via the SLLOD equations of motion. The velocity
profile setup in the fluid by this procedure can be monitored by the
<A HREF = "fix_ave_spatial.html">fix ave/spatial</A> command, which determines
grad(Vstream) in the equation above. E.g. the derivative in the
y-direction of the Vx component of fluid motion or grad(Vstream) =
dVx/dy. In this case, the Pxy off-diagonal component of the pressure
or stress tensor, as calculated by the <A HREF = "compute_pressure.html">compute
pressure</A> command, can also be monitored, which
is the J term in the equation above. See <A HREF = "Section_howto.html#howto_13">this
section</A> of the manual for details on NEMD
simulations.
</P>
<P>The second method is to perform a reverse non-equilibrium MD
simulation using the <A HREF = "fix_viscosity.html">fix viscosity</A> command which
implements the rNEMD algorithm of Muller-Plathe. Momentum in one
dimension is swapped between atoms in two different layers of the
simulation box in a different dimension. This induces a velocity
gradient which can be monitored with the <A HREF = "fix_ave_spatial.html">fix
ave/spatial</A> command. The fix tallies the
cummulative momentum transfer that it performs. See the <A HREF = "fix_viscosity.html">fix
viscosity</A> command for details.
</P>
<P>The third method is based on the Green-Kubo (GK) formula which relates
the ensemble average of the auto-correlation of the stress/pressure
tensor to eta. This can be done in a steady-state equilibrated
simulation which is in contrast to the two preceding non-equilibrium
methods, where momentum flows continuously through the simulation box.
</P>
<P>Here is an example input script that calculates the viscosity of
liquid Ar via the GK formalism:
</P>
<PRE># Sample LAMMPS input script for viscosity of liquid Ar
</PRE>
<PRE>units real
variable T equal 86.4956
variable V equal vol
variable dt equal 4.0
variable p equal 400 # correlation length
variable s equal 5 # sample interval
variable d equal $p*$s # dump interval
</PRE>
<PRE># convert from LAMMPS real units to SI
</PRE>
<PRE>variable kB equal 1.3806504e-23 # [J/K/</B> Boltzmann
variable atm2Pa equal 101325.0
variable A2m equal 1.0e-10
variable fs2s equal 1.0e-15
variable convert equal ${atm2Pa}*${atm2Pa}*${fs2s}*${A2m}*${A2m}*${A2m}
</PRE>
<PRE># setup problem
</PRE>
<PRE>dimension 3
boundary p p p
lattice fcc 5.376 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
region box block 0 4 0 4 0 4
create_box 1 box
create_atoms 1 box
mass 1 39.948
pair_style lj/cut 13.0
pair_coeff * * 0.2381 3.405
timestep ${dt}
thermo $d
</PRE>
<PRE># equilibration and thermalization
</PRE>
<PRE>velocity all create $T 102486 mom yes rot yes dist gaussian
fix NVT all nvt temp $T $T 10 drag 0.2
run 8000
</PRE>
<PRE># viscosity calculation, switch to NVE if desired
</PRE>
<PRE>#unfix NVT
#fix NVE all nve
</PRE>
<PRE>reset_timestep 0
variable pxy equal pxy
variable pxz equal pxz
variable pyz equal pyz
fix SS all ave/correlate $s $p $d &
v_pxy v_pxz v_pyz type auto file S0St.dat ave running
variable scale equal ${convert}/(${kB}*$T)*$V*$s*${dt}
variable v11 equal trap(f_SS[3/</B>)*${scale}
variable v22 equal trap(f_SS[4/</B>)*${scale}
variable v33 equal trap(f_SS[5/</B>)*${scale}
thermo_style custom step temp press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33
run 100000
variable v equal (v_v11+v_v22+v_v33)/3.0
variable ndens equal count(all)/vol
print "average viscosity: $v [Pa.s/</B> @ $T K, ${ndens} /A^3"
</PRE>
<HR>
<HR>
<A NAME = "Berendsen"></A>
<P><B>(Berendsen)</B> Berendsen, Grigera, Straatsma, J Phys Chem, 91,
6269-6271 (1987).
</P>
<A NAME = "Cornell"></A>
<P><B>(Cornell)</B> Cornell, Cieplak, Bayly, Gould, Merz, Ferguson,
Spellmeyer, Fox, Caldwell, Kollman, JACS 117, 5179-5197 (1995).
</P>
<A NAME = "Horn"></A>
<P><B>(Horn)</B> Horn, Swope, Pitera, Madura, Dick, Hura, and Head-Gordon,
J Chem Phys, 120, 9665 (2004).
</P>
<A NAME = "Ikeshoji"></A>
<P><B>(Ikeshoji)</B> Ikeshoji and Hafskjold, Molecular Physics, 81, 251-261
(1994).
</P>
<A NAME = "MacKerell"></A>
<P><B>(MacKerell)</B> MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field,
Fischer, Gao, Guo, Ha, et al, J Phys Chem, 102, 3586 (1998).
</P>
<A NAME = "Mayo"></A>
<P><B>(Mayo)</B> Mayo, Olfason, Goddard III, J Phys Chem, 94, 8897-8909
(1990).
</P>
<A NAME = "Jorgensen"></A>
<P><B>(Jorgensen)</B> Jorgensen, Chandrasekhar, Madura, Impey, Klein, J Chem
Phys, 79, 926 (1983).
</P>
<A NAME = "Price"></A>
<P><B>(Price)</B> Price and Brooks, J Chem Phys, 121, 10096 (2004).
</P>
<A NAME = "Shinoda"></A>
<P><B>(Shinoda)</B> Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004).
</P>
</HTML>
diff --git a/doc/Section_howto.txt b/doc/Section_howto.txt
index 50c24d8b8..923153d81 100644
--- a/doc/Section_howto.txt
+++ b/doc/Section_howto.txt
@@ -1,1938 +1,1942 @@
"Previous Section"_Section_accelerate.html - "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc - "Next Section"_Section_example.html :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
6. How-to discussions :h3
This section describes how to perform common tasks using LAMMPS.
6.1 "Restarting a simulation"_#howto_1
6.2 "2d simulations"_#howto_2
6.3 "CHARMM, AMBER, and DREIDING force fields"_#howto_3
6.4 "Running multiple simulations from one input script"_#howto_4
6.5 "Multi-replica simulations"_#howto_5
6.6 "Granular models"_#howto_6
6.7 "TIP3P water model"_#howto_7
6.8 "TIP4P water model"_#howto_8
6.9 "SPC water model"_#howto_9
6.10 "Coupling LAMMPS to other codes"_#howto_10
6.11 "Visualizing LAMMPS snapshots"_#howto_11
6.12 "Triclinic (non-orthogonal) simulation boxes"_#howto_12
6.13 "NEMD simulations"_#howto_13
6.14 "Extended spherical and aspherical particles"_#howto_14
6.15 "Output from LAMMPS (thermo, dumps, computes, fixes, variables)"_#howto_15
6.16 "Thermostatting, barostatting and computing temperature"_#howto_16
6.17 "Walls"_#howto_17
6.18 "Elastic constants"_#howto_18
6.19 "Library interface to LAMMPS"_#howto_19
6.20 "Calculating thermal conductivity"_#howto_20
6.21 "Calculating viscosity"_#howto_21 :all(b)
The example input scripts included in the LAMMPS distribution and
highlighted in "Section_example"_Section_example.html also show how to
setup and run various kinds of simulations.
:line
:line
6.1 Restarting a simulation :link(howto_1),h4
There are 3 ways to continue a long LAMMPS simulation. Multiple
"run"_run.html 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"_restart.html
command. At a later time, these binary files can be read via a
"read_restart"_read_restart.html command in a new script. Or they can
be converted to text data files and read by a
"read_data"_read_data.html command in a new script. "This
section"_Section_tools.html 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"_read_restart.html and
"read_data"_read_data.html 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 :pre
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 :pre
neighbor 0.4 bin
neigh_modify every 1 delay 1 :pre
fix 1 all nve
fix 2 all langevin 1.0 1.0 10.0 904297 :pre
timestep 0.012 :pre
run 50 :pre
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"_fix_langevin.html 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 :pre
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 :pre
read_data tmp.restart.data :pre
neighbor 0.4 bin
neigh_modify every 1 delay 1 :pre
fix 1 all nve
fix 2 all langevin 1.0 1.0 10.0 904297 :pre
timestep 0.012 :pre
reset_timestep 50
run 50 :pre
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"_reset_timestep.html command is used to tell
LAMMPS the current timestep. This value is stored in restart files,
but not in data files.
:line
6.2 2d simulations :link(howto_2),h4
Use the "dimension"_dimension.html command to specify a 2d simulation.
Make the simulation box periodic in z via the "boundary"_boundary.html
command. This is the default.
If using the "create box"_create_box.html 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"_create_box.html 1 -10 10 -10 10 -0.25 0.25 :pre
If using the "read data"_read_data.html 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"_fix_enforce2d.html 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.
IMPORTANT NOTE: Some models in LAMMPS treat particles as extended
spheres, as opposed to point particles. In 2d, the particles will
still be spheres, not disks, meaning their moment of inertia will be
the same as in 3d.
:line
6.3 CHARMM, AMBER, and DREIDING force fields :link(howto_3),h4
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 that correspond to formulas commonly
used in the CHARMM, AMBER, and DREIDING force fields. Setting
coefficients is done in the input data file via the
"read_data"_read_data.html command or in the input script with
commands like "pair_coeff"_pair_coeff.html or
"bond_coeff"_bond_coeff.html. See "Section_tools"_Section_tools.html
for additional tools that can use CHARMM or AMBER to assign force
field coefficients and convert their output into LAMMPS input.
See "(MacKerell)"_#MacKerell for a description of the CHARMM force
field. See "(Cornell)"_#Cornell for a description of the AMBER force
field.
:link(charmm,http://www.scripps.edu/brooks)
:link(amber,http://amber.scripps.edu)
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"_bond_harmonic.html harmonic
"angle_style"_angle_charmm.html charmm
"dihedral_style"_dihedral_charmm.html charmm
"pair_style"_pair_charmm.html lj/charmm/coul/charmm
"pair_style"_pair_charmm.html lj/charmm/coul/charmm/implicit
"pair_style"_pair_charmm.html lj/charmm/coul/long :ul
"special_bonds"_special_bonds.html charmm
"special_bonds"_special_bonds.html amber :ul
DREIDING is a generic force field developed by the "Goddard
group"_http://www.wag.caltech.edu at Caltech and is useful for
predicting structures and dynamics of organic, biological and
main-group inorganic molecules. The philosophy in DREIDING is to use
general force constants and geometry parameters based on simple
hybridization considerations, rather than individual force constants
and geometric parameters that depend on the particular combinations of
atoms involved in the bond, angle, or torsion terms. DREIDING has an
"explicit hydrogen bond term"_pair_hbond_dreiding.html to describe
interactions involving a hydrogen atom on very electronegative atoms
(N, O, F).
See "(Mayo)"_#Mayo for a description of the DREIDING force field
These style choices compute force field formulas that are consistent
with the DREIDING force field. See each command's
documentation for the formula it computes.
"bond_style"_bond_harmonic.html harmonic
"bond_style"_bond_morse.html morse :ul
"angle_style"_angle_harmonic.html harmonic
"angle_style"_angle_cosine.html cosine
"angle_style"_angle_cosine_periodic.html cosine/periodic :ul
"dihedral_style"_dihedral_charmm.html charmm
"improper_style"_improper_umbrella.html umbrella :ul
"pair_style"_pair_buck.html buck
"pair_style"_pair_buck.html buck/coul/cut
"pair_style"_pair_buck.html buck/coul/long
"pair_style"_pair_lj.html lj/cut
"pair_style"_pair_lj.html lj/cut/coul/cut
"pair_style"_pair_lj.html lj/cut/coul/long :ul
"pair_style"_pair_hbond_dreiding.html hbond/dreiding/lj
"pair_style"_pair_hbond_dreiding.html hbond/dreiding/morse :ul
"special_bonds"_special_bonds.html dreiding :ul
:line
6.4 Running multiple simulations from one input script :link(howto_4),h4
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"_run.html 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 :pre
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"_clear.html 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 :pre
would run 2 independent simulations, one after the other.
For large numbers of independent simulations, you can use
"variables"_variable.html and the "next"_next.html and
"jump"_jump.html 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
shell cd $d
read_data data.polymer
run 10000
shell cd ..
clear
next d
jump in.polymer :pre
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 :pre
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"_Section_start.html#start_6 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"_variable.html 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.
:line
6.5 Multi-replica simulations :link(howto_5),h4
Several commands in LAMMPS run mutli-replica simulations, meaning
that multiple instances (replicas) of your simulation are run
simultaneously, with small amounts of data exchanged between replicas
periodically.
These are the relevant commands:
"neb"_neb.html for nudged elastic band calculations
"prd"_prd.html for parallel replica dynamics
"tad"_tad.html for temperature accelerated dynamics
"temper"_temper.html for parallel tempering :ul
NEB is a method for finding transition states and barrier energies.
PRD and TAD are methods for performing accelerated dynamics to find
and perform infrequent events. Parallel tempering or replica exchange
runs different replicas at a series of temperature to facilitate
rare-event sampling.
These command can only be used if LAMMPS was built with the "replica"
package. See the "Making LAMMPS"_Section_start.html#start_3 section
for more info on packages.
In all these cases, you must run with one or more processors per
replica. The processors assigned to each replica are determined at
run-time by using the "-partition command-line
switch"_Section_start.html#start_6 to launch LAMMPS on multiple
partitions, which in this context are the same as replicas. E.g.
these commands:
mpirun -np 16 lmp_linux -partition 8x2 -in in.temper
mpirun -np 8 lmp_linux -partition 8x1 -in in.neb :pre
would each run 8 replicas, on either 16 or 8 processors. Note the use
of the "-in command-line switch"_Section_start.html#start_6 to specify
the input script which is required when running in multi-replica mode.
Also note that with MPI installed on a machine (e.g. your desktop),
you can run on more (virtual) processors than you have physical
processors. Thus the above commands could be run on a
single-processor (or few-processor) desktop so that you can run
a multi-replica simulation on more replicas than you have
physical processors.
:line
6.6 Granular models :link(howto_6),h4
Granular system are composed of spherical particles with a diameter,
as opposed to point particles. This means they have an angular
velocity and torque can be imparted to them to cause them to rotate.
To run a simulation of a granular model, you will want to use
the following commands:
"atom_style sphere"_atom_style.html
"fix nve/sphere"_fix_nve_sphere.html
"fix gravity"_fix_gravity.html :ul
This compute
"compute erotate/sphere"_compute_erotate_sphere.html :ul
calculates rotational kinetic energy which can be "output with
thermodynamic info"_Section_howto.html#howto_15.
Use one of these 3 pair potentials, which compute forces and torques
between interacting pairs of particles:
"pair_style"_pair_style.html gran/history
"pair_style"_pair_style.html gran/no_history
"pair_style"_pair_style.html gran/hertzian :ul
These commands implement fix options specific to granular systems:
"fix freeze"_fix_freeze.html
"fix pour"_fix_pour.html
"fix viscous"_fix_viscous.html
"fix wall/gran"_fix_wall_gran.html :ul
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"_neigh_modify.html exclude :ul
:line
6.7 TIP3P water model :link(howto_7),h4
The TIP3P water model as implemented in CHARMM
"(MacKerell)"_#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"_fix_shake.html 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)"_#Jorgensen.
O mass = 15.9994
H mass = 1.008 :all(b),p
O charge = -0.834
H charge = 0.417 :all(b),p
LJ epsilon of OO = 0.1521
LJ sigma of OO = 3.1507
LJ epsilon of HH = 0.0460
LJ sigma of HH = 0.4000
LJ epsilon of OH = 0.0836
LJ sigma of OH = 1.7753 :all(b),p
K of OH bond = 450
r0 of OH bond = 0.9572 :all(b),p
K of HOH angle = 55
theta of HOH angle = 104.52 :all(b),p
These are the parameters to use for TIP3P with a long-range Coulombic
solver (Ewald or PPPM in LAMMPS), see "(Price)"_#Price for details:
O mass = 15.9994
H mass = 1.008 :all(b),p
O charge = -0.830
H charge = 0.415 :all(b),p
LJ epsilon of OO = 0.102
LJ sigma of OO = 3.188
LJ epsilon, sigma of OH, HH = 0.0 :all(b),p
K of OH bond = 450
r0 of OH bond = 0.9572 :all(b),p
K of HOH angle = 55
theta of HOH angle = 104.52 :all(b),p
Wikipedia also has a nice article on "water
models"_http://en.wikipedia.org/wiki/Water_model.
:line
6.8 TIP4P water model :link(howto_8),h4
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"_pair_lj.html. A cutoff version may be added
-the future. For both models, the bond lengths and bond angles should
-be held fixed using the "fix shake"_fix_shake.html command.
+A TIP4P model is run with LAMMPS using two commands:
+
+"pair_style lj/cut/coul/long/tip4p"_pair_lj.html
+"kspace_style pppm/tip4p"_kspace_style.html :ul
+
+Note that only a TIP4P model with long-range Coulombics is currently
+implemented. A cutoff version may be added in the future. for both
+models, the bond lengths and bond angles should be held fixed using
+the "fix shake"_fix_shake.html 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)"_#Jorgensen. Note that the OM distance is specified in
the "pair_style"_pair_style.html command, not as part of the pair
coefficients.
O mass = 15.9994
H mass = 1.008 :all(b),p
O charge = -1.040
H charge = 0.520 :all(b),p
r0 of OH bond = 0.9572
theta of HOH angle = 104.52 :all(b),p
OM distance = 0.15 :all(b),p
LJ epsilon of O-O = 0.1550
LJ sigma of O-O = 3.1536
LJ epsilon, sigma of OH, HH = 0.0 :all(b),p
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 :all(b),p
O charge = -1.0484
H charge = 0.5242 :all(b),p
r0 of OH bond = 0.9572
theta of HOH angle = 104.52 :all(b),p
OM distance = 0.1250 :all(b),p
LJ epsilon of O-O = 0.16275
LJ sigma of O-O = 3.16435
LJ epsilon, sigma of OH, HH = 0.0 :all(b),p
Wikipedia also has a nice article on "water
models"_http://en.wikipedia.org/wiki/Water_model.
:line
6.9 SPC water model :link(howto_9),h4
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"_fix_shake.html 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.
O mass = 15.9994
H mass = 1.008 :all(b),p
O charge = -0.820
H charge = 0.410 :all(b),p
LJ epsilon of OO = 0.1553
LJ sigma of OO = 3.166
LJ epsilon, sigma of OH, HH = 0.0 :all(b),p
r0 of OH bond = 1.0
theta of HOH angle = 109.47 :all(b),p
Note that as originally proposed, the SPC model was run with a 9
Angstrom cutoff for both LJ and Coulommbic terms. It can also be used
with long-range Coulombics (Ewald or PPPM in LAMMPS), without changing
any of the parameters above, though it becomes a different model in
that mode of usage.
The SPC/E (extended) water model is the same, except
the partial charge assignemnts change:
O charge = -0.8476
H charge = 0.4238 :all(b),p
See the "(Berendsen)"_#Berendsen reference for more details on both
the SPC and SPC/E models.
Wikipedia also has a nice article on "water
models"_http://en.wikipedia.org/wiki/Water_model.
:line
6.10 Coupling LAMMPS to other codes :link(howto_10),h4
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"_fix.html 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"_poems package that performs constrained rigid-body motion on
groups of atoms is hooked to LAMMPS. See the
"fix_poems"_fix_poems.html command for more details. See "this
section"_Section_modify.html of the documentation for info on how to add
a new fix to LAMMPS.
:link(poems,http://www.rpi.edu/~anderk5/lab)
(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"_run.html 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 "Section_modify"_Section_modify.html 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"_run.html 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.
Examples of driver codes that call LAMMPS as a library are included in
the "couple" directory of the LAMMPS distribution; see couple/README
for more details:
simple: simple driver programs in C++ and C which invoke LAMMPS as a
library :ulb,l
lammps_quest: coupling of LAMMPS and "Quest"_quest, to run classical
MD with quantum forces calculated by a density functional code :l
lammps_spparks: coupling of LAMMPS and "SPPARKS"_spparks, to couple
a kinetic Monte Carlo model for grain growth using MD to calculate
strain induced across grain boundaries :l,ule
:link(quest,http://dft.sandia.gov/Quest)
:link(spparks,http://www.sandia.gov/~sjplimp/spparks.html)
"This section"_Section_start.html#start_4 of the documentation
describes how to build LAMMPS as a library. Once this is done, you
can interface with LAMMPS either via C++, C, Fortran, or Python (or
any other language that supports a vanilla C-like interface). 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. See
"Section_python"_Section_python.html of the manual for a description
of the Python wrapper provided with LAMMPS that operates through the
LAMMPS library interface.
The files src/library.cpp and library.h contain the C-style interface
to LAMMPS. See "Section_howto 19"_Section_howto.html#howto_19 of the
manual for a description of the interface and how to extend it for
your needs.
Note that the lammps_open() function that creates an instance of
LAMMPS takes an MPI communicator as an argument. This means that
instance of LAMMPS will run on the set of processors in the
communicator. Thus 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. Or it might instantiate multiple
instances of LAMMPS to perform different calculations.
:line
6.11 Visualizing LAMMPS snapshots :link(howto_11),h4
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"_dump.html 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"_Section_tools.html#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 geometry
and atom trajectories.
Several programs included with LAMMPS as auxiliary tools can convert
native LAMMPS dump files to other formats. See the
"Section_tools"_Section_tools.html doc page for details. The first is
the "ch2lmp tool"_Section_tools.html#charmm, which contains a
lammps2pdb Perl script which converts LAMMPS dump files into PDB
files. The second is the "lmp2arc tool"_Section_tools.html#arc which
converts LAMMPS dump files into Accelrys' Insight MD program files.
The third is the "lmp2cfg tool"_Section_tools.html#cfg which converts
LAMMPS dump files into CFG files which can be read into the
"AtomEye"_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"_pizza for details. Specifically, Pizza.py can convert
LAMMPS dump files into PDB, XYZ, "Ensight"_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"_vmd.
LAMMPS can create DCD files directly (via "dump dcd") which can be
read by "VMD"_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"_dump.html 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"_vmd for visualization.
See the "dump"_dump.html command for more information on XTC files.
:link(pizza,http://www.sandia.gov/~sjplimp/pizza.html)
:link(vmd,http://www.ks.uiuc.edu/Research/vmd)
:link(ensight,http://www.ensight.com)
:link(atomeye,http://mt.seas.upenn.edu/Archive/Graphics/A)
:line
6.12 Triclinic (non-orthogonal) simulation boxes :link(howto_12),h4
By default, LAMMPS uses an orthogonal simulation box to encompass the
particles. The "boundary"_boundary.html command sets the boundary
conditions of the box (periodic, non-periodic, etc). The orthogonal
box has its "origin" at (xlo,ylo,zlo) and is defined by 3 edge vectors
starting from the origin given by [a] = (xhi-xlo,0,0); [b] =
(0,yhi-ylo,0); [c] = (0,0,zhi-zlo). The 6 parameters
(xlo,xhi,ylo,yhi,zlo,zhi) are defined at the time the simluation box
is created, e.g. by the "create_box"_create_box.html or
"read_data"_read_data.html or "read_restart"_read_restart.html
commands. Additionally, LAMMPS defines box size parameters lx,ly,lz
where lx = xhi-xlo, and similarly in the y and z dimensions. The 6
parameters, as well as lx,ly,lz, can be output via the "thermo_style
custom"_thermo_style.html command.
LAMMPS also allows simulations to be perfored in non-orthogonal
simulation boxes shaped as a parallelepiped with triclinic symmetry.
The parallelepiped has its "origin" at (xlo,ylo,zlo) and is defined by
3 edge vectors starting from the origin given by [a] = (xhi-xlo,0,0); [b]
= (xy,yhi-ylo,0); [c] = (xz,yz,zhi-zlo). {Xy,xz,yz} can be 0.0 or
positive or negative values and are called "tilt factors" because they
are the amount of displacement applied to faces of an originally
orthogonal box to transform it into the parallelepiped. Note that in
LAMMPS the triclinic simulation box edge vectors [a], [b], and [c] cannot be
arbitrary vectors. As indicated, [a] must be aligned with the x axis, [b]
must be in the xy plane, and [c] is arbitrary. However, this is not a
restriction since it is possible to rotate any set of 3 crystal basis
vectors so that they meet this restriction.
The 9 parameters (xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz) are defined at the
time the simluation box is created. This happens in one of 3 ways.
If the "create_box"_create_box.html command is used with a region of
style {prism}, then a triclinic box is setup. See the
"region"_region.html command for details. If the
"read_data"_read_data.html 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 triclinic box is setup. See the
"read_data"_read_data.html command for details. Finally, if the
"read_restart"_read_restart.html command reads a restart file which
was written from a simulation using a triclinic box, then a triclinic
box will be setup for the restarted simulation.
Note that you can define a triclinic box with all 3 tilt factors =
0.0, so that it is initially orthogonal. This is necessary if the box
will become non-orthogonal, e.g. due to the "fix npt"_fix_nh.html or
"fix deform"_fix_deform.html commands. Alternatively, you can use the
"change_box"_change_box.html command to convert a simulation box from
orthogonal to triclinic and vice versa.
As with orthogonal boxes, LAMMPS defines triclinic box size parameters
lx,ly,lz where lx = xhi-xlo, and similarly in the y and z dimensions.
The 9 parameters, as well as lx,ly,lz, can be output via the
"thermo_style custom"_thermo_style.html command.
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 geometrically all
equivalent.
Triclinic crystal structures are often defined using three lattice
constants {a}, {b}, and {c}, and three angles {alpha}, {beta} and
{gamma}. Note that in this nomenclature, the a, b, and c lattice constants
are the scalar lengths of the edge vectors [a], [b], and [c] defined
above. The
relationship between these 6 quantities (a,b,c,alpha,beta,gamma) and
the LAMMPS box sizes (lx,ly,lz) = (xhi-xlo,yhi-ylo,zhi-zlo) and tilt
factors (xy,xz,yz) is as follows:
:c,image(Eqs/box.jpg)
The inverse relationship can be written as follows:
:c,image(Eqs/box_inverse.jpg)
The values of {a}, {b}, {c} , {alpha}, {beta} , and {gamma} can be printed
out or accessed by computes using the
"thermo_style custom"_thermo_style.html keywords
{cella}, {cellb}, {cellc}, {cellalpha}, {cellbeta}, {cellgamma},
respectively.
As discussed on the "dump"_dump.html command doc page, when the BOX
BOUNDS for a snapshot is written to a dump file for a triclinic box,
an orthogonal bounding box which encloses the triclinic simulation box
is output, along with the 3 tilt factors (xy, xz, yz) of the triclinic
box, formatted as follows:
ITEM: BOX BOUNDS xy xz yz
xlo_bound xhi_bound xy
ylo_bound yhi_bound xz
zlo_bound zhi_bound yz :pre
This bounding box is convenient for many visualization programs and is
calculated from the 9 triclinic box parameters
(xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz) as follows:
xlo_bound = xlo + MIN(0.0,xy,xz,xy+xz)
xhi_bound = xhi + MAX(0.0,xy,xz,xy+xz)
ylo_bound = ylo + MIN(0.0,yz)
yhi_bound = yhi + MAX(0.0,yz)
zlo_bound = zlo
zhi_bound = zhi :pre
These formulas can be inverted if you need to convert the bounding box
back into the triclinic box parameters, e.g. xlo = xlo_bound -
MIN(0.0,xy,xz,xy+xz).
One use of triclinic simulation boxes is to model solid-state crystals
with triclinic symmetry. The "lattice"_lattice.html command can be
used with non-orthogonal basis vectors to define a lattice that will
tile a triclinic simulation box via the
"create_atoms"_create_atoms.html command.
A second use is to run Parinello-Rahman dyanamics via the "fix
npt"_fix_nh.html command, which will adjust the xy, xz, yz tilt
factors to compensate for off-diagonal components of the pressure
tensor. The analalog for an "energy minimization"_minimize.html is
the "fix box/relax"_fix_box_relax.html command.
A third use is to shear a bulk solid to study the response of the
material. The "fix deform"_fix_deform.html command can be used for
this purpose. It allows dynamic control of the xy, xz, yz tilt
factors as a simulation runs. This is discussed in the next section
on non-equilibrium MD (NEMD) simulations.
:line
6.13 NEMD simulations :link(howto_13),h4
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 preceding Howto section).
A shear strain can be applied to the simulation box at a desired
strain rate by using the "fix deform"_fix_deform.html command. The
"fix nvt/sllod"_fix_nvt_sllod.html 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"_compute_temp_deform.html 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"_fix_ave_spatial.html 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"_fix_deform.html can continuously strain
a box by an arbitrary amount. As discussed in the "fix
deform"_fix_deform.html 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"_fix_deform.html should be set to "remap v", since that is what
"fix nvt/sllod"_fix_nvt_sllod.html 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"_fix_viscosity.html command.
:line
6.14 Extended spherical and aspherical particles :link(howto_14),h4
Typical MD models treat atoms or particles as point masses.
Sometimes, however, it is desirable to have a model with finite-size
particles such as spheres or aspherical ellipsoids. The difference is
that such particles have a moment of inertia, rotational energy, and
angular momentum. Rotation is induced by torque from interactions
with other particles.
LAMMPS has several options for running simulations with these kinds of
particles. The following aspects are discussed in turn:
atom styles
pair potentials
time integration
computes, thermodynamics, and dump output
rigid bodies composed of extended particles :ul
Atom styles :h5
There are 2 "atom styles"_atom_style.html that allow for definition of
finite-size particles: sphere and ellipsoid. The peri atom style also
treats particles as having a volume, but that is internal to the
"pair_style peri"_pair_peri.html potentials. The dipole atom style is
most often used in conjunction with finite-size particles.
The sphere style defines particles that are spheriods and each
particle can have a unique diameter and mass (or density). These
particles store an angular velocity (omega) and can be acted upon by
torque. The "set" command can be used to modify the diameter and mass
of individual particles, after then are created.
The ellipsoid style defines particles that are ellipsoids and thus can
be aspherical. Each particle has a shape, specified by 3 diameters,
and mass (or density). These particles store an angular momentum and
their orientation (quaternion), and can be acted upon by torque. They
do not store an angular velocity (omega), which can be in a different
direction than angular momentum, rather they compute it as needed.
The "set" command can be used to modify the diameter, orientation, and
mass of individual particles, after then are created. It also has a
brief explanation of what quaternions are.
The dipole style does not define extended particles, but is often
used in conjunction with spherical particles, via a command like
atom_style hybrid sphere dipole :pre
This is because when dipoles interact with each other, they induce
torques, and a particle must be extended (i.e. have a moment of
inertia) in order to respond and rotate. See the "atom_style
dipole"_atom_style.html command for details. The "set" command can be
used to modify the orientation and length of the dipole moment of
individual particles, after then are created.
Note that if one of these atom styles is used (or multiple styles via
the "atom_style hybrid"_atom_style.html command), not all particles in
the system are required to be finite-size or aspherical. For example,
if the 3 shape parameters are set to the same value, the particle will
be a sphere rather than an ellipsoid. If the 3 shape parameters are
all set to 0.0 or if the diameter is set to 0.0, it will be a point
particle. If the length of the dipole moment is set to zero, the
particle will not have a point dipole associated with it. The pair
styles used to compute pairwise interactions will typically compute
the correct interaction in these simplified (cheaper) cases.
"Pair_style hybrid"_pair_hybrid.html can be used to insure the correct
interactions are computed for the appropriate style of interactions.
Likewise, using groups to partition particles (ellipsoids versus
spheres versus point particles) will allow you to use the appropriate
time integrators and temperature computations for each class of
particles. See the doc pages for various commands for details.
Also note that for "2d simulations"_dimension.html, finite-size
spheres and ellipsoids are still treated as 3d particles, rather than
as circular disks or ellipses. This means they have the same moment
of inertia for a 3d extended object. When their temperature is
coomputed, the correct degrees of freedom are used for rotation in a
2d versus 3d system.
Pair potentials :h5
When a system with extended particles is defined, the particles will
only rotate and experience torque if the force field computes such
interactions. These are the various "pair styles"_pair_style.html
that generate torque:
"pair_style gran/history"_pair_gran.html
"pair_style gran/hertzian"_pair_gran.html
"pair_style gran/no_history"_pair_gran.html
"pair_style dipole/cut"_pair_dipole.html
"pair_style gayberne"_pair_gayberne.html
"pair_style resquared"_pair_resquared.html
"pair_style lubricate"_pair_lubricate.html :ul
The "granular pair styles"_pair_gran.html are used with spherical
particles. The "dipole pair style"_pair_dipole.html is used with
"atom_style dipole"_atom_style.html, which could be applied to
spherical or ellipsoidal particles. The "GayBerne"_pair_gayberne.html
and "REsquared"_pair_resquared.html potentials require ellipsoidal
particles, though they will also work if the 3 shape parameters are
the same (a sphere). The "lubrication potential"_pair_lubricate.html
works with spherical particles.
Time integration :h5
There are 3 fixes that perform time integration on extended spherical
particles, meaning the integrators update the rotational orientation
and angular velocity or angular momentum of the particles:
"fix nve/sphere"_fix_nve_sphere.html
"fix nvt/sphere"_fix_nvt_sphere.html
"fix npt/sphere"_fix_npt_sphere.html :ul
Likewise, there are 3 fixes that perform time integration on
ellipsoids as extended aspherical particles:
"fix nve/asphere"_fix_nve_asphere.html
"fix nvt/asphere"_fix_nvt_asphere.html
"fix npt/asphere"_fix_npt_asphere.html :ul
The advantage of these fixes is that those which thermostat the
particles include the rotational degrees of freedom in the temperature
calculation and thermostatting. Other thermostats can be used with
fix nve/sphere or fix nve/asphere, such as fix langevin or fix
temp/berendsen, but those thermostats only operate on the
translational kinetic energy of the extended particles.
Note that for mixtures of point and extended particles, you should
only use these integration fixes on "groups"_group.html which contain
extended particles.
Computes, thermodynamics, and dump output :h5
There are 4 computes that calculate the temperature or rotational energy
of extended spherical or aspherical particles (ellipsoids):
"compute temp/sphere"_compute_temp_sphere.html
"compute temp/asphere"_compute_temp_asphere.html
"compute erotate/sphere"_compute_erotate_sphere.html
"compute erotate/asphere"_compute_erotate_asphere.html :ul
These include rotational degrees of freedom in their computation. If
you wish the thermodynamic output of temperature or pressure to use
one of these computes (e.g. for a system entirely composed of extended
particles), then the compute can be defined and the
"thermo_modify"_thermo_modify.html command used. Note that by
default thermodynamic quantities will be calculated with a temperature
that only includes translational degrees of freedom. See the
"thermo_style"_thermo_style.html command for details.
The "dump custom"_dump.html command can output various attributes of
extended particles, including the dipole moment (mu), the angular
velocity (omega), the angular momentum (angmom), the quaternion
(quat), and the torque (tq) on the particle.
Rigid bodies composed of extended particles :h5
The "fix rigid"_fix_rigid.html command treats a collection of
particles as a rigid body, computes its inertia tensor, sums the total
force and torque on the rigid body each timestep due to forces on its
constituent particles, and integrates the motion of the rigid body.
If any of the constituent particles of a rigid body are extended
particles (spheres or ellipsoids), then their contribution to the
inertia tensor of the body is different than if they were point
particles. This means the rotational dynamics of the rigid body will
be different. Thus a model of a dimer is different if the dimer
consists of two point masses versus two extended sphereoids, even if
the two particles have the same mass. Extended particles that
experience torque due to their interaction with other particles will
also impart that torque to a rigid body they are part of.
See the "fix rigid" command for example of complex rigid-body models
it is possible to define in LAMMPS.
Note that the "fix shake"_fix_shake.html command can also be used to
treat 2, 3, or 4 particles as a rigid body, but it always assumes the
particles are point masses.
:line
6.15 Output from LAMMPS (thermo, dumps, computes, fixes, variables) :link(howto_15),h4
There are four basic kinds of LAMMPS output:
"Thermodynamic output"_thermo_style.html, which is a list
of quantities printed every few timesteps to the screen and logfile. :ulb,l
"Dump files"_dump.html, which contain snapshots of atoms and various
per-atom values and are written at a specified frequency. :l
Certain fixes can output user-specified quantities to files: "fix
ave/time"_fix_ave_time.html for time averaging, "fix
ave/spatial"_fix_ave_spatial.html for spatial averaging, and "fix
print"_fix_print.html for single-line output of
"variables"_variable.html. Fix print can also output to the
screen. :l
"Restart files"_restart.html. :l,ule
A simulation prints one set of thermodynamic output and (optionally)
restart files. It can generate any number of dump files and fix
output files, depending on what "dump"_dump.html and "fix"_fix.html
commands you specify.
As discussed below, LAMMPS gives you a variety of ways to determine
what quantities are computed and printed when the thermodynamics,
dump, or fix commands listed above perform output. Throughout this
discussion, note that users can also "add their own computes and fixes
to LAMMPS"_Section_modify.html which can then generate values that can
then be output with these commands.
The following sub-sections discuss different LAMMPS command related
to output and the kind of data they operate on and produce:
"Global/per-atom/local data"_#global
"Scalar/vector/array data"_#scalar
"Thermodynamic output"_#thermo
"Dump file output"_#dump
"Fixes that write output files"_#fixoutput
"Computes that process output quantities"_#computeoutput
"Fixes that process output quantities"_#fixoutput
"Computes that generate values to output"_#compute
"Fixes that generate values to output"_#fix
"Variables that generate values to output"_#variable
"Summary table of output options and data flow between commands"_#table :ul
Global/per-atom/local data :h5,link(global)
Various output-related commands work with three different styles of
data: global, per-atom, or local. A global datum is one or more
system-wide values, e.g. the temperature of the system. A per-atom
datum is one or more values per atom, e.g. the kinetic energy of each
atom. Local datums are calculated by each processor based on the
atoms it owns, but there may be zero or more per atom, e.g. a list of
bond distances.
Scalar/vector/array data :h5,link(scalar)
Global, per-atom, and local datums can each come in three kinds: a
single scalar value, a vector of values, or a 2d array of values. The
doc page for a "compute" or "fix" or "variable" that generates data
will specify both the style and kind of data it produces, e.g. a
per-atom vector.
When a quantity is accessed, as in many of the output commands
discussed below, it can be referenced via the following bracket
notation, where ID in this case is the ID of a compute. The leading
"c_" would be replaced by "f_" for a fix, or "v_" for a variable:
c_ID | entire scalar, vector, or array
c_ID\[I\] | one element of vector, one column of array
c_ID\[I\]\[J\] | one element of array :tb(s=|)
In other words, using one bracket reduces the dimension of the data
once (vector -> scalar, array -> vector). Using two brackets reduces
the dimension twice (array -> scalar). Thus a command that uses
scalar values as input can typically also process elements of a vector
or array.
Thermodynamic output :h5,link(thermo)
The frequency and format of thermodynamic output is set by the
"thermo"_thermo.html, "thermo_style"_thermo_style.html, and
"thermo_modify"_thermo_modify.html commands. The
"thermo_style"_thermo_style.html command also specifies what values
are calculated and written out. Pre-defined keywords can be specified
(e.g. press, etotal, etc). Three additional kinds of keywords can
also be specified (c_ID, f_ID, v_name), where a "compute"_compute.html
or "fix"_fix.html or "variable"_variable.html provides the value to be
output. In each case, the compute, fix, or variable must generate
global values for input to the "thermo_style custom"_dump.html
command.
Dump file output :h5,link(dump)
Dump file output is specified by the "dump"_dump.html and
"dump_modify"_dump_modify.html commands. There are several
pre-defined formats (dump atom, dump xtc, etc).
There is also a "dump custom"_dump.html format where the user
specifies what values are output with each atom. Pre-defined atom
attributes can be specified (id, x, fx, etc). Three additional kinds
of keywords can also be specified (c_ID, f_ID, v_name), where a
"compute"_compute.html or "fix"_fix.html or "variable"_variable.html
provides the values to be output. In each case, the compute, fix, or
variable must generate per-atom values for input to the "dump
custom"_dump.html command.
There is also a "dump local"_dump.html format where the user specifies
what local values to output. A pre-defined index keyword can be
specified to enumuerate the local values. Two additional kinds of
keywords can also be specified (c_ID, f_ID), where a
"compute"_compute.html or "fix"_fix.html or "variable"_variable.html
provides the values to be output. In each case, the compute or fix
must generate local values for input to the "dump local"_dump.html
command.
Fixes that write output files :h5,link(fixoutput)
Sevarl fixes take various quantities as input and can write output
files: "fix ave/time"_fix_ave_time.html, "fix
ave/spatial"_fix_ave_spatial.html, "fix ave/histo"_fix_ave_histo.html,
"fix ave/correlate"_fix_ave_correlate.html, and "fix
print"_fix_print.html.
The "fix ave/time"_fix_ave_time.html command enables direct output to
a file and/or time-averaging of global scalars or vectors. The user
specifies one or more quantities as input. These can be global
"compute"_compute.html values, global "fix"_fix.html values, or
"variables"_variable.html of any style except the atom style which
produces per-atom values. Since a variable can refer to keywords used
by the "thermo_style custom"_thermo_style.html command (like temp or
press) and individual per-atom values, a wide variety of quantities
can be time averaged and/or output in this way. If the inputs are one
or more scalar values, then the fix generate a global scalar or vector
of output. If the inputs are one or more vector values, then the fix
generates a global vector or array of output. The time-averaged
output of this fix can also be used as input to other output commands.
The "fix ave/spatial"_fix_ave_spatial.html command enables direct
output to a file of spatial-averaged per-atom quantities like those
output in dump files, within 1d layers of the simulation box. The
per-atom quantities can be atom density (mass or number) or atom
attributes such as position, velocity, force. They can also be
per-atom quantities calculated by a "compute"_compute.html, by a
"fix"_fix.html, or by an atom-style "variable"_variable.html. The
spatial-averaged output of this fix can also be used as input to other
output commands.
The "fix ave/histo"_fix_ave_histo.html command enables direct output
to a file of histogrammed quantities, which can be global or per-atom
or local quantities. The histogram output of this fix can also be
used as input to other output commands.
The "fix ave/correlate"_fix_ave_histo.html command enables direct
output to a file of time-correlated quantities, which can be global
scalars. The correlation matrix output of this fix can also be used
as input to other output commands.
The "fix print"_fix_print.html command can generate a line of output
written to the screen and log file or to a separate file, periodically
during a running simulation. The line can contain one or more
"variable"_variable.html values for any style variable except the atom
style). As explained above, variables themselves can contain
references to global values generated by "thermodynamic
keywords"_thermo_style.html, "computes"_compute.html,
"fixes"_fix.html, or other "variables"_variable.html, or to per-atom
values for a specific atom. Thus the "fix print"_fix_print.html
command is a means to output a wide variety of quantities separate
from normal thermodynamic or dump file output.
Computes that process output quantities :h5,link(computeoutput)
The "compute reduce"_compute_reduce.html and "compute
reduce/region"_compute_reduce.html commands take one or more per-atom
or local vector quantities as inputs and "reduce" them (sum, min, max,
ave) to scalar quantities. These are produced as output values which
can be used as input to other output commands.
The "compute slice"_compute_slice.html command take one or more global
vector or array quantities as inputs and extracts a subset of their
values to create a new vector or array. These are produced as output
values which can be used as input to other output commands.
The "compute property/atom"_compute_property_atom.html command takes a
list of one or more pre-defined atom attributes (id, x, fx, etc) and
stores the values in a per-atom vector or array. These are produced
as output values which can be used as input to other output commands.
The list of atom attributes is the same as for the "dump
custom"_dump.html command.
The "compute property/local"_compute_property_local.html command takes
a list of one or more pre-defined local attributes (bond info, angle
info, etc) and stores the values in a local vector or array. These
are produced as output values which can be used as input to other
output commands.
The "compute atom/molecule"_compute_atom_molecule.html command takes a
list of one or more per-atom quantities (from a compute, fix, per-atom
variable) and sums the quantities on a per-molecule basis. It
produces a global vector or array as output values which can be used
as input to other output commands.
Fixes that process output quantities :h5,link(fixoutput)
The "fix ave/atom"_fix_ave_atom.html command performs time-averaging
of per-atom vectors. The per-atom quantities can be atom attributes
such as position, velocity, force. They can also be per-atom
quantities calculated by a "compute"_compute.html, by a
"fix"_fix.html, or by an atom-style "variable"_variable.html. The
time-averaged per-atom output of this fix can be used as input to
other output commands.
The "fix store/state"_fix_store_state.html command can archive one or
more per-atom attributes at a particular time, so that the old values
can be used in a future calculation or output. The list of atom
attributes is the same as for the "dump custom"_dump.html command,
including per-atom quantities calculated by a "compute"_compute.html,
by a "fix"_fix.html, or by an atom-style "variable"_variable.html.
The output of this fix can be used as input to other output commands.
Computes that generate values to output :h5,link(compute)
Every "compute"_compute.html in LAMMPS produces either global or
per-atom or local values. The values can be scalars or vectors or
arrays of data. These values can be output using the other commands
described in this section. The doc page for each compute command
describes what it produces. Computes that produce per-atom or local
values have the word "atom" or "local" in their style name. Computes
without the word "atom" or "local" produce global values.
Fixes that generate values to output :h5,link(fix)
Some "fixes"_fix.html in LAMMPS produces either global or per-atom or
local values which can be accessed by other commands. The values can
be scalars or vectors or arrays of data. These values can be output
using the other commands described in this section. The doc page for
each fix command tells whether it produces any output quantities and
describes them.
Variables that generate values to output :h5,link(variable)
Every "variables"_variable.html defined in an input script generates
either a global scalar value or a per-atom vector (only atom-style
variables) when it is accessed. The formulas used to define equal-
and atom-style variables can contain references to the thermodynamic
keywords and to global and per-atom data generated by computes, fixes,
and other variables. The values generated by variables can be output
using the other commands described in this section.
Summary table of output options and data flow between commands :h5,link(table)
This table summarizes the various commands that can be used for
generating output from LAMMPS. Each command produces output data of
some kind and/or writes data to a file. Most of the commands can take
data from other commands as input. Thus you can link many of these
commands together in pipeline form, where data produced by one command
is used as input to another command and eventually written to the
screen or to a file. Note that to hook two commands together the
output and input data types must match, e.g. global/per-atom/local
data and scalar/vector/array data.
Also note that, as described above, when a command takes a scalar as
input, that could be an element of a vector or array. Likewise a
vector input could be a column of an array.
Command: Input: Output:
"thermo_style custom"_thermo_style.html: global scalars: screen, log file:
"dump custom"_dump.html: per-atom vectors: dump file:
"dump local"_dump.html: local vectors: dump file:
"fix print"_fix_print.html: global scalar from variable: screen, file:
"print"_print.html: global scalar from variable: screen:
"computes"_compute.html: N/A: global/per-atom/local scalar/vector/array:
"fixes"_fix.html: N/A: global/per-atom/local scalar/vector/array:
"variables"_variable.html: global scalars, per-atom vectors: global scalar, per-atom vector:
"compute reduce"_compute_reduce.html: per-atom/local vectors: global scalar/vector:
"compute slice"_compute_slice.html: global vectors/arrays: global vector/array:
"compute property/atom"_compute_property_atom.html: per-atom vectors: per-atom vector/array:
"compute property/local"_compute_property_local.html: local vectors: local vector/array:
"compute atom/molecule"_compute_atom_molecule.html: per-atom vectors: global vector/array:
"fix ave/atom"_fix_ave_atom.html: per-atom vectors: per-atom vector/array:
"fix ave/time"_fix_ave_time.html: global scalars/vectors: global scalar/vector/array, file:
"fix ave/spatial"_fix_ave_spatial.html: per-atom vectors: global array, file:
"fix ave/histo"_fix_ave_histo.html: global/per-atom/local scalars and vectors: global array, file:
"fix ave/correlate"_fix_ave_correlate.html: global scalars: global array, file:
"fix store/state"_fix_store_state.html: per-atom vectors: per-atom vector/array:
:tb(s=:)
:line
6.16 Thermostatting, barostatting, and computing temperature :link(howto_16),h4
Thermostatting means controlling the temperature of particles in an MD
simulation. Barostatting means controlling the pressure. Since the
pressure includes a kinetic component due to particle velocities, both
these operations require calculation of the temperature. Typically a
target temperature (T) and/or pressure (P) is specified by the user,
and the thermostat or barostat attempts to equilibrate the system to
the requested T and/or P.
Temperature is computed as kinetic energy divided by some number of
degrees of freedom (and the Boltzmann constant). Since kinetic energy
is a function of particle velocity, there is often a need to
distinguish between a particle's advection velocity (due to some
aggregate motiion of particles) and its thermal velocity. The sum of
the two is the particle's total velocity, but the latter is often what
is wanted to compute a temperature.
LAMMPS has several options for computing temperatures, any of which
can be used in thermostatting and barostatting. These "compute
commands"_compute.html calculate temperature, and the "compute
pressure"_compute_pressure.html command calculates pressure.
"compute temp"_compute_temp.html
"compute temp/sphere"_compute_temp_sphere.html
"compute temp/asphere"_compute_temp_asphere.html
"compute temp/com"_compute_temp_com.html
"compute temp/deform"_compute_temp_deform.html
"compute temp/partial"_compute_temp_partial.html
"compute temp/profile"_compute_temp_profile.html
"compute temp/ramp"_compute_temp_ramp.html
"compute temp/region"_compute_temp_region.html :ul
All but the first 3 calculate velocity biases (i.e. advection
velocities) that are removed when computing the thermal temperature.
"Compute temp/sphere"_compute_temp_sphere.html and "compute
temp/asphere"_compute_temp_asphere.html compute kinetic energy for
extended particles that includes rotational degrees of freedom. They
both allow, as an extra argument, which is another temperature compute
that subtracts a velocity bias. This allows the translational
velocity of extended spherical or aspherical particles to be adjusted
in prescribed ways.
Thermostatting in LAMMPS is performed by "fixes"_fix.html, or in one
case by a pair style. Four thermostatting fixes are currently
available: Nose-Hoover (nvt), Berendsen, Langevin, and direct
rescaling (temp/rescale). Dissipative particle dynamics (DPD)
thermostatting can be invoked via the {dpd/tstat} pair style:
"fix nvt"_fix_nh.html
"fix nvt/sphere"_fix_nvt_sphere.html
"fix nvt/asphere"_fix_nvt_asphere.html
"fix nvt/sllod"_fix_nvt_sllod.html
"fix temp/berendsen"_fix_temp_berendsen.html
"fix langevin"_fix_langevin.html
"fix temp/rescale"_fix_temp_rescale.html
"pair_style dpd/tstat"_pair_dpd.html :ul
"Fix nvt"_fix_nh.html only thermostats the translational velocity of
particles. "Fix nvt/sllod"_fix_nvt_sllod.html also does this, except
that it subtracts out a velocity bias due to a deforming box and
integrates the SLLOD equations of motion. See the "NEMD
simulations"_#howto_13 section of this page for further details. "Fix
nvt/sphere"_fix_nvt_sphere.html and "fix
nvt/asphere"_fix_nvt_asphere.html thermostat not only translation
velocities but also rotational velocities for spherical and aspherical
particles.
DPD thermostatting alters pairwise interactions in a manner analagous
to the per-particle thermostatting of "fix
langevin"_fix_langevin.html.
Any of the thermostatting fixes can use temperature computes that
remove bias for two purposes: (a) computing the current temperature to
compare to the requested target temperature, and (b) adjusting only
the thermal temperature component of the particle's velocities. See
the doc pages for the individual fixes and for the
"fix_modify"_fix_modify.html command for instructions on how to assign
a temperature compute to a thermostatting fix. For example, you can
apply a thermostat to only the x and z components of velocity by using
it in conjunction with "compute
temp/partial"_compute_temp_partial.html.
IMPORTANT NOTE: Only the nvt fixes perform time integration, meaning
they update the velocities and positions of particles due to forces
and velocities respectively. The other thermostat fixes only adjust
velocities; they do NOT perform time integration updates. Thus they
should be used in conjunction with a constant NVE integration fix such
as these:
"fix nve"_fix_nve.html
"fix nve/sphere"_fix_nve_sphere.html
"fix nve/asphere"_fix_nve_asphere.html :ul
Barostatting in LAMMPS is also performed by "fixes"_fix.html. Two
barosttating methods are currently available: Nose-Hoover (npt and
nph) and Berendsen:
"fix npt"_fix_nh.html
"fix npt/sphere"_fix_npt_sphere.html
"fix npt/asphere"_fix_npt_asphere.html
"fix nph"_fix_nh.html
"fix press/berendsen"_fix_press_berendsen.html :ul
The "fix npt"_fix_nh.html commands include a Nose-Hoover thermostat
and barostat. "Fix nph"_fix_nh.html is just a Nose/Hoover barostat;
it does no thermostatting. Both "fix nph"_fix_nh.html and "fix
press/bernendsen"_fix_press_berendsen.html can be used in conjunction
with any of the thermostatting fixes.
As with the thermostats, "fix npt"_fix_nh.html and "fix
nph"_fix_nh.html only use translational motion of the particles in
computing T and P and performing thermo/barostatting. "Fix
npt/sphere"_fix_npt_sphere.html and "fix
npt/asphere"_fix_npt_asphere.html thermo/barostat using not only
translation velocities but also rotational velocities for spherical
and aspherical particles.
All of the barostatting fixes use the "compute
pressure"_compute_pressure.html compute to calculate a current
pressure. By default, this compute is created with a simple "compute
temp"_compute_temp.html (see the last argument of the "compute
pressure"_compute_pressure.html command), which is used to calculated
the kinetic componenet of the pressure. The barostatting fixes can
also use temperature computes that remove bias for the purpose of
computing the kinetic componenet which contributes to the current
pressure. See the doc pages for the individual fixes and for the
"fix_modify"_fix_modify.html command for instructions on how to assign
a temperature or pressure compute to a barostatting fix.
IMPORTANT NOTE: As with the thermostats, the Nose/Hoover methods ("fix
npt"_fix_nh.html and "fix nph"_fix_nh.html) perform time
integration. "Fix press/berendsen"_fix_press_berendsen.html does NOT,
so it should be used with one of the constant NVE fixes or with one of
the NVT fixes.
Finally, thermodynamic output, which can be setup via the
"thermo_style"_thermo_style.html command, often includes temperature
and pressure values. As explained on the doc page for the
"thermo_style"_thermo_style.html command, the default T and P are
setup by the thermo command itself. They are NOT the ones associated
with any thermostatting or barostatting fix you have defined or with
any compute that calculates a temperature or pressure. Thus if you
want to view these values of T and P, you need to specify them
explicitly via a "thermo_style custom"_thermo_style.html command. Or
you can use the "thermo_modify"_thermo_modify.html command to
re-define what temperature or pressure compute is used for default
thermodynamic output.
:line
6.17 Walls :link(howto_17),h4
Walls in an MD simulation are typically used to bound particle motion,
i.e. to serve as a boundary condition.
Walls in LAMMPS can be of rough (made of particles) or idealized
surfaces. Ideal walls can be smooth, generating forces only in the
normal direction, or frictional, generating forces also in the
tangential direction.
Rough walls, built of particles, can be created in various ways. The
particles themselves can be generated like any other particle, via the
"lattice"_lattice.html and "create_atoms"_create_atoms.html commands,
or read in via the "read_data"_read_data.html command.
Their motion can be constrained by many different commands, so that
they do not move at all, move together as a group at constant velocity
or in response to a net force acting on them, move in a prescribed
fashion (e.g. rotate around a point), etc. Note that if a time
integration fix like "fix nve"_fix_nve.html or "fix nvt"_fix_nh.html
is not used with the group that contains wall particles, their
positions and velocities will not be updated.
"fix aveforce"_fix_aveforce.html - set force on particles to average value, so they move together
"fix setforce"_fix_setforce.html - set force on particles to a value, e.g. 0.0
"fix freeze"_fix_freeze.html - freeze particles for use as granular walls
"fix nve/noforce"_fix_nve_noforce.html - advect particles by their velocity, but without force
"fix move"_fix_move.html - prescribe motion of particles by a linear velocity, oscillation, rotation, variable :ul
The "fix move"_fix_move.html command offers the most generality, since
the motion of individual particles can be specified with
"variable"_variable.html formula which depends on time and/or the
particle position.
For rough walls, it may be useful to turn off pairwise interactions
between wall particles via the "neigh_modify
exclude"_neigh_modify.html command.
Rough walls can also be created by specifying frozen particles that do
not move and do not interact with mobile particles, and then tethering
other particles to the fixed particles, via a "bond"_bond_style.html.
The bonded particles do interact with other mobile particles.
Idealized walls can be specified via several fix commands. "Fix
wall/gran"_fix_wall_gran.html creates frictional walls for use with
granular particles; all the other commands create smooth walls.
"fix wall/reflect"_fix_wall_reflect.html - reflective flat walls
"fix wall/lj93"_fix_wall.html - flat walls, with Lennard-Jones 9/3 potential
"fix wall/lj126"_fix_wall.html - flat walls, with Lennard-Jones 12/6 potential
"fix wall/colloid"_fix_wall.html - flat walls, with "pair_style colloid"_pair_colloid.html potential
"fix wall/harmonic"_fix_wall.html - flat walls, with repulsive harmonic spring potential
"fix wall/region"_fix_wall_region.html - use region surface as wall
"fix wall/gran"_fix_wall_gran.html - flat or curved walls with "pair_style granular"_pair_gran.html potential :ul
The {lj93}, {lj126}, {colloid}, and {harmonic} styles all allow the
flat walls to move with a constant velocity, or oscillate in time.
The "fix wall/region"_fix_wall_region.html command offers the most
generality, since the region surface is treated as a wall, and the
geometry of the region can be a simple primitive volume (e.g. a
sphere, or cube, or plane), or a complex volume made from the union
and intersection of primitive volumes. "Regions"_region.html can also
specify a volume "interior" or "exterior" to the specified primitive
shape or {union} or {intersection}. "Regions"_region.html can also be
"dynamic" meaning they move with constant velocity, oscillate, or
rotate.
The only frictional idealized walls currently in LAMMPS are flat or
curved surfaces specified by the "fix wall/gran"_fix_wall_gran.html
command. At some point we plan to allow regoin surfaces to be used as
frictional walls, as well as triangulated surfaces.
:line
6.18 Elastic constants :link(howto_18),h4
Elastic constants characterize the stiffness of a material. The formal
definition is provided by the linear relation that holds between the
stress and strain tensors in the limit of infinitesimal deformation.
In tensor notation, this is expressed as s_ij = C_ijkl * e_kl, where
the repeated indices imply summation. s_ij are the elements of the
symmetric stress tensor. e_kl are the elements of the symmetric strain
tensor. C_ijkl are the elements of the fourth rank tensor of elastic
constants. In three dimensions, this tensor has 3^4=81 elements. Using
Voigt notation, the tensor can be written as a 6x6 matrix, where C_ij
is now the derivative of s_i w.r.t. e_j. Because s_i is itself a
derivative w.r.t. e_i, it follows that C_ij is also symmetric, with at
most 7*6/2 = 21 distinct elements.
At zero temperature, it is easy to estimate these derivatives by
deforming the cell in one of the six directions using the command
"displace_box"_displace_box.html and measuring the change in the
stress tensor. A general-purpose script that does this is given in the
examples/elastic directory described in "this
section"_Section_example.html.
Calculating elastic constants at finite temperature is more
challenging, because it is necessary to run a simulation that perfoms
time averages of differential properties. One way to do this is to
measure the change in average stress tensor in an NVT simulations when
the cell volume undergoes a finite deformation. In order to balance
the systematic and statistical errors in this method, the magnitude of
the deformation must be chosen judiciously, and care must be taken to
fully equilibrate the deformed cell before sampling the stress
tensor. Another approach is to sample the triclinic cell fluctuations
that occur in an NPT simulation. This method can also be slow to
converge and requires careful post-processing "(Shinoda)"_#Shinoda
:line
6.19 Library interface to LAMMPS :link(howto_19),h4
As described in "Section_start 4"_Section_start.html#start_4, LAMMPS
can be built as a library, so that it can be called by another code,
used in a "coupled manner"_Section_howto.html#howto_10 with other
codes, or driven through a "Python interface"_Section_python.html.
All of these methodologies use a C-style interface to LAMMPS that is
provided in the files src/library.cpp and src/library.h. The
functions therein have a C-style argument list, but contain C++ code
you could write yourself in a C++ application that was invoking LAMMPS
directly. The C++ code in the functions illustrates how to invoke
internal LAMMPS operations. Note that LAMMPS classes are defined
within a LAMMPS namespace (LAMMPS_NS) if you use them from another C++
application.
Library.cpp contains these 4 functions:
void lammps_open(int, char **, MPI_Comm, void **);
void lammps_close(void *);
void lammps_file(void *, char *);
char *lammps_command(void *, char *); :pre
The lammps_open() function is used to initialize LAMMPS, passing in a
list of strings as if they were "command-line
arguments"_Section_start.html#start_6 when LAMMPS is run in
stand-alone mode from the command line, and a MPI communicator for
LAMMPS to run under. It returns a ptr to the LAMMPS object that is
created, and which is used in subsequent library calls. The
lammps_open() function can be called multiple 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. Or it might instantiate multiple
instances of LAMMPS to perform different calculations.
The lammps_close() function is used to shut down an instance of LAMMPS
and free all its memory.
The lammps_file() and lammps_command() functions are used to pass a
file or string to LAMMPS as if it were an input script or single
command in an input script. Thus the calling code can read or
generate a series of LAMMPS commands one line at a time and pass it
thru the library interface to setup a problem and then run it,
interleaving the lammps_command() calls with other calls to extract
information from LAMMPS, perform its own operations, or call another
code's library.
Other useful functions are also included in library.cpp. For example:
void *lammps_extract_global(void *, char *)
void *lammps_extract_atom(void *, char *)
void *lammps_extract_compute(void *, char *, int, int)
void *lammps_extract_fix(void *, char *, int, int, int, int)
void *lammps_extract_variable(void *, char *, char *)
int lammps_get_natoms(void *)
void lammps_get_coords(void *, double *)
void lammps_put_coords(void *, double *) :pre
These can extract various global or per-atom quantities from LAMMPS as
well as values calculated by a compute, fix, or variable. The "get"
and "put" operations can retrieve and reset atom coordinates.
See the library.cpp file and its associated header file library.h for
details.
The key idea of the library interface is that you can write any
functions you wish to define how your code talks to LAMMPS and add
them to src/library.cpp and src/library.h, as well as to the "Python
interface"_Section_python.html. The routines you add can access
or change any LAMMPS data you wish. The couple and python directories
have example C++ and C and Python codes which show how a driver code
can link to LAMMPS as a library, run LAMMPS on a subset of processors,
grab data from LAMMPS, change it, and put it back into LAMMPS.
:line
6.20 Calculating thermal conductivity :link(howto_20),h4
The thermal conductivity kappa of a material can be measured in at
least 3 ways using various options in LAMMPS. (See "this
section"_Section_howto.html#howto_21 of the manual for an analogous
discussion for viscosity). The thermal conducitivity tensor kappa is
a measure of the propensity of a material to transmit heat energy in a
diffusive manner as given by Fourier's law
J = -kappa grad(T)
where J is the heat flux in units of energy per area per time and
grad(T) is the spatial gradient of temperature. The thermal
conductivity thus has units of energy per distance per time per degree
K and is often approximated as an isotropic quantity, i.e. as a
scalar.
The first method is to setup two thermostatted regions at opposite
ends of a simulation box, or one in the middle and one at the end of a
periodic box. By holding the two regions at different temperatures
with a "thermostatting fix"_Section_howto.html#howto_13, the energy added
to the hot region should equal the energy subtracted from the cold
region and be proportional to the heat flux moving between the
regions. See the paper by "Ikeshoji and Hafskjold"_#Ikeshoji for
details of this idea. Note that thermostatting fixes such as "fix
nvt"_fix_nh.html, "fix langevin"_fix_langevin.html, and "fix
temp/rescale"_fix_temp_rescale.html store the cumulative energy they
add/subtract. Alternatively, the "fix heat"_fix_heat.html command can
used in place of thermostats on each of two regions, and the resulting
temperatures of the two regions monitored with the "compute
temp/region" command or the temperature profile of the intermediate
region monitored with the "fix ave/spatial"_fix_ave_spatial.html and
"compute ke/atom"_compute_ke_atom.html commands.
The second method is to perform a reverse non-equilibrium MD
simulation using the "fix
thermal/conductivity"_fix_thermal_conductivity.html command which
implements the rNEMD algorithm of Muller-Plathe. Kinetic energy is
swapped between atoms in two different layers of the simulation box.
This induces a temperature gradient between the two layers which can
be monitored with the "fix ave/spatial"_fix_ave_spatial.html and
"compute ke/atom"_compute_ke_atom.html commands. The fix tallies the
cumulative energy transfer that it performs. See the "fix
thermal/conductivity"_fix_thermal_conductivity.html command for
details.
The third method is based on the Green-Kubo (GK) formula which relates
the ensemble average of the auto-correlation of the heat flux to
kappa. The heat flux can be calculated from the fluctuations of
per-atom potential and kinetic energies and per-atom stress tensor in
a steady-state equilibrated simulation. This is in contrast to the
two preceding non-equilibrium methods, where energy flows continuously
between hot and cold regions of the simulation box.
The "compute heat/flux"_compute_heat_flux.html command can calculate
the needed heat flux and describes how to implement the Green_Kubo
formalism using additional LAMMPS commands, such as the "fix
ave/correlate"_fix_ave_correlate.html command to calculate the needed
auto-correlation. See the doc page for the "compute
heat/flux"_compute_heat_flux.html command for an example input script
that calculates the thermal conductivity of solid Ar via the GK
formalism.
:line
6.21 Calculating viscosity :link(howto_21),h4
The shear viscosity eta of a fluid can be measured in at least 3 ways
using various options in LAMMPS. (See "this
section"_Section_howto.html#howto_20 of the manual for an analogous
discussion for thermal conductivity). Eta is a measure of the
propensity of a fluid to transmit momentum in a direction
perpendicular to the direction of velocity or momentum flow.
Alternatively it is the resistance the fluid has to being sheared. It
is given by
J = -eta grad(Vstream)
where J is the momentum flux in units of momentum per area per time.
and grad(Vstream) is the spatial gradient of the velocity of the fluid
moving in another direction, normal to the area through which the
momentum flows. Viscosity thus has units of pressure-time.
The first method is to perform a non-equlibrium MD (NEMD) simulation
by shearing the simulation box via the "fix deform"_fix_deform.html
command, and using the "fix nvt/sllod"_fix_nvt_sllod.html command to
thermostat the fluid via the SLLOD equations of motion. The velocity
profile setup in the fluid by this procedure can be monitored by the
"fix ave/spatial"_fix_ave_spatial.html command, which determines
grad(Vstream) in the equation above. E.g. the derivative in the
y-direction of the Vx component of fluid motion or grad(Vstream) =
dVx/dy. In this case, the Pxy off-diagonal component of the pressure
or stress tensor, as calculated by the "compute
pressure"_compute_pressure.html command, can also be monitored, which
is the J term in the equation above. See "this
section"_Section_howto.html#howto_13 of the manual for details on NEMD
simulations.
The second method is to perform a reverse non-equilibrium MD
simulation using the "fix viscosity"_fix_viscosity.html command which
implements the rNEMD algorithm of Muller-Plathe. Momentum in one
dimension is swapped between atoms in two different layers of the
simulation box in a different dimension. This induces a velocity
gradient which can be monitored with the "fix
ave/spatial"_fix_ave_spatial.html command. The fix tallies the
cummulative momentum transfer that it performs. See the "fix
viscosity"_fix_viscosity.html command for details.
The third method is based on the Green-Kubo (GK) formula which relates
the ensemble average of the auto-correlation of the stress/pressure
tensor to eta. This can be done in a steady-state equilibrated
simulation which is in contrast to the two preceding non-equilibrium
methods, where momentum flows continuously through the simulation box.
Here is an example input script that calculates the viscosity of
liquid Ar via the GK formalism:
# Sample LAMMPS input script for viscosity of liquid Ar :pre
units real
variable T equal 86.4956
variable V equal vol
variable dt equal 4.0
variable p equal 400 # correlation length
variable s equal 5 # sample interval
variable d equal $p*$s # dump interval :pre
# convert from LAMMPS real units to SI :pre
variable kB equal 1.3806504e-23 # \[J/K/] Boltzmann
variable atm2Pa equal 101325.0
variable A2m equal 1.0e-10
variable fs2s equal 1.0e-15
variable convert equal $\{atm2Pa\}*$\{atm2Pa\}*$\{fs2s\}*$\{A2m\}*$\{A2m\}*$\{A2m\} :pre
# setup problem :pre
dimension 3
boundary p p p
lattice fcc 5.376 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
region box block 0 4 0 4 0 4
create_box 1 box
create_atoms 1 box
mass 1 39.948
pair_style lj/cut 13.0
pair_coeff * * 0.2381 3.405
timestep $\{dt\}
thermo $d :pre
# equilibration and thermalization :pre
velocity all create $T 102486 mom yes rot yes dist gaussian
fix NVT all nvt temp $T $T 10 drag 0.2
run 8000 :pre
# viscosity calculation, switch to NVE if desired :pre
#unfix NVT
#fix NVE all nve :pre
reset_timestep 0
variable pxy equal pxy
variable pxz equal pxz
variable pyz equal pyz
fix SS all ave/correlate $s $p $d &
v_pxy v_pxz v_pyz type auto file S0St.dat ave running
variable scale equal $\{convert\}/($\{kB\}*$T)*$V*$s*$\{dt\}
variable v11 equal trap(f_SS\[3/])*$\{scale\}
variable v22 equal trap(f_SS\[4/])*$\{scale\}
variable v33 equal trap(f_SS\[5/])*$\{scale\}
thermo_style custom step temp press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33
run 100000
variable v equal (v_v11+v_v22+v_v33)/3.0
variable ndens equal count(all)/vol
print "average viscosity: $v \[Pa.s/] @ $T K, $\{ndens\} /A^3" :pre
:line
:line
:link(Berendsen)
[(Berendsen)] Berendsen, Grigera, Straatsma, J Phys Chem, 91,
6269-6271 (1987).
:link(Cornell)
[(Cornell)] Cornell, Cieplak, Bayly, Gould, Merz, Ferguson,
Spellmeyer, Fox, Caldwell, Kollman, JACS 117, 5179-5197 (1995).
:link(Horn)
[(Horn)] Horn, Swope, Pitera, Madura, Dick, Hura, and Head-Gordon,
J Chem Phys, 120, 9665 (2004).
:link(Ikeshoji)
[(Ikeshoji)] Ikeshoji and Hafskjold, Molecular Physics, 81, 251-261
(1994).
:link(MacKerell)
[(MacKerell)] MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field,
Fischer, Gao, Guo, Ha, et al, J Phys Chem, 102, 3586 (1998).
:link(Mayo)
[(Mayo)] Mayo, Olfason, Goddard III, J Phys Chem, 94, 8897-8909
(1990).
:link(Jorgensen)
[(Jorgensen)] Jorgensen, Chandrasekhar, Madura, Impey, Klein, J Chem
Phys, 79, 926 (1983).
:link(Price)
[(Price)] Price and Brooks, J Chem Phys, 121, 10096 (2004).
:link(Shinoda)
[(Shinoda)] Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004).