Section_howto.html
1d3f94fbb20clammps-icms

Section_howto.html
View Options



			<!DOCTYPE html>
			<!--[if IE 8]><html class="no-js lt-ie9" lang="en" > <![endif]-->
			<!--[if gt IE 8]><!--> <html class="no-js" lang="en" > <!--<![endif]-->
			<head>
			<meta charset="utf-8">

			<meta name="viewport" content="width=device-width, initial-scale=1.0">

			<title>6. How-to discussions — LAMMPS documentation</title>















			<link rel="stylesheet" href="_static/css/theme.css" type="text/css" />



			<link rel="stylesheet" href="_static/sphinxcontrib-images/LightBox2/lightbox2/css/lightbox.css" type="text/css" />



			<link rel="top" title="LAMMPS documentation" href="index.html"/>
			<link rel="next" title="7. Example problems" href="Section_example.html"/>
			<link rel="prev" title="5. Accelerating LAMMPS performance" href="Section_accelerate.html"/>


			<script src="_static/js/modernizr.min.js"></script>

			</head>

			<body class="wy-body-for-nav" role="document">

			<div class="wy-grid-for-nav">


			<nav data-toggle="wy-nav-shift" class="wy-nav-side">
			<div class="wy-side-nav-search">



			<a href="Manual.html" class="icon icon-home"> LAMMPS



			</a>


			<div role="search">
			<form id="rtd-search-form" class="wy-form" action="search.html" method="get">
			<input type="text" name="q" placeholder="Search docs" />
			<input type="hidden" name="check_keywords" value="yes" />
			<input type="hidden" name="area" value="default" />
			</form>
			</div>


			</div>

			<div class="wy-menu wy-menu-vertical" data-spy="affix" role="navigation" aria-label="main navigation">



			<ul class="current">
			<li class="toctree-l1"><a class="reference internal" href="Section_intro.html">1. Introduction</a></li>
			<li class="toctree-l1"><a class="reference internal" href="Section_start.html">2. Getting Started</a></li>
			<li class="toctree-l1"><a class="reference internal" href="Section_commands.html">3. Commands</a></li>
			<li class="toctree-l1"><a class="reference internal" href="Section_packages.html">4. Packages</a></li>
			<li class="toctree-l1"><a class="reference internal" href="Section_accelerate.html">5. Accelerating LAMMPS performance</a></li>
			<li class="toctree-l1 current"><a class="current reference internal" href="#">6. How-to discussions</a><ul>
			<li class="toctree-l2"><a class="reference internal" href="#restarting-a-simulation">6.1. Restarting a simulation</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#d-simulations">6.2. 2d simulations</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#charmm-amber-and-dreiding-force-fields">6.3. CHARMM, AMBER, and DREIDING force fields</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#running-multiple-simulations-from-one-input-script">6.4. Running multiple simulations from one input script</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#multi-replica-simulations">6.5. Multi-replica simulations</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#granular-models">6.6. Granular models</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#tip3p-water-model">6.7. TIP3P water model</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#tip4p-water-model">6.8. TIP4P water model</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#spc-water-model">6.9. SPC water model</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#coupling-lammps-to-other-codes">6.10. Coupling LAMMPS to other codes</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#visualizing-lammps-snapshots">6.11. Visualizing LAMMPS snapshots</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#triclinic-non-orthogonal-simulation-boxes">6.12. Triclinic (non-orthogonal) simulation boxes</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#nemd-simulations">6.13. NEMD simulations</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#finite-size-spherical-and-aspherical-particles">6.14. Finite-size spherical and aspherical particles</a><ul>
			<li class="toctree-l3"><a class="reference internal" href="#atom-styles">6.14.1. Atom styles</a></li>
			<li class="toctree-l3"><a class="reference internal" href="#pair-potentials">6.14.2. Pair potentials</a></li>
			<li class="toctree-l3"><a class="reference internal" href="#time-integration">6.14.3. Time integration</a></li>
			<li class="toctree-l3"><a class="reference internal" href="#computes-thermodynamics-and-dump-output">6.14.4. Computes, thermodynamics, and dump output</a></li>
			<li class="toctree-l3"><a class="reference internal" href="#rigid-bodies-composed-of-finite-size-particles">6.14.5. Rigid bodies composed of finite-size particles</a></li>
			</ul>
			</li>
			<li class="toctree-l2"><a class="reference internal" href="#output-from-lammps-thermo-dumps-computes-fixes-variables">6.15. Output from LAMMPS (thermo, dumps, computes, fixes, variables)</a><ul>
			<li class="toctree-l3"><a class="reference internal" href="#global-per-atom-local-data">6.15.1. Global/per-atom/local data</a></li>
			<li class="toctree-l3"><a class="reference internal" href="#scalar-vector-array-data">6.15.2. Scalar/vector/array data</a></li>
			<li class="toctree-l3"><a class="reference internal" href="#thermodynamic-output">6.15.3. Thermodynamic output</a></li>
			<li class="toctree-l3"><a class="reference internal" href="#dump-file-output">6.15.4. Dump file output</a></li>
			<li class="toctree-l3"><a class="reference internal" href="#fixes-that-write-output-files">6.15.5. Fixes that write output files</a></li>
			<li class="toctree-l3"><a class="reference internal" href="#computes-that-process-output-quantities">6.15.6. Computes that process output quantities</a></li>
			<li class="toctree-l3"><a class="reference internal" href="#fixes-that-process-output-quantities">6.15.7. Fixes that process output quantities</a></li>
			<li class="toctree-l3"><a class="reference internal" href="#computes-that-generate-values-to-output">6.15.8. Computes that generate values to output</a></li>
			<li class="toctree-l3"><a class="reference internal" href="#fixes-that-generate-values-to-output">6.15.9. Fixes that generate values to output</a></li>
			<li class="toctree-l3"><a class="reference internal" href="#variables-that-generate-values-to-output">6.15.10. Variables that generate values to output</a></li>
			<li class="toctree-l3"><a class="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>
			<li class="toctree-l2"><a class="reference internal" href="#thermostatting-barostatting-and-computing-temperature">6.16. Thermostatting, barostatting, and computing temperature</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#walls">6.17. Walls</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#elastic-constants">6.18. Elastic constants</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#library-interface-to-lammps">6.19. Library interface to LAMMPS</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#calculating-thermal-conductivity">6.20. Calculating thermal conductivity</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#calculating-viscosity">6.21. Calculating viscosity</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#calculating-a-diffusion-coefficient">6.22. Calculating a diffusion coefficient</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#using-chunks-to-calculate-system-properties">6.23. Using chunks to calculate system properties</a><ul>
			<li class="toctree-l3"><a class="reference internal" href="#compute-chunk-atom-command">6.23.1. Compute chunk/atom command:</a></li>
			<li class="toctree-l3"><a class="reference internal" href="#fix-ave-chunk-command">6.23.2. Fix ave/chunk command:</a></li>
			<li class="toctree-l3"><a class="reference internal" href="#compute-chunk-commands">6.23.3. Compute */chunk commands:</a></li>
			<li class="toctree-l3"><a class="reference internal" href="#example-calculations-with-chunks">6.23.4. Example calculations with chunks</a></li>
			</ul>
			</li>
			<li class="toctree-l2"><a class="reference internal" href="#setting-parameters-for-the-kspace-style-pppm-disp-command">6.24. Setting parameters for the <code class="docutils literal"><span class="pre">kspace_style</span> <span class="pre">pppm/disp</span></code> command</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#polarizable-models">6.25. Polarizable models</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#adiabatic-core-shell-model">6.26. Adiabatic core/shell model</a></li>
			<li class="toctree-l2"><a class="reference internal" href="#drude-induced-dipoles">6.27. Drude induced dipoles</a></li>
			</ul>
			</li>
			<li class="toctree-l1"><a class="reference internal" href="Section_example.html">7. Example problems</a></li>
			<li class="toctree-l1"><a class="reference internal" href="Section_perf.html">8. Performance & scalability</a></li>
			<li class="toctree-l1"><a class="reference internal" href="Section_tools.html">9. Additional tools</a></li>
			<li class="toctree-l1"><a class="reference internal" href="Section_modify.html">10. Modifying & extending LAMMPS</a></li>
			<li class="toctree-l1"><a class="reference internal" href="Section_python.html">11. Python interface to LAMMPS</a></li>
			<li class="toctree-l1"><a class="reference internal" href="Section_errors.html">12. Errors</a></li>
			<li class="toctree-l1"><a class="reference internal" href="Section_history.html">13. Future and history</a></li>
			</ul>



			</div>

			</nav>

			<section data-toggle="wy-nav-shift" class="wy-nav-content-wrap">


			<nav class="wy-nav-top" role="navigation" aria-label="top navigation">
			<i data-toggle="wy-nav-top" class="fa fa-bars"></i>
			<a href="Manual.html">LAMMPS</a>
			</nav>



			<div class="wy-nav-content">
			<div class="rst-content">
			<div role="navigation" aria-label="breadcrumbs navigation">
			<ul class="wy-breadcrumbs">
			<li><a href="Manual.html">Docs</a> »</li>

			<li>6. How-to discussions</li>
			<li class="wy-breadcrumbs-aside">


			<a href="http://lammps.sandia.gov">Website</a>
			<a href="Section_commands.html#comm">Commands</a>

			</li>
			</ul>
			<hr/>

			<div class="rst-footer-buttons" style="margin-bottom: 1em" role="navigation" aria-label="footer navigation">

			<a href="Section_example.html" class="btn btn-neutral float-right" title="7. Example problems" accesskey="n">Next <span class="fa fa-arrow-circle-right"></span></a>


			<a href="Section_accelerate.html" class="btn btn-neutral" title="5. Accelerating LAMMPS performance" accesskey="p"><span class="fa fa-arrow-circle-left"></span> Previous</a>

			</div>

			</div>
			<div role="main" class="document" itemscope="itemscope" itemtype="http://schema.org/Article">
			<div itemprop="articleBody">

			<div class="section" id="how-to-discussions">
			<h1>6. How-to discussions</h1>
			<p>This section describes how to perform common tasks using LAMMPS.</p>
			<div class="line-block">
			<div class="line">6.1 <a class="reference internal" href="#howto-1"><span class="std std-ref">Restarting a simulation</span></a></div>
			<div class="line">6.2 <a class="reference internal" href="#howto-2"><span class="std std-ref">2d simulations</span></a></div>
			<div class="line">6.3 <a class="reference internal" href="#howto-3"><span class="std std-ref">CHARMM, AMBER, and DREIDING force fields</span></a></div>
			<div class="line">6.4 <a class="reference internal" href="#howto-4"><span class="std std-ref">Running multiple simulations from one input script</span></a></div>
			<div class="line">6.5 <a class="reference internal" href="#howto-5"><span class="std std-ref">Multi-replica simulations</span></a></div>
			<div class="line">6.6 <a class="reference internal" href="#howto-6"><span class="std std-ref">Granular models</span></a></div>
			<div class="line">6.7 <a class="reference internal" href="#howto-7"><span class="std std-ref">TIP3P water model</span></a></div>
			<div class="line">6.8 <a class="reference internal" href="#howto-8"><span class="std std-ref">TIP4P water model</span></a></div>
			<div class="line">6.9 <a class="reference internal" href="#howto-9"><span class="std std-ref">SPC water model</span></a></div>
			<div class="line">6.10 <a class="reference internal" href="#howto-10"><span class="std std-ref">Coupling LAMMPS to other codes</span></a></div>
			<div class="line">6.11 <a class="reference internal" href="#howto-11"><span class="std std-ref">Visualizing LAMMPS snapshots</span></a></div>
			<div class="line">6.12 <a class="reference internal" href="#howto-12"><span class="std std-ref">Triclinic (non-orthogonal) simulation boxes</span></a></div>
			<div class="line">6.13 <a class="reference internal" href="#howto-13"><span class="std std-ref">NEMD simulations</span></a></div>
			<div class="line">6.14 <a class="reference internal" href="#howto-14"><span class="std std-ref">Finite-size spherical and aspherical particles</span></a></div>
			<div class="line">6.15 <a class="reference internal" href="#howto-15"><span class="std std-ref">Output from LAMMPS (thermo, dumps, computes, fixes, variables)</span></a></div>
			<div class="line">6.16 <a class="reference internal" href="#howto-16"><span class="std std-ref">Thermostatting, barostatting and computing temperature</span></a></div>
			<div class="line">6.17 <a class="reference internal" href="#howto-17"><span class="std std-ref">Walls</span></a></div>
			<div class="line">6.18 <a class="reference internal" href="#howto-18"><span class="std std-ref">Elastic constants</span></a></div>
			<div class="line">6.19 <a class="reference internal" href="#howto-19"><span class="std std-ref">Library interface to LAMMPS</span></a></div>
			<div class="line">6.20 <a class="reference internal" href="#howto-20"><span class="std std-ref">Calculating thermal conductivity</span></a></div>
			<div class="line">6.21 <a class="reference internal" href="#howto-21"><span class="std std-ref">Calculating viscosity</span></a></div>
			<div class="line">6.22 <a class="reference internal" href="#howto-22"><span class="std std-ref">Calculating a diffusion coefficient</span></a></div>
			<div class="line">6.23 <a class="reference internal" href="#howto-23"><span class="std std-ref">Using chunks to calculate system properties</span></a></div>
			<div class="line">6.24 <a class="reference internal" href="#howto-24"><span class="std std-ref">Setting parameters for the kspace_style pppm/disp command</span></a></div>
			<div class="line">6.25 <a class="reference internal" href="#howto-25"><span class="std std-ref">Polarizable models</span></a></div>
			<div class="line">6.26 <a class="reference internal" href="#howto-26"><span class="std std-ref">Adiabatic core/shell model</span></a></div>
			<div class="line">6.27 <a class="reference internal" href="#howto-27"><span class="std std-ref">Drude induced dipoles</span></a></div>
			<div class="line"><br /></div>
			</div>
			<p>The example input scripts included in the LAMMPS distribution and
			highlighted in <a class="reference internal" href="Section_example.html"><span class="doc">Section_example</span></a> also show how to
			setup and run various kinds of simulations.</p>
			<div class="section" id="restarting-a-simulation">
			<span id="howto-1"></span><h2>6.1. Restarting a simulation</h2>
			<p>There are 3 ways to continue a long LAMMPS simulation. Multiple
			<a class="reference internal" href="run.html"><span class="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 <a class="reference internal" href="restart.html"><span class="doc">restart</span></a>
			command. At a later time, these binary files can be read via a
			<a class="reference internal" href="read_restart.html"><span class="doc">read_restart</span></a> command in a new script. Or they can
			be converted to text data files using the <a class="reference internal" href="Section_start.html#start-7"><span class="std std-ref">-r command-line switch</span></a> and read by a
			<a class="reference internal" href="read_data.html"><span class="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 <a class="reference internal" href="read_restart.html"><span class="doc">read_restart</span></a> and
			<a class="reference internal" href="read_data.html"><span class="doc">read_data</span></a> commands.</p>
			<p>Look at the <em>in.chain</em> input script provided in the <em>bench</em> directory
			of the LAMMPS distribution to see the original script that these 2
			scripts are based on. If that script had the line</p>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">restart</span> <span class="mi">50</span> <span class="n">tmp</span><span class="o">.</span><span class="n">restart</span>
			</pre></div>
			</div>
			<p>added to it, it would produce 2 binary restart files (tmp.restart.50
			and tmp.restart.100) as it ran.</p>
			<p>This script could be used to read the 1st restart file and re-run the
			last 50 timesteps:</p>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">read_restart</span> <span class="n">tmp</span><span class="o">.</span><span class="n">restart</span><span class="o">.</span><span class="mi">50</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">neighbor</span> <span class="mf">0.4</span> <span class="nb">bin</span>
			<span class="n">neigh_modify</span> <span class="n">every</span> <span class="mi">1</span> <span class="n">delay</span> <span class="mi">1</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">nve</span>
			<span class="n">fix</span> <span class="mi">2</span> <span class="nb">all</span> <span class="n">langevin</span> <span class="mf">1.0</span> <span class="mf">1.0</span> <span class="mf">10.0</span> <span class="mi">904297</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">timestep</span> <span class="mf">0.012</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">run</span> <span class="mi">50</span>
			</pre></div>
			</div>
			<p>Note that the following commands do not need to be repeated because
			their settings are included in the restart file: <em>units, atom_style,
			special_bonds, pair_style, bond_style</em>. However these commands do
			need to be used, since their settings are not in the restart file:
			<em>neighbor, fix, timestep</em>.</p>
			<p>If you actually use this script to perform a restarted run, you will
			notice that the thermodynamic data match at step 50 (if you also put a
			“thermo 50” command in the original script), but do not match at step
			100. This is because the <a class="reference internal" href="fix_langevin.html"><span class="doc">fix langevin</span></a> command
			uses random numbers in a way that does not allow for perfect restarts.</p>
			<p>As an alternate approach, the restart file could be converted to a data
			file as follows:</p>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">lmp_g</span><span class="o">++</span> <span class="o">-</span><span class="n">r</span> <span class="n">tmp</span><span class="o">.</span><span class="n">restart</span><span class="o">.</span><span class="mi">50</span> <span class="n">tmp</span><span class="o">.</span><span class="n">restart</span><span class="o">.</span><span class="n">data</span>
			</pre></div>
			</div>
			<p>Then, this script could be used to re-run the last 50 steps:</p>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">units</span> <span class="n">lj</span>
			<span class="n">atom_style</span> <span class="n">bond</span>
			<span class="n">pair_style</span> <span class="n">lj</span><span class="o">/</span><span class="n">cut</span> <span class="mf">1.12</span>
			<span class="n">pair_modify</span> <span class="n">shift</span> <span class="n">yes</span>
			<span class="n">bond_style</span> <span class="n">fene</span>
			<span class="n">special_bonds</span> <span class="mf">0.0</span> <span class="mf">1.0</span> <span class="mf">1.0</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">read_data</span> <span class="n">tmp</span><span class="o">.</span><span class="n">restart</span><span class="o">.</span><span class="n">data</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">neighbor</span> <span class="mf">0.4</span> <span class="nb">bin</span>
			<span class="n">neigh_modify</span> <span class="n">every</span> <span class="mi">1</span> <span class="n">delay</span> <span class="mi">1</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">nve</span>
			<span class="n">fix</span> <span class="mi">2</span> <span class="nb">all</span> <span class="n">langevin</span> <span class="mf">1.0</span> <span class="mf">1.0</span> <span class="mf">10.0</span> <span class="mi">904297</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">timestep</span> <span class="mf">0.012</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">reset_timestep</span> <span class="mi">50</span>
			<span class="n">run</span> <span class="mi">50</span>
			</pre></div>
			</div>
			<p>Note that nearly all the settings specified in the original <em>in.chain</em>
			script must be repeated, except the <em>pair_coeff</em> and <em>bond_coeff</em>
			commands since the new data file lists the force field coefficients.
			Also, the <a class="reference internal" href="reset_timestep.html"><span class="doc">reset_timestep</span></a> command is used to tell
			LAMMPS the current timestep. This value is stored in restart files,
			but not in data files.</p>
			<hr class="docutils" />
			</div>
			<div class="section" id="d-simulations">
			<span id="howto-2"></span><h2>6.2. 2d simulations</h2>
			<p>Use the <a class="reference internal" href="dimension.html"><span class="doc">dimension</span></a> command to specify a 2d simulation.</p>
			<p>Make the simulation box periodic in z via the <a class="reference internal" href="boundary.html"><span class="doc">boundary</span></a>
			command. This is the default.</p>
			<p>If using the <a class="reference internal" href="create_box.html"><span class="doc">create box</span></a> command to define a
			simulation box, set the z dimensions narrow, but finite, so that the
			create_atoms command will tile the 3d simulation box with a single z
			plane of atoms - e.g.</p>
			<pre class="literal-block">
			<a class="reference internal" href="create_box.html"><span class="doc">create box</span></a> 1 -10 10 -10 10 -0.25 0.25
			</pre>
			<p>If using the <a class="reference internal" href="read_data.html"><span class="doc">read data</span></a> command to read in a file of
			atom coordinates, set the “zlo zhi” values to be finite but narrow,
			similar to the create_box command settings just described. For each
			atom in the file, assign a z coordinate so it falls inside the
			z-boundaries of the box - e.g. 0.0.</p>
			<p>Use the <a class="reference internal" href="fix_enforce2d.html"><span class="doc">fix enforce2d</span></a> command as the last
			defined fix to insure that the z-components of velocities and forces
			are zeroed out every timestep. The reason to make it the last fix is
			so that any forces induced by other fixes will be zeroed out.</p>
			<p>Many of the example input scripts included in the LAMMPS distribution
			are for 2d models.</p>
			<div class="admonition note">
			<p class="first admonition-title">Note</p>
			<p class="last">Some models in LAMMPS treat particles as finite-size spheres, as
			opposed to point particles. In 2d, the particles will still be
			spheres, not disks, meaning their moment of inertia will be the same
			as in 3d.</p>
			</div>
			<hr class="docutils" />
			</div>
			<div class="section" id="charmm-amber-and-dreiding-force-fields">
			<span id="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
			<a class="reference internal" href="read_data.html"><span class="doc">read_data</span></a> command or in the input script with
			commands like <a class="reference internal" href="pair_coeff.html"><span class="doc">pair_coeff</span></a> or
			<a class="reference internal" href="bond_coeff.html"><span class="doc">bond_coeff</span></a>. See <a class="reference internal" href="Section_tools.html"><span class="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 <a class="reference internal" href="pair_charmm.html#mackerell"><span class="std std-ref">(MacKerell)</span></a> for a description of the CHARMM force
			field. See <a class="reference internal" href="dihedral_charmm.html#cornell"><span class="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
			documentation for the formula it computes.</p>
			<ul class="simple">
			<li><a class="reference internal" href="bond_harmonic.html"><span class="doc">bond_style</span></a> harmonic</li>
			<li><a class="reference internal" href="angle_charmm.html"><span class="doc">angle_style</span></a> charmm</li>
			<li><a class="reference internal" href="dihedral_charmm.html"><span class="doc">dihedral_style</span></a> charmm</li>
			<li><a class="reference internal" href="pair_charmm.html"><span class="doc">pair_style</span></a> lj/charmm/coul/charmm</li>
			<li><a class="reference internal" href="pair_charmm.html"><span class="doc">pair_style</span></a> lj/charmm/coul/charmm/implicit</li>
			<li><a class="reference internal" href="pair_charmm.html"><span class="doc">pair_style</span></a> lj/charmm/coul/long</li>
			<li><a class="reference internal" href="special_bonds.html"><span class="doc">special_bonds</span></a> charmm</li>
			<li><a class="reference internal" href="special_bonds.html"><span class="doc">special_bonds</span></a> amber</li>
			</ul>
			<p>DREIDING is a generic force field developed by the <a class="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
			<a class="reference internal" href="pair_hbond_dreiding.html"><span class="doc">explicit hydrogen bond term</span></a> to describe
			interactions involving a hydrogen atom on very electronegative atoms
			(N, O, F).</p>
			<p>See <a class="reference internal" href="pair_hbond_dreiding.html#mayo"><span class="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
			documentation for the formula it computes.</p>
			<ul class="simple">
			<li><a class="reference internal" href="bond_harmonic.html"><span class="doc">bond_style</span></a> harmonic</li>
			<li><a class="reference internal" href="bond_morse.html"><span class="doc">bond_style</span></a> morse</li>
			<li><a class="reference internal" href="angle_harmonic.html"><span class="doc">angle_style</span></a> harmonic</li>
			<li><a class="reference internal" href="angle_cosine.html"><span class="doc">angle_style</span></a> cosine</li>
			<li><a class="reference internal" href="angle_cosine_periodic.html"><span class="doc">angle_style</span></a> cosine/periodic</li>
			<li><a class="reference internal" href="dihedral_charmm.html"><span class="doc">dihedral_style</span></a> charmm</li>
			<li><a class="reference internal" href="improper_umbrella.html"><span class="doc">improper_style</span></a> umbrella</li>
			<li><a class="reference internal" href="pair_buck.html"><span class="doc">pair_style</span></a> buck</li>
			<li><a class="reference internal" href="pair_buck.html"><span class="doc">pair_style</span></a> buck/coul/cut</li>
			<li><a class="reference internal" href="pair_buck.html"><span class="doc">pair_style</span></a> buck/coul/long</li>
			<li><a class="reference internal" href="pair_lj.html"><span class="doc">pair_style</span></a> lj/cut</li>
			<li><a class="reference internal" href="pair_lj.html"><span class="doc">pair_style</span></a> lj/cut/coul/cut</li>
			<li><a class="reference internal" href="pair_lj.html"><span class="doc">pair_style</span></a> lj/cut/coul/long</li>
			<li><a class="reference internal" href="pair_hbond_dreiding.html"><span class="doc">pair_style</span></a> hbond/dreiding/lj</li>
			<li><a class="reference internal" href="pair_hbond_dreiding.html"><span class="doc">pair_style</span></a> hbond/dreiding/morse</li>
			<li><a class="reference internal" href="special_bonds.html"><span class="doc">special_bonds</span></a> dreiding</li>
			</ul>
			<hr class="docutils" />
			</div>
			<div class="section" id="running-multiple-simulations-from-one-input-script">
			<span id="howto-4"></span><h2>6.4. Running multiple simulations from one input script</h2>
			<p>This can be done in several ways. See the documentation for
			individual commands for more details on how these examples work.</p>
			<p>If “multiple simulations” means continue a previous simulation for
			more timesteps, then you simply use the <a class="reference internal" href="run.html"><span class="doc">run</span></a> command
			multiple times. For example, this script</p>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">units</span> <span class="n">lj</span>
			<span class="n">atom_style</span> <span class="n">atomic</span>
			<span class="n">read_data</span> <span class="n">data</span><span class="o">.</span><span class="n">lj</span>
			<span class="n">run</span> <span class="mi">10000</span>
			<span class="n">run</span> <span class="mi">10000</span>
			<span class="n">run</span> <span class="mi">10000</span>
			<span class="n">run</span> <span class="mi">10000</span>
			<span class="n">run</span> <span class="mi">10000</span>
			</pre></div>
			</div>
			<p>would run 5 successive simulations of the same system for a total of
			50,000 timesteps.</p>
			<p>If you wish to run totally different simulations, one after the other,
			the <a class="reference internal" href="clear.html"><span class="doc">clear</span></a> command can be used in between them to
			re-initialize LAMMPS. For example, this script</p>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">units</span> <span class="n">lj</span>
			<span class="n">atom_style</span> <span class="n">atomic</span>
			<span class="n">read_data</span> <span class="n">data</span><span class="o">.</span><span class="n">lj</span>
			<span class="n">run</span> <span class="mi">10000</span>
			<span class="n">clear</span>
			<span class="n">units</span> <span class="n">lj</span>
			<span class="n">atom_style</span> <span class="n">atomic</span>
			<span class="n">read_data</span> <span class="n">data</span><span class="o">.</span><span class="n">lj</span><span class="o">.</span><span class="n">new</span>
			<span class="n">run</span> <span class="mi">10000</span>
			</pre></div>
			</div>
			<p>would run 2 independent simulations, one after the other.</p>
			<p>For large numbers of independent simulations, you can use
			<a class="reference internal" href="variable.html"><span class="doc">variables</span></a> and the <a class="reference internal" href="next.html"><span class="doc">next</span></a> and
			<a class="reference internal" href="jump.html"><span class="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>
			<div class="highlight-default"><div class="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>
			<div class="highlight-default"><div class="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 <a class="reference internal" href="Section_start.html#start-7"><span class="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
			<a class="reference internal" href="variable.html"><span class="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>
			<hr class="docutils" />
			</div>
			<div class="section" id="multi-replica-simulations">
			<span id="howto-5"></span><h2>6.5. Multi-replica simulations</h2>
			<p>Several commands in LAMMPS run mutli-replica simulations, meaning
			that multiple instances (replicas) of your simulation are run
			simultaneously, with small amounts of data exchanged between replicas
			periodically.</p>
			<p>These are the relevant commands:</p>
			<ul class="simple">
			<li><a class="reference internal" href="neb.html"><span class="doc">neb</span></a> for nudged elastic band calculations</li>
			<li><a class="reference internal" href="prd.html"><span class="doc">prd</span></a> for parallel replica dynamics</li>
			<li><a class="reference internal" href="tad.html"><span class="doc">tad</span></a> for temperature accelerated dynamics</li>
			<li><a class="reference internal" href="temper.html"><span class="doc">temper</span></a> for parallel tempering</li>
			<li><a class="reference internal" href="fix_pimd.html"><span class="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 <a class="reference internal" href="Section_start.html#start-3"><span class="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 <a class="reference internal" href="Section_start.html#start-3"><span class="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 <a class="reference internal" href="Section_start.html#start-7"><span class="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.
			these commands:</p>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">mpirun</span> <span class="o">-</span><span class="n">np</span> <span class="mi">16</span> <span class="n">lmp_linux</span> <span class="o">-</span><span class="n">partition</span> <span class="mi">8</span><span class="n">x2</span> <span class="o">-</span><span class="ow">in</span> <span class="ow">in</span><span class="o">.</span><span class="n">temper</span>
			<span class="n">mpirun</span> <span class="o">-</span><span class="n">np</span> <span class="mi">8</span> <span class="n">lmp_linux</span> <span class="o">-</span><span class="n">partition</span> <span class="mi">8</span><span class="n">x1</span> <span class="o">-</span><span class="ow">in</span> <span class="ow">in</span><span class="o">.</span><span class="n">neb</span>
			</pre></div>
			</div>
			<p>would each run 8 replicas, on either 16 or 8 processors. Note the use
			of the <a class="reference internal" href="Section_start.html#start-7"><span class="std std-ref">-in command-line switch</span></a> to specify
			the input script which is required when running in multi-replica mode.</p>
			<p>Also note that with MPI installed on a machine (e.g. your desktop),
			you can run on more (virtual) processors than you have physical
			processors. Thus the above commands could be run on a
			single-processor (or few-processor) desktop so that you can run
			a multi-replica simulation on more replicas than you have
			physical processors.</p>
			<hr class="docutils" />
			</div>
			<div class="section" id="granular-models">
			<span id="howto-6"></span><h2>6.6. Granular models</h2>
			<p>Granular system are composed of spherical particles with a diameter,
			as opposed to point particles. This means they have an angular
			velocity and torque can be imparted to them to cause them to rotate.</p>
			<p>To run a simulation of a granular model, you will want to use
			the following commands:</p>
			<ul class="simple">
			<li><a class="reference internal" href="atom_style.html"><span class="doc">atom_style sphere</span></a></li>
			<li><a class="reference internal" href="fix_nve_sphere.html"><span class="doc">fix nve/sphere</span></a></li>
			<li><a class="reference internal" href="fix_gravity.html"><span class="doc">fix gravity</span></a></li>
			</ul>
			<p>This compute</p>
			<ul class="simple">
			<li><a class="reference internal" href="compute_erotate_sphere.html"><span class="doc">compute erotate/sphere</span></a></li>
			</ul>
			<p>calculates rotational kinetic energy which can be <a class="reference internal" href="#howto-15"><span class="std std-ref">output with thermodynamic info</span></a>.</p>
			<p>Use one of these 3 pair potentials, which compute forces and torques
			between interacting pairs of particles:</p>
			<ul class="simple">
			<li><a class="reference internal" href="pair_style.html"><span class="doc">pair_style</span></a> gran/history</li>
			<li><a class="reference internal" href="pair_style.html"><span class="doc">pair_style</span></a> gran/no_history</li>
			<li><a class="reference internal" href="pair_style.html"><span class="doc">pair_style</span></a> gran/hertzian</li>
			</ul>
			<p>These commands implement fix options specific to granular systems:</p>
			<ul class="simple">
			<li><a class="reference internal" href="fix_freeze.html"><span class="doc">fix freeze</span></a></li>
			<li><a class="reference internal" href="fix_pour.html"><span class="doc">fix pour</span></a></li>
			<li><a class="reference internal" href="fix_viscous.html"><span class="doc">fix viscous</span></a></li>
			<li><a class="reference internal" href="fix_wall_gran.html"><span class="doc">fix wall/gran</span></a></li>
			</ul>
			<p>The fix style <em>freeze</em> zeroes both the force and torque of frozen
			atoms, and should be used for granular system instead of the fix style
			<em>setforce</em>.</p>
			<p>For computational efficiency, you can eliminate needless pairwise
			computations between frozen atoms by using this command:</p>
			<ul class="simple">
			<li><a class="reference internal" href="neigh_modify.html"><span class="doc">neigh_modify</span></a> exclude</li>
			</ul>
			<hr class="docutils" />
			</div>
			<div class="section" id="tip3p-water-model">
			<span id="howto-7"></span><h2>6.7. TIP3P water model</h2>
			<p>The TIP3P water model as implemented in CHARMM
			<a class="reference internal" href="pair_charmm.html#mackerell"><span class="std std-ref">(MacKerell)</span></a> specifies a 3-site rigid water molecule with
			charges and Lennard-Jones parameters assigned to each of the 3 atoms.
			In LAMMPS the <a class="reference internal" href="fix_shake.html"><span class="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 TIP3P-CHARMM model with a
			cutoff. The K values can be used if a flexible TIP3P model (without
			fix shake) is desired. If the LJ epsilon and sigma for HH and OH are
			set to 0.0, it corresponds to the original 1983 TIP3P model
			<a class="reference internal" href="pair_lj.html#jorgensen"><span class="std std-ref">(Jorgensen)</span></a>.</p>
			<div class="line-block">
			<div class="line">O mass = 15.9994</div>
			<div class="line">H mass = 1.008</div>
			<div class="line">O charge = -0.834</div>
			<div class="line">H charge = 0.417</div>
			<div class="line">LJ epsilon of OO = 0.1521</div>
			<div class="line">LJ sigma of OO = 3.1507</div>
			<div class="line">LJ epsilon of HH = 0.0460</div>
			<div class="line">LJ sigma of HH = 0.4000</div>
			<div class="line">LJ epsilon of OH = 0.0836</div>
			<div class="line">LJ sigma of OH = 1.7753</div>
			<div class="line">K of OH bond = 450</div>
			<div class="line">r0 of OH bond = 0.9572</div>
			<div class="line">K of HOH angle = 55</div>
			<div class="line">theta of HOH angle = 104.52</div>
			<div class="line"><br /></div>
			</div>
			<p>These are the parameters to use for TIP3P with a long-range Coulombic
			solver (e.g. Ewald or PPPM in LAMMPS), see <a class="reference internal" href="pair_dipole.html#price"><span class="std std-ref">(Price)</span></a> for
			details:</p>
			<div class="line-block">
			<div class="line">O mass = 15.9994</div>
			<div class="line">H mass = 1.008</div>
			<div class="line">O charge = -0.830</div>
			<div class="line">H charge = 0.415</div>
			<div class="line">LJ epsilon of OO = 0.102</div>
			<div class="line">LJ sigma of OO = 3.188</div>
			<div class="line">LJ epsilon, sigma of OH, HH = 0.0</div>
			<div class="line">K of OH bond = 450</div>
			<div class="line">r0 of OH bond = 0.9572</div>
			<div class="line">K of HOH angle = 55</div>
			<div class="line">theta of HOH angle = 104.52</div>
			<div class="line"><br /></div>
			</div>
			<p>Wikipedia also has a nice article on <a class="reference external" href="http://en.wikipedia.org/wiki/Water_model">water models</a>.</p>
			<hr class="docutils" />
			</div>
			<div class="section" id="tip4p-water-model">
			<span id="howto-8"></span><h2>6.8. TIP4P water model</h2>
			<p>The four-point TIP4P rigid water model extends the traditional
			three-point TIP3P model by adding an additional site, usually
			massless, where the charge associated with the oxygen atom is placed.
			This site M is located at a fixed distance away from the oxygen along
			the bisector of the HOH bond angle. A bond style of <em>harmonic</em> and an
			angle style of <em>harmonic</em> or <em>charmm</em> should also be used.</p>
			<p>A TIP4P model is run with LAMMPS using either this command
			for a cutoff model:</p>
			<p><a class="reference internal" href="pair_lj.html"><span class="doc">pair_style lj/cut/tip4p/cut</span></a></p>
			<p>or these two commands for a long-range model:</p>
			<ul class="simple">
			<li><a class="reference internal" href="pair_lj.html"><span class="doc">pair_style lj/cut/tip4p/long</span></a></li>
			<li><a class="reference internal" href="kspace_style.html"><span class="doc">kspace_style pppm/tip4p</span></a></li>
			</ul>
			<p>For both models, the bond lengths and bond angles should be held fixed
			using the <a class="reference internal" href="fix_shake.html"><span class="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
			<a class="reference internal" href="pair_lj.html#jorgensen"><span class="std std-ref">(Jorgensen)</span></a>. Note that the OM distance is specified in
			the <a class="reference internal" href="pair_style.html"><span class="doc">pair_style</span></a> command, not as part of the pair
			coefficients.</p>
			<div class="line-block">
			<div class="line">O mass = 15.9994</div>
			<div class="line">H mass = 1.008</div>
			<div class="line">O charge = -1.040</div>
			<div class="line">H charge = 0.520</div>
			<div class="line">r0 of OH bond = 0.9572</div>
			<div class="line">theta of HOH angle = 104.52</div>
			<div class="line">OM distance = 0.15</div>
			<div class="line">LJ epsilon of O-O = 0.1550</div>
			<div class="line">LJ sigma of O-O = 3.1536</div>
			<div class="line">LJ epsilon, sigma of OH, HH = 0.0</div>
			<div class="line">Coulombic cutoff = 8.5</div>
			<div class="line"><br /></div>
			</div>
			<p>For the TIP4/Ice model (J Chem Phys, 122, 234511 (2005);
			<a class="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>
			<div class="line-block">
			<div class="line">O mass = 15.9994</div>
			<div class="line">H mass = 1.008</div>
			<div class="line">O charge = -1.1794</div>
			<div class="line">H charge = 0.5897</div>
			<div class="line">r0 of OH bond = 0.9572</div>
			<div class="line">theta of HOH angle = 104.52</div>
			<div class="line">OM distance = 0.1577</div>
			<div class="line">LJ epsilon of O-O = 0.21084</div>
			<div class="line">LJ sigma of O-O = 3.1668</div>
			<div class="line">LJ epsilon, sigma of OH, HH = 0.0</div>
			<div class="line">Coulombic cutoff = 8.5</div>
			<div class="line"><br /></div>
			</div>
			<p>For the TIP4P/2005 model (J Chem Phys, 123, 234505 (2005);
			<a class="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>
			<div class="line-block">
			<div class="line">O mass = 15.9994</div>
			<div class="line">H mass = 1.008</div>
			<div class="line">O charge = -1.1128</div>
			<div class="line">H charge = 0.5564</div>
			<div class="line">r0 of OH bond = 0.9572</div>
			<div class="line">theta of HOH angle = 104.52</div>
			<div class="line">OM distance = 0.1546</div>
			<div class="line">LJ epsilon of O-O = 0.1852</div>
			<div class="line">LJ sigma of O-O = 3.1589</div>
			<div class="line">LJ epsilon, sigma of OH, HH = 0.0</div>
			<div class="line">Coulombic cutoff = 8.5</div>
			<div class="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>
			<div class="line-block">
			<div class="line">O mass = 15.9994</div>
			<div class="line">H mass = 1.008</div>
			<div class="line">O charge = -1.0484</div>
			<div class="line">H charge = 0.5242</div>
			<div class="line">r0 of OH bond = 0.9572</div>
			<div class="line">theta of HOH angle = 104.52</div>
			<div class="line">OM distance = 0.1250</div>
			<div class="line">LJ epsilon of O-O = 0.16275</div>
			<div class="line">LJ sigma of O-O = 3.16435</div>
			<div class="line">LJ epsilon, sigma of OH, HH = 0.0</div>
			<div class="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 <a class="reference internal" href="pair_lj.html"><span class="doc">pair_style lj/cut/tip4p/long</span></a> command.</p>
			<p>Wikipedia also has a nice article on <a class="reference external" href="http://en.wikipedia.org/wiki/Water_model">water models</a>.</p>
			<hr class="docutils" />
			</div>
			<div class="section" id="spc-water-model">
			<span id="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 <a class="reference internal" href="fix_shake.html"><span class="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>
			<div class="line-block">
			<div class="line">O mass = 15.9994</div>
			<div class="line">H mass = 1.008</div>
			<div class="line">O charge = -0.820</div>
			<div class="line">H charge = 0.410</div>
			<div class="line">LJ epsilon of OO = 0.1553</div>
			<div class="line">LJ sigma of OO = 3.166</div>
			<div class="line">LJ epsilon, sigma of OH, HH = 0.0</div>
			<div class="line">r0 of OH bond = 1.0</div>
			<div class="line">theta of HOH angle = 109.47</div>
			<div class="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>
			<div class="line-block">
			<div class="line">O charge = -0.8476</div>
			<div class="line">H charge = 0.4238</div>
			<div class="line"><br /></div>
			</div>
			<p>See the <a class="reference internal" href="#berendsen"><span class="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 <a class="reference external" href="http://en.wikipedia.org/wiki/Water_model">water models</a>.</p>
			<hr class="docutils" />
			</div>
			<div class="section" id="coupling-lammps-to-other-codes">
			<span id="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 <a class="reference internal" href="fix.html"><span class="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
			<a class="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
			<a class="reference internal" href="fix_poems.html"><span class="doc">fix poems</span></a> command for more details. See <a class="reference internal" href="Section_modify.html"><span class="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
			<a class="reference internal" href="run.html"><span class="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 <a class="reference internal" href="Section_modify.html"><span class="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 <a class="reference internal" href="run.html"><span class="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>
			<ul class="simple">
			<li>simple: simple driver programs in C++ and C which invoke LAMMPS as a
			library</li>
			<li>lammps_quest: coupling of LAMMPS and <a class="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 <a class="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><a class="reference internal" href="Section_start.html#start-5"><span class="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
			<a class="reference internal" href="Section_python.html"><span class="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 <a class="reference internal" href="#howto-19"><span class="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>
			<hr class="docutils" />
			</div>
			<div class="section" id="visualizing-lammps-snapshots">
			<span id="howto-11"></span><h2>6.11. Visualizing LAMMPS snapshots</h2>
			<p>LAMMPS itself does not do visualization, but snapshots from LAMMPS
			simulations can be visualized (and analyzed) in a variety of ways.</p>
			<p>LAMMPS snapshots are created by the <a class="reference internal" href="dump.html"><span class="doc">dump</span></a> command which can
			create files in several formats. The native LAMMPS dump format is a
			text file (see “dump atom” or “dump custom”) which can be visualized
			by the <a class="reference internal" href="Section_tools.html#xmovie"><span class="std std-ref">xmovie</span></a> program, included with the
			LAMMPS package. This produces simple, fast 2d projections of 3d
			systems, and can be useful for rapid debugging of simulation geometry
			and atom trajectories.</p>
			<p>Several programs included with LAMMPS as auxiliary tools can convert
			native LAMMPS dump files to other formats. See the
			<a class="reference internal" href="Section_tools.html"><span class="doc">Section_tools</span></a> doc page for details. The first is
			the <a class="reference internal" href="Section_tools.html#charmm"><span class="std std-ref">ch2lmp tool</span></a>, which contains a
			lammps2pdb Perl script which converts LAMMPS dump files into PDB
			files. The second is the <a class="reference internal" href="Section_tools.html#arc"><span class="std std-ref">lmp2arc tool</span></a> which
			converts LAMMPS dump files into Accelrys’ Insight MD program files.
			The third is the <a class="reference internal" href="Section_tools.html#cfg"><span class="std std-ref">lmp2cfg tool</span></a> which converts
			LAMMPS dump files into CFG files which can be read into the
			<a class="reference external" href="http://mt.seas.upenn.edu/Archive/Graphics/A">AtomEye</a> visualizer.</p>
			<p>A Python-based toolkit distributed by our group can read native LAMMPS
			dump files, including custom dump files with additional columns of
			user-specified atom information, and convert them to various formats
			or pipe them into visualization software directly. See the <a class="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, <a class="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 <a class="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 <a class="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 <a class="reference internal" href="dump.html"><span class="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 <a class="reference external" href="http://www.ks.uiuc.edu/Research/vmd">VMD</a> for visualization.
			See the <a class="reference internal" href="dump.html"><span class="doc">dump</span></a> command for more information on XTC files.</p>
			<hr class="docutils" />
			</div>
			<div class="section" id="triclinic-non-orthogonal-simulation-boxes">
			<span id="howto-12"></span><h2>6.12. Triclinic (non-orthogonal) simulation boxes</h2>
			<p>By default, LAMMPS uses an orthogonal simulation box to encompass the
			particles. The <a class="reference internal" href="boundary.html"><span class="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 <a class="reference internal" href="create_box.html"><span class="doc">create_box</span></a> or
			<a class="reference internal" href="read_data.html"><span class="doc">read_data</span></a> or <a class="reference internal" href="read_restart.html"><span class="doc">read_restart</span></a>
			commands. Additionally, LAMMPS defines box size parameters lx,ly,lz
			where lx = xhi-xlo, and similarly in the y and z dimensions. The 6
			parameters, as well as lx,ly,lz, can be output via the <a class="reference internal" href="thermo_style.html"><span class="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>
			<img alt="_images/transform.jpg" class="align-center" src="_images/transform.jpg" />
			<p>where A = <a href="#id3"><span class="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>
			<img alt="_images/rotate.jpg" class="align-center" src="_images/rotate.jpg" />
			<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 <a class="reference internal" href="fix_deform.html"><span class="doc">fix deform</span></a> and <a class="reference internal" href="fix_nh.html"><span class="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 <a class="reference internal" href="create_box.html"><span class="doc">create_box</span></a> command is used with a region of
			style <em>prism</em>, then a triclinic box is setup. See the
			<a class="reference internal" href="region.html"><span class="doc">region</span></a> command for details. If the
			<a class="reference internal" href="read_data.html"><span class="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
			<a class="reference internal" href="read_data.html"><span class="doc">read_data</span></a> command for details. Finally, if the
			<a class="reference internal" href="read_restart.html"><span class="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 <a class="reference internal" href="fix_nh.html"><span class="doc">fix npt</span></a> or
			<a class="reference internal" href="fix_deform.html"><span class="doc">fix deform</span></a> commands. Alternatively, you can use the
			<a class="reference internal" href="change_box.html"><span class="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
			<a class="reference internal" href="thermo_style.html"><span class="doc">thermo_style custom</span></a> command.</p>
			<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
			<a class="reference internal" href="create_box.html"><span class="doc">create_box</span></a> or <a class="reference internal" href="read_data.html"><span class="doc">read_data</span></a> commands,
			as well as when the box shape changes dynamically during a simulation,
			e.g. via the <a class="reference internal" href="fix_deform.html"><span class="doc">fix deform</span></a> or <a class="reference internal" href="fix_nh.html"><span class="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 <a class="reference internal" href="fix_deform.html"><span class="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 <a class="reference internal" href="fix_deform.html"><span class="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 <a class="reference internal" href="box.html"><span class="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 <a class="reference internal" href="fix_deform.html"><span class="doc">fix deform</span></a> or
			<a class="reference internal" href="fix_nh.html"><span class="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>
			<img alt="_images/box.jpg" class="align-center" src="_images/box.jpg" />
			<p>The inverse relationship can be written as follows:</p>
			<img alt="_images/box_inverse.jpg" class="align-center" src="_images/box_inverse.jpg" />
			<p>The values of <em>a</em>, <em>b</em>, <em>c</em> , <em>alpha</em>, <em>beta</em> , and <em>gamma</em> can be printed
			out or accessed by computes using the
			<a class="reference internal" href="thermo_style.html"><span class="doc">thermo_style custom</span></a> keywords
			<em>cella</em>, <em>cellb</em>, <em>cellc</em>, <em>cellalpha</em>, <em>cellbeta</em>, <em>cellgamma</em>,
			respectively.</p>
			<p>As discussed on the <a class="reference internal" href="dump.html"><span class="doc">dump</span></a> command doc page, when the BOX
			BOUNDS for a snapshot is written to a dump file for a triclinic box,
			an orthogonal bounding box which encloses the triclinic simulation box
			is output, along with the 3 tilt factors (xy, xz, yz) of the triclinic
			box, formatted as follows:</p>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">ITEM</span><span class="p">:</span> <span class="n">BOX</span> <span class="n">BOUNDS</span> <span class="n">xy</span> <span class="n">xz</span> <span class="n">yz</span>
			<span class="n">xlo_bound</span> <span class="n">xhi_bound</span> <span class="n">xy</span>
			<span class="n">ylo_bound</span> <span class="n">yhi_bound</span> <span class="n">xz</span>
			<span class="n">zlo_bound</span> <span class="n">zhi_bound</span> <span class="n">yz</span>
			</pre></div>
			</div>
			<p>This bounding box is convenient for many visualization programs and is
			calculated from the 9 triclinic box parameters
			(xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz) as follows:</p>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">xlo_bound</span> <span class="o">=</span> <span class="n">xlo</span> <span class="o">+</span> <span class="n">MIN</span><span class="p">(</span><span class="mf">0.0</span><span class="p">,</span><span class="n">xy</span><span class="p">,</span><span class="n">xz</span><span class="p">,</span><span class="n">xy</span><span class="o">+</span><span class="n">xz</span><span class="p">)</span>
			<span class="n">xhi_bound</span> <span class="o">=</span> <span class="n">xhi</span> <span class="o">+</span> <span class="n">MAX</span><span class="p">(</span><span class="mf">0.0</span><span class="p">,</span><span class="n">xy</span><span class="p">,</span><span class="n">xz</span><span class="p">,</span><span class="n">xy</span><span class="o">+</span><span class="n">xz</span><span class="p">)</span>
			<span class="n">ylo_bound</span> <span class="o">=</span> <span class="n">ylo</span> <span class="o">+</span> <span class="n">MIN</span><span class="p">(</span><span class="mf">0.0</span><span class="p">,</span><span class="n">yz</span><span class="p">)</span>
			<span class="n">yhi_bound</span> <span class="o">=</span> <span class="n">yhi</span> <span class="o">+</span> <span class="n">MAX</span><span class="p">(</span><span class="mf">0.0</span><span class="p">,</span><span class="n">yz</span><span class="p">)</span>
			<span class="n">zlo_bound</span> <span class="o">=</span> <span class="n">zlo</span>
			<span class="n">zhi_bound</span> <span class="o">=</span> <span class="n">zhi</span>
			</pre></div>
			</div>
			<p>These formulas can be inverted if you need to convert the bounding box
			back into the triclinic box parameters, e.g. xlo = xlo_bound -
			MIN(0.0,xy,xz,xy+xz).</p>
			<p>One use of triclinic simulation boxes is to model solid-state crystals
			with triclinic symmetry. The <a class="reference internal" href="lattice.html"><span class="doc">lattice</span></a> command can be
			used with non-orthogonal basis vectors to define a lattice that will
			tile a triclinic simulation box via the
			<a class="reference internal" href="create_atoms.html"><span class="doc">create_atoms</span></a> command.</p>
			<p>A second use is to run Parinello-Rahman dyanamics via the <a class="reference internal" href="fix_nh.html"><span class="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 <a class="reference internal" href="minimize.html"><span class="doc">energy minimization</span></a> is
			the <a class="reference internal" href="fix_box_relax.html"><span class="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 <a class="reference internal" href="fix_deform.html"><span class="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
			on non-equilibrium MD (NEMD) simulations.</p>
			<hr class="docutils" />
			</div>
			<div class="section" id="nemd-simulations">
			<span id="howto-13"></span><h2>6.13. NEMD simulations</h2>
			<p>Non-equilibrium molecular dynamics or NEMD simulations are typically
			used to measure a fluid’s rheological properties such as viscosity.
			In LAMMPS, such simulations can be performed by first setting up a
			non-orthogonal simulation box (see the preceding Howto section).</p>
			<p>A shear strain can be applied to the simulation box at a desired
			strain rate by using the <a class="reference internal" href="fix_deform.html"><span class="doc">fix deform</span></a> command. The
			<a class="reference internal" href="fix_nvt_sllod.html"><span class="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 <a class="reference internal" href="compute_temp_deform.html"><span class="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 <a class="reference internal" href="fix_ave_spatial.html"><span class="doc">fix ave/spatial</span></a> command.</p>
			<p>As discussed in the previous section on non-orthogonal simulation
			boxes, the amount of tilt or skew that can be applied is limited by
			LAMMPS for computational efficiency to be 1/2 of the parallel box
			length. However, <a class="reference internal" href="fix_deform.html"><span class="doc">fix deform</span></a> can continuously strain
			a box by an arbitrary amount. As discussed in the <a class="reference internal" href="fix_deform.html"><span class="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 <a class="reference internal" href="fix_deform.html"><span class="doc">fix deform</span></a> should be set to “remap v”, since that is what
			<a class="reference internal" href="fix_nvt_sllod.html"><span class="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
			<a class="reference internal" href="fix_viscosity.html"><span class="doc">fix viscosity</span></a> command.</p>
			<hr class="docutils" />
			</div>
			<div class="section" id="finite-size-spherical-and-aspherical-particles">
			<span id="howto-14"></span><h2>6.14. Finite-size spherical and aspherical particles</h2>
			<p>Typical MD models treat atoms or particles as point masses. Sometimes
			it is desirable to have a model with finite-size particles such as
			spheroids or ellipsoids or generalized aspherical bodies. The
			difference is that such particles have a moment of inertia, rotational
			energy, and angular momentum. Rotation is induced by torque coming
			from interactions with other particles.</p>
			<p>LAMMPS has several options for running simulations with these kinds of
			particles. The following aspects are discussed in turn:</p>
			<ul class="simple">
			<li>atom styles</li>
			<li>pair potentials</li>
			<li>time integration</li>
			<li>computes, thermodynamics, and dump output</li>
			<li>rigid bodies composed of finite-size particles</li>
			</ul>
			<p>Example input scripts for these kinds of models are in the body,
			colloid, dipole, ellipse, line, peri, pour, and tri directories of the
			<a class="reference internal" href="Section_example.html"><span class="doc">examples directory</span></a> in the LAMMPS distribution.</p>
			<div class="section" id="atom-styles">
			<h3>6.14.1. Atom styles</h3>
			<p>There are several <a class="reference internal" href="atom_style.html"><span class="doc">atom styles</span></a> that allow for
			definition of finite-size particles: sphere, dipole, ellipsoid, line,
			tri, peri, and body.</p>
			<p>The sphere style defines particles that are spheriods and each
			particle can have a unique diameter and mass (or density). These
			particles store an angular velocity (omega) and can be acted upon by
			torque. The “set” command can be used to modify the diameter and mass
			of individual particles, after then are created.</p>
			<p>The dipole style does not actually define finite-size particles, but
			is often used in conjunction with spherical particles, via a command
			like</p>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">atom_style</span> <span class="n">hybrid</span> <span class="n">sphere</span> <span class="n">dipole</span>
			</pre></div>
			</div>
			<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 <a class="reference internal" href="atom_style.html"><span class="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 <a class="reference internal" href="pair_peri.html"><span class="doc">Peridynamic models</span></a> and
			defines particles as having a volume, that is used internally in the
			<a class="reference internal" href="pair_peri.html"><span class="doc">pair_style peri</span></a> potentials.</p>
			<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 <a class="reference internal" href="body.html"><span class="doc">body</span></a> doc page.</p>
			<p>Note that if one of these atom styles is used (or multiple styles via
			the <a class="reference internal" href="atom_style.html"><span class="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,
			<a class="reference internal" href="pair_hybrid.html"><span class="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 <a class="reference internal" href="dimension.html"><span class="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>
			<div class="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 <a class="reference internal" href="pair_style.html"><span class="doc">pair styles</span></a> that generate torque:</p>
			<ul class="simple">
			<li><a class="reference internal" href="pair_gran.html"><span class="doc">pair_style gran/history</span></a></li>
			<li><a class="reference internal" href="pair_gran.html"><span class="doc">pair_style gran/hertzian</span></a></li>
			<li><a class="reference internal" href="pair_gran.html"><span class="doc">pair_style gran/no_history</span></a></li>
			<li><a class="reference internal" href="pair_dipole.html"><span class="doc">pair_style dipole/cut</span></a></li>
			<li><a class="reference internal" href="pair_gayberne.html"><span class="doc">pair_style gayberne</span></a></li>
			<li><a class="reference internal" href="pair_resquared.html"><span class="doc">pair_style resquared</span></a></li>
			<li><a class="reference internal" href="pair_brownian.html"><span class="doc">pair_style brownian</span></a></li>
			<li><a class="reference internal" href="pair_lubricate.html"><span class="doc">pair_style lubricate</span></a></li>
			<li><a class="reference internal" href="pair_line_lj.html"><span class="doc">pair_style line/lj</span></a></li>
			<li><a class="reference internal" href="pair_tri_lj.html"><span class="doc">pair_style tri/lj</span></a></li>
			<li><a class="reference internal" href="pair_body.html"><span class="doc">pair_style body</span></a></li>
			</ul>
			<p>The granular pair styles are used with spherical particles. The
			dipole pair style is used with the dipole atom style, which could be
			applied to spherical or ellipsoidal particles. The GayBerne and
			REsquared potentials require ellipsoidal particles, though they will
			also work if the 3 shape parameters are the same (a sphere). The
			Brownian and lubrication potentials are used with spherical particles.
			The line, tri, and body potentials are used with line segment,
			triangular, and body particles respectively.</p>
			</div>
			<div class="section" id="time-integration">
			<h3>6.14.3. Time integration</h3>
			<p>There are several fixes that perform time integration on finite-size
			spherical particles, meaning the integrators update the rotational
			orientation and angular velocity or angular momentum of the particles:</p>
			<ul class="simple">
			<li><a class="reference internal" href="fix_nve_sphere.html"><span class="doc">fix nve/sphere</span></a></li>
			<li><a class="reference internal" href="fix_nvt_sphere.html"><span class="doc">fix nvt/sphere</span></a></li>
			<li><a class="reference internal" href="fix_npt_sphere.html"><span class="doc">fix npt/sphere</span></a></li>
			</ul>
			<p>Likewise, there are 3 fixes that perform time integration on
			ellipsoidal particles:</p>
			<ul class="simple">
			<li><a class="reference internal" href="fix_nve_asphere.html"><span class="doc">fix nve/asphere</span></a></li>
			<li><a class="reference internal" href="fix_nvt_asphere.html"><span class="doc">fix nvt/asphere</span></a></li>
			<li><a class="reference internal" href="fix_npt_asphere.html"><span class="doc">fix npt/asphere</span></a></li>
			</ul>
			<p>The advantage of these fixes is that those which thermostat the
			particles include the rotational degrees of freedom in the temperature
			calculation and thermostatting. The <a class="reference external" href="fix_langevin">fix langevin</a>
			command can also be used with its <em>omgea</em> or <em>angmom</em> options to
			thermostat the rotational degrees of freedom for spherical or
			ellipsoidal particles. Other thermostatting fixes only operate on the
			translational kinetic energy of finite-size particles.</p>
			<p>These fixes perform constant NVE time integration on line segment,
			triangular, and body particles:</p>
			<ul class="simple">
			<li><a class="reference internal" href="fix_nve_line.html"><span class="doc">fix nve/line</span></a></li>
			<li><a class="reference internal" href="fix_nve_tri.html"><span class="doc">fix nve/tri</span></a></li>
			<li><a class="reference internal" href="fix_nve_body.html"><span class="doc">fix nve/body</span></a></li>
			</ul>
			<p>Note that for mixtures of point and finite-size particles, these
			integration fixes can only be used with <a class="reference internal" href="group.html"><span class="doc">groups</span></a> which
			contain finite-size particles.</p>
			</div>
			<div class="section" id="computes-thermodynamics-and-dump-output">
			<h3>6.14.4. Computes, thermodynamics, and dump output</h3>
			<p>There are several computes that calculate the temperature or
			rotational energy of spherical or ellipsoidal particles:</p>
			<ul class="simple">
			<li><a class="reference internal" href="compute_temp_sphere.html"><span class="doc">compute temp/sphere</span></a></li>
			<li><a class="reference internal" href="compute_temp_asphere.html"><span class="doc">compute temp/asphere</span></a></li>
			<li><a class="reference internal" href="compute_erotate_sphere.html"><span class="doc">compute erotate/sphere</span></a></li>
			<li><a class="reference internal" href="compute_erotate_asphere.html"><span class="doc">compute erotate/asphere</span></a></li>
			</ul>
			<p>These include rotational degrees of freedom in their computation. If
			you wish the thermodynamic output of temperature or pressure to use
			one of these computes (e.g. for a system entirely composed of
			finite-size particles), then the compute can be defined and the
			<a class="reference internal" href="thermo_modify.html"><span class="doc">thermo_modify</span></a> command used. Note that by default
			thermodynamic quantities will be calculated with a temperature that
			only includes translational degrees of freedom. See the
			<a class="reference internal" href="thermo_style.html"><span class="doc">thermo_style</span></a> command for details.</p>
			<p>These commands can be used to output various attributes of finite-size
			particles:</p>
			<ul class="simple">
			<li><a class="reference internal" href="dump.html"><span class="doc">dump custom</span></a></li>
			<li><a class="reference internal" href="compute_property_atom.html"><span class="doc">compute property/atom</span></a></li>
			<li><a class="reference internal" href="dump.html"><span class="doc">dump local</span></a></li>
			<li><a class="reference internal" href="compute_body_local.html"><span class="doc">compute body/local</span></a></li>
			</ul>
			<p>Attributes include the dipole moment, the angular velocity, the
			angular momentum, the quaternion, the torque, the end-point and
			corner-point coordinates (for line and tri particles), and
			sub-particle attributes of body particles.</p>
			</div>
			<div class="section" id="rigid-bodies-composed-of-finite-size-particles">
			<h3>6.14.5. Rigid bodies composed of finite-size particles</h3>
			<p>The <a class="reference internal" href="fix_rigid.html"><span class="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 <a class="reference internal" href="fix_shake.html"><span class="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 <a class="reference internal" href="fix_rigid.html"><span class="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 <a class="reference internal" href="fix_nve_body.html"><span class="doc">fix nve/body</span></a>. Interactions between pairs of body
			particles are computed via a command like <a class="reference internal" href="pair_body.html"><span class="doc">pair_style body</span></a>.</p>
			<hr class="docutils" />
			</div>
			</div>
			<div class="section" id="output-from-lammps-thermo-dumps-computes-fixes-variables">
			<span id="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>
			<ul class="simple">
			<li><a class="reference internal" href="thermo_style.html"><span class="doc">Thermodynamic output</span></a>, which is a list
			of quantities printed every few timesteps to the screen and logfile.</li>
			<li><a class="reference internal" href="dump.html"><span class="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: <a class="reference internal" href="fix_ave_time.html"><span class="doc">fix ave/time</span></a> for time averaging, <a class="reference internal" href="fix_ave_chunk.html"><span class="doc">fix ave/chunk</span></a> for spatial or other averaging, and <a class="reference internal" href="fix_print.html"><span class="doc">fix print</span></a> for single-line output of
			<a class="reference internal" href="variable.html"><span class="doc">variables</span></a>. Fix print can also output to the
			screen.</li>
			<li><a class="reference internal" href="restart.html"><span class="doc">Restart files</span></a>.</li>
			</ul>
			<p>A simulation prints one set of thermodynamic output and (optionally)
			restart files. It can generate any number of dump files and fix
			output files, depending on what <a class="reference internal" href="dump.html"><span class="doc">dump</span></a> and <a class="reference internal" href="fix.html"><span class="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 <a class="reference internal" href="Section_modify.html"><span class="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>
			<ul class="simple">
			<li><a class="reference internal" href="#global"><span class="std std-ref">Global/per-atom/local data</span></a></li>
			<li><a class="reference internal" href="#scalar"><span class="std std-ref">Scalar/vector/array data</span></a></li>
			<li><a class="reference internal" href="#thermo"><span class="std std-ref">Thermodynamic output</span></a></li>
			<li><a class="reference internal" href="#dump"><span class="std std-ref">Dump file output</span></a></li>
			<li><span class="xref std std-ref">Fixes that write output files</span></li>
			<li><a class="reference internal" href="#computeoutput"><span class="std std-ref">Computes that process output quantities</span></a></li>
			<li><span class="xref std std-ref">Fixes that process output quantities</span></li>
			<li><a class="reference internal" href="#compute"><span class="std std-ref">Computes that generate values to output</span></a></li>
			<li><a class="reference internal" href="#fix"><span class="std std-ref">Fixes that generate values to output</span></a></li>
			<li><a class="reference internal" href="#variable"><span class="std std-ref">Variables that generate values to output</span></a></li>
			<li><a class="reference internal" href="#table"><span class="std std-ref">Summary table of output options and data flow between commands</span></a></li>
			</ul>
			<div class="section" id="global-per-atom-local-data">
			<span id="global"></span><h3>6.15.1. Global/per-atom/local data</h3>
			<p>Various output-related commands work with three different styles of
			data: global, per-atom, or local. A global datum is one or more
			system-wide values, e.g. the temperature of the system. A per-atom
			datum is one or more values per atom, e.g. the kinetic energy of each
			atom. Local datums are calculated by each processor based on the
			atoms it owns, but there may be zero or more per atom, e.g. a list of
			bond distances.</p>
			</div>
			<div class="section" id="scalar-vector-array-data">
			<span id="scalar"></span><h3>6.15.2. Scalar/vector/array data</h3>
			<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
			“<a href="#id75"><span class="problematic" id="id76">c_</span></a>” would be replaced by “<a href="#id77"><span class="problematic" id="id78">f_</span></a>” for a fix, or “<a href="#id79"><span class="problematic" id="id80">v_</span></a>” for a variable:</p>
			<table border="1" class="docutils">
			<colgroup>
			<col width="21%" />
			<col width="79%" />
			</colgroup>
			<tbody valign="top">
			<tr class="row-odd"><td>c_ID</td>
			<td>entire scalar, vector, or array</td>
			</tr>
			<tr class="row-even"><td>c_ID[I]</td>
			<td>one element of vector, one column of array</td>
			</tr>
			<tr class="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
			or array.</p>
			</div>
			<div class="section" id="thermodynamic-output">
			<span id="thermo"></span><h3>6.15.3. Thermodynamic output</h3>
			<p>The frequency and format of thermodynamic output is set by the
			<a class="reference internal" href="thermo.html"><span class="doc">thermo</span></a>, <a class="reference internal" href="thermo_style.html"><span class="doc">thermo_style</span></a>, and
			<a class="reference internal" href="thermo_modify.html"><span class="doc">thermo_modify</span></a> commands. The
			<a class="reference internal" href="thermo_style.html"><span class="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 <a class="reference internal" href="compute.html"><span class="doc">compute</span></a>
			or <a class="reference internal" href="fix.html"><span class="doc">fix</span></a> or <a class="reference internal" href="variable.html"><span class="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 <a class="reference internal" href="dump.html"><span class="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 <a class="reference internal" href="thermo_modify.html"><span class="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. <a class="reference internal" href="variable.html"><span class="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
			intensive result.</p>
			</div>
			<div class="section" id="dump-file-output">
			<span id="dump"></span><h3>6.15.4. Dump file output</h3>
			<p>Dump file output is specified by the <a class="reference internal" href="dump.html"><span class="doc">dump</span></a> and
			<a class="reference internal" href="dump_modify.html"><span class="doc">dump_modify</span></a> commands. There are several
			pre-defined formats (dump atom, dump xtc, etc).</p>
			<p>There is also a <a class="reference internal" href="dump.html"><span class="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
			<a class="reference internal" href="compute.html"><span class="doc">compute</span></a> or <a class="reference internal" href="fix.html"><span class="doc">fix</span></a> or <a class="reference internal" href="variable.html"><span class="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 <a class="reference internal" href="dump.html"><span class="doc">dump custom</span></a> command.</p>
			<p>There is also a <a class="reference internal" href="dump.html"><span class="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
			<a class="reference internal" href="compute.html"><span class="doc">compute</span></a> or <a class="reference internal" href="fix.html"><span class="doc">fix</span></a> or <a class="reference internal" href="variable.html"><span class="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 <a class="reference internal" href="dump.html"><span class="doc">dump local</span></a>
			command.</p>
			</div>
			<div class="section" id="fixes-that-write-output-files">
			<span id="fixoutput"></span><h3>6.15.5. Fixes that write output files</h3>
			<p>Several fixes take various quantities as input and can write output
			files: <a class="reference internal" href="fix_ave_time.html"><span class="doc">fix ave/time</span></a>, <a class="reference internal" href="fix_ave_chunk.html"><span class="doc">fix ave/chunk</span></a>, <a class="reference internal" href="fix_ave_histo.html"><span class="doc">fix ave/histo</span></a>,
			<a class="reference internal" href="fix_ave_correlate.html"><span class="doc">fix ave/correlate</span></a>, and <a class="reference internal" href="fix_print.html"><span class="doc">fix print</span></a>.</p>
			<p>The <a class="reference internal" href="fix_ave_time.html"><span class="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
			<a class="reference internal" href="compute.html"><span class="doc">compute</span></a> values, global <a class="reference internal" href="fix.html"><span class="doc">fix</span></a> values, or
			<a class="reference internal" href="variable.html"><span class="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 <a class="reference internal" href="thermo_style.html"><span class="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 <a class="reference internal" href="fix_ave_chunk.html"><span class="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
			<a class="reference internal" href="compute.html"><span class="doc">compute</span></a>, by a <a class="reference internal" href="fix.html"><span class="doc">fix</span></a>, or by an atom-style
			<a class="reference internal" href="variable.html"><span class="doc">variable</span></a>. The chunk-averaged output of this fix can
			also be used as input to other output commands.</p>
			<p>The <a class="reference internal" href="fix_ave_histo.html"><span class="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 <a class="reference internal" href="fix_ave_correlate.html"><span class="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 <a class="reference internal" href="fix_print.html"><span class="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
			<a class="reference internal" href="variable.html"><span class="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 <a class="reference internal" href="thermo_style.html"><span class="doc">thermodynamic keywords</span></a>, <a class="reference internal" href="compute.html"><span class="doc">computes</span></a>,
			<a class="reference internal" href="fix.html"><span class="doc">fixes</span></a>, or other <a class="reference internal" href="variable.html"><span class="doc">variables</span></a>, or to per-atom
			values for a specific atom. Thus the <a class="reference internal" href="fix_print.html"><span class="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>
			</div>
			<div class="section" id="computes-that-process-output-quantities">
			<span id="computeoutput"></span><h3>6.15.6. Computes that process output quantities</h3>
			<p>The <a class="reference internal" href="compute_reduce.html"><span class="doc">compute reduce</span></a> and <a class="reference internal" href="compute_reduce.html"><span class="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 <a class="reference internal" href="compute_slice.html"><span class="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 <a class="reference internal" href="compute_property_atom.html"><span class="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 <a class="reference internal" href="dump.html"><span class="doc">dump custom</span></a> command.</p>
			<p>The <a class="reference internal" href="compute_property_local.html"><span class="doc">compute property/local</span></a> command takes
			a list of one or more pre-defined local attributes (bond info, angle
			info, etc) and stores the values in a local vector or array. These
			are produced as output values which can be used as input to other
			output commands.</p>
			</div>
			<div class="section" id="fixes-that-process-output-quantities">
			<span id="id5"></span><h3>6.15.7. Fixes that process output quantities</h3>
			<p>The <a class="reference internal" href="fix_vector.html"><span class="doc">fix vector</span></a> command can create global vectors as
			output from global scalars as input, accumulating them one element at
			a time.</p>
			<p>The <a class="reference internal" href="fix_ave_atom.html"><span class="doc">fix ave/atom</span></a> command performs time-averaging
			of per-atom vectors. The per-atom quantities can be atom attributes
			such as position, velocity, force. They can also be per-atom
			quantities calculated by a <a class="reference internal" href="compute.html"><span class="doc">compute</span></a>, by a
			<a class="reference internal" href="fix.html"><span class="doc">fix</span></a>, or by an atom-style <a class="reference internal" href="variable.html"><span class="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 <a class="reference internal" href="fix_store_state.html"><span class="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 <a class="reference internal" href="dump.html"><span class="doc">dump custom</span></a> command,
			including per-atom quantities calculated by a <a class="reference internal" href="compute.html"><span class="doc">compute</span></a>,
			by a <a class="reference internal" href="fix.html"><span class="doc">fix</span></a>, or by an atom-style <a class="reference internal" href="variable.html"><span class="doc">variable</span></a>.
			The output of this fix can be used as input to other output commands.</p>
			</div>
			<div class="section" id="computes-that-generate-values-to-output">
			<span id="compute"></span><h3>6.15.8. Computes that generate values to output</h3>
			<p>Every <a class="reference internal" href="compute.html"><span class="doc">compute</span></a> in LAMMPS produces either global or
			per-atom or local values. The values can be scalars or vectors or
			arrays of data. These values can be output using the other commands
			described in this section. The doc page for each compute command
			describes what it produces. Computes that produce per-atom or local
			values have the word “atom” or “local” in their style name. Computes
			without the word “atom” or “local” produce global values.</p>
			</div>
			<div class="section" id="fixes-that-generate-values-to-output">
			<span id="fix"></span><h3>6.15.9. Fixes that generate values to output</h3>
			<p>Some <a class="reference internal" href="fix.html"><span class="doc">fixes</span></a> in LAMMPS produces either global or per-atom or
			local values which can be accessed by other commands. The values can
			be scalars or vectors or arrays of data. These values can be output
			using the other commands described in this section. The doc page for
			each fix command tells whether it produces any output quantities and
			describes them.</p>
			</div>
			<div class="section" id="variables-that-generate-values-to-output">
			<span id="variable"></span><h3>6.15.10. Variables that generate values to output</h3>
			<p><a class="reference internal" href="variable.html"><span class="doc">Variables</span></a> defined in an input script can store one or
			more strings. But equal-style, vector-style, and atom-style or
			atomfile-style variables generate a global scalar value, global vector
			or values, or a per-atom vector, resepctively, when accessed. The
			formulas used to define these variables can contain references to the
			thermodynamic keywords and to global and per-atom data generated by
			computes, fixes, and other variables. The values generated by
			variables can be used as input to and thus output by the other
			commands described in this section.</p>
			</div>
			<div class="section" id="summary-table-of-output-options-and-data-flow-between-commands">
			<span id="table"></span><h3>6.15.11. Summary table of output options and data flow between commands</h3>
			<p>This table summarizes the various commands that can be used for
			generating output from LAMMPS. Each command produces output data of
			some kind and/or writes data to a file. Most of the commands can take
			data from other commands as input. Thus you can link many of these
			commands together in pipeline form, where data produced by one command
			is used as input to another command and eventually written to the
			screen or to a file. Note that to hook two commands together the
			output and input data types must match, e.g. global/per-atom/local
			data and scalar/vector/array data.</p>
			<p>Also note that, as described above, when a command takes a scalar as
			input, that could be an element of a vector or array. Likewise a
			vector input could be a column of an array.</p>
			<table border="1" class="docutils">
			<colgroup>
			<col width="39%" />
			<col width="32%" />
			<col width="30%" />
			</colgroup>
			<tbody valign="top">
			<tr class="row-odd"><td>Command</td>
			<td>Input</td>
			<td>Output</td>
			</tr>
			<tr class="row-even"><td><a class="reference internal" href="thermo_style.html"><span class="doc">thermo_style custom</span></a></td>
			<td>global scalars</td>
			<td>screen, log file</td>
			</tr>
			<tr class="row-odd"><td><a class="reference internal" href="dump.html"><span class="doc">dump custom</span></a></td>
			<td>per-atom vectors</td>
			<td>dump file</td>
			</tr>
			<tr class="row-even"><td><a class="reference internal" href="dump.html"><span class="doc">dump local</span></a></td>
			<td>local vectors</td>
			<td>dump file</td>
			</tr>
			<tr class="row-odd"><td><a class="reference internal" href="fix_print.html"><span class="doc">fix print</span></a></td>
			<td>global scalar from variable</td>
			<td>screen, file</td>
			</tr>
			<tr class="row-even"><td><a class="reference internal" href="print.html"><span class="doc">print</span></a></td>
			<td>global scalar from variable</td>
			<td>screen</td>
			</tr>
			<tr class="row-odd"><td><a class="reference internal" href="compute.html"><span class="doc">computes</span></a></td>
			<td>N/A</td>
			<td>global/per-atom/local scalar/vector/array</td>
			</tr>
			<tr class="row-even"><td><a class="reference internal" href="fix.html"><span class="doc">fixes</span></a></td>
			<td>N/A</td>
			<td>global/per-atom/local scalar/vector/array</td>
			</tr>
			<tr class="row-odd"><td><a class="reference internal" href="variable.html"><span class="doc">variables</span></a></td>
			<td>global scalars and vectors, per-atom vectors</td>
			<td>global scalar and vector, per-atom vector</td>
			</tr>
			<tr class="row-even"><td><a class="reference internal" href="compute_reduce.html"><span class="doc">compute reduce</span></a></td>
			<td>per-atom/local vectors</td>
			<td>global scalar/vector</td>
			</tr>
			<tr class="row-odd"><td><a class="reference internal" href="compute_slice.html"><span class="doc">compute slice</span></a></td>
			<td>global vectors/arrays</td>
			<td>global vector/array</td>
			</tr>
			<tr class="row-even"><td><a class="reference internal" href="compute_property_atom.html"><span class="doc">compute property/atom</span></a></td>
			<td>per-atom vectors</td>
			<td>per-atom vector/array</td>
			</tr>
			<tr class="row-odd"><td><a class="reference internal" href="compute_property_local.html"><span class="doc">compute property/local</span></a></td>
			<td>local vectors</td>
			<td>local vector/array</td>
			</tr>
			<tr class="row-even"><td><a class="reference internal" href="fix_vector.html"><span class="doc">fix vector</span></a></td>
			<td>global scalars</td>
			<td>global vector</td>
			</tr>
			<tr class="row-odd"><td><a class="reference internal" href="fix_ave_atom.html"><span class="doc">fix ave/atom</span></a></td>
			<td>per-atom vectors</td>
			<td>per-atom vector/array</td>
			</tr>
			<tr class="row-even"><td><a class="reference internal" href="fix_ave_time.html"><span class="doc">fix ave/time</span></a></td>
			<td>global scalars/vectors</td>
			<td>global scalar/vector/array, file</td>
			</tr>
			<tr class="row-odd"><td><a class="reference internal" href="fix_ave_chunk.html"><span class="doc">fix ave/chunk</span></a></td>
			<td>per-atom vectors</td>
			<td>global array, file</td>
			</tr>
			<tr class="row-even"><td><a class="reference internal" href="fix_ave_histo.html"><span class="doc">fix ave/histo</span></a></td>
			<td>global/per-atom/local scalars and vectors</td>
			<td>global array, file</td>
			</tr>
			<tr class="row-odd"><td><a class="reference internal" href="fix_ave_correlate.html"><span class="doc">fix ave/correlate</span></a></td>
			<td>global scalars</td>
			<td>global array, file</td>
			</tr>
			<tr class="row-even"><td><a class="reference internal" href="fix_store_state.html"><span class="doc">fix store/state</span></a></td>
			<td>per-atom vectors</td>
			<td>per-atom vector/array</td>
			</tr>
			</tbody>
			</table>
			<hr class="docutils" />
			</div>
			</div>
			<div class="section" id="thermostatting-barostatting-and-computing-temperature">
			<span id="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 <a class="reference internal" href="compute.html"><span class="doc">compute commands</span></a> calculate temperature, and the <a class="reference internal" href="compute_pressure.html"><span class="doc">compute pressure</span></a> command calculates pressure.</p>
			<ul class="simple">
			<li><a class="reference internal" href="compute_temp.html"><span class="doc">compute temp</span></a></li>
			<li><a class="reference internal" href="compute_temp_sphere.html"><span class="doc">compute temp/sphere</span></a></li>
			<li><a class="reference internal" href="compute_temp_asphere.html"><span class="doc">compute temp/asphere</span></a></li>
			<li><a class="reference internal" href="compute_temp_com.html"><span class="doc">compute temp/com</span></a></li>
			<li><a class="reference internal" href="compute_temp_deform.html"><span class="doc">compute temp/deform</span></a></li>
			<li><a class="reference internal" href="compute_temp_partial.html"><span class="doc">compute temp/partial</span></a></li>
			<li><a class="reference internal" href="compute_temp_profile.html"><span class="doc">compute temp/profile</span></a></li>
			<li><a class="reference internal" href="compute_temp_ramp.html"><span class="doc">compute temp/ramp</span></a></li>
			<li><a class="reference internal" href="compute_temp_region.html"><span class="doc">compute temp/region</span></a></li>
			</ul>
			<p>All but the first 3 calculate velocity biases directly (e.g. advection
			velocities) that are removed when computing the thermal temperature.
			<a class="reference internal" href="compute_temp_sphere.html"><span class="doc">Compute temp/sphere</span></a> and <a class="reference internal" href="compute_temp_asphere.html"><span class="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 <a class="reference internal" href="fix.html"><span class="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
			(temp/rescale). Dissipative particle dynamics (DPD) thermostatting
			can be invoked via the <em>dpd/tstat</em> pair style:</p>
			<ul class="simple">
			<li><a class="reference internal" href="fix_nh.html"><span class="doc">fix nvt</span></a></li>
			<li><a class="reference internal" href="fix_nvt_sphere.html"><span class="doc">fix nvt/sphere</span></a></li>
			<li><a class="reference internal" href="fix_nvt_asphere.html"><span class="doc">fix nvt/asphere</span></a></li>
			<li><a class="reference internal" href="fix_nvt_sllod.html"><span class="doc">fix nvt/sllod</span></a></li>
			<li><a class="reference internal" href="fix_temp_berendsen.html"><span class="doc">fix temp/berendsen</span></a></li>
			<li><a class="reference internal" href="fix_temp_csvr.html"><span class="doc">fix temp/csvr</span></a></li>
			<li><a class="reference internal" href="fix_langevin.html"><span class="doc">fix langevin</span></a></li>
			<li><a class="reference internal" href="fix_temp_rescale.html"><span class="doc">fix temp/rescale</span></a></li>
			<li><a class="reference internal" href="pair_dpd.html"><span class="doc">pair_style dpd/tstat</span></a></li>
			</ul>
			<p><a class="reference internal" href="fix_nh.html"><span class="doc">Fix nvt</span></a> only thermostats the translational velocity of
			particles. <a class="reference internal" href="fix_nvt_sllod.html"><span class="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 <a class="reference internal" href="#howto-13"><span class="std std-ref">NEMD simulations</span></a> section of this page for further details. <a class="reference internal" href="fix_nvt_sphere.html"><span class="doc">Fix nvt/sphere</span></a> and <a class="reference internal" href="fix_nvt_asphere.html"><span class="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 <a class="reference internal" href="fix_langevin.html"><span class="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
			<a class="reference internal" href="fix_modify.html"><span class="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 <a class="reference internal" href="compute_temp_partial.html"><span class="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 <a class="reference internal" href="compute_temp_profile.html"><span class="doc">compute temp/profile</span></a>.</p>
			<div class="admonition note">
			<p class="first admonition-title">Note</p>
			<p class="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>
			</div>
			<ul class="simple">
			<li><a class="reference internal" href="fix_nve.html"><span class="doc">fix nve</span></a></li>
			<li><a class="reference internal" href="fix_nve_sphere.html"><span class="doc">fix nve/sphere</span></a></li>
			<li><a class="reference internal" href="fix_nve_asphere.html"><span class="doc">fix nve/asphere</span></a></li>
			</ul>
			<p>Barostatting in LAMMPS is also performed by <a class="reference internal" href="fix.html"><span class="doc">fixes</span></a>. Two
			barosttating methods are currently available: Nose-Hoover (npt and
			nph) and Berendsen:</p>
			<ul class="simple">
			<li><a class="reference internal" href="fix_nh.html"><span class="doc">fix npt</span></a></li>
			<li><a class="reference internal" href="fix_npt_sphere.html"><span class="doc">fix npt/sphere</span></a></li>
			<li><a class="reference internal" href="fix_npt_asphere.html"><span class="doc">fix npt/asphere</span></a></li>
			<li><a class="reference internal" href="fix_nh.html"><span class="doc">fix nph</span></a></li>
			<li><a class="reference internal" href="fix_press_berendsen.html"><span class="doc">fix press/berendsen</span></a></li>
			</ul>
			<p>The <a class="reference internal" href="fix_nh.html"><span class="doc">fix npt</span></a> commands include a Nose-Hoover thermostat
			and barostat. <a class="reference internal" href="fix_nh.html"><span class="doc">Fix nph</span></a> is just a Nose/Hoover barostat;
			it does no thermostatting. Both <a class="reference internal" href="fix_nh.html"><span class="doc">fix nph</span></a> and <a class="reference internal" href="fix_press_berendsen.html"><span class="doc">fix press/bernendsen</span></a> can be used in conjunction
			with any of the thermostatting fixes.</p>
			<p>As with the thermostats, <a class="reference internal" href="fix_nh.html"><span class="doc">fix npt</span></a> and <a class="reference internal" href="fix_nh.html"><span class="doc">fix nph</span></a> only use translational motion of the particles in
			computing T and P and performing thermo/barostatting. <a class="reference internal" href="fix_npt_sphere.html"><span class="doc">Fix npt/sphere</span></a> and <a class="reference internal" href="fix_npt_asphere.html"><span class="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 <a class="reference internal" href="compute_pressure.html"><span class="doc">compute pressure</span></a> compute to calculate a current
			pressure. By default, this compute is created with a simple <a class="reference internal" href="compute_temp.html"><span class="doc">compute temp</span></a> (see the last argument of the <a class="reference internal" href="compute_pressure.html"><span class="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
			<a class="reference internal" href="fix_modify.html"><span class="doc">fix_modify</span></a> command for instructions on how to assign
			a temperature or pressure compute to a barostatting fix.</p>
			<div class="admonition note">
			<p class="first admonition-title">Note</p>
			<p class="last">As with the thermostats, the Nose/Hoover methods (<a class="reference internal" href="fix_nh.html"><span class="doc">fix npt</span></a> and <a class="reference internal" href="fix_nh.html"><span class="doc">fix nph</span></a>) perform time integration.
			<a class="reference internal" href="fix_press_berendsen.html"><span class="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
			<a class="reference internal" href="thermo_style.html"><span class="doc">thermo_style</span></a> command, often includes temperature
			and pressure values. As explained on the doc page for the
			<a class="reference internal" href="thermo_style.html"><span class="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 <a class="reference internal" href="thermo_style.html"><span class="doc">thermo_style custom</span></a> command. Or
			you can use the <a class="reference internal" href="thermo_modify.html"><span class="doc">thermo_modify</span></a> command to
			re-define what temperature or pressure compute is used for default
			thermodynamic output.</p>
			<hr class="docutils" />
			</div>
			<div class="section" id="walls">
			<span id="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
			<a class="reference internal" href="lattice.html"><span class="doc">lattice</span></a> and <a class="reference internal" href="create_atoms.html"><span class="doc">create_atoms</span></a> commands,
			or read in via the <a class="reference internal" href="read_data.html"><span class="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 <a class="reference internal" href="fix_nve.html"><span class="doc">fix nve</span></a> or <a class="reference internal" href="fix_nh.html"><span class="doc">fix nvt</span></a>
			is not used with the group that contains wall particles, their
			positions and velocities will not be updated.</p>
			<ul class="simple">
			<li><a class="reference internal" href="fix_aveforce.html"><span class="doc">fix aveforce</span></a> - set force on particles to average value, so they move together</li>
			<li><a class="reference internal" href="fix_setforce.html"><span class="doc">fix setforce</span></a> - set force on particles to a value, e.g. 0.0</li>
			<li><a class="reference internal" href="fix_freeze.html"><span class="doc">fix freeze</span></a> - freeze particles for use as granular walls</li>
			<li><a class="reference internal" href="fix_nve_noforce.html"><span class="doc">fix nve/noforce</span></a> - advect particles by their velocity, but without force</li>
			<li><a class="reference internal" href="fix_move.html"><span class="doc">fix move</span></a> - prescribe motion of particles by a linear velocity, oscillation, rotation, variable</li>
			</ul>
			<p>The <a class="reference internal" href="fix_move.html"><span class="doc">fix move</span></a> command offers the most generality, since
			the motion of individual particles can be specified with
			<a class="reference internal" href="variable.html"><span class="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 <a class="reference internal" href="neigh_modify.html"><span class="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 <a class="reference internal" href="bond_style.html"><span class="doc">bond</span></a>.
			The bonded particles do interact with other mobile particles.</p>
			<p>Idealized walls can be specified via several fix commands. <a class="reference internal" href="fix_wall_gran.html"><span class="doc">Fix wall/gran</span></a> creates frictional walls for use with
			granular particles; all the other commands create smooth walls.</p>
			<ul class="simple">
			<li><a class="reference internal" href="fix_wall_reflect.html"><span class="doc">fix wall/reflect</span></a> - reflective flat walls</li>
			<li><a class="reference internal" href="fix_wall.html"><span class="doc">fix wall/lj93</span></a> - flat walls, with Lennard-Jones 9/3 potential</li>
			<li><a class="reference internal" href="fix_wall.html"><span class="doc">fix wall/lj126</span></a> - flat walls, with Lennard-Jones 12/6 potential</li>
			<li><a class="reference internal" href="fix_wall.html"><span class="doc">fix wall/colloid</span></a> - flat walls, with <a class="reference internal" href="pair_colloid.html"><span class="doc">pair_style colloid</span></a> potential</li>
			<li><a class="reference internal" href="fix_wall.html"><span class="doc">fix wall/harmonic</span></a> - flat walls, with repulsive harmonic spring potential</li>
			<li><a class="reference internal" href="fix_wall_region.html"><span class="doc">fix wall/region</span></a> - use region surface as wall</li>
			<li><a class="reference internal" href="fix_wall_gran.html"><span class="doc">fix wall/gran</span></a> - flat or curved walls with <a class="reference internal" href="pair_gran.html"><span class="doc">pair_style granular</span></a> potential</li>
			</ul>
			<p>The <em>lj93</em>, <em>lj126</em>, <em>colloid</em>, and <em>harmonic</em> styles all allow the
			flat walls to move with a constant velocity, or oscillate in time.
			The <a class="reference internal" href="fix_wall_region.html"><span class="doc">fix wall/region</span></a> command offers the most
			generality, since the region surface is treated as a wall, and the
			geometry of the region can be a simple primitive volume (e.g. a
			sphere, or cube, or plane), or a complex volume made from the union
			and intersection of primitive volumes. <a class="reference internal" href="region.html"><span class="doc">Regions</span></a> can also
			specify a volume “interior” or “exterior” to the specified primitive
			shape or <em>union</em> or <em>intersection</em>. <a class="reference internal" href="region.html"><span class="doc">Regions</span></a> can also be
			“dynamic” meaning they move with constant velocity, oscillate, or
			rotate.</p>
			<p>The only frictional idealized walls currently in LAMMPS are flat or
			curved surfaces specified by the <a class="reference internal" href="fix_wall_gran.html"><span class="doc">fix wall/gran</span></a>
			command. At some point we plan to allow regoin surfaces to be used as
			frictional walls, as well as triangulated surfaces.</p>
			<hr class="docutils" />
			</div>
			<div class="section" id="elastic-constants">
			<span id="howto-18"></span><h2>6.18. Elastic constants</h2>
			<p>Elastic constants characterize the stiffness of a material. The formal
			definition is provided by the linear relation that holds between the
			stress and strain tensors in the limit of infinitesimal deformation.
			In tensor notation, this is expressed as s_ij = C_ijkl * e_kl, where
			the repeated indices imply summation. s_ij are the elements of the
			symmetric stress tensor. e_kl are the elements of the symmetric strain
			tensor. C_ijkl are the elements of the fourth rank tensor of elastic
			constants. In three dimensions, this tensor has 3^4=81 elements. Using
			Voigt notation, the tensor can be written as a 6x6 matrix, where C_ij
			is now the derivative of s_i w.r.t. e_j. Because s_i is itself a
			derivative w.r.t. e_i, it follows that C_ij is also symmetric, with at
			most 7*6/2 = 21 distinct elements.</p>
			<p>At zero temperature, it is easy to estimate these derivatives by
			deforming the simulation box in one of the six directions using the
			<a class="reference internal" href="change_box.html"><span class="doc">change_box</span></a> command and measuring the change in the
			stress tensor. A general-purpose script that does this is given in the
			examples/elastic directory described in <a class="reference internal" href="Section_example.html"><span class="doc">this section</span></a>.</p>
			<p>Calculating elastic constants at finite temperature is more
			challenging, because it is necessary to run a simulation that perfoms
			time averages of differential properties. One way to do this is to
			measure the change in average stress tensor in an NVT simulations when
			the cell volume undergoes a finite deformation. In order to balance
			the systematic and statistical errors in this method, the magnitude of
			the deformation must be chosen judiciously, and care must be taken to
			fully equilibrate the deformed cell before sampling the stress
			tensor. Another approach is to sample the triclinic cell fluctuations
			that occur in an NPT simulation. This method can also be slow to
			converge and requires careful post-processing <a class="reference internal" href="pair_sdk.html#shinoda"><span class="std std-ref">(Shinoda)</span></a></p>
			<hr class="docutils" />
			</div>
			<div class="section" id="library-interface-to-lammps">
			<span id="howto-19"></span><h2>6.19. Library interface to LAMMPS</h2>
			<p>As described in <a class="reference internal" href="Section_start.html#start-5"><span class="std std-ref">Section_start 5</span></a>, LAMMPS
			can be built as a library, so that it can be called by another code,
			used in a <a class="reference internal" href="#howto-10"><span class="std std-ref">coupled manner</span></a> with other
			codes, or driven through a <a class="reference internal" href="Section_python.html"><span class="doc">Python interface</span></a>.</p>
			<p>All of these methodologies use a C-style interface to LAMMPS that is
			provided in the files src/library.cpp and src/library.h. The
			functions therein have a C-style argument list, but contain C++ code
			you could write yourself in a C++ application that was invoking LAMMPS
			directly. The C++ code in the functions illustrates how to invoke
			internal LAMMPS operations. Note that LAMMPS classes are defined
			within a LAMMPS namespace (LAMMPS_NS) if you use them from another C++
			application.</p>
			<p>Library.cpp contains these 5 basic functions:</p>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">void</span> <span class="n">lammps_open</span><span class="p">(</span><span class="nb">int</span><span class="p">,</span> <span class="n">char</span> <span class="o"></span><span class="p">,</span> <span class="n">MPI_Comm</span><span class="p">,</span> <span class="n">void</span> <span class="o"></span><span class="p">)</span>
			<span class="n">void</span> <span class="n">lammps_close</span><span class="p">(</span><span class="n">void</span> <span class="o">*</span><span class="p">)</span>
			<span class="nb">int</span> <span class="n">lammps_version</span><span class="p">(</span><span class="n">void</span> <span class="o">*</span><span class="p">)</span>
			<span class="n">void</span> <span class="n">lammps_file</span><span class="p">(</span><span class="n">void</span> <span class="o"></span><span class="p">,</span> <span class="n">char</span> <span class="o"></span><span class="p">)</span>
			<span class="n">char</span> <span class="o"></span><span class="n">lammps_command</span><span class="p">(</span><span class="n">void</span> <span class="o"></span><span class="p">,</span> <span class="n">char</span> <span class="o">*</span><span class="p">)</span>
			</pre></div>
			</div>
			<p>The lammps_open() function is used to initialize LAMMPS, passing in a
			list of strings as if they were <a class="reference internal" href="Section_start.html#start-7"><span class="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>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">void</span> <span class="o"></span><span class="n">lammps_extract_global</span><span class="p">(</span><span class="n">void</span> <span class="o"></span><span class="p">,</span> <span class="n">char</span> <span class="o">*</span><span class="p">)</span>
			<span class="n">void</span> <span class="o"></span><span class="n">lammps_extract_atom</span><span class="p">(</span><span class="n">void</span> <span class="o"></span><span class="p">,</span> <span class="n">char</span> <span class="o">*</span><span class="p">)</span>
			<span class="n">void</span> <span class="o"></span><span class="n">lammps_extract_compute</span><span class="p">(</span><span class="n">void</span> <span class="o"></span><span class="p">,</span> <span class="n">char</span> <span class="o">*</span><span class="p">,</span> <span class="nb">int</span><span class="p">,</span> <span class="nb">int</span><span class="p">)</span>
			<span class="n">void</span> <span class="o"></span><span class="n">lammps_extract_fix</span><span class="p">(</span><span class="n">void</span> <span class="o"></span><span class="p">,</span> <span class="n">char</span> <span class="o">*</span><span class="p">,</span> <span class="nb">int</span><span class="p">,</span> <span class="nb">int</span><span class="p">,</span> <span class="nb">int</span><span class="p">,</span> <span class="nb">int</span><span class="p">)</span>
			<span class="n">void</span> <span class="o"></span><span class="n">lammps_extract_variable</span><span class="p">(</span><span class="n">void</span> <span class="o"></span><span class="p">,</span> <span class="n">char</span> <span class="o"></span><span class="p">,</span> <span class="n">char</span> <span class="o"></span><span class="p">)</span>
			<span class="nb">int</span> <span class="n">lammps_set_variable</span><span class="p">(</span><span class="n">void</span> <span class="o"></span><span class="p">,</span> <span class="n">char</span> <span class="o"></span><span class="p">,</span> <span class="n">char</span> <span class="o">*</span><span class="p">)</span>
			<span class="nb">int</span> <span class="n">lammps_get_natoms</span><span class="p">(</span><span class="n">void</span> <span class="o">*</span><span class="p">)</span>
			<span class="n">void</span> <span class="n">lammps_get_coords</span><span class="p">(</span><span class="n">void</span> <span class="o"></span><span class="p">,</span> <span class="n">double</span> <span class="o"></span><span class="p">)</span>
			<span class="n">void</span> <span class="n">lammps_put_coords</span><span class="p">(</span><span class="n">void</span> <span class="o"></span><span class="p">,</span> <span class="n">double</span> <span class="o"></span><span class="p">)</span>
			</pre></div>
			</div>
			<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 <a class="reference internal" href="Section_python.html"><span class="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
			LAMMPS.</p>
			<hr class="docutils" />
			</div>
			<div class="section" id="calculating-thermal-conductivity">
			<span id="howto-20"></span><h2>6.20. Calculating thermal conductivity</h2>
			<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 <a class="reference internal" href="#howto-21"><span class="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 <a class="reference internal" href="#howto-13"><span class="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 paper by <a class="reference internal" href="#ikeshoji"><span class="std std-ref">Ikeshoji and Hafskjold</span></a> for
			details of this idea. Note that thermostatting fixes such as <a class="reference internal" href="fix_nh.html"><span class="doc">fix nvt</span></a>, <a class="reference internal" href="fix_langevin.html"><span class="doc">fix langevin</span></a>, and <a class="reference internal" href="fix_temp_rescale.html"><span class="doc">fix temp/rescale</span></a> store the cumulative energy they
			add/subtract.</p>
			<p>Alternatively, as a second method, the <a class="reference internal" href="fix_heat.html"><span class="doc">fix heat</span></a>
			command can 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 <a class="reference internal" href="fix_ave_spatial.html"><span class="doc">fix ave/spatial</span></a> and <a class="reference internal" href="compute_ke_atom.html"><span class="doc">compute ke/atom</span></a> commands.</p>
			<p>The third method is to perform a reverse non-equilibrium MD simulation
			using the <a class="reference internal" href="fix_thermal_conductivity.html"><span class="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 <a class="reference internal" href="fix_ave_spatial.html"><span class="doc">fix ave/spatial</span></a> and <a class="reference internal" href="compute_ke_atom.html"><span class="doc">compute ke/atom</span></a> commands. The fix tallies the
			cumulative energy transfer that it performs. See the <a class="reference internal" href="fix_thermal_conductivity.html"><span class="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 <a class="reference internal" href="compute_heat_flux.html"><span class="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 <a class="reference internal" href="fix_ave_correlate.html"><span class="doc">fix ave/correlate</span></a> command to calculate the needed
			auto-correlation. See the doc page for the <a class="reference internal" href="compute_heat_flux.html"><span class="doc">compute heat/flux</span></a> command for an example input script
			that calculates the thermal conductivity of solid Ar via the GK
			formalism.</p>
			<hr class="docutils" />
			</div>
			<div class="section" id="calculating-viscosity">
			<span id="howto-21"></span><h2>6.21. Calculating viscosity</h2>
			<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 <a class="reference internal" href="#howto-20"><span class="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 <a class="reference internal" href="fix_deform.html"><span class="doc">fix deform</span></a>
			command, and using the <a class="reference internal" href="fix_nvt_sllod.html"><span class="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 <a class="reference internal" href="fix_ave_spatial.html"><span class="doc">fix ave/spatial</span></a> command, which determines
			grad(Vstream) in the equation above. E.g. the derivative in the
			y-direction of the Vx component of fluid motion or grad(Vstream) =
			dVx/dy. The Pxy off-diagonal component of the pressure or stress
			tensor, as calculated by the <a class="reference internal" href="compute_pressure.html"><span class="doc">compute pressure</span></a>
			command, can also be monitored, which is the J term in the equation
			above. See <a class="reference internal" href="#howto-13"><span class="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 <a class="reference internal" href="fix_viscosity.html"><span class="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 <a class="reference internal" href="fix_ave_spatial.html"><span class="doc">fix ave/spatial</span></a> command.
			The fix tallies the cummulative momentum transfer that it performs.
			See the <a class="reference internal" href="fix_viscosity.html"><span class="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>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="c1"># Sample LAMMPS input script for viscosity of liquid Ar</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="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>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="c1"># convert from LAMMPS real units to SI</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span>variable kB equal 1.3806504e-23 # [J/K/** Boltzmann
			variable atm2Pa equal 101325.0
			variable A2m equal 1.0e-10
			variable fs2s equal 1.0e-15
			variable convert equal ${atm2Pa}${atm2Pa}${fs2s}${A2m}${A2m}*${A2m}
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="c1"># setup problem</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span>dimension 3
			boundary p p p
			lattice fcc 5.376 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
			region box block 0 4 0 4 0 4
			create_box 1 box
			create_atoms 1 box
			mass 1 39.948
			pair_style lj/cut 13.0
			pair_coeff * * 0.2381 3.405
			timestep ${dt}
			thermo $d
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="c1"># equilibration and thermalization</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span>velocity all create $T 102486 mom yes rot yes dist gaussian
			fix NVT all nvt temp $T $T 10 drag 0.2
			run 8000
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="c1"># viscosity calculation, switch to NVE if desired</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="c1">#unfix NVT</span>
			<span class="c1">#fix NVE all nve</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span>reset_timestep 0
			variable pxy equal pxy
			variable pxz equal pxz
			variable pyz equal pyz
			fix SS all ave/correlate $s $p $d &
			v_pxy v_pxz v_pyz type auto file S0St.dat ave running
			variable scale equal ${convert}/(${kB}$T)$V$s${dt}
			variable v11 equal trap(f_SS[3])*${scale}
			variable v22 equal trap(f_SS[4])*${scale}
			variable v33 equal trap(f_SS[5])*${scale}
			thermo_style custom step temp press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33
			run 100000
			variable v equal (v_v11+v_v22+v_v33)/3.0
			variable ndens equal count(all)/vol
			print "average viscosity: $v [Pa.s/** @ $T K, ${ndens} /A^3"
			</pre></div>
			</div>
			<p>The fifth method is related to the above Green-Kubo method,
			but uses the Einstein formulation, analogous to the Einstein
			mean-square-displacement formulation for self-diffusivity. The
			time-integrated momentum fluxes play the role of Cartesian
			coordinates, whose mean-square displacement increases linearly
			with time at sufficiently long times.</p>
			<hr class="docutils" />
			</div>
			<div class="section" id="calculating-a-diffusion-coefficient">
			<span id="howto-22"></span><h2>6.22. Calculating a diffusion coefficient</h2>
			<p>The diffusion coefficient D of a material can be measured in at least
			2 ways using various options in LAMMPS. See the examples/DIFFUSE
			directory for scripts that implement the 2 methods discussed here for
			a simple Lennard-Jones fluid model.</p>
			<p>The first method is to measure the mean-squared displacement (MSD) of
			the system, via the <a class="reference internal" href="compute_msd.html"><span class="doc">compute msd</span></a> command. The slope
			of the MSD versus time is proportional to the diffusion coefficient.
			The instantaneous MSD values can be accumulated in a vector via the
			<a class="reference internal" href="fix_vector.html"><span class="doc">fix vector</span></a> command, and a line fit to the vector to
			compute its slope via the <a class="reference internal" href="variable.html"><span class="doc">variable slope</span></a> function, and
			thus extract D.</p>
			<p>The second method is to measure the velocity auto-correlation function
			(VACF) of the system, via the <a class="reference internal" href="compute_vacf.html"><span class="doc">compute vacf</span></a>
			command. The time-integral of the VACF is proportional to the
			diffusion coefficient. The instantaneous VACF values can be
			accumulated in a vector via the <a class="reference internal" href="fix_vector.html"><span class="doc">fix vector</span></a> command,
			and time integrated via the <a class="reference internal" href="variable.html"><span class="doc">variable trap</span></a> function,
			and thus extract D.</p>
			<hr class="docutils" />
			</div>
			<div class="section" id="using-chunks-to-calculate-system-properties">
			<span id="howto-23"></span><h2>6.23. Using chunks to calculate system properties</h2>
			<p>In LAMMS, “chunks” are collections of atoms, as defined by the
			<a class="reference internal" href="compute_chunk_atom.html"><span class="doc">compute chunk/atom</span></a> command, which assigns
			each atom to a chunk ID (or to no chunk at all). The number of chunks
			and the assignment of chunk IDs to atoms can be static or change over
			time. Examples of “chunks” are molecules or spatial bins or atoms
			with similar values (e.g. coordination number or potential energy).</p>
			<p>The per-atom chunk IDs can be used as input to two other kinds of
			commands, to calculate various properties of a system:</p>
			<ul class="simple">
			<li><a class="reference internal" href="fix_ave_chunk.html"><span class="doc">fix ave/chunk</span></a></li>
			<li>any of the <a class="reference internal" href="compute.html"><span class="doc">compute */chunk</span></a> commands</li>
			</ul>
			<p>Here, each of the 3 kinds of chunk-related commands is briefly
			overviewed. Then some examples are given of how to compute different
			properties with chunk commands.</p>
			<div class="section" id="compute-chunk-atom-command">
			<h3>6.23.1. Compute chunk/atom command:</h3>
			<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>
			<table border="1" class="docutils">
			<colgroup>
			<col width="31%" />
			<col width="69%" />
			</colgroup>
			<tbody valign="top">
			<tr class="row-odd"><td>atoms in same molecule</td>
			<td>chunk ID = molecule ID</td>
			</tr>
			<tr class="row-even"><td>atoms of same atom type</td>
			<td>chunk ID = atom type</td>
			</tr>
			<tr class="row-odd"><td>all atoms with same atom property (charge, radius, etc)</td>
			<td>chunk ID = output of compute property/atom</td>
			</tr>
			<tr class="row-even"><td>atoms in same cluster</td>
			<td>chunk ID = output of <a class="reference internal" href="compute_cluster_atom.html"><span class="doc">compute cluster/atom</span></a> command</td>
			</tr>
			<tr class="row-odd"><td>atoms in same spatial bin</td>
			<td>chunk ID = bin ID</td>
			</tr>
			<tr class="row-even"><td>atoms in same rigid body</td>
			<td>chunk ID = molecule ID used to define rigid bodies</td>
			</tr>
			<tr class="row-odd"><td>atoms with similar potential energy</td>
			<td>chunk ID = output of <a class="reference internal" href="compute_pe_atom.html"><span class="doc">compute pe/atom</span></a></td>
			</tr>
			<tr class="row-even"><td>atoms with same local defect structure</td>
			<td>chunk ID = output of <a class="reference internal" href="compute_centro_atom.html"><span class="doc">compute centro/atom</span></a> or <a class="reference internal" href="compute_coord_atom.html"><span class="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
			<a class="reference internal" href="compute.html"><span class="doc">computes</span></a>, <a class="reference internal" href="fix.html"><span class="doc">fixes</span></a>, and
			<a class="reference internal" href="variable.html"><span class="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
			<a class="reference internal" href="Section_modify.html"><span class="doc">Section_modify</span></a> of the documentation for how to
			do this. You can also define a <a class="reference internal" href="variable.html"><span class="doc">per-atom variable</span></a> in
			the input script that uses a formula to generate a chunk ID for each
			atom.</p>
			</div>
			<div class="section" id="fix-ave-chunk-command">
			<h3>6.23.2. Fix ave/chunk command:</h3>
			<p>This fix takes the ID of a <a class="reference internal" href="compute_chunk_atom.html"><span class="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
			velocity, force, charge, potential energy, kinetic energy, stress,
			etc. Additional keywords are defined for per-chunk properties like
			density and temperature. More generally any per-atom value generated
			by other <a class="reference internal" href="compute.html"><span class="doc">computes</span></a>, <a class="reference internal" href="fix.html"><span class="doc">fixes</span></a>, and <a class="reference internal" href="variable.html"><span class="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
			chunk.</p>
			</div>
			<div class="section" id="compute-chunk-commands">
			<h3>6.23.3. Compute <a href="#id72"><span class="problematic" id="id73">*</span></a>/chunk commands:</h3>
			<p>Currently the following computes operate on chunks of atoms to produce
			per-chunk values.</p>
			<ul class="simple">
			<li><a class="reference internal" href="compute_com_chunk.html"><span class="doc">compute com/chunk</span></a></li>
			<li><a class="reference internal" href="compute_gyration_chunk.html"><span class="doc">compute gyration/chunk</span></a></li>
			<li><a class="reference internal" href="compute_inertia_chunk.html"><span class="doc">compute inertia/chunk</span></a></li>
			<li><a class="reference internal" href="compute_msd_chunk.html"><span class="doc">compute msd/chunk</span></a></li>
			<li><a class="reference internal" href="compute_property_chunk.html"><span class="doc">compute property/chunk</span></a></li>
			<li><a class="reference internal" href="compute_temp_chunk.html"><span class="doc">compute temp/chunk</span></a></li>
			<li><a class="reference internal" href="compute_vcm_chunk.html"><span class="doc">compute torque/chunk</span></a></li>
			<li><a class="reference internal" href="compute_vcm_chunk.html"><span class="doc">compute vcm/chunk</span></a></li>
			</ul>
			<p>They each take the ID of a <a class="reference internal" href="compute_chunk_atom.html"><span class="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 <a class="reference internal" href="compute_property_chunk.html"><span class="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 <a class="reference internal" href="fix_ave_chunk.html"><span class="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>
			<ul class="simple">
			<li>As input to the <a class="reference internal" href="fix_ave_time.html"><span class="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 <a class="reference internal" href="fix_ave_histo.html"><span class="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 <a class="reference internal" href="variable.html"><span class="doc">equal-style variables</span></a>, like sum() and max(). E.g. to find the
			largest cluster or fastest diffusing molecule.</li>
			</ul>
			</div>
			<div class="section" id="example-calculations-with-chunks">
			<h3>6.23.4. Example calculations with chunks</h3>
			<p>Here are eaxmples using chunk commands to calculate various
			properties:</p>
			<ol class="arabic simple">
			<li>Average velocity in each of 1000 2d spatial bins:</li>
			</ol>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">compute</span> <span class="n">cc1</span> <span class="nb">all</span> <span class="n">chunk</span><span class="o">/</span><span class="n">atom</span> <span class="nb">bin</span><span class="o">/</span><span class="mi">2</span><span class="n">d</span> <span class="n">x</span> <span class="mf">0.0</span> <span class="mf">0.1</span> <span class="n">y</span> <span class="n">lower</span> <span class="mf">0.01</span> <span class="n">units</span> <span class="n">reduced</span>
			<span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">ave</span><span class="o">/</span><span class="n">chunk</span> <span class="mi">100</span> <span class="mi">10</span> <span class="mi">1000</span> <span class="n">cc1</span> <span class="n">vx</span> <span class="n">vy</span> <span class="n">file</span> <span class="n">tmp</span><span class="o">.</span><span class="n">out</span>
			</pre></div>
			</div>
			<p>(2) Temperature in each spatial bin, after subtracting a flow
			velocity:</p>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">compute</span> <span class="n">cc1</span> <span class="nb">all</span> <span class="n">chunk</span><span class="o">/</span><span class="n">atom</span> <span class="nb">bin</span><span class="o">/</span><span class="mi">2</span><span class="n">d</span> <span class="n">x</span> <span class="mf">0.0</span> <span class="mf">0.1</span> <span class="n">y</span> <span class="n">lower</span> <span class="mf">0.1</span> <span class="n">units</span> <span class="n">reduced</span>
			<span class="n">compute</span> <span class="n">vbias</span> <span class="nb">all</span> <span class="n">temp</span><span class="o">/</span><span class="n">profile</span> <span class="mi">1</span> <span class="mi">0</span> <span class="mi">0</span> <span class="n">y</span> <span class="mi">10</span>
			<span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">ave</span><span class="o">/</span><span class="n">chunk</span> <span class="mi">100</span> <span class="mi">10</span> <span class="mi">1000</span> <span class="n">cc1</span> <span class="n">temp</span> <span class="n">bias</span> <span class="n">vbias</span> <span class="n">file</span> <span class="n">tmp</span><span class="o">.</span><span class="n">out</span>
			</pre></div>
			</div>
			<ol class="arabic simple" start="3">
			<li>Center of mass of each molecule:</li>
			</ol>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">compute</span> <span class="n">cc1</span> <span class="nb">all</span> <span class="n">chunk</span><span class="o">/</span><span class="n">atom</span> <span class="n">molecule</span>
			<span class="n">compute</span> <span class="n">myChunk</span> <span class="nb">all</span> <span class="n">com</span><span class="o">/</span><span class="n">chunk</span> <span class="n">cc1</span>
			<span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">ave</span><span class="o">/</span><span class="n">time</span> <span class="mi">100</span> <span class="mi">1</span> <span class="mi">100</span> <span class="n">c_myChunk</span> <span class="n">file</span> <span class="n">tmp</span><span class="o">.</span><span class="n">out</span> <span class="n">mode</span> <span class="n">vector</span>
			</pre></div>
			</div>
			<ol class="arabic simple" start="4">
			<li>Total force on each molecule and ave/max across all molecules:</li>
			</ol>
			<pre class="literal-block">
			compute cc1 all chunk/atom molecule
			fix 1 all ave/chunk 1000 1 1000 cc1 fx fy fz file tmp.out
			variable xave equal ave(f_1<strong>2</strong>)
			variable xmax equal max(f_1<strong>2</strong>)
			thermo 1000
			thermo_style custom step temp v_xave v_xmax
			</pre>
			<ol class="arabic simple" start="5">
			<li>Histogram of cluster sizes:</li>
			</ol>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">compute</span> <span class="n">cluster</span> <span class="nb">all</span> <span class="n">cluster</span><span class="o">/</span><span class="n">atom</span> <span class="mf">1.0</span>
			<span class="n">compute</span> <span class="n">cc1</span> <span class="nb">all</span> <span class="n">chunk</span><span class="o">/</span><span class="n">atom</span> <span class="n">c_cluster</span> <span class="n">compress</span> <span class="n">yes</span>
			<span class="n">compute</span> <span class="n">size</span> <span class="nb">all</span> <span class="nb">property</span><span class="o">/</span><span class="n">chunk</span> <span class="n">cc1</span> <span class="n">count</span>
			<span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">ave</span><span class="o">/</span><span class="n">histo</span> <span class="mi">100</span> <span class="mi">1</span> <span class="mi">100</span> <span class="mi">0</span> <span class="mi">20</span> <span class="mi">20</span> <span class="n">c_size</span> <span class="n">mode</span> <span class="n">vector</span> <span class="n">ave</span> <span class="n">running</span> <span class="n">beyond</span> <span class="n">ignore</span> <span class="n">file</span> <span class="n">tmp</span><span class="o">.</span><span class="n">histo</span>
			</pre></div>
			</div>
			<hr class="docutils" />
			</div>
			</div>
			<div class="section" id="setting-parameters-for-the-kspace-style-pppm-disp-command">
			<span id="howto-24"></span><h2>6.24. Setting parameters for the <a class="reference internal" href="kspace_style.html"><span class="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
			to specify this command explicitly.</p>
			<hr class="docutils" />
			</div>
			<div class="section" id="polarizable-models">
			<span id="howto-25"></span><h2>6.25. Polarizable models</h2>
			<p>In polarizable force fields the charge distributions in molecules and
			materials respond to their electrostatic environements. Polarizable
			systems can be simulated in LAMMPS using three methods:</p>
			<ul class="simple">
			<li>the fluctuating charge method, implemented in the <a class="reference internal" href="fix_qeq.html"><span class="doc">QEQ</span></a>
			package,</li>
			<li>the adiabatic core-shell method, implemented in the
			<a class="reference internal" href="#howto-26"><span class="std std-ref">CORESHELL</span></a> package,</li>
			<li>the thermalized Drude dipole method, implemented in the
			<a class="reference internal" href="#howto-27"><span class="std std-ref">USER-DRUDE</span></a> package.</li>
			</ul>
			<p>The fluctuating charge method calculates instantaneous charges on
			interacting atoms based on the electronegativity equalization
			principle. It is implemented in the <a class="reference internal" href="fix_qeq.html"><span class="doc">fix qeq</span></a> which is
			available in several variants. It is a relatively efficient technique
			since no additional particles are introduced. This method allows for
			charge transfer between molecules or atom groups. However, because the
			charges are located at the interaction sites, off-plane components of
			polarization cannot be represented in planar molecules or atom groups.</p>
			<p>The two other methods share the same basic idea: polarizable atoms are
			split into one core atom and one satellite particle (called shell or
			Drude particle) attached to it by a harmonic spring. Both atoms bear
			a charge and they represent collectively an induced electric dipole.
			These techniques are computationally more expensive than the QEq
			method because of additional particles and bonds. These two
			charge-on-spring methods differ in certain features, with the
			core-shell model being normally used for ionic/crystalline materials,
			whereas the so-called Drude model is normally used for molecular
			systems and fluid states.</p>
			<p>The core-shell model is applicable to crystalline materials where the
			high symmetry around each site leads to stable trajectories of the
			core-shell pairs. However, bonded atoms in molecules can be so close
			that a core would interact too strongly or even capture the Drude
			particle of a neighbor. The Drude dipole model is relatively more
			complex in order to remediate this and other issues. Specifically, the
			Drude model includes specific thermostating of the core-Drude pairs
			and short-range damping of the induced dipoles.</p>
			<p>The three polarization methods can be implemented through a
			self-consistent calculation of charges or induced dipoles at each
			timestep. In the fluctuating charge scheme this is done by the matrix
			inversion method in <a class="reference internal" href="fix_qeq.html"><span class="doc">fix qeq/point</span></a>, but for core-shell
			or Drude-dipoles the relaxed-dipoles technique would require an slow
			iterative procedure. These self-consistent solutions yield accurate
			trajectories since the additional degrees of freedom representing
			polarization are massless. An alternative is to attribute a mass to
			the additional degrees of freedom and perform time integration using
			an extended Lagrangian technique. For the fluctuating charge scheme
			this is done by <a class="reference internal" href="fix_qeq.html"><span class="doc">fix qeq/dynamic</span></a>, and for the
			charge-on-spring models by the methods outlined in the next two
			sections. The assignment of masses to the additional degrees of
			freedom can lead to unphysical trajectories if care is not exerted in
			choosing the parameters of the poarizable models and the simulation
			conditions.</p>
			<p>In the core-shell model the vibration of the shells is kept faster
			than the ionic vibrations to mimic the fast response of the
			polarizable electrons. But in molecular systems thermalizing the
			core-Drude pairs at temperatures comparable to the rest of the
			simulation leads to several problems (kinetic energy transfer, too
			short a timestep, etc.) In order to avoid these problems the relative
			motion of the Drude particles with respect to their cores is kept
			“cold” so the vibration of the core-Drude pairs is very slow,
			approaching the self-consistent regime. In both models the
			temperature is regulated using the velocities of the center of mass of
			core+shell (or Drude) pairs, but in the Drude model the actual
			relative core-Drude particle motion is thermostated separately as
			well.</p>
			<hr class="docutils" />
			</div>
			<div class="section" id="adiabatic-core-shell-model">
			<span id="howto-26"></span><h2>6.26. Adiabatic core/shell model</h2>
			<p>The adiabatic core-shell model by <a class="reference internal" href="pair_cs.html#mitchellfinchham"><span class="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, <a class="reference internal" href="atom_style.html"><span class="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>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="mi">432</span> <span class="n">atoms</span> <span class="c1"># core and shell atoms</span>
			<span class="mi">216</span> <span class="n">bonds</span> <span class="c1"># number of core/shell springs</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="mi">4</span> <span class="n">atom</span> <span class="n">types</span> <span class="c1"># 2 cores and 2 shells for Na and Cl</span>
			<span class="mi">2</span> <span class="n">bond</span> <span class="n">types</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="mf">0.0</span> <span class="mf">24.09597</span> <span class="n">xlo</span> <span class="n">xhi</span>
			<span class="mf">0.0</span> <span class="mf">24.09597</span> <span class="n">ylo</span> <span class="n">yhi</span>
			<span class="mf">0.0</span> <span class="mf">24.09597</span> <span class="n">zlo</span> <span class="n">zhi</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">Masses</span> <span class="c1"># core/shell mass ratio = 0.1</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="mi">1</span> <span class="mf">20.690784</span> <span class="c1"># Na core</span>
			<span class="mi">2</span> <span class="mf">31.90500</span> <span class="c1"># Cl core</span>
			<span class="mi">3</span> <span class="mf">2.298976</span> <span class="c1"># Na shell</span>
			<span class="mi">4</span> <span class="mf">3.54500</span> <span class="c1"># Cl shell</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">Atoms</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="mi">1</span> <span class="mi">1</span> <span class="mi">2</span> <span class="mf">1.5005</span> <span class="mf">0.00000000</span> <span class="mf">0.00000000</span> <span class="mf">0.00000000</span> <span class="c1"># core of core/shell pair 1</span>
			<span class="mi">2</span> <span class="mi">1</span> <span class="mi">4</span> <span class="o">-</span><span class="mf">2.5005</span> <span class="mf">0.00000000</span> <span class="mf">0.00000000</span> <span class="mf">0.00000000</span> <span class="c1"># shell of core/shell pair 1</span>
			<span class="mi">3</span> <span class="mi">2</span> <span class="mi">1</span> <span class="mf">1.5056</span> <span class="mf">4.01599500</span> <span class="mf">4.01599500</span> <span class="mf">4.01599500</span> <span class="c1"># core of core/shell pair 2</span>
			<span class="mi">4</span> <span class="mi">2</span> <span class="mi">3</span> <span class="o">-</span><span class="mf">0.5056</span> <span class="mf">4.01599500</span> <span class="mf">4.01599500</span> <span class="mf">4.01599500</span> <span class="c1"># shell of core/shell pair 2</span>
			<span class="p">(</span><span class="o">...</span><span class="p">)</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">Bonds</span> <span class="c1"># Bond topology for spring forces</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="mi">1</span> <span class="mi">2</span> <span class="mi">1</span> <span class="mi">2</span> <span class="c1"># spring for core/shell pair 1</span>
			<span class="mi">2</span> <span class="mi">2</span> <span class="mi">3</span> <span class="mi">4</span> <span class="c1"># spring for core/shell pair 2</span>
			<span class="p">(</span><span class="o">...</span><span class="p">)</span>
			</pre></div>
			</div>
			<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 <a class="reference internal" href="special_bonds.html"><span class="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
			<a class="reference internal" href="special_bonds.html"><span class="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 <a class="reference internal" href="special_bonds.html"><span class="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 <a class="reference internal" href="pair_cs.html"><span class="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 <a class="reference internal" href="kspace_style.html"><span class="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
			are used:</p>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">pair_style</span> <span class="n">born</span><span class="o">/</span><span class="n">coul</span><span class="o">/</span><span class="n">long</span><span class="o">/</span><span class="n">cs</span> <span class="mf">20.0</span> <span class="mf">20.0</span>
			<span class="n">pair_coeff</span> <span class="o"></span> <span class="o"></span> <span class="mf">0.0</span> <span class="mf">1.000</span> <span class="mf">0.00</span> <span class="mf">0.00</span> <span class="mf">0.00</span>
			<span class="n">pair_coeff</span> <span class="mi">3</span> <span class="mi">3</span> <span class="mf">487.0</span> <span class="mf">0.23768</span> <span class="mf">0.00</span> <span class="mf">1.05</span> <span class="mf">0.50</span> <span class="c1">#Na-Na</span>
			<span class="n">pair_coeff</span> <span class="mi">3</span> <span class="mi">4</span> <span class="mf">145134.0</span> <span class="mf">0.23768</span> <span class="mf">0.00</span> <span class="mf">6.99</span> <span class="mf">8.70</span> <span class="c1">#Na-Cl</span>
			<span class="n">pair_coeff</span> <span class="mi">4</span> <span class="mi">4</span> <span class="mf">405774.0</span> <span class="mf">0.23768</span> <span class="mf">0.00</span> <span class="mf">72.40</span> <span class="mf">145.40</span> <span class="c1">#Cl-Cl</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">bond_style</span> <span class="n">harmonic</span>
			<span class="n">bond_coeff</span> <span class="mi">1</span> <span class="mf">63.014</span> <span class="mf">0.0</span>
			<span class="n">bond_coeff</span> <span class="mi">2</span> <span class="mf">25.724</span> <span class="mf">0.0</span>
			</pre></div>
			</div>
			<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 <a class="reference internal" href="compute_temp_cs.html"><span class="doc">compute temp/cs</span></a> command can be used, in conjunction with
			any of the thermostat fixes, such as <a class="reference internal" href="fix_nh.html"><span class="doc">fix nvt</span></a> or <a class="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
			<a class="reference internal" href="compute_temp_cs.html"><span class="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 <a class="reference internal" href="group.html"><span class="doc">group type</span></a> command.
			Note that to perform thermostatting using this definition of
			temperature, the <a class="reference internal" href="fix_modify.html"><span class="doc">fix modify temp</span></a> command should be
			used to assign the compute to the thermostat fix. Likewise the
			<a class="reference internal" href="thermo_modify.html"><span class="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>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">group</span> <span class="n">cores</span> <span class="nb">type</span> <span class="mi">1</span> <span class="mi">2</span>
			<span class="n">group</span> <span class="n">shells</span> <span class="nb">type</span> <span class="mi">3</span> <span class="mi">4</span>
			<span class="n">compute</span> <span class="n">CSequ</span> <span class="nb">all</span> <span class="n">temp</span><span class="o">/</span><span class="n">cs</span> <span class="n">cores</span> <span class="n">shells</span>
			<span class="n">fix</span> <span class="n">thermoberendsen</span> <span class="nb">all</span> <span class="n">temp</span><span class="o">/</span><span class="n">berendsen</span> <span class="mi">1427</span> <span class="mi">1427</span> <span class="mf">0.4</span> <span class="c1"># thermostat for the true physical system</span>
			<span class="n">fix</span> <span class="n">thermostatequ</span> <span class="nb">all</span> <span class="n">nve</span> <span class="c1"># integrator as needed for the berendsen thermostat</span>
			<span class="n">fix_modify</span> <span class="n">thermoberendsen</span> <span class="n">temp</span> <span class="n">CSequ</span>
			<span class="n">thermo_modify</span> <span class="n">temp</span> <span class="n">CSequ</span> <span class="c1"># output of center-of-mass derived temperature</span>
			</pre></div>
			</div>
			<p>If <a class="reference internal" href="compute_temp_cs.html"><span class="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 <a class="reference internal" href="fix_momentum.html"><span class="doc">fix momentum</span></a> command in combination with <a class="reference internal" href="compute_temp_cs.html"><span class="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
			<a class="reference internal" href="velocity.html"><span class="doc">velocity create</span></a> command and assigning the <a class="reference internal" href="compute_temp_cs.html"><span class="doc">compute temp/cs</span></a> command to the <em>temp</em> keyword of the
			<a class="reference internal" href="velocity.html"><span class="doc">velocity</span></a> commmand, e.g.</p>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">velocity</span> <span class="nb">all</span> <span class="n">create</span> <span class="mi">1427</span> <span class="mi">134</span> <span class="n">bias</span> <span class="n">yes</span> <span class="n">temp</span> <span class="n">CSequ</span>
			<span class="n">velocity</span> <span class="nb">all</span> <span class="n">scale</span> <span class="mi">1427</span> <span class="n">temp</span> <span class="n">CSequ</span>
			</pre></div>
			</div>
			<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 <a class="reference internal" href="compute_temp_cs.html"><span class="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 <a class="reference internal" href="compute_chunk_atom.html"><span class="doc">compute chunk/atom</span></a> and <a class="reference internal" href="compute_temp_chunk.html"><span class="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 <a class="reference internal" href="variable.html"><span class="doc">variable</span></a> command. Or they can
			be time/averaged and output using the <a class="reference internal" href="fix_ave_time.html"><span class="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 <a class="reference internal" href="fix_property_atom.html"><span class="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 <a class="reference internal" href="read_data.html"><span class="doc">read_data</span></a> command. This field can then be
			accessed by the <a class="reference internal" href="compute_property_atom.html"><span class="doc">compute property/atom</span></a>
			command, to use as input to the <a class="reference internal" href="compute_chunk_atom.html"><span class="doc">compute chunk/atom</span></a> command to define the core/shell
			pairs as chunks.</p>
			<p>For example,</p>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">fix</span> <span class="n">csinfo</span> <span class="nb">all</span> <span class="nb">property</span><span class="o">/</span><span class="n">atom</span> <span class="n">i_CSID</span> <span class="c1"># property/atom command</span>
			<span class="n">read_data</span> <span class="n">NaCl_CS_x0</span><span class="o">.</span><span class="mi">1</span><span class="n">_prop</span><span class="o">.</span><span class="n">data</span> <span class="n">fix</span> <span class="n">csinfo</span> <span class="n">NULL</span> <span class="n">CS</span><span class="o">-</span><span class="n">Info</span> <span class="c1"># atom property added in the data-file</span>
			<span class="n">compute</span> <span class="n">prop</span> <span class="nb">all</span> <span class="nb">property</span><span class="o">/</span><span class="n">atom</span> <span class="n">i_CSID</span>
			<span class="n">compute</span> <span class="n">cs_chunk</span> <span class="nb">all</span> <span class="n">chunk</span><span class="o">/</span><span class="n">atom</span> <span class="n">c_prop</span>
			<span class="n">compute</span> <span class="n">cstherm</span> <span class="nb">all</span> <span class="n">temp</span><span class="o">/</span><span class="n">chunk</span> <span class="n">cs_chunk</span> <span class="n">temp</span> <span class="n">internal</span> <span class="n">com</span> <span class="n">yes</span> <span class="n">cdof</span> <span class="mf">3.0</span> <span class="c1"># note the chosen degrees of freedom for the core/shell pairs</span>
			<span class="n">fix</span> <span class="n">ave_chunk</span> <span class="nb">all</span> <span class="n">ave</span><span class="o">/</span><span class="n">time</span> <span class="mi">10</span> <span class="mi">1</span> <span class="mi">10</span> <span class="n">c_cstherm</span> <span class="n">file</span> <span class="n">chunk</span><span class="o">.</span><span class="n">dump</span> <span class="n">mode</span> <span class="n">vector</span>
			</pre></div>
			</div>
			<p>The additional section in the date file would be formatted like this:</p>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">CS</span><span class="o">-</span><span class="n">Info</span> <span class="c1"># header of additional section</span>
			</pre></div>
			</div>
			<div class="highlight-default"><div class="highlight"><pre><span></span><span class="mi">1</span> <span class="mi">1</span> <span class="c1"># column 1 = atom ID, column 2 = core/shell ID</span>
			<span class="mi">2</span> <span class="mi">1</span>
			<span class="mi">3</span> <span class="mi">2</span>
			<span class="mi">4</span> <span class="mi">2</span>
			<span class="mi">5</span> <span class="mi">3</span>
			<span class="mi">6</span> <span class="mi">3</span>
			<span class="mi">7</span> <span class="mi">4</span>
			<span class="mi">8</span> <span class="mi">4</span>
			<span class="p">(</span><span class="o">...</span><span class="p">)</span>
			</pre></div>
			</div>
			<hr class="docutils" />
			</div>
			<div class="section" id="drude-induced-dipoles">
			<span id="howto-27"></span><h2>6.27. Drude induced dipoles</h2>
			<p>The thermalized Drude model, similarly to the <a class="reference internal" href="#howto-26"><span class="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
			(<a class="reference internal" href="#lamoureux"><span class="std std-ref">Lamoureux and Roux</span></a>):</p>
			<ul class="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 <a class="reference internal" href="tutorial_drude.html"><span class="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 <a class="reference internal" href="fix_drude.html"><span class="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
			<a class="reference internal" href="special_bonds.html"><span class="doc">special_bonds</span></a> command. If using <a class="reference internal" href="fix_shake.html"><span class="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 <a class="reference internal" href="fix_langevin_drude.html"><span class="doc">fix langevin/drude</span></a>
			for a Langevin thermostat, or <a class="reference internal" href="fix_drude_transform.html"><span class="doc">fix drude/transform/*</span></a> for a Nose-Hoover
			thermostat. The former requires use of the command <a class="reference internal" href="comm_modify.html"><span class="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 <a class="reference internal" href="compute_temp_drude.html"><span class="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 <a class="reference internal" href="pair_thole.html"><span class="doc">pair style thole</span></a> in <a class="reference internal" href="pair_hybrid.html"><span class="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
			too close, which can cause numerical issues.</p>
			<p id="berendsen"><strong>(Berendsen)</strong> Berendsen, Grigera, Straatsma, J Phys Chem, 91,
			6269-6271 (1987).</p>
			<p id="cornell"><strong>(Cornell)</strong> Cornell, Cieplak, Bayly, Gould, Merz, Ferguson,
			Spellmeyer, Fox, Caldwell, Kollman, JACS 117, 5179-5197 (1995).</p>
			<p id="horn"><strong>(Horn)</strong> Horn, Swope, Pitera, Madura, Dick, Hura, and Head-Gordon,
			J Chem Phys, 120, 9665 (2004).</p>
			<p id="ikeshoji"><strong>(Ikeshoji)</strong> Ikeshoji and Hafskjold, Molecular Physics, 81, 251-261
			(1994).</p>
			<p id="mackerell"><strong>(MacKerell)</strong> MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field,
			Fischer, Gao, Guo, Ha, et al, J Phys Chem, 102, 3586 (1998).</p>
			<p id="mayo"><strong>(Mayo)</strong> Mayo, Olfason, Goddard III, J Phys Chem, 94, 8897-8909
			(1990).</p>
			<p id="jorgensen"><strong>(Jorgensen)</strong> Jorgensen, Chandrasekhar, Madura, Impey, Klein, J Chem
			Phys, 79, 926 (1983).</p>
			<p id="price"><strong>(Price)</strong> Price and Brooks, J Chem Phys, 121, 10096 (2004).</p>
			<p id="shinoda"><strong>(Shinoda)</strong> Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004).</p>
			<p id="mitchellfinchham"><strong>(Mitchell and Finchham)</strong> Mitchell, Finchham, J Phys Condensed Matter,
			5, 1031-1038 (1993).</p>
			<p id="lamoureux"><strong>(Lamoureux and Roux)</strong> G. Lamoureux, B. Roux, J. Chem. Phys 119, 3025 (2003)</p>
			</div>
			</div>


			</div>
			</div>
			<footer>

			<div class="rst-footer-buttons" role="navigation" aria-label="footer navigation">

			<a href="Section_example.html" class="btn btn-neutral float-right" title="7. Example problems" accesskey="n">Next <span class="fa fa-arrow-circle-right"></span></a>


			<a href="Section_accelerate.html" class="btn btn-neutral" title="5. Accelerating LAMMPS performance" accesskey="p"><span class="fa fa-arrow-circle-left"></span> Previous</a>

			</div>


			<hr/>

			<div role="contentinfo">
			<p>
			© Copyright 2013 Sandia Corporation.
			</p>
			</div>
			Built with <a href="http://sphinx-doc.org/">Sphinx</a> using a <a href="https://github.com/snide/sphinx_rtd_theme">theme</a> provided by <a href="https://readthedocs.org">Read the Docs</a>.

			</footer>

			</div>
			</div>

			</section>

			</div>





			<script type="text/javascript">
			var DOCUMENTATION_OPTIONS = {
			URL_ROOT:'./',
			VERSION:'',
			COLLAPSE_INDEX:false,
			FILE_SUFFIX:'.html',
			HAS_SOURCE: true
			};
			</script>
			<script type="text/javascript" src="_static/jquery.js"></script>
			<script type="text/javascript" src="_static/underscore.js"></script>
			<script type="text/javascript" src="_static/doctools.js"></script>
			<script type="text/javascript" src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
			<script type="text/javascript" src="_static/sphinxcontrib-images/LightBox2/lightbox2/js/jquery-1.11.0.min.js"></script>
			<script type="text/javascript" src="_static/sphinxcontrib-images/LightBox2/lightbox2/js/lightbox.min.js"></script>
			<script type="text/javascript" src="_static/sphinxcontrib-images/LightBox2/lightbox2-customize/jquery-noconflict.js"></script>





			<script type="text/javascript" src="_static/js/theme.js"></script>




			<script type="text/javascript">
			jQuery(function () {
			SphinxRtdTheme.StickyNav.enable();
			});
			</script>


			</body>
			</html>

lammps/doc/html/Section_howto.html1d3f94fbb20clammps-icms

Section_howto.htmlView Options

lammps/doc/html/Section_howto.html
1d3f94fbb20clammps-icms

Section_howto.html
View Options