<liclass="toctree-l2"><aclass="reference internal"href="#charmm-amber-and-dreiding-force-fields">6.3. CHARMM, AMBER, and DREIDING force fields</a></li>
<liclass="toctree-l2"><aclass="reference internal"href="#running-multiple-simulations-from-one-input-script">6.4. Running multiple simulations from one input script</a></li>
<liclass="toctree-l3"><aclass="reference internal"href="#fixes-that-write-output-files">6.15.5. Fixes that write output files</a></li>
<liclass="toctree-l3"><aclass="reference internal"href="#computes-that-process-output-quantities">6.15.6. Computes that process output quantities</a></li>
<liclass="toctree-l3"><aclass="reference internal"href="#fixes-that-process-output-quantities">6.15.7. Fixes that process output quantities</a></li>
<liclass="toctree-l3"><aclass="reference internal"href="#computes-that-generate-values-to-output">6.15.8. Computes that generate values to output</a></li>
<liclass="toctree-l3"><aclass="reference internal"href="#fixes-that-generate-values-to-output">6.15.9. Fixes that generate values to output</a></li>
<liclass="toctree-l3"><aclass="reference internal"href="#variables-that-generate-values-to-output">6.15.10. Variables that generate values to output</a></li>
<liclass="toctree-l3"><aclass="reference internal"href="#summary-table-of-output-options-and-data-flow-between-commands">6.15.11. Summary table of output options and data flow between commands</a></li>
</ul>
</li>
<liclass="toctree-l2"><aclass="reference internal"href="#thermostatting-barostatting-and-computing-temperature">6.16. Thermostatting, barostatting, and computing temperature</a></li>
<liclass="toctree-l2"><aclass="reference internal"href="#calculating-a-diffusion-coefficient">6.22. Calculating a diffusion coefficient</a></li>
<liclass="toctree-l2"><aclass="reference internal"href="#using-chunks-to-calculate-system-properties">6.23. Using chunks to calculate system properties</a><ul>
<liclass="toctree-l3"><aclass="reference internal"href="#example-calculations-with-chunks">6.23.4. Example calculations with chunks</a></li>
</ul>
</li>
<liclass="toctree-l2"><aclass="reference internal"href="#setting-parameters-for-the-kspace-style-pppm-disp-command">6.24. Setting parameters for the <codeclass="docutils literal"><spanclass="pre">kspace_style</span><spanclass="pre">pppm/disp</span></code> command</a></li>
<divclass="line">6.3 <aclass="reference internal"href="#howto-3"><spanclass="std std-ref">CHARMM, AMBER, and DREIDING force fields</span></a></div>
<divclass="line">6.4 <aclass="reference internal"href="#howto-4"><spanclass="std std-ref">Running multiple simulations from one input script</span></a></div>
<divclass="line">6.22 <aclass="reference internal"href="#howto-22"><spanclass="std std-ref">Calculating a diffusion coefficient</span></a></div>
<divclass="line">6.23 <aclass="reference internal"href="#howto-23"><spanclass="std std-ref">Using chunks to calculate system properties</span></a></div>
<divclass="line">6.24 <aclass="reference internal"href="#howto-24"><spanclass="std std-ref">Setting parameters for the kspace_style pppm/disp command</span></a></div>
<p>The example input scripts included in the LAMMPS distribution and
highlighted in <aclass="reference internal"href="Section_example.html"><spanclass="doc">Section_example</span></a> also show how to
setup and run various kinds of simulations.</p>
<divclass="section"id="restarting-a-simulation">
<spanid="howto-1"></span><h2>6.1. Restarting a simulation</h2>
<p>There are 3 ways to continue a long LAMMPS simulation. Multiple
<aclass="reference internal"href="run.html"><spanclass="doc">run</span></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 <aclass="reference internal"href="restart.html"><spanclass="doc">restart</span></a>
command. At a later time, these binary files can be read via a
<aclass="reference internal"href="read_restart.html"><spanclass="doc">read_restart</span></a> command in a new script. Or they can
be converted to text data files using the <aclass="reference internal"href="Section_start.html#start-7"><spanclass="std std-ref">-r command-line switch</span></a> and read by a
<aclass="reference internal"href="read_data.html"><spanclass="doc">read_data</span></a> command in a new script.</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 <aclass="reference internal"href="read_restart.html"><spanclass="doc">read_restart</span></a> and
<spanid="howto-3"></span><h2>6.3. CHARMM, AMBER, and DREIDING force fields</h2>
<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
<aclass="reference internal"href="read_data.html"><spanclass="doc">read_data</span></a> command or in the input script with
commands like <aclass="reference internal"href="pair_coeff.html"><spanclass="doc">pair_coeff</span></a> or
<aclass="reference internal"href="bond_coeff.html"><spanclass="doc">bond_coeff</span></a>. See <aclass="reference internal"href="Section_tools.html"><spanclass="doc">Section_tools</span></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 <aclass="reference internal"href="#howto-mackerell"><spanclass="std std-ref">(MacKerell)</span></a> for a description of the CHARMM force
field. See <aclass="reference internal"href="#howto-cornell"><spanclass="std std-ref">(Cornell)</span></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
<p>DREIDING is a generic force field developed by the <aclass="reference external"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
<aclass="reference internal"href="pair_hbond_dreiding.html"><spanclass="doc">explicit hydrogen bond term</span></a> to describe
interactions involving a hydrogen atom on very electronegative atoms
(N, O, F).</p>
<p>See <aclass="reference internal"href="#howto-mayo"><spanclass="std std-ref">(Mayo)</span></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
<p>would run 2 independent simulations, one after the other.</p>
<p>For large numbers of independent simulations, you can use
<aclass="reference internal"href="variable.html"><spanclass="doc">variables</span></a> and the <aclass="reference internal"href="next.html"><spanclass="doc">next</span></a> and
<aclass="reference internal"href="jump.html"><spanclass="doc">jump</span></a> commands to loop over the same input script
multiple times with different settings. For example, this
script, named in.polymer</p>
<divclass="highlight-default"><divclass="highlight"><pre><span></span>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></div>
</div>
<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>
<divclass="highlight-default"><divclass="highlight"><pre><span></span>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
clear
next t
next a
jump in.polymer
</pre></div>
</div>
<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 <aclass="reference internal"href="Section_start.html#start-7"><spanclass="std std-ref">this section</span></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 <em>universe</em>-style variables, as described in the
<aclass="reference internal"href="variable.html"><spanclass="doc">variable</span></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>
<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>
<ulclass="simple">
<li><aclass="reference internal"href="neb.html"><spanclass="doc">neb</span></a> for nudged elastic band calculations</li>
<li><aclass="reference internal"href="prd.html"><spanclass="doc">prd</span></a> for parallel replica dynamics</li>
<li><aclass="reference internal"href="tad.html"><spanclass="doc">tad</span></a> for temperature accelerated dynamics</li>
<li><aclass="reference internal"href="temper.html"><spanclass="doc">temper</span></a> for parallel tempering</li>
<li><aclass="reference internal"href="fix_pimd.html"><spanclass="doc">fix pimd</span></a> for path-integral molecular dynamics (PIMD)</li>
</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 commands can only be used if LAMMPS was built with the REPLICA
package. See the <aclass="reference internal"href="Section_start.html#start-3"><spanclass="std std-ref">Making LAMMPS</span></a> section
for more info on packages.</p>
<p>PIMD runs different replicas whose individual particles are coupled
together by springs to model a system or ring-polymers.</p>
<p>This commands can only be used if LAMMPS was built with the USER-MISC
package. See the <aclass="reference internal"href="Section_start.html#start-3"><spanclass="std std-ref">Making LAMMPS</span></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 <aclass="reference internal"href="Section_start.html#start-7"><spanclass="std std-ref">-partition command-line switch</span></a> to launch LAMMPS on multiple
partitions, which in this context are the same as replicas. E.g.
<p>calculates rotational kinetic energy which can be <aclass="reference internal"href="#howto-15"><spanclass="std std-ref">output with thermodynamic info</span></a>.</p>
<p>Use one of these 3 pair potentials, which compute forces and torques
<p>For both models, the bond lengths and bond angles should be held fixed
using the <aclass="reference internal"href="fix_shake.html"><spanclass="doc">fix shake</span></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
<aclass="reference internal"href="pair_lj.html#jorgensen"><spanclass="std std-ref">(Jorgensen)</span></a>. Note that the OM distance is specified in
the <aclass="reference internal"href="pair_style.html"><spanclass="doc">pair_style</span></a> command, not as part of the pair
coefficients.</p>
<divclass="line-block">
<divclass="line">O mass = 15.9994</div>
<divclass="line">H mass = 1.008</div>
<divclass="line">O charge = -1.040</div>
<divclass="line">H charge = 0.520</div>
<divclass="line">r0 of OH bond = 0.9572</div>
<divclass="line">theta of HOH angle = 104.52</div>
<divclass="line">OM distance = 0.15</div>
<divclass="line">LJ epsilon of O-O = 0.1550</div>
<divclass="line">LJ sigma of O-O = 3.1536</div>
<divclass="line">LJ epsilon, sigma of OH, HH = 0.0</div>
<divclass="line">Coulombic cutoff = 8.5</div>
<divclass="line"><br/></div>
</div>
<p>For the TIP4/Ice model (J Chem Phys, 122, 234511 (2005);
<aclass="reference external"href="http://dx.doi.org/10.1063/1.1931662">http://dx.doi.org/10.1063/1.1931662</a>) these values can be used:</p>
<divclass="line-block">
<divclass="line">O mass = 15.9994</div>
<divclass="line">H mass = 1.008</div>
<divclass="line">O charge = -1.1794</div>
<divclass="line">H charge = 0.5897</div>
<divclass="line">r0 of OH bond = 0.9572</div>
<divclass="line">theta of HOH angle = 104.52</div>
<divclass="line">OM distance = 0.1577</div>
<divclass="line">LJ epsilon of O-O = 0.21084</div>
<divclass="line">LJ sigma of O-O = 3.1668</div>
<divclass="line">LJ epsilon, sigma of OH, HH = 0.0</div>
<divclass="line">Coulombic cutoff = 8.5</div>
<divclass="line"><br/></div>
</div>
<p>For the TIP4P/2005 model (J Chem Phys, 123, 234505 (2005);
<aclass="reference external"href="http://dx.doi.org/10.1063/1.2121687">http://dx.doi.org/10.1063/1.2121687</a>), these values can be used:</p>
<divclass="line-block">
<divclass="line">O mass = 15.9994</div>
<divclass="line">H mass = 1.008</div>
<divclass="line">O charge = -1.1128</div>
<divclass="line">H charge = 0.5564</div>
<divclass="line">r0 of OH bond = 0.9572</div>
<divclass="line">theta of HOH angle = 104.52</div>
<divclass="line">OM distance = 0.1546</div>
<divclass="line">LJ epsilon of O-O = 0.1852</div>
<divclass="line">LJ sigma of O-O = 3.1589</div>
<divclass="line">LJ epsilon, sigma of OH, HH = 0.0</div>
<divclass="line">Coulombic cutoff = 8.5</div>
<divclass="line"><br/></div>
</div>
<p>These are the parameters to use for TIP4P with a long-range Coulombic
solver (e.g. Ewald or PPPM in LAMMPS):</p>
<divclass="line-block">
<divclass="line">O mass = 15.9994</div>
<divclass="line">H mass = 1.008</div>
<divclass="line">O charge = -1.0484</div>
<divclass="line">H charge = 0.5242</div>
<divclass="line">r0 of OH bond = 0.9572</div>
<divclass="line">theta of HOH angle = 104.52</div>
<divclass="line">OM distance = 0.1250</div>
<divclass="line">LJ epsilon of O-O = 0.16275</div>
<divclass="line">LJ sigma of O-O = 3.16435</div>
<divclass="line">LJ epsilon, sigma of OH, HH = 0.0</div>
<divclass="line"><br/></div>
</div>
<p>Note that the when using the TIP4P pair style, the neighobr list
cutoff for Coulomb interactions is effectively extended by a distance
2 * (OM distance), to account for the offset distance of the
fictitious charges on O atoms in water molecules. Thus it is
typically best in an efficiency sense to use a LJ cutoff >= Coulomb
cutoff + 2*(OM distance), to shrink the size of the neighbor list.
This leads to slightly larger cost for the long-range calculation, so
you can test the trade-off for your model. The OM distance and the LJ
and Coulombic cutoffs are set in the <aclass="reference internal"href="pair_lj.html"><spanclass="doc">pair_style lj/cut/tip4p/long</span></a> command.</p>
<p>Wikipedia also has a nice article on <aclass="reference external"href="http://en.wikipedia.org/wiki/Water_model">water models</a>.</p>
<hrclass="docutils"/>
</div>
<divclass="section"id="spc-water-model">
<spanid="howto-9"></span><h2>6.9. SPC water model</h2>
<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 <aclass="reference internal"href="fix_shake.html"><spanclass="doc">fix shake</span></a> command can be used to hold
the two O-H bonds and the H-O-H angle rigid. A bond style of
<em>harmonic</em> and an angle style of <em>harmonic</em> or <em>charmm</em> 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>
<divclass="line-block">
<divclass="line">O mass = 15.9994</div>
<divclass="line">H mass = 1.008</div>
<divclass="line">O charge = -0.820</div>
<divclass="line">H charge = 0.410</div>
<divclass="line">LJ epsilon of OO = 0.1553</div>
<divclass="line">LJ sigma of OO = 3.166</div>
<divclass="line">LJ epsilon, sigma of OH, HH = 0.0</div>
<divclass="line">r0 of OH bond = 1.0</div>
<divclass="line">theta of HOH angle = 109.47</div>
<divclass="line"><br/></div>
</div>
<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>
<divclass="line-block">
<divclass="line">O charge = -0.8476</div>
<divclass="line">H charge = 0.4238</div>
<divclass="line"><br/></div>
</div>
<p>See the <aclass="reference internal"href="#howto-berendsen"><spanclass="std std-ref">(Berendsen)</span></a> reference for more details on both
the SPC and SPC/E models.</p>
<p>Wikipedia also has a nice article on <aclass="reference external"href="http://en.wikipedia.org/wiki/Water_model">water models</a>.</p>
<spanid="howto-10"></span><h2>6.10. Coupling LAMMPS to other codes</h2>
<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 <aclass="reference internal"href="fix.html"><spanclass="doc">fix</span></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
<aclass="reference external"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
<aclass="reference internal"href="fix_poems.html"><spanclass="doc">fix poems</span></a> command for more details. See <aclass="reference internal"href="Section_modify.html"><spanclass="doc">this section</span></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
<aclass="reference internal"href="run.html"><spanclass="doc">run</span></a> command facilitates this with its <em>every</em> 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 <aclass="reference internal"href="Section_modify.html"><spanclass="doc">Section_modify</span></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 <aclass="reference internal"href="run.html"><spanclass="doc">run</span></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 examples/COUPLE directory of the LAMMPS distribution; see
examples/COUPLE/README for more details:</p>
<ulclass="simple">
<li>simple: simple driver programs in C++ and C which invoke LAMMPS as a
library</li>
<li>lammps_quest: coupling of LAMMPS and <aclass="reference external"href="http://dft.sandia.gov/Quest">Quest</a>, to run classical
MD with quantum forces calculated by a density functional code</li>
<li>lammps_spparks: coupling of LAMMPS and <aclass="reference external"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</li>
</ul>
<p><aclass="reference internal"href="Section_start.html#start-5"><spanclass="std std-ref">This section</span></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
<aclass="reference internal"href="Section_python.html"><spanclass="doc">Section_python</span></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 <aclass="reference internal"href="#howto-19"><spanclass="std std-ref">Section_howto 19</span></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>
<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 <aclass="reference external"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, <aclass="reference external"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 <aclass="reference external"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 <aclass="reference external"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 <aclass="reference internal"href="dump.html"><spanclass="doc">dump</span></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 <aclass="reference external"href="http://www.ks.uiuc.edu/Research/vmd">VMD</a> for visualization.
See the <aclass="reference internal"href="dump.html"><spanclass="doc">dump</span></a> command for more information on XTC files.</p>
<p>By default, LAMMPS uses an orthogonal simulation box to encompass the
particles. The <aclass="reference internal"href="boundary.html"><spanclass="doc">boundary</span></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 <strong>a</strong> = (xhi-xlo,0,0); <strong>b</strong> =
(0,yhi-ylo,0); <strong>c</strong> = (0,0,zhi-zlo). The 6 parameters
(xlo,xhi,ylo,yhi,zlo,zhi) are defined at the time the simulation box
is created, e.g. by the <aclass="reference internal"href="create_box.html"><spanclass="doc">create_box</span></a> or
<aclass="reference internal"href="read_data.html"><spanclass="doc">read_data</span></a> or <aclass="reference internal"href="read_restart.html"><spanclass="doc">read_restart</span></a>
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 <aclass="reference internal"href="thermo_style.html"><spanclass="doc">thermo_style custom</span></a> command.</p>
<p>LAMMPS also allows simulations to be performed in triclinic
(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 <strong>a</strong> = (xhi-xlo,0,0); <strong>b</strong> = (xy,yhi-ylo,0); <strong>c</strong> =
(xz,yz,zhi-zlo). <em>xy,xz,yz</em> 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. In LAMMPS the triclinic
simulation box edge vectors <strong>a</strong>, <strong>b</strong>, and <strong>c</strong> cannot be arbitrary
vectors. As indicated, <strong>a</strong> must lie on the positive x axis. <strong>b</strong> must
lie in the xy plane, with strictly positive y component. <strong>c</strong> may have
any orientation with strictly positive z component. The requirement
that <strong>a</strong>, <strong>b</strong>, and <strong>c</strong> have strictly positive x, y, and z components,
respectively, ensures that <strong>a</strong>, <strong>b</strong>, and <strong>c</strong> form a complete
right-handed basis. These restrictions impose no loss of generality,
since it is possible to rotate/invert any set of 3 crystal basis
vectors so that they conform to the restrictions.</p>
<p>For example, assume that the 3 vectors <strong>A</strong>,<strong>B</strong>,<strong>C</strong> are the edge
vectors of a general parallelepiped, where there is no restriction on
<strong>A</strong>,<strong>B</strong>,<strong>C</strong> other than they form a complete right-handed basis i.e.
<strong>A</strong> x <strong>B</strong> . <strong>C</strong>> 0. The equivalent LAMMPS <strong>a</strong>,<strong>b</strong>,<strong>c</strong> are a linear
rotation of <strong>A</strong>, <strong>B</strong>, and <strong>C</strong> and can be computed as follows:</p>
<p>where A = <ahref="#id3"><spanclass="problematic"id="id4">|</span></a><strong>A</strong>| indicates the scalar length of <strong>A</strong>. The ^ hat symbol
indicates the corresponding unit vector. <em>beta</em> and <em>gamma</em> are angles
between the vectors described below. Note that by construction,
<strong>a</strong>, <strong>b</strong>, and <strong>c</strong> have strictly positive x, y, and z components, respectively.
If it should happen that
<strong>A</strong>, <strong>B</strong>, and <strong>C</strong> form a left-handed basis, then the above equations
are not valid for <strong>c</strong>. In this case, it is necessary
to first apply an inversion. This can be achieved
by interchanging two basis vectors or by changing the sign of one of them.</p>
<p>For consistency, the same rotation/inversion applied to the basis vectors
must also be applied to atom positions, velocities,
and any other vector quantities.
This can be conveniently achieved by first converting to
fractional coordinates in the
old basis and then converting to distance coordinates in the new basis.
The transformation is given by the following equation:</p>
<p>where <em>V</em> is the volume of the box, <strong>X</strong> is the original vector quantity and
<strong>x</strong> is the vector in the LAMMPS basis.</p>
<p>There is no requirement that a triclinic box be periodic in any
dimension, though it typically should be in at least the 2nd dimension
of the tilt (y in xy) if you want to enforce a shift in periodic
boundary conditions across that boundary. Some commands that work
with triclinic boxes, e.g. the <aclass="reference internal"href="fix_deform.html"><spanclass="doc">fix deform</span></a> and <aclass="reference internal"href="fix_nh.html"><spanclass="doc">fix npt</span></a> commands, require periodicity or non-shrink-wrap
boundary conditions in specific dimensions. See the command doc pages
for details.</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 <aclass="reference internal"href="create_box.html"><spanclass="doc">create_box</span></a> command is used with a region of
style <em>prism</em>, then a triclinic box is setup. See the
<aclass="reference internal"href="region.html"><spanclass="doc">region</span></a> command for details. If the
<aclass="reference internal"href="read_data.html"><spanclass="doc">read_data</span></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
<aclass="reference internal"href="read_data.html"><spanclass="doc">read_data</span></a> command for details. Finally, if the
<aclass="reference internal"href="read_restart.html"><spanclass="doc">read_restart</span></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 <aclass="reference internal"href="fix_nh.html"><spanclass="doc">fix npt</span></a> or
<aclass="reference internal"href="fix_deform.html"><spanclass="doc">fix deform</span></a> commands. Alternatively, you can use the
<aclass="reference internal"href="change_box.html"><spanclass="doc">change_box</span></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
<p>To avoid extremely tilted boxes (which would be computationally
inefficient), LAMMPS normally requires that 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). This is required
both when the simulation box is created, e.g. via the
<aclass="reference internal"href="create_box.html"><spanclass="doc">create_box</span></a> or <aclass="reference internal"href="read_data.html"><spanclass="doc">read_data</span></a> commands,
as well as when the box shape changes dynamically during a simulation,
e.g. via the <aclass="reference internal"href="fix_deform.html"><spanclass="doc">fix deform</span></a> or <aclass="reference internal"href="fix_nh.html"><spanclass="doc">fix npt</span></a>
commands.</p>
<p>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. If the box tilt exceeds this
limit during a dynamics run (e.g. via the <aclass="reference internal"href="fix_deform.html"><spanclass="doc">fix deform</span></a>
command), then the box is “flipped” to an equivalent shape with a tilt
factor within the bounds, so the run can continue. See the <aclass="reference internal"href="fix_deform.html"><spanclass="doc">fix deform</span></a> doc page for further details.</p>
<p>One exception to this rule is if the 1st dimension in the tilt
factor (x for xy) is non-periodic. In that case, the limits on the
tilt factor are not enforced, since flipping the box in that dimension
does not change the atom positions due to non-periodicity. In this
mode, if you tilt the system to extreme angles, the simulation will
simply become inefficient, due to the highly skewed simulation box.</p>
<p>The limitation on not creating a simulation box with a tilt factor
skewing the box more than half the distance of the parallel box length
can be overridden via the <aclass="reference internal"href="box.html"><spanclass="doc">box</span></a> command. Setting the <em>tilt</em>
keyword to <em>large</em> allows any tilt factors to be specified.</p>
<p>Box flips that may occur using the <aclass="reference internal"href="fix_deform.html"><spanclass="doc">fix deform</span></a> or
<aclass="reference internal"href="fix_nh.html"><spanclass="doc">fix npt</span></a> commands can be turned off using the <em>flip no</em>
option with either of the commands.</p>
<p>Note that if a simulation box has a large tilt factor, LAMMPS will run
less efficiently, due to the large volume of communication needed to
acquire ghost atoms around a processor’s irregular-shaped sub-domain.
For extreme values of tilt, LAMMPS may also lose atoms and generate an
error.</p>
<p>Triclinic crystal structures are often defined using three lattice
constants <em>a</em>, <em>b</em>, and <em>c</em>, and three angles <em>alpha</em>, <em>beta</em> and
<em>gamma</em>. Note that in this nomenclature, the a, b, and c lattice
constants are the scalar lengths of the edge vectors <strong>a</strong>, <strong>b</strong>, and <strong>c</strong>
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>
<p>A second use is to run Parinello-Rahman dyanamics via the <aclass="reference internal"href="fix_nh.html"><spanclass="doc">fix npt</span></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 <aclass="reference internal"href="minimize.html"><spanclass="doc">energy minimization</span></a> is
the <aclass="reference internal"href="fix_box_relax.html"><spanclass="doc">fix box/relax</span></a> command.</p>
<p>A third use is to shear a bulk solid to study the response of the
material. The <aclass="reference internal"href="fix_deform.html"><spanclass="doc">fix deform</span></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
<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 <aclass="reference internal"href="fix_deform.html"><spanclass="doc">fix deform</span></a> command. The
<aclass="reference internal"href="fix_nvt_sllod.html"><spanclass="doc">fix nvt/sllod</span></a> command can be used to thermostat
the sheared fluid and integrate the SLLOD equations of motion for the
system. Fix nvt/sllod uses <aclass="reference internal"href="compute_temp_deform.html"><spanclass="doc">compute temp/deform</span></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 <spanclass="xref doc">fix ave/spatial</span> 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, <aclass="reference internal"href="fix_deform.html"><spanclass="doc">fix deform</span></a> can continuously strain
a box by an arbitrary amount. As discussed in the <aclass="reference internal"href="fix_deform.html"><spanclass="doc">fix deform</span></a> command, when the tilt value reaches a limit,
the box is flipped 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 <aclass="reference internal"href="fix_deform.html"><spanclass="doc">fix deform</span></a> should be set to “remap v”, since that is what
<aclass="reference internal"href="fix_nvt_sllod.html"><spanclass="doc">fix nvt/sllod</span></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
<p>This is because when dipoles interact with each other, they induce
torques, and a particle must be finite-size (i.e. have a moment of
inertia) in order to respond and rotate. See the <aclass="reference internal"href="atom_style.html"><spanclass="doc">atom_style dipole</span></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>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 line style defines line segment particles with two end points and
a mass (or density). They can be used in 2d simulations, and they can
be joined together to form rigid bodies which represent arbitrary
polygons.</p>
<p>The tri style defines triangular particles with three corner points
and a mass (or density). They can be used in 3d simulations, and they
can be joined together to form rigid bodies which represent arbitrary
particles with a triangulated surface.</p>
<p>The peri style is used with <aclass="reference internal"href="pair_peri.html"><spanclass="doc">Peridynamic models</span></a> and
defines particles as having a volume, that is used internally in the
<p>The body style allows for definition of particles which can represent
complex entities, such as surface meshes of discrete points,
collections of sub-particles, deformable objects, etc. The body style
is discussed in more detail on the <aclass="reference internal"href="body.html"><spanclass="doc">body</span></a> doc page.</p>
<p>Note that if one of these atom styles is used (or multiple styles via
the <aclass="reference internal"href="atom_style.html"><spanclass="doc">atom_style hybrid</span></a> command), not all particles in
the system are required to be finite-size or aspherical.</p>
<p>For example, in the ellipsoid style, 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. In the line or
tri style, if the lineflag or triflag is specified as 0, then it
will be a point particle.</p>
<p>Some of the pair styles used to compute pairwise interactions between
finite-size particles also compute the correct interaction with point
particles as well, e.g. the interaction between a point particle and a
finite-size particle or between two point particles. If necessary,
<aclass="reference internal"href="pair_hybrid.html"><spanclass="doc">pair_style hybrid</span></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 <aclass="reference internal"href="dimension.html"><spanclass="doc">2d simulations</span></a>, atom styles sphere
and ellipsoid still use 3d particles, rather than as circular disks or
ellipses. This means they have the same moment of inertia as the 3d
object. When temperature is computed, the correct degrees of freedom
are used for rotation in a 2d versus 3d system.</p>
</div>
<divclass="section"id="pair-potentials">
<h3>6.14.2. Pair potentials</h3>
<p>When a system with finite-size particles is defined, the particles
will only rotate and experience torque if the force field computes
such interactions. These are the various <aclass="reference internal"href="pair_style.html"><spanclass="doc">pair styles</span></a> that generate torque:</p>
<h3>6.14.5. Rigid bodies composed of finite-size particles</h3>
<p>The <aclass="reference internal"href="fix_rigid.html"><spanclass="doc">fix rigid</span></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 finite-size
particles (spheres or ellipsoids or line segments or triangles), 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
spheroids, even if the two particles have the same mass. Finite-size
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 <aclass="reference internal"href="fix_shake.html"><spanclass="doc">fix shake</span></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>
<p>Also note that body particles cannot be modeled with the <aclass="reference internal"href="fix_rigid.html"><spanclass="doc">fix rigid</span></a> command. Body particles are treated by LAMMPS
as single particles, though they can store internal state, such as a
list of sub-particles. Individual body partices are typically treated
as rigid bodies, and their motion integrated with a command like <aclass="reference internal"href="fix_nve_body.html"><spanclass="doc">fix nve/body</span></a>. Interactions between pairs of body
particles are computed via a command like <aclass="reference internal"href="pair_body.html"><spanclass="doc">pair_style body</span></a>.</p>
<spanid="howto-15"></span><h2>6.15. Output from LAMMPS (thermo, dumps, computes, fixes, variables)</h2>
<p>There are four basic kinds of LAMMPS output:</p>
<ulclass="simple">
<li><aclass="reference internal"href="thermo_style.html"><spanclass="doc">Thermodynamic output</span></a>, which is a list
of quantities printed every few timesteps to the screen and logfile.</li>
<li><aclass="reference internal"href="dump.html"><spanclass="doc">Dump files</span></a>, which contain snapshots of atoms and various
per-atom values and are written at a specified frequency.</li>
<li>Certain fixes can output user-specified quantities to files: <aclass="reference internal"href="fix_ave_time.html"><spanclass="doc">fix ave/time</span></a> for time averaging, <aclass="reference internal"href="fix_ave_chunk.html"><spanclass="doc">fix ave/chunk</span></a> for spatial or other averaging, and <aclass="reference internal"href="fix_print.html"><spanclass="doc">fix print</span></a> for single-line output of
<aclass="reference internal"href="variable.html"><spanclass="doc">variables</span></a>. Fix print can also output to the
<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 <aclass="reference internal"href="dump.html"><spanclass="doc">dump</span></a> and <aclass="reference internal"href="fix.html"><spanclass="doc">fix</span></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 <aclass="reference internal"href="Section_modify.html"><spanclass="doc">add their own computes and fixes to LAMMPS</span></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>
<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
“<ahref="#id74"><spanclass="problematic"id="id75">c_</span></a>” would be replaced by “<ahref="#id76"><spanclass="problematic"id="id77">f_</span></a>” for a fix, or “<ahref="#id78"><spanclass="problematic"id="id79">v_</span></a>” for a variable:</p>
<tableborder="1"class="docutils">
<colgroup>
<colwidth="21%"/>
<colwidth="79%"/>
</colgroup>
<tbodyvalign="top">
<trclass="row-odd"><td>c_ID</td>
<td>entire scalar, vector, or array</td>
</tr>
<trclass="row-even"><td>c_ID[I]</td>
<td>one element of vector, one column of array</td>
</tr>
<trclass="row-odd"><td>c_ID[I][J]</td>
<td>one element of array</td>
</tr>
</tbody>
</table>
<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
<p>The frequency and format of thermodynamic output is set by the
<aclass="reference internal"href="thermo.html"><spanclass="doc">thermo</span></a>, <aclass="reference internal"href="thermo_style.html"><spanclass="doc">thermo_style</span></a>, and
<aclass="reference internal"href="thermo_modify.html"><spanclass="doc">thermo_modify</span></a> commands. The
<aclass="reference internal"href="thermo_style.html"><spanclass="doc">thermo_style</span></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 <aclass="reference internal"href="compute.html"><spanclass="doc">compute</span></a>
or <aclass="reference internal"href="fix.html"><spanclass="doc">fix</span></a> or <aclass="reference internal"href="variable.html"><spanclass="doc">variable</span></a> provides the value to be
output. In each case, the compute, fix, or variable must generate
global values for input to the <aclass="reference internal"href="dump.html"><spanclass="doc">thermo_style custom</span></a>
command.</p>
<p>Note that thermodynamic output values can be “extensive” or
“intensive”. The former scale with the number of atoms in the system
(e.g. total energy), the latter do not (e.g. temperature). The
setting for <aclass="reference internal"href="thermo_modify.html"><spanclass="doc">thermo_modify norm</span></a> determines whether
extensive quantities are normalized or not. Computes and fixes
produce either extensive or intensive values; see their individual doc
pages for details. <aclass="reference internal"href="variable.html"><spanclass="doc">Equal-style variables</span></a> produce only
intensive values; you can include a division by “natoms” in the
formula if desired, to make an extensive calculation produce an
<p>There is also a <aclass="reference internal"href="dump.html"><spanclass="doc">dump custom</span></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
<aclass="reference internal"href="compute.html"><spanclass="doc">compute</span></a> or <aclass="reference internal"href="fix.html"><spanclass="doc">fix</span></a> or <aclass="reference internal"href="variable.html"><spanclass="doc">variable</span></a>
provides the values to be output. In each case, the compute, fix, or
variable must generate per-atom values for input to the <aclass="reference internal"href="dump.html"><spanclass="doc">dump custom</span></a> command.</p>
<p>There is also a <aclass="reference internal"href="dump.html"><spanclass="doc">dump local</span></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
<aclass="reference internal"href="compute.html"><spanclass="doc">compute</span></a> or <aclass="reference internal"href="fix.html"><spanclass="doc">fix</span></a> or <aclass="reference internal"href="variable.html"><spanclass="doc">variable</span></a>
provides the values to be output. In each case, the compute or fix
must generate local values for input to the <aclass="reference internal"href="dump.html"><spanclass="doc">dump local</span></a>
<aclass="reference internal"href="fix_ave_correlate.html"><spanclass="doc">fix ave/correlate</span></a>, and <aclass="reference internal"href="fix_print.html"><spanclass="doc">fix print</span></a>.</p>
<p>The <aclass="reference internal"href="fix_ave_time.html"><spanclass="doc">fix ave/time</span></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
<aclass="reference internal"href="compute.html"><spanclass="doc">compute</span></a> values, global <aclass="reference internal"href="fix.html"><spanclass="doc">fix</span></a> values, or
<aclass="reference internal"href="variable.html"><spanclass="doc">variables</span></a> of any style except the atom style which
produces per-atom values. Since a variable can refer to keywords used
by the <aclass="reference internal"href="thermo_style.html"><spanclass="doc">thermo_style custom</span></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 <aclass="reference internal"href="fix_ave_chunk.html"><spanclass="doc">fix ave/chunk</span></a> command enables direct output
to a file of chunk-averaged per-atom quantities like those output in
dump files. Chunks can represent spatial bins or other collections of
atoms, e.g. individual molecules. 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
<aclass="reference internal"href="compute.html"><spanclass="doc">compute</span></a>, by a <aclass="reference internal"href="fix.html"><spanclass="doc">fix</span></a>, or by an atom-style
<aclass="reference internal"href="variable.html"><spanclass="doc">variable</span></a>. The chunk-averaged output of this fix can
also be used as input to other output commands.</p>
<p>The <aclass="reference internal"href="fix_ave_histo.html"><spanclass="doc">fix ave/histo</span></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 <aclass="reference internal"href="fix_ave_correlate.html"><spanclass="doc">fix ave/correlate</span></a> command enables direct
output to a file of time-correlated quantities, which can be global
values. The correlation matrix output of this fix can also be used as
input to other output commands.</p>
<p>The <aclass="reference internal"href="fix_print.html"><spanclass="doc">fix print</span></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
<aclass="reference internal"href="variable.html"><spanclass="doc">variable</span></a> values for any style variable except the
vector or atom styles). As explained above, variables themselves can
contain references to global values generated by <aclass="reference internal"href="thermo_style.html"><spanclass="doc">thermodynamic keywords</span></a>, <aclass="reference internal"href="compute.html"><spanclass="doc">computes</span></a>,
<aclass="reference internal"href="fix.html"><spanclass="doc">fixes</span></a>, or other <aclass="reference internal"href="variable.html"><spanclass="doc">variables</span></a>, or to per-atom
values for a specific atom. Thus the <aclass="reference internal"href="fix_print.html"><spanclass="doc">fix print</span></a>
command is a means to output a wide variety of quantities separate
from normal thermodynamic or dump file output.</p>
<spanid="computeoutput"></span><h3>6.15.6. Computes that process output quantities</h3>
<p>The <aclass="reference internal"href="compute_reduce.html"><spanclass="doc">compute reduce</span></a> and <aclass="reference internal"href="compute_reduce.html"><spanclass="doc">compute reduce/region</span></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 <aclass="reference internal"href="compute_slice.html"><spanclass="doc">compute slice</span></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 <aclass="reference internal"href="compute_property_atom.html"><spanclass="doc">compute property/atom</span></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 <aclass="reference internal"href="dump.html"><spanclass="doc">dump custom</span></a> command.</p>
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 <aclass="reference internal"href="compute.html"><spanclass="doc">compute</span></a>, by a
<aclass="reference internal"href="fix.html"><spanclass="doc">fix</span></a>, or by an atom-style <aclass="reference internal"href="variable.html"><spanclass="doc">variable</span></a>. The
time-averaged per-atom output of this fix can be used as input to
other output commands.</p>
<p>The <aclass="reference internal"href="fix_store_state.html"><spanclass="doc">fix store/state</span></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 <aclass="reference internal"href="dump.html"><spanclass="doc">dump custom</span></a> command,
including per-atom quantities calculated by a <aclass="reference internal"href="compute.html"><spanclass="doc">compute</span></a>,
by a <aclass="reference internal"href="fix.html"><spanclass="doc">fix</span></a>, or by an atom-style <aclass="reference internal"href="variable.html"><spanclass="doc">variable</span></a>.
The output of this fix can be used as input to other output commands.</p>
<spanid="howto-16"></span><h2>6.16. Thermostatting, barostatting, and computing temperature</h2>
<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 <aclass="reference internal"href="compute.html"><spanclass="doc">compute commands</span></a> calculate temperature, and the <aclass="reference internal"href="compute_pressure.html"><spanclass="doc">compute pressure</span></a> command calculates pressure.</p>
<p>All but the first 3 calculate velocity biases directly (e.g. advection
velocities) that are removed when computing the thermal temperature.
<aclass="reference internal"href="compute_temp_sphere.html"><spanclass="doc">Compute temp/sphere</span></a> and <aclass="reference internal"href="compute_temp_asphere.html"><spanclass="doc">compute temp/asphere</span></a> compute kinetic energy for
finite-size particles that includes rotational degrees of freedom.
They both allow for velocity biases indirectly, via an optional extra
argument, another temperature compute that subtracts a velocity bias.
This allows the translational velocity of spherical or aspherical
particles to be adjusted in prescribed ways.</p>
<p>Thermostatting in LAMMPS is performed by <aclass="reference internal"href="fix.html"><spanclass="doc">fixes</span></a>, or in one
case by a pair style. Several thermostatting fixes are available:
Nose-Hoover (nvt), Berendsen, CSVR, Langevin, and direct rescaling
<p><aclass="reference internal"href="fix_nh.html"><spanclass="doc">Fix nvt</span></a> only thermostats the translational velocity of
particles. <aclass="reference internal"href="fix_nvt_sllod.html"><spanclass="doc">Fix nvt/sllod</span></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 <aclass="reference internal"href="#howto-13"><spanclass="std std-ref">NEMD simulations</span></a> section of this page for further details. <aclass="reference internal"href="fix_nvt_sphere.html"><spanclass="doc">Fix nvt/sphere</span></a> and <aclass="reference internal"href="fix_nvt_asphere.html"><spanclass="doc">fix nvt/asphere</span></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 <aclass="reference internal"href="fix_langevin.html"><spanclass="doc">fix langevin</span></a>.</p>
<p>Any of the thermostatting fixes can use temperature computes that
remove bias which has two effects. First, the current calculated
temperature, which is compared to the requested target temperature, is
caluclated with the velocity bias removed. Second, the thermostat
adjusts only the thermal temperature component of the particle’s
velocities, which are the velocities with the bias removed. The
removed bias is then added back to the adjusted velocities. See the
doc pages for the individual fixes and for the
<aclass="reference internal"href="fix_modify.html"><spanclass="doc">fix_modify</span></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 <aclass="reference internal"href="compute_temp_partial.html"><spanclass="doc">compute temp/partial</span></a>. Of you could thermostat only
the thermal temperature of a streaming flow of particles without
affecting the streaming velocity, by using <aclass="reference internal"href="compute_temp_profile.html"><spanclass="doc">compute temp/profile</span></a>.</p>
<divclass="admonition note">
<pclass="first admonition-title">Note</p>
<pclass="last">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>
<p>The <aclass="reference internal"href="fix_nh.html"><spanclass="doc">fix npt</span></a> commands include a Nose-Hoover thermostat
and barostat. <aclass="reference internal"href="fix_nh.html"><spanclass="doc">Fix nph</span></a> is just a Nose/Hoover barostat;
it does no thermostatting. Both <aclass="reference internal"href="fix_nh.html"><spanclass="doc">fix nph</span></a> and <aclass="reference internal"href="fix_press_berendsen.html"><spanclass="doc">fix press/bernendsen</span></a> can be used in conjunction
with any of the thermostatting fixes.</p>
<p>As with the thermostats, <aclass="reference internal"href="fix_nh.html"><spanclass="doc">fix npt</span></a> and <aclass="reference internal"href="fix_nh.html"><spanclass="doc">fix nph</span></a> only use translational motion of the particles in
computing T and P and performing thermo/barostatting. <aclass="reference internal"href="fix_npt_sphere.html"><spanclass="doc">Fix npt/sphere</span></a> and <aclass="reference internal"href="fix_npt_asphere.html"><spanclass="doc">fix npt/asphere</span></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 <aclass="reference internal"href="compute_pressure.html"><spanclass="doc">compute pressure</span></a> compute to calculate a current
pressure. By default, this compute is created with a simple <aclass="reference internal"href="compute_temp.html"><spanclass="doc">compute temp</span></a> (see the last argument of the <aclass="reference internal"href="compute_pressure.html"><spanclass="doc">compute pressure</span></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
<aclass="reference internal"href="fix_modify.html"><spanclass="doc">fix_modify</span></a> command for instructions on how to assign
a temperature or pressure compute to a barostatting fix.</p>
<divclass="admonition note">
<pclass="first admonition-title">Note</p>
<pclass="last">As with the thermostats, the Nose/Hoover methods (<aclass="reference internal"href="fix_nh.html"><spanclass="doc">fix npt</span></a> and <aclass="reference internal"href="fix_nh.html"><spanclass="doc">fix nph</span></a>) perform time integration.
<aclass="reference internal"href="fix_press_berendsen.html"><spanclass="doc">Fix press/berendsen</span></a> does NOT, so it should
be used with one of the constant NVE fixes or with one of the NVT
fixes.</p>
</div>
<p>Finally, thermodynamic output, which can be setup via the
<aclass="reference internal"href="thermo_style.html"><spanclass="doc">thermo_style</span></a> command, often includes temperature
and pressure values. As explained on the doc page for the
<aclass="reference internal"href="thermo_style.html"><spanclass="doc">thermo_style</span></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 <aclass="reference internal"href="thermo_style.html"><spanclass="doc">thermo_style custom</span></a> command. Or
you can use the <aclass="reference internal"href="thermo_modify.html"><spanclass="doc">thermo_modify</span></a> command to
re-define what temperature or pressure compute is used for default
thermodynamic output.</p>
<hrclass="docutils"/>
</div>
<divclass="section"id="walls">
<spanid="howto-17"></span><h2>6.17. Walls</h2>
<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
<aclass="reference internal"href="lattice.html"><spanclass="doc">lattice</span></a> and <aclass="reference internal"href="create_atoms.html"><spanclass="doc">create_atoms</span></a> commands,
or read in via the <aclass="reference internal"href="read_data.html"><spanclass="doc">read_data</span></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 <aclass="reference internal"href="fix_nve.html"><spanclass="doc">fix nve</span></a> or <aclass="reference internal"href="fix_nh.html"><spanclass="doc">fix nvt</span></a>
is not used with the group that contains wall particles, their
positions and velocities will not be updated.</p>
<ulclass="simple">
<li><aclass="reference internal"href="fix_aveforce.html"><spanclass="doc">fix aveforce</span></a> - set force on particles to average value, so they move together</li>
<li><aclass="reference internal"href="fix_setforce.html"><spanclass="doc">fix setforce</span></a> - set force on particles to a value, e.g. 0.0</li>
<li><aclass="reference internal"href="fix_freeze.html"><spanclass="doc">fix freeze</span></a> - freeze particles for use as granular walls</li>
<li><aclass="reference internal"href="fix_nve_noforce.html"><spanclass="doc">fix nve/noforce</span></a> - advect particles by their velocity, but without force</li>
<li><aclass="reference internal"href="fix_move.html"><spanclass="doc">fix move</span></a> - prescribe motion of particles by a linear velocity, oscillation, rotation, variable</li>
</ul>
<p>The <aclass="reference internal"href="fix_move.html"><spanclass="doc">fix move</span></a> command offers the most generality, since
the motion of individual particles can be specified with
<aclass="reference internal"href="variable.html"><spanclass="doc">variable</span></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 <aclass="reference internal"href="neigh_modify.html"><spanclass="doc">neigh_modify exclude</span></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 <aclass="reference internal"href="bond_style.html"><spanclass="doc">bond</span></a>.
The bonded particles do interact with other mobile particles.</p>
<p>Idealized walls can be specified via several fix commands. <aclass="reference internal"href="fix_wall_gran.html"><spanclass="doc">Fix wall/gran</span></a> creates frictional walls for use with
granular particles; all the other commands create smooth walls.</p>
<p>The lammps_open() function is used to initialize LAMMPS, passing in a
list of strings as if they were <aclass="reference internal"href="Section_start.html#start-7"><spanclass="std std-ref">command-line arguments</span></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_version() function can be used to determined the specific
version of the underlying LAMMPS code. This is particularly useful
when loading LAMMPS as a shared library via dlopen(). The code using
the library interface can than use this information to adapt to
changes to the LAMMPS command syntax between versions. The returned
LAMMPS version code is an integer (e.g. 2 Sep 2015 results in
20150902) that grows with every new LAMMPS version.</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>
<p>These can extract various global or per-atom quantities from LAMMPS as
well as values calculated by a compute, fix, or variable. The
“set_variable” function can set an existing string-style variable to a
new value, so that subsequent LAMMPS commands can access the 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 <aclass="reference internal"href="Section_python.html"><spanclass="doc">Python interface</span></a>. The routines you add can access or
change any LAMMPS data you wish. The examples/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
<p>The thermal conductivity kappa of a material can be measured in at
least 4 ways using various options in LAMMPS. See the examples/KAPPA
directory for scripts that implement the 4 methods discussed here for
a simple Lennard-Jones fluid model. Also, see <aclass="reference internal"href="#howto-21"><spanclass="std std-ref">this section</span></a> of the manual for an analogous
discussion for viscosity.</p>
<p>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 <aclass="reference internal"href="#howto-13"><spanclass="std std-ref">thermostatting fix</span></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 papers by <aclass="reference internal"href="#howto-ikeshoji"><spanclass="std std-ref">Ikeshoji and Hafskjold</span></a>
and <aclass="reference internal"href="#howto-wirnsberger"><spanclass="std std-ref">Wirnsberger et al</span></a> for details of this idea.
Note that thermostatting fixes such as <aclass="reference internal"href="fix_nh.html"><spanclass="doc">fix nvt</span></a>, <aclass="reference internal"href="fix_langevin.html"><spanclass="doc">fix langevin</span></a>, and <aclass="reference internal"href="fix_temp_rescale.html"><spanclass="doc">fix temp/rescale</span></a> store the cumulative energy they
add/subtract.</p>
<p>Alternatively, as a second method, the <aclass="reference internal"href="fix_heat.html"><spanclass="doc">fix heat</span></a> or
<aclass="reference internal"href="fix_ehex.html"><spanclass="doc">fix ehex</span></a> commands can be used in place of thermostats
on each of two regions to add/subtract specified amounts of energy to
both regions. In both cases, the resulting temperatures of the two
regions can be monitored with the “compute temp/region” command and
the temperature profile of the intermediate region can be monitored
with the <spanclass="xref doc">fix ave/spatial</span> and <aclass="reference internal"href="compute_ke_atom.html"><spanclass="doc">compute ke/atom</span></a> commands.</p>
<p>The third method is to perform a reverse non-equilibrium MD simulation
using the <aclass="reference internal"href="fix_thermal_conductivity.html"><spanclass="doc">fix thermal/conductivity</span></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 <spanclass="xref doc">fix ave/spatial</span> and <aclass="reference internal"href="compute_ke_atom.html"><spanclass="doc">compute ke/atom</span></a> commands. The fix tallies the
cumulative energy transfer that it performs. See the <aclass="reference internal"href="fix_thermal_conductivity.html"><spanclass="doc">fix thermal/conductivity</span></a> command for
details.</p>
<p>The fourth 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 <aclass="reference internal"href="compute_heat_flux.html"><spanclass="doc">compute heat/flux</span></a> command can calculate
the needed heat flux and describes how to implement the Green_Kubo
formalism using additional LAMMPS commands, such as the <aclass="reference internal"href="fix_ave_correlate.html"><spanclass="doc">fix ave/correlate</span></a> command to calculate the needed
auto-correlation. See the doc page for the <aclass="reference internal"href="compute_heat_flux.html"><spanclass="doc">compute heat/flux</span></a> command for an example input script
that calculates the thermal conductivity of solid Ar via the GK
<p>The shear viscosity eta of a fluid can be measured in at least 5 ways
using various options in LAMMPS. See the examples/VISCOSITY directory
for scripts that implement the 5 methods discussed here for a simple
Lennard-Jones fluid model. Also, see <aclass="reference internal"href="#howto-20"><spanclass="std std-ref">this section</span></a> of the manual for an analogous
discussion for thermal conductivity.</p>
<p>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 <aclass="reference internal"href="fix_deform.html"><spanclass="doc">fix deform</span></a>
command, and using the <aclass="reference internal"href="fix_nvt_sllod.html"><spanclass="doc">fix nvt/sllod</span></a> command to
thermostat the fluid via the SLLOD equations of motion.
Alternatively, as a second method, one or more moving walls can be
used to shear the fluid in between them, again with some kind of
thermostat that modifies only the thermal (non-shearing) components of
velocity to prevent the fluid from heating up.</p>
<p>In both cases, the velocity profile setup in the fluid by this
procedure can be monitored by the <spanclass="xref doc">fix ave/spatial</span> 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. The Pxy off-diagonal component of the pressure or stress
tensor, as calculated by the <aclass="reference internal"href="compute_pressure.html"><spanclass="doc">compute pressure</span></a>
command, can also be monitored, which is the J term in the equation
above. See <aclass="reference internal"href="#howto-13"><spanclass="std std-ref">this section</span></a> of the manual
for details on NEMD simulations.</p>
<p>The third method is to perform a reverse non-equilibrium MD simulation
using the <aclass="reference internal"href="fix_viscosity.html"><spanclass="doc">fix viscosity</span></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 <spanclass="xref doc">fix ave/spatial</span> command.
The fix tallies the cummulative momentum transfer that it performs.
See the <aclass="reference internal"href="fix_viscosity.html"><spanclass="doc">fix viscosity</span></a> command for details.</p>
<p>The fourth 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 fully
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>
<divclass="highlight-default"><divclass="highlight"><pre><span></span><spanclass="c1"># Sample LAMMPS input script for viscosity of liquid Ar</span>
</pre></div>
</div>
<divclass="highlight-default"><divclass="highlight"><pre><span></span>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></div>
</div>
<divclass="highlight-default"><divclass="highlight"><pre><span></span><spanclass="c1"># convert from LAMMPS real units to SI</span>
<p>This compute can assign atoms to chunks of various styles. Only atoms
in the specified group and optional specified region are assigned to a
chunk. Here are some possible chunk definitions:</p>
<tableborder="1"class="docutils">
<colgroup>
<colwidth="31%"/>
<colwidth="69%"/>
</colgroup>
<tbodyvalign="top">
<trclass="row-odd"><td>atoms in same molecule</td>
<td>chunk ID = molecule ID</td>
</tr>
<trclass="row-even"><td>atoms of same atom type</td>
<td>chunk ID = atom type</td>
</tr>
<trclass="row-odd"><td>all atoms with same atom property (charge, radius, etc)</td>
<td>chunk ID = output of compute property/atom</td>
</tr>
<trclass="row-even"><td>atoms in same cluster</td>
<td>chunk ID = output of <aclass="reference internal"href="compute_cluster_atom.html"><spanclass="doc">compute cluster/atom</span></a> command</td>
</tr>
<trclass="row-odd"><td>atoms in same spatial bin</td>
<td>chunk ID = bin ID</td>
</tr>
<trclass="row-even"><td>atoms in same rigid body</td>
<td>chunk ID = molecule ID used to define rigid bodies</td>
</tr>
<trclass="row-odd"><td>atoms with similar potential energy</td>
<td>chunk ID = output of <aclass="reference internal"href="compute_pe_atom.html"><spanclass="doc">compute pe/atom</span></a></td>
</tr>
<trclass="row-even"><td>atoms with same local defect structure</td>
<td>chunk ID = output of <aclass="reference internal"href="compute_centro_atom.html"><spanclass="doc">compute centro/atom</span></a> or <aclass="reference internal"href="compute_coord_atom.html"><spanclass="doc">compute coord/atom</span></a> command</td>
</tr>
</tbody>
</table>
<p>Note that chunk IDs are integer values, so for atom properties or
computes that produce a floating point value, they will be truncated
to an integer. You could also use the compute in a variable that
scales the floating point value to spread it across multiple intergers.</p>
<p>Spatial bins can be of various kinds, e.g. 1d bins = slabs, 2d bins =
pencils, 3d bins = boxes, spherical bins, cylindrical bins.</p>
<p>This compute also calculates the number of chunks <em>Nchunk</em>, which is
used by other commands to tally per-chunk data. <em>Nchunk</em> can be a
static value or change over time (e.g. the number of clusters). The
chunk ID for an individual atom can also be static (e.g. a molecule
ID), or dynamic (e.g. what spatial bin an atom is in as it moves).</p>
<p>Note that this compute allows the per-atom output of other
<aclass="reference internal"href="compute.html"><spanclass="doc">computes</span></a>, <aclass="reference internal"href="fix.html"><spanclass="doc">fixes</span></a>, and
<aclass="reference internal"href="variable.html"><spanclass="doc">variables</span></a> to be used to define chunk IDs for each
atom. This means you can write your own compute or fix to output a
per-atom quantity to use as chunk ID. See
<aclass="reference internal"href="Section_modify.html"><spanclass="doc">Section_modify</span></a> of the documentation for how to
do this. You can also define a <aclass="reference internal"href="variable.html"><spanclass="doc">per-atom variable</span></a> in
the input script that uses a formula to generate a chunk ID for each
atom.</p>
</div>
<divclass="section"id="fix-ave-chunk-command">
<h3>6.23.2. Fix ave/chunk command:</h3>
<p>This fix takes the ID of a <aclass="reference internal"href="compute_chunk_atom.html"><spanclass="doc">compute chunk/atom</span></a> command as input. For each chunk,
it then sums one or more specified per-atom values over the atoms in
each chunk. The per-atom values can be any atom property, such as
etc. Additional keywords are defined for per-chunk properties like
density and temperature. More generally any per-atom value generated
by other <aclass="reference internal"href="compute.html"><spanclass="doc">computes</span></a>, <aclass="reference internal"href="fix.html"><spanclass="doc">fixes</span></a>, and <aclass="reference internal"href="variable.html"><spanclass="doc">per-atom variables</span></a>, can be summed over atoms in each chunk.</p>
<p>Similar to other averaging fixes, this fix allows the summed per-chunk
values to be time-averaged in various ways, and output to a file. The
fix produces a global array as output with one row of values per
<p>They each take the ID of a <aclass="reference internal"href="compute_chunk_atom.html"><spanclass="doc">compute chunk/atom</span></a> command as input. As their names
indicate, they calculate the center-of-mass, radius of gyration,
moments of inertia, mean-squared displacement, temperature, torque,
and velocity of center-of-mass for each chunk of atoms. The <aclass="reference internal"href="compute_property_chunk.html"><spanclass="doc">compute property/chunk</span></a> command can tally the
count of atoms in each chunk and extract other per-chunk properties.</p>
<p>The reason these various calculations are not part of the <aclass="reference internal"href="fix_ave_chunk.html"><spanclass="doc">fix ave/chunk command</span></a>, is that each requires a more
complicated operation than simply summing and averaging over per-atom
values in each chunk. For example, many of them require calculation
of a center of mass, which requires summing mass*position over the
atoms and then dividing by summed mass.</p>
<p>All of these computes produce a global vector or global array as
output, wih one or more values per chunk. They can be used
in various ways:</p>
<ulclass="simple">
<li>As input to the <aclass="reference internal"href="fix_ave_time.html"><spanclass="doc">fix ave/time</span></a> command, which can
write the values to a file and optionally time average them.</li>
<li>As input to the <aclass="reference internal"href="fix_ave_histo.html"><spanclass="doc">fix ave/histo</span></a> command to
histogram values across chunks. E.g. a histogram of cluster sizes or
molecule diffusion rates.</li>
<li>As input to special functions of <aclass="reference internal"href="variable.html"><spanclass="doc">equal-style variables</span></a>, like sum() and max(). E.g. to find the
largest cluster or fastest diffusing molecule.</li>
<spanid="howto-24"></span><h2>6.24. Setting parameters for the <aclass="reference internal"href="kspace_style.html"><spanclass="doc">kspace_style pppm/disp</span></a> command</h2>
<p>The PPPM method computes interactions by splitting the pair potential
into two parts, one of which is computed in a normal pairwise fashion,
the so-called real-space part, and one of which is computed using the
Fourier transform, the so called reciprocal-space or kspace part. For
both parts, the potential is not computed exactly but is approximated.
Thus, there is an error in both parts of the computation, the
real-space and the kspace error. The just mentioned facts are true
both for the PPPM for Coulomb as well as dispersion interactions. The
deciding difference - and also the reason why the parameters for
pppm/disp have to be selected with more care - is the impact of the
errors on the results: The kspace error of the PPPM for Coulomb and
dispersion interaction and the real-space error of the PPPM for
Coulomb interaction have the character of noise. In contrast, the
real-space error of the PPPM for dispersion has a clear physical
interpretation: the underprediction of cohesion. As a consequence, the
real-space error has a much stronger effect than the kspace error on
simulation results for pppm/disp. Parameters must thus be chosen in a
way that this error is much smaller than the kspace error.</p>
<p>When using pppm/disp and not making any specifications on the PPPM
parameters via the kspace modify command, parameters will be tuned
such that the real-space error and the kspace error are equal. This
will result in simulations that are either inaccurate or slow, both of
which is not desirable. For selecting parameters for the pppm/disp
that provide fast and accurate simulations, there are two approaches,
which both have their up- and downsides.</p>
<p>The first approach is to set desired real-space an kspace accuracies
via the <em>kspace_modify force/disp/real</em> and <em>kspace_modify
force/disp/kspace</em> commands. Note that the accuracies have to be
specified in force units and are thus dependend on the chosen unit
settings. For real units, 0.0001 and 0.002 seem to provide reasonable
accurate and efficient computations for the real-space and kspace
accuracies. 0.002 and 0.05 work well for most systems using lj
units. PPPM parameters will be generated based on the desired
accuracies. The upside of this approach is that it usually provides a
good set of parameters and will work for both the <em>kspace_modify diff
ad</em> and <em>kspace_modify diff ik</em> options. The downside of the method
is that setting the PPPM parameters will take some time during the
initialization of the simulation.</p>
<p>The second approach is to set the parameters for the pppm/disp
explicitly using the <em>kspace_modify mesh/disp</em>, <em>kspace_modify
order/disp</em>, and <em>kspace_modify gewald/disp</em> commands. This approach
requires a more experienced user who understands well the impact of
the choice of parameters on the simulation accuracy and
performance. This approach provides a fast initialization of the
simulation. However, it is sensitive to errors: A combination of
parameters that will perform well for one system might result in
far-from-optimal conditions for other simulations. For example,
parametes that provide accurate and fast computations for
all-atomistic force fields can provide insufficient accuracy or
united-atomistic force fields (which is related to that the latter
typically have larger dispersion coefficients).</p>
<p>To avoid inaccurate or inefficient simulations, the pppm/disp stops
simulations with an error message if no action is taken to control the
PPPM parameters. If the automatic parameter generation is desired and
real-space and kspace accuracies are desired to be equal, this error
message can be suppressed using the <em>kspace_modify disp/auto yes</em>
command.</p>
<p>A reasonable approach that combines the upsides of both methods is to
make the first run using the <em>kspace_modify force/disp/real</em> and
<em>kspace_modify force/disp/kspace</em> commands, write down the PPPM
parameters from the outut, and specify these parameters using the
second approach in subsequent runs (which have the same composition,
force field, and approximately the same volume).</p>
<p>Concerning the performance of the pppm/disp there are two more things
to consider. The first is that when using the pppm/disp, the cutoff
parameter does no longer affect the accuracy of the simulation
(subject to that gewald/disp is adjusted when changing the cutoff).
The performance can thus be increased by examining different values
for the cutoff parameter. A lower bound for the cutoff is only set by
the truncation error of the repulsive term of pair potentials.</p>
<p>The second is that the mixing rule of the pair style has an impact on
the computation time when using the pppm/disp. Fastest computations
are achieved when using the geometric mixing rule. Using the
arithmetic mixing rule substantially increases the computational cost.
The computational overhead can be reduced using the <em>kspace_modify
mix/disp geom</em> and <em>kspace_modify splittol</em> commands. The first
command simply enforces geometric mixing of the dispersion
coeffiecients in kspace computations. This introduces some error in
the computations but will also significantly speed-up the
simulations. The second keyword sets the accuracy with which the
dispersion coefficients are approximated using a matrix factorization
approach. This may result in better accuracy then using the first
command, but will usually also not provide an equally good increase of
efficiency.</p>
<p>Finally, pppm/disp can also be used when no mixing rules apply.
This can be achieved using the <em>kspace_modify mix/disp none</em> command.
Note that the code does not check automatically whether any mixing
rule is fulfilled. If mixing rules do not apply, the user will have
<p>The adiabatic core-shell model by <aclass="reference internal"href="pair_cs.html#mitchellfinchham"><spanclass="std std-ref">Mitchell and Finchham</span></a> is a simple method for adding
polarizability to a system. In order to mimic the electron shell of
an ion, a satellite particle is attached to it. This way the ions are
split into a core and a shell where the latter is meant to react to
the electrostatic environment inducing polarizability.</p>
<p>Technically, shells are attached to the cores by a spring force f =
k*r where k is a parametrized spring constant and r is the distance
between the core and the shell. The charges of the core and the shell
add up to the ion charge, thus q(ion) = q(core) + q(shell). This
setup introduces the ion polarizability (alpha) given by
alpha = q(shell)^2 / k. In a
similar fashion the mass of the ion is distributed on the core and the
shell with the core having the larger mass.</p>
<p>To run this model in LAMMPS, <aclass="reference internal"href="atom_style.html"><spanclass="doc">atom_style</span></a><em>full</em> can
be used since atom charge and bonds are needed. Each kind of
core/shell pair requires two atom types and a bond type. The core and
shell of a core/shell pair should be bonded to each other with a
harmonic bond that provides the spring force. For example, a data file
for NaCl, as found in examples/coreshell, has this format:</p>
<divclass="highlight-default"><divclass="highlight"><pre><span></span><spanclass="mi">432</span><spanclass="n">atoms</span><spanclass="c1"># core and shell atoms</span>
<spanclass="mi">216</span><spanclass="n">bonds</span><spanclass="c1"># number of core/shell springs</span>
</pre></div>
</div>
<divclass="highlight-default"><divclass="highlight"><pre><span></span><spanclass="mi">4</span><spanclass="n">atom</span><spanclass="n">types</span><spanclass="c1"># 2 cores and 2 shells for Na and Cl</span>
<divclass="highlight-default"><divclass="highlight"><pre><span></span><spanclass="n">Masses</span><spanclass="c1"># core/shell mass ratio = 0.1</span>
</pre></div>
</div>
<divclass="highlight-default"><divclass="highlight"><pre><span></span><spanclass="mi">1</span><spanclass="mf">20.690784</span><spanclass="c1"># Na core</span>
<divclass="highlight-default"><divclass="highlight"><pre><span></span><spanclass="n">Bonds</span><spanclass="c1"># Bond topology for spring forces</span>
</pre></div>
</div>
<divclass="highlight-default"><divclass="highlight"><pre><span></span><spanclass="mi">1</span><spanclass="mi">2</span><spanclass="mi">1</span><spanclass="mi">2</span><spanclass="c1"># spring for core/shell pair 1</span>
<spanclass="mi">2</span><spanclass="mi">2</span><spanclass="mi">3</span><spanclass="mi">4</span><spanclass="c1"># spring for core/shell pair 2</span>
<p>Non-Coulombic (e.g. Lennard-Jones) pairwise interactions are only
defined between the shells. Coulombic interactions are defined
between all cores and shells. If desired, additional bonds can be
specified between cores.</p>
<p>The <aclass="reference internal"href="special_bonds.html"><spanclass="doc">special_bonds</span></a> command should be used to
turn-off the Coulombic interaction within core/shell pairs, since that
interaction is set by the bond spring. This is done using the
<aclass="reference internal"href="special_bonds.html"><spanclass="doc">special_bonds</span></a> command with a 1-2 weight = 0.0,
which is the default value. It needs to be considered whether one has
to adjust the <aclass="reference internal"href="special_bonds.html"><spanclass="doc">special_bonds</span></a> weighting according
to the molecular topology since the interactions of the shells are
bypassed over an extra bond.</p>
<p>Note that this core/shell implementation does not require all ions to
be polarized. One can mix core/shell pairs and ions without a
satellite particle if desired.</p>
<p>Since the core/shell model permits distances of r = 0.0 between the
core and shell, a pair style with a “cs” suffix needs to be used to
implement a valid long-range Coulombic correction. Several such pair
styles are provided in the CORESHELL package. See <aclass="reference internal"href="pair_cs.html"><spanclass="doc">this doc page</span></a> for details. All of the core/shell enabled pair
styles require the use of a long-range Coulombic solver, as specified
by the <aclass="reference internal"href="kspace_style.html"><spanclass="doc">kspace_style</span></a> command. Either the PPPM or
Ewald solvers can be used.</p>
<p>For the NaCL example problem, these pair style and bond style settings
<p>When running dynamics with the adiabatic core/shell model, the
following issues should be considered. Since the relative motion of
the core and shell particles corresponds to the polarization, typical
thermostats can alter the polarization behaviour, meaning the shell
will not react freely to its electrostatic environment. This is
critical during the equilibration of the system. Therefore
it’s typically desirable to decouple the relative motion of the
core/shell pair, which is an imaginary degree of freedom, from the
real physical system. To do that, the <aclass="reference internal"href="compute_temp_cs.html"><spanclass="doc">compute temp/cs</span></a> command can be used, in conjunction with
any of the thermostat fixes, such as <aclass="reference internal"href="fix_nh.html"><spanclass="doc">fix nvt</span></a> or <aclass="reference external"href="fix_langevin">fix langevin</a>. This compute uses the center-of-mass velocity
of the core/shell pairs to calculate a temperature, and insures that
velocity is what is rescaled for thermostatting purposes. This
compute also works for a system with both core/shell pairs and
non-polarized ions (ions without an attached satellite particle). The
<aclass="reference internal"href="compute_temp_cs.html"><spanclass="doc">compute temp/cs</span></a> command requires input of two
groups, one for the core atoms, another for the shell atoms.
Non-polarized ions which might also be included in the treated system
should not be included into either of these groups, they are taken
into account by the <em>group-ID</em> (2nd argument) of the compute. The
groups can be defined using the <aclass="reference internal"href="group.html"><spanclass="doc">group *type*</span></a> command.
Note that to perform thermostatting using this definition of
temperature, the <aclass="reference internal"href="fix_modify.html"><spanclass="doc">fix modify temp</span></a> command should be
used to assign the compute to the thermostat fix. Likewise the
<aclass="reference internal"href="thermo_modify.html"><spanclass="doc">thermo_modify temp</span></a> command can be used to make
this temperature be output for the overall system.</p>
<p>For the NaCl example, this can be done as follows:</p>
<spanclass="n">fix</span><spanclass="n">thermoberendsen</span><spanclass="nb">all</span><spanclass="n">temp</span><spanclass="o">/</span><spanclass="n">berendsen</span><spanclass="mi">1427</span><spanclass="mi">1427</span><spanclass="mf">0.4</span><spanclass="c1"># thermostat for the true physical system</span>
<spanclass="n">fix</span><spanclass="n">thermostatequ</span><spanclass="nb">all</span><spanclass="n">nve</span><spanclass="c1"># integrator as needed for the berendsen thermostat</span>
<spanclass="n">thermo_modify</span><spanclass="n">temp</span><spanclass="n">CSequ</span><spanclass="c1"># output of center-of-mass derived temperature</span>
</pre></div>
</div>
<p>If <aclass="reference internal"href="compute_temp_cs.html"><spanclass="doc">compute temp/cs</span></a> is used, the decoupled
relative motion of the core and the shell should in theory be
stable. However numerical fluctuation can introduce a small
momentum to the system, which is noticable over long trajectories.
Therefore it is recomendable to use the <aclass="reference internal"href="fix_momentum.html"><spanclass="doc">fix momentum</span></a> command in combination with <aclass="reference internal"href="compute_temp_cs.html"><spanclass="doc">compute temp/cs</span></a> when equilibrating the system to
prevent any drift.</p>
<p>When intializing the velocities of a system with core/shell pairs, it
is also desirable to not introduce energy into the relative motion of
the core/shell particles, but only assign a center-of-mass velocity to
the pairs. This can be done by using the <em>bias</em> keyword of the
<aclass="reference internal"href="velocity.html"><spanclass="doc">velocity create</span></a> command and assigning the <aclass="reference internal"href="compute_temp_cs.html"><spanclass="doc">compute temp/cs</span></a> command to the <em>temp</em> keyword of the
<p>It is important to note that the polarizability of the core/shell
pairs is based on their relative motion. Therefore the choice of
spring force and mass ratio need to ensure much faster relative motion
of the 2 atoms within the core/shell pair than their center-of-mass
velocity. This allow the shells to effectively react instantaneously
to the electrostatic environment. This fast movement also limits the
timestep size that can be used.</p>
<p>The primary literature of the adiabatic core/shell model suggests that
the fast relative motion of the core/shell pairs only allows negligible
energy transfer to the environment. Therefore it is not intended to
decouple the core/shell degree of freedom from the physical system
during production runs. In other words, the <aclass="reference internal"href="compute_temp_cs.html"><spanclass="doc">compute temp/cs</span></a> command should not be used during
production runs and is only required during equilibration. This way one
is consistent with literature (based on the code packages DL_POLY or
GULP for instance).</p>
<p>The mentioned energy transfer will typically lead to a a small drift
in total energy over time. This internal energy can be monitored
using the <aclass="reference internal"href="compute_chunk_atom.html"><spanclass="doc">compute chunk/atom</span></a> and <aclass="reference internal"href="compute_temp_chunk.html"><spanclass="doc">compute temp/chunk</span></a> commands. The internal kinetic
energies of each core/shell pair can then be summed using the sum()
special function of the <aclass="reference internal"href="variable.html"><spanclass="doc">variable</span></a> command. Or they can
be time/averaged and output using the <aclass="reference internal"href="fix_ave_time.html"><spanclass="doc">fix ave/time</span></a>
command. To use these commands, each core/shell pair must be defined
as a “chunk”. If each core/shell pair is defined as its own molecule,
the molecule ID can be used to define the chunks. If cores are bonded
to each other to form larger molecules, the chunks can be identified
by the <aclass="reference internal"href="fix_property_atom.html"><spanclass="doc">fix property/atom</span></a> via assigning a
core/shell ID to each atom using a special field in the data file read
by the <aclass="reference internal"href="read_data.html"><spanclass="doc">read_data</span></a> command. This field can then be
accessed by the <aclass="reference internal"href="compute_property_atom.html"><spanclass="doc">compute property/atom</span></a>
command, to use as input to the <aclass="reference internal"href="compute_chunk_atom.html"><spanclass="doc">compute chunk/atom</span></a> command to define the core/shell
<spanclass="n">read_data</span><spanclass="n">NaCl_CS_x0</span><spanclass="o">.</span><spanclass="mi">1</span><spanclass="n">_prop</span><spanclass="o">.</span><spanclass="n">data</span><spanclass="n">fix</span><spanclass="n">csinfo</span><spanclass="n">NULL</span><spanclass="n">CS</span><spanclass="o">-</span><spanclass="n">Info</span><spanclass="c1"># atom property added in the data-file</span>
<spanclass="n">compute</span><spanclass="n">cstherm</span><spanclass="nb">all</span><spanclass="n">temp</span><spanclass="o">/</span><spanclass="n">chunk</span><spanclass="n">cs_chunk</span><spanclass="n">temp</span><spanclass="n">internal</span><spanclass="n">com</span><spanclass="n">yes</span><spanclass="n">cdof</span><spanclass="mf">3.0</span><spanclass="c1"># note the chosen degrees of freedom for the core/shell pairs</span>
<p>The additional section in the date file would be formatted like this:</p>
<divclass="highlight-default"><divclass="highlight"><pre><span></span><spanclass="n">CS</span><spanclass="o">-</span><spanclass="n">Info</span><spanclass="c1"># header of additional section</span>
<p>The thermalized Drude model, similarly to the <aclass="reference internal"href="#howto-26"><spanclass="std std-ref">core-shell</span></a>
model, representes induced dipoles by a pair of charges (the core atom
and the Drude particle) connected by a harmonic spring. The Drude
model has a number of features aimed at its use in molecular systems
(<aclass="reference internal"href="#howto-lamoureux"><spanclass="std std-ref">Lamoureux and Roux</span></a>):</p>
<ulclass="simple">
<li>Thermostating of the additional degrees of freedom associated with the
induced dipoles at very low temperature, in terms of the reduced
coordinates of the Drude particles with respect to their cores. This
makes the trajectory close to that of relaxed induced dipoles.</li>
<li>Consistent definition of 1-2 to 1-4 neighbors. A core-Drude particle
pair represents a single (polarizable) atom, so the special screening
factors in a covalent structure should be the same for the core and
the Drude particle. Drude particles have to inherit the 1-2, 1-3, 1-4
special neighbor relations from their respective cores.</li>
<li>Stabilization of the interactions between induced dipoles. Drude
dipoles on covalently bonded atoms interact too strongly due to the
short distances, so an atom may capture the Drude particle of a
neighbor, or the induced dipoles within the same molecule may align
too much. To avoid this, damping at short range can be done by Thole
functions (for which there are physical grounds). This Thole damping
is applied to the point charges composing the induced dipole (the
charge of the Drude particle and the opposite charge on the core, not
to the total charge of the core atom).</li>
</ul>
<p>A detailed tutorial covering the usage of Drude induced dipoles in
LAMMPS is <aclass="reference internal"href="tutorial_drude.html"><spanclass="doc">available here</span></a>.</p>
<p>As with the core-shell model, the cores and Drude particles should
appear in the data file as standard atoms. The same holds for the
springs between them, which are described by standard harmonic bonds.
The nature of the atoms (core, Drude particle or non-polarizable) is
specified via the <aclass="reference internal"href="fix_drude.html"><spanclass="doc">fix drude</span></a> command. The special
list of neighbors is automatically refactored to account for the
equivalence of core and Drude particles as regards special 1-2 to 1-4
screening. It may be necessary to use the <em>extra</em> keyword of the
<aclass="reference internal"href="special_bonds.html"><spanclass="doc">special_bonds</span></a> command. If using <aclass="reference internal"href="fix_shake.html"><spanclass="doc">fix shake</span></a>, make sure no Drude particle is in this fix
group.</p>
<p>There are two ways to thermostat the Drude particles at a low
temperature: use either <aclass="reference internal"href="fix_langevin_drude.html"><spanclass="doc">fix langevin/drude</span></a>
for a Langevin thermostat, or <aclass="reference internal"href="fix_drude_transform.html"><spanclass="doc">fix drude/transform/*</span></a> for a Nose-Hoover
thermostat. The former requires use of the command <aclass="reference internal"href="comm_modify.html"><spanclass="doc">comm_modify vel yes</span></a>. The latter requires two separate integration
fixes like <em>nvt</em> or <em>npt</em>. The correct temperatures of the reduced
degrees of freedom can be calculated using the <aclass="reference internal"href="compute_temp_drude.html"><spanclass="doc">compute temp/drude</span></a>. This requires also to use the
command <em>comm_modify vel yes</em>.</p>
<p>Short-range damping of the induced dipole interactions can be achieved
using Thole functions through the the <aclass="reference internal"href="pair_thole.html"><spanclass="doc">pair style thole</span></a> in <aclass="reference internal"href="pair_hybrid.html"><spanclass="doc">pair_style hybrid/overlay</span></a>
with a Coulomb pair style. It may be useful to use <em>coul/long/cs</em> or
similar from the CORESHELL package if the core and Drude particle come
Built with <ahref="http://sphinx-doc.org/">Sphinx</a> using a <ahref="https://github.com/snide/sphinx_rtd_theme">theme</a> provided by <ahref="https://readthedocs.org">Read the Docs</a>.