forked from lijiext/lammps
2890 lines
189 KiB
HTML
2890 lines
189 KiB
HTML
|
|
|
|
<!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<a class="headerlink" href="#how-to-discussions" title="Permalink to this headline">¶</a></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>Restarting a simulation</span></a></div>
|
|
<div class="line">6.2 <a class="reference internal" href="#howto-2"><span>2d simulations</span></a></div>
|
|
<div class="line">6.3 <a class="reference internal" href="#howto-3"><span>CHARMM, AMBER, and DREIDING force fields</span></a></div>
|
|
<div class="line">6.4 <a class="reference internal" href="#howto-4"><span>Running multiple simulations from one input script</span></a></div>
|
|
<div class="line">6.5 <a class="reference internal" href="#howto-5"><span>Multi-replica simulations</span></a></div>
|
|
<div class="line">6.6 <a class="reference internal" href="#howto-6"><span>Granular models</span></a></div>
|
|
<div class="line">6.7 <a class="reference internal" href="#howto-7"><span>TIP3P water model</span></a></div>
|
|
<div class="line">6.8 <a class="reference internal" href="#howto-8"><span>TIP4P water model</span></a></div>
|
|
<div class="line">6.9 <a class="reference internal" href="#howto-9"><span>SPC water model</span></a></div>
|
|
<div class="line">6.10 <a class="reference internal" href="#howto-10"><span>Coupling LAMMPS to other codes</span></a></div>
|
|
<div class="line">6.11 <a class="reference internal" href="#howto-11"><span>Visualizing LAMMPS snapshots</span></a></div>
|
|
<div class="line">6.12 <a class="reference internal" href="#howto-12"><span>Triclinic (non-orthogonal) simulation boxes</span></a></div>
|
|
<div class="line">6.13 <a class="reference internal" href="#howto-13"><span>NEMD simulations</span></a></div>
|
|
<div class="line">6.14 <a class="reference internal" href="#howto-14"><span>Finite-size spherical and aspherical particles</span></a></div>
|
|
<div class="line">6.15 <a class="reference internal" href="#howto-15"><span>Output from LAMMPS (thermo, dumps, computes, fixes, variables)</span></a></div>
|
|
<div class="line">6.16 <a class="reference internal" href="#howto-16"><span>Thermostatting, barostatting and computing temperature</span></a></div>
|
|
<div class="line">6.17 <a class="reference internal" href="#howto-17"><span>Walls</span></a></div>
|
|
<div class="line">6.18 <a class="reference internal" href="#howto-18"><span>Elastic constants</span></a></div>
|
|
<div class="line">6.19 <a class="reference internal" href="#howto-19"><span>Library interface to LAMMPS</span></a></div>
|
|
<div class="line">6.20 <a class="reference internal" href="#howto-20"><span>Calculating thermal conductivity</span></a></div>
|
|
<div class="line">6.21 <a class="reference internal" href="#howto-21"><span>Calculating viscosity</span></a></div>
|
|
<div class="line">6.22 <a class="reference internal" href="#howto-22"><span>Calculating a diffusion coefficient</span></a></div>
|
|
<div class="line">6.23 <a class="reference internal" href="#howto-23"><span>Using chunks to calculate system properties</span></a></div>
|
|
<div class="line">6.24 <a class="reference internal" href="#howto-24"><span>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>Polarizable models</span></a></div>
|
|
<div class="line">6.26 <a class="reference internal" href="#howto-26"><span>Adiabatic core/shell model</span></a></div>
|
|
<div class="line">6.27 <a class="reference internal" href="#howto-27"><span>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"><em>Section_example</em></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<a class="headerlink" href="#restarting-a-simulation" title="Permalink to this headline">¶</a></h2>
|
|
<p>There are 3 ways to continue a long LAMMPS simulation. Multiple
|
|
<a class="reference internal" href="run.html"><em>run</em></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"><em>restart</em></a>
|
|
command. At a later time, these binary files can be read via a
|
|
<a class="reference internal" href="read_restart.html"><em>read_restart</em></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>-r command-line switch</span></a> and read by a
|
|
<a class="reference internal" href="read_data.html"><em>read_data</em></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"><em>read_restart</em></a> and
|
|
<a class="reference internal" href="read_data.html"><em>read_data</em></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-python"><div class="highlight"><pre>restart 50 tmp.restart
|
|
</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-python"><div class="highlight"><pre>read_restart tmp.restart.50
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>neighbor 0.4 bin
|
|
neigh_modify every 1 delay 1
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>fix 1 all nve
|
|
fix 2 all langevin 1.0 1.0 10.0 904297
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>timestep 0.012
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>run 50
|
|
</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"><em>fix langevin</em></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-python"><div class="highlight"><pre>lmp_g++ -r tmp.restart.50 tmp.restart.data
|
|
</pre></div>
|
|
</div>
|
|
<p>Then, this script could be used to re-run the last 50 steps:</p>
|
|
<div class="highlight-python"><div class="highlight"><pre>units lj
|
|
atom_style bond
|
|
pair_style lj/cut 1.12
|
|
pair_modify shift yes
|
|
bond_style fene
|
|
special_bonds 0.0 1.0 1.0
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>read_data tmp.restart.data
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>neighbor 0.4 bin
|
|
neigh_modify every 1 delay 1
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>fix 1 all nve
|
|
fix 2 all langevin 1.0 1.0 10.0 904297
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>timestep 0.012
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>reset_timestep 50
|
|
run 50
|
|
</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"><em>reset_timestep</em></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<a class="headerlink" href="#d-simulations" title="Permalink to this headline">¶</a></h2>
|
|
<p>Use the <a class="reference internal" href="dimension.html"><em>dimension</em></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"><em>boundary</em></a>
|
|
command. This is the default.</p>
|
|
<p>If using the <a class="reference internal" href="create_box.html"><em>create box</em></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"><em>create box</em></a> 1 -10 10 -10 10 -0.25 0.25
|
|
</pre>
|
|
<p>If using the <a class="reference internal" href="read_data.html"><em>read data</em></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"><em>fix enforce2d</em></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<a class="headerlink" href="#charmm-amber-and-dreiding-force-fields" title="Permalink to this headline">¶</a></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"><em>read_data</em></a> command or in the input script with
|
|
commands like <a class="reference internal" href="pair_coeff.html"><em>pair_coeff</em></a> or
|
|
<a class="reference internal" href="bond_coeff.html"><em>bond_coeff</em></a>. See <a class="reference internal" href="Section_tools.html"><em>Section_tools</em></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="special_bonds.html#mackerell"><span>(MacKerell)</span></a> for a description of the CHARMM force
|
|
field. See <a class="reference internal" href="special_bonds.html#cornell"><span>(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"><em>bond_style</em></a> harmonic</li>
|
|
<li><a class="reference internal" href="angle_charmm.html"><em>angle_style</em></a> charmm</li>
|
|
<li><a class="reference internal" href="dihedral_charmm.html"><em>dihedral_style</em></a> charmm</li>
|
|
<li><a class="reference internal" href="pair_charmm.html"><em>pair_style</em></a> lj/charmm/coul/charmm</li>
|
|
<li><a class="reference internal" href="pair_charmm.html"><em>pair_style</em></a> lj/charmm/coul/charmm/implicit</li>
|
|
<li><a class="reference internal" href="pair_charmm.html"><em>pair_style</em></a> lj/charmm/coul/long</li>
|
|
<li><a class="reference internal" href="special_bonds.html"><em>special_bonds</em></a> charmm</li>
|
|
<li><a class="reference internal" href="special_bonds.html"><em>special_bonds</em></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"><em>explicit hydrogen bond term</em></a> to describe
|
|
interactions involving a hydrogen atom on very electronegative atoms
|
|
(N, O, F).</p>
|
|
<p>See <a class="reference internal" href="special_bonds.html#mayo"><span>(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"><em>bond_style</em></a> harmonic</li>
|
|
<li><a class="reference internal" href="bond_morse.html"><em>bond_style</em></a> morse</li>
|
|
<li><a class="reference internal" href="angle_harmonic.html"><em>angle_style</em></a> harmonic</li>
|
|
<li><a class="reference internal" href="angle_cosine.html"><em>angle_style</em></a> cosine</li>
|
|
<li><a class="reference internal" href="angle_cosine_periodic.html"><em>angle_style</em></a> cosine/periodic</li>
|
|
<li><a class="reference internal" href="dihedral_charmm.html"><em>dihedral_style</em></a> charmm</li>
|
|
<li><a class="reference internal" href="improper_umbrella.html"><em>improper_style</em></a> umbrella</li>
|
|
<li><a class="reference internal" href="pair_buck.html"><em>pair_style</em></a> buck</li>
|
|
<li><a class="reference internal" href="pair_buck.html"><em>pair_style</em></a> buck/coul/cut</li>
|
|
<li><a class="reference internal" href="pair_buck.html"><em>pair_style</em></a> buck/coul/long</li>
|
|
<li><a class="reference internal" href="pair_lj.html"><em>pair_style</em></a> lj/cut</li>
|
|
<li><a class="reference internal" href="pair_lj.html"><em>pair_style</em></a> lj/cut/coul/cut</li>
|
|
<li><a class="reference internal" href="pair_lj.html"><em>pair_style</em></a> lj/cut/coul/long</li>
|
|
<li><a class="reference internal" href="pair_hbond_dreiding.html"><em>pair_style</em></a> hbond/dreiding/lj</li>
|
|
<li><a class="reference internal" href="pair_hbond_dreiding.html"><em>pair_style</em></a> hbond/dreiding/morse</li>
|
|
<li><a class="reference internal" href="special_bonds.html"><em>special_bonds</em></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<a class="headerlink" href="#running-multiple-simulations-from-one-input-script" title="Permalink to this headline">¶</a></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"><em>run</em></a> command
|
|
multiple times. For example, this script</p>
|
|
<div class="highlight-python"><div class="highlight"><pre>units lj
|
|
atom_style atomic
|
|
read_data data.lj
|
|
run 10000
|
|
run 10000
|
|
run 10000
|
|
run 10000
|
|
run 10000
|
|
</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"><em>clear</em></a> command can be used in between them to
|
|
re-initialize LAMMPS. For example, this script</p>
|
|
<div class="highlight-python"><div class="highlight"><pre>units lj
|
|
atom_style atomic
|
|
read_data data.lj
|
|
run 10000
|
|
clear
|
|
units lj
|
|
atom_style atomic
|
|
read_data data.lj.new
|
|
run 10000
|
|
</pre></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"><em>variables</em></a> and the <a class="reference internal" href="next.html"><em>next</em></a> and
|
|
<a class="reference internal" href="jump.html"><em>jump</em></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-python"><div class="highlight"><pre>variable d index run1 run2 run3 run4 run5 run6 run7 run8
|
|
shell cd $d
|
|
read_data data.polymer
|
|
run 10000
|
|
shell cd ..
|
|
clear
|
|
next d
|
|
jump in.polymer
|
|
</pre></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-python"><div class="highlight"><pre>variable a loop 8
|
|
variable t index 0.8 0.85 0.9 0.95 1.0 1.05 1.1 1.15
|
|
log log.$a
|
|
read data.polymer
|
|
velocity all create $t 352839
|
|
fix 1 all nvt $t $t 100.0
|
|
dump 1 all atom 1000 dump.$a
|
|
run 100000
|
|
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>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"><em>variable</em></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<a class="headerlink" href="#multi-replica-simulations" title="Permalink to this headline">¶</a></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"><em>neb</em></a> for nudged elastic band calculations</li>
|
|
<li><a class="reference internal" href="prd.html"><em>prd</em></a> for parallel replica dynamics</li>
|
|
<li><a class="reference internal" href="tad.html"><em>tad</em></a> for temperature accelerated dynamics</li>
|
|
<li><a class="reference internal" href="temper.html"><em>temper</em></a> for parallel tempering</li>
|
|
<li><a class="reference internal" href="fix_pimd.html"><em>fix pimd</em></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>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>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>-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-python"><div class="highlight"><pre>mpirun -np 16 lmp_linux -partition 8x2 -in in.temper
|
|
mpirun -np 8 lmp_linux -partition 8x1 -in in.neb
|
|
</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>-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<a class="headerlink" href="#granular-models" title="Permalink to this headline">¶</a></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"><em>atom_style sphere</em></a></li>
|
|
<li><a class="reference internal" href="fix_nve_sphere.html"><em>fix nve/sphere</em></a></li>
|
|
<li><a class="reference internal" href="fix_gravity.html"><em>fix gravity</em></a></li>
|
|
</ul>
|
|
<p>This compute</p>
|
|
<ul class="simple">
|
|
<li><a class="reference internal" href="compute_erotate_sphere.html"><em>compute erotate/sphere</em></a></li>
|
|
</ul>
|
|
<p>calculates rotational kinetic energy which can be <a class="reference internal" href="#howto-15"><span>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"><em>pair_style</em></a> gran/history</li>
|
|
<li><a class="reference internal" href="pair_style.html"><em>pair_style</em></a> gran/no_history</li>
|
|
<li><a class="reference internal" href="pair_style.html"><em>pair_style</em></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"><em>fix freeze</em></a></li>
|
|
<li><a class="reference internal" href="fix_pour.html"><em>fix pour</em></a></li>
|
|
<li><a class="reference internal" href="fix_viscous.html"><em>fix viscous</em></a></li>
|
|
<li><a class="reference internal" href="fix_wall_gran.html"><em>fix wall/gran</em></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"><em>neigh_modify</em></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<a class="headerlink" href="#tip3p-water-model" title="Permalink to this headline">¶</a></h2>
|
|
<p>The TIP3P water model as implemented in CHARMM
|
|
<a class="reference internal" href="special_bonds.html#mackerell"><span>(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"><em>fix shake</em></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>(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>(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<a class="headerlink" href="#tip4p-water-model" title="Permalink to this headline">¶</a></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"><em>pair_style lj/cut/tip4p/cut</em></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"><em>pair_style lj/cut/tip4p/long</em></a></li>
|
|
<li><a class="reference internal" href="kspace_style.html"><em>kspace_style pppm/tip4p</em></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"><em>fix shake</em></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>(Jorgensen)</span></a>. Note that the OM distance is specified in
|
|
the <a class="reference internal" href="pair_style.html"><em>pair_style</em></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"><em>pair_style lj/cut/tip4p/long</em></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<a class="headerlink" href="#spc-water-model" title="Permalink to this headline">¶</a></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"><em>fix shake</em></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="fix_temp_berendsen.html#berendsen"><span>(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<a class="headerlink" href="#coupling-lammps-to-other-codes" title="Permalink to this headline">¶</a></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"><em>fix</em></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"><em>fix_poems</em></a> command for more details. See <a class="reference internal" href="Section_modify.html"><em>this section</em></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"><em>run</em></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"><em>Section_modify</em></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"><em>run</em></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>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"><em>Section_python</em></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>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<a class="headerlink" href="#visualizing-lammps-snapshots" title="Permalink to this headline">¶</a></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"><em>dump</em></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>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"><em>Section_tools</em></a> doc page for details. The first is
|
|
the <a class="reference internal" href="Section_tools.html#charmm"><span>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>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>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"><em>dump</em></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"><em>dump</em></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<a class="headerlink" href="#triclinic-non-orthogonal-simulation-boxes" title="Permalink to this headline">¶</a></h2>
|
|
<p>By default, LAMMPS uses an orthogonal simulation box to encompass the
|
|
particles. The <a class="reference internal" href="boundary.html"><em>boundary</em></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"><em>create_box</em></a> or
|
|
<a class="reference internal" href="read_data.html"><em>read_data</em></a> or <a class="reference internal" href="read_restart.html"><em>read_restart</em></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"><em>thermo_style custom</em></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>,**B**,**C** are the edge
|
|
vectors of a general parallelepiped, where there is no restriction on
|
|
<strong>A</strong>,**B**,**C** 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>,**b**,**c** 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="#id73"><span class="problematic" id="id74">|**A**|</span></a> 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"><em>fix deform</em></a> and <a class="reference internal" href="fix_nh.html"><em>fix npt</em></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"><em>create_box</em></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"><em>region</em></a> command for details. If the
|
|
<a class="reference internal" href="read_data.html"><em>read_data</em></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"><em>read_data</em></a> command for details. Finally, if the
|
|
<a class="reference internal" href="read_restart.html"><em>read_restart</em></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"><em>fix npt</em></a> or
|
|
<a class="reference internal" href="fix_deform.html"><em>fix deform</em></a> commands. Alternatively, you can use the
|
|
<a class="reference internal" href="change_box.html"><em>change_box</em></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"><em>thermo_style custom</em></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"><em>create_box</em></a> or <a class="reference internal" href="read_data.html"><em>read_data</em></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"><em>fix deform</em></a> or <a class="reference internal" href="fix_nh.html"><em>fix npt</em></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"><em>fix deform</em></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"><em>fix deform</em></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"><em>box</em></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"><em>fix deform</em></a> or
|
|
<a class="reference internal" href="fix_nh.html"><em>fix npt</em></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"><em>thermo_style custom</em></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"><em>dump</em></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-python"><div class="highlight"><pre>ITEM: BOX BOUNDS xy xz yz
|
|
xlo_bound xhi_bound xy
|
|
ylo_bound yhi_bound xz
|
|
zlo_bound zhi_bound yz
|
|
</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-python"><div class="highlight"><pre><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"><em>lattice</em></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"><em>create_atoms</em></a> command.</p>
|
|
<p>A second use is to run Parinello-Rahman dyanamics via the <a class="reference internal" href="fix_nh.html"><em>fix npt</em></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"><em>energy minimization</em></a> is
|
|
the <a class="reference internal" href="fix_box_relax.html"><em>fix box/relax</em></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"><em>fix deform</em></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<a class="headerlink" href="#nemd-simulations" title="Permalink to this headline">¶</a></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"><em>fix deform</em></a> command. The
|
|
<a class="reference internal" href="fix_nvt_sllod.html"><em>fix nvt/sllod</em></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"><em>compute temp/deform</em></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"><em>fix ave/spatial</em></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"><em>fix deform</em></a> can continuously strain
|
|
a box by an arbitrary amount. As discussed in the <a class="reference internal" href="fix_deform.html"><em>fix deform</em></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"><em>fix deform</em></a> should be set to “remap v”, since that is what
|
|
<a class="reference internal" href="fix_nvt_sllod.html"><em>fix nvt/sllod</em></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"><em>fix viscosity</em></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<a class="headerlink" href="#finite-size-spherical-and-aspherical-particles" title="Permalink to this headline">¶</a></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"><em>examples directory</em></a> in the LAMMPS distribution.</p>
|
|
<div class="section" id="atom-styles">
|
|
<h3>6.14.1. Atom styles<a class="headerlink" href="#atom-styles" title="Permalink to this headline">¶</a></h3>
|
|
<p>There are several <a class="reference internal" href="atom_style.html"><em>atom styles</em></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-python"><div class="highlight"><pre>atom_style hybrid sphere dipole
|
|
</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"><em>atom_style dipole</em></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"><em>Peridynamic models</em></a> and
|
|
defines particles as having a volume, that is used internally in the
|
|
<a class="reference internal" href="pair_peri.html"><em>pair_style peri</em></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"><em>body</em></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"><em>atom_style hybrid</em></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"><em>pair_style hybrid</em></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"><em>2d simulations</em></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<a class="headerlink" href="#pair-potentials" title="Permalink to this headline">¶</a></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"><em>pair styles</em></a> that generate torque:</p>
|
|
<ul class="simple">
|
|
<li><a class="reference internal" href="pair_gran.html"><em>pair_style gran/history</em></a></li>
|
|
<li><a class="reference internal" href="pair_gran.html"><em>pair_style gran/hertzian</em></a></li>
|
|
<li><a class="reference internal" href="pair_gran.html"><em>pair_style gran/no_history</em></a></li>
|
|
<li><a class="reference internal" href="pair_dipole.html"><em>pair_style dipole/cut</em></a></li>
|
|
<li><a class="reference internal" href="pair_gayberne.html"><em>pair_style gayberne</em></a></li>
|
|
<li><a class="reference internal" href="pair_resquared.html"><em>pair_style resquared</em></a></li>
|
|
<li><a class="reference internal" href="pair_brownian.html"><em>pair_style brownian</em></a></li>
|
|
<li><a class="reference internal" href="pair_lubricate.html"><em>pair_style lubricate</em></a></li>
|
|
<li><a class="reference internal" href="pair_line_lj.html"><em>pair_style line/lj</em></a></li>
|
|
<li><a class="reference internal" href="pair_tri_lj.html"><em>pair_style tri/lj</em></a></li>
|
|
<li><a class="reference internal" href="pair_body.html"><em>pair_style body</em></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<a class="headerlink" href="#time-integration" title="Permalink to this headline">¶</a></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"><em>fix nve/sphere</em></a></li>
|
|
<li><a class="reference internal" href="fix_nvt_sphere.html"><em>fix nvt/sphere</em></a></li>
|
|
<li><a class="reference internal" href="fix_npt_sphere.html"><em>fix npt/sphere</em></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"><em>fix nve/asphere</em></a></li>
|
|
<li><a class="reference internal" href="fix_nvt_asphere.html"><em>fix nvt/asphere</em></a></li>
|
|
<li><a class="reference internal" href="fix_npt_asphere.html"><em>fix npt/asphere</em></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"><em>fix nve/line</em></a></li>
|
|
<li><a class="reference internal" href="fix_nve_tri.html"><em>fix nve/tri</em></a></li>
|
|
<li><a class="reference internal" href="fix_nve_body.html"><em>fix nve/body</em></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"><em>groups</em></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<a class="headerlink" href="#computes-thermodynamics-and-dump-output" title="Permalink to this headline">¶</a></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"><em>compute temp/sphere</em></a></li>
|
|
<li><a class="reference internal" href="compute_temp_asphere.html"><em>compute temp/asphere</em></a></li>
|
|
<li><a class="reference internal" href="compute_erotate_sphere.html"><em>compute erotate/sphere</em></a></li>
|
|
<li><a class="reference internal" href="compute_erotate_asphere.html"><em>compute erotate/asphere</em></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"><em>thermo_modify</em></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"><em>thermo_style</em></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"><em>dump custom</em></a></li>
|
|
<li><a class="reference internal" href="compute_property_atom.html"><em>compute property/atom</em></a></li>
|
|
<li><a class="reference internal" href="dump.html"><em>dump local</em></a></li>
|
|
<li><a class="reference internal" href="compute_body_local.html"><em>compute body/local</em></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<a class="headerlink" href="#rigid-bodies-composed-of-finite-size-particles" title="Permalink to this headline">¶</a></h3>
|
|
<p>The <a class="reference internal" href="fix_rigid.html"><em>fix rigid</em></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"><em>fix shake</em></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"><em>fix rigid</em></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"><em>fix nve/body</em></a>. Interactions between pairs of body
|
|
particles are computed via a command like <a class="reference internal" href="pair_body.html"><em>pair_style body</em></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)<a class="headerlink" href="#output-from-lammps-thermo-dumps-computes-fixes-variables" title="Permalink to this headline">¶</a></h2>
|
|
<p>There are four basic kinds of LAMMPS output:</p>
|
|
<ul class="simple">
|
|
<li><a class="reference internal" href="thermo_style.html"><em>Thermodynamic output</em></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"><em>Dump files</em></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"><em>fix ave/time</em></a> for time averaging, <a class="reference internal" href="fix_ave_spatial.html"><em>fix ave/spatial</em></a> for spatial averaging, and <a class="reference internal" href="fix_print.html"><em>fix print</em></a> for single-line output of
|
|
<a class="reference internal" href="variable.html"><em>variables</em></a>. Fix print can also output to the
|
|
screen.</li>
|
|
<li><a class="reference internal" href="restart.html"><em>Restart files</em></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"><em>dump</em></a> and <a class="reference internal" href="fix.html"><em>fix</em></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"><em>add their own computes and fixes to LAMMPS</em></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>Global/per-atom/local data</span></a></li>
|
|
<li><a class="reference internal" href="#scalar"><span>Scalar/vector/array data</span></a></li>
|
|
<li><a class="reference internal" href="#thermo"><span>Thermodynamic output</span></a></li>
|
|
<li><a class="reference internal" href="#dump"><span>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>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>Computes that generate values to output</span></a></li>
|
|
<li><a class="reference internal" href="#fix"><span>Fixes that generate values to output</span></a></li>
|
|
<li><a class="reference internal" href="#variable"><span>Variables that generate values to output</span></a></li>
|
|
<li><a class="reference internal" href="#table"><span>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<a class="headerlink" href="#global-per-atom-local-data" title="Permalink to this headline">¶</a></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<a class="headerlink" href="#scalar-vector-array-data" title="Permalink to this headline">¶</a></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<a class="headerlink" href="#thermodynamic-output" title="Permalink to this headline">¶</a></h3>
|
|
<p>The frequency and format of thermodynamic output is set by the
|
|
<a class="reference internal" href="thermo.html"><em>thermo</em></a>, <a class="reference internal" href="thermo_style.html"><em>thermo_style</em></a>, and
|
|
<a class="reference internal" href="thermo_modify.html"><em>thermo_modify</em></a> commands. The
|
|
<a class="reference internal" href="thermo_style.html"><em>thermo_style</em></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"><em>compute</em></a>
|
|
or <a class="reference internal" href="fix.html"><em>fix</em></a> or <a class="reference internal" href="variable.html"><em>variable</em></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"><em>thermo_style custom</em></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"><em>thermo_modify norm</em></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"><em>Equal-style variables</em></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<a class="headerlink" href="#dump-file-output" title="Permalink to this headline">¶</a></h3>
|
|
<p>Dump file output is specified by the <a class="reference internal" href="dump.html"><em>dump</em></a> and
|
|
<a class="reference internal" href="dump_modify.html"><em>dump_modify</em></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"><em>dump custom</em></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"><em>compute</em></a> or <a class="reference internal" href="fix.html"><em>fix</em></a> or <a class="reference internal" href="variable.html"><em>variable</em></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"><em>dump custom</em></a> command.</p>
|
|
<p>There is also a <a class="reference internal" href="dump.html"><em>dump local</em></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"><em>compute</em></a> or <a class="reference internal" href="fix.html"><em>fix</em></a> or <a class="reference internal" href="variable.html"><em>variable</em></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"><em>dump local</em></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<a class="headerlink" href="#fixes-that-write-output-files" title="Permalink to this headline">¶</a></h3>
|
|
<p>Several fixes take various quantities as input and can write output
|
|
files: <a class="reference internal" href="fix_ave_time.html"><em>fix ave/time</em></a>, <a class="reference internal" href="fix_ave_spatial.html"><em>fix ave/spatial</em></a>, <a class="reference internal" href="fix_ave_histo.html"><em>fix ave/histo</em></a>,
|
|
<a class="reference internal" href="fix_ave_correlate.html"><em>fix ave/correlate</em></a>, and <a class="reference internal" href="fix_print.html"><em>fix print</em></a>.</p>
|
|
<p>The <a class="reference internal" href="fix_ave_time.html"><em>fix ave/time</em></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"><em>compute</em></a> values, global <a class="reference internal" href="fix.html"><em>fix</em></a> values, or
|
|
<a class="reference internal" href="variable.html"><em>variables</em></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"><em>thermo_style custom</em></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_spatial.html"><em>fix ave/spatial</em></a> command enables direct
|
|
output to a file of spatial-averaged per-atom quantities like those
|
|
output in dump files, within 1d layers of the simulation box. The
|
|
per-atom quantities can be atom density (mass or number) or atom
|
|
attributes such as position, velocity, force. They can also be
|
|
per-atom quantities calculated by a <a class="reference internal" href="compute.html"><em>compute</em></a>, by a
|
|
<a class="reference internal" href="fix.html"><em>fix</em></a>, or by an atom-style <a class="reference internal" href="variable.html"><em>variable</em></a>. The
|
|
spatial-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"><em>fix ave/histo</em></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"><em>fix ave/correlate</em></a> command enables direct
|
|
output to a file of time-correlated quantities, which can be global
|
|
scalars. The correlation matrix output of this fix can also be used
|
|
as input to other output commands.</p>
|
|
<p>The <a class="reference internal" href="fix_print.html"><em>fix print</em></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"><em>variable</em></a> values for any style variable except the atom
|
|
style). As explained above, variables themselves can contain
|
|
references to global values generated by <a class="reference internal" href="thermo_style.html"><em>thermodynamic keywords</em></a>, <a class="reference internal" href="compute.html"><em>computes</em></a>,
|
|
<a class="reference internal" href="fix.html"><em>fixes</em></a>, or other <a class="reference internal" href="variable.html"><em>variables</em></a>, or to per-atom
|
|
values for a specific atom. Thus the <a class="reference internal" href="fix_print.html"><em>fix print</em></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<a class="headerlink" href="#computes-that-process-output-quantities" title="Permalink to this headline">¶</a></h3>
|
|
<p>The <a class="reference internal" href="compute_reduce.html"><em>compute reduce</em></a> and <a class="reference internal" href="compute_reduce.html"><em>compute reduce/region</em></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"><em>compute slice</em></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"><em>compute property/atom</em></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"><em>dump custom</em></a> command.</p>
|
|
<p>The <a class="reference internal" href="compute_property_local.html"><em>compute property/local</em></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="id3"></span><h3>6.15.7. Fixes that process output quantities<a class="headerlink" href="#fixes-that-process-output-quantities" title="Permalink to this headline">¶</a></h3>
|
|
<p>The <a class="reference internal" href="fix_vector.html"><em>fix vector</em></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"><em>fix ave/atom</em></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"><em>compute</em></a>, by a
|
|
<a class="reference internal" href="fix.html"><em>fix</em></a>, or by an atom-style <a class="reference internal" href="variable.html"><em>variable</em></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"><em>fix store/state</em></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"><em>dump custom</em></a> command,
|
|
including per-atom quantities calculated by a <a class="reference internal" href="compute.html"><em>compute</em></a>,
|
|
by a <a class="reference internal" href="fix.html"><em>fix</em></a>, or by an atom-style <a class="reference internal" href="variable.html"><em>variable</em></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<a class="headerlink" href="#computes-that-generate-values-to-output" title="Permalink to this headline">¶</a></h3>
|
|
<p>Every <a class="reference internal" href="compute.html"><em>compute</em></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<a class="headerlink" href="#fixes-that-generate-values-to-output" title="Permalink to this headline">¶</a></h3>
|
|
<p>Some <a class="reference internal" href="fix.html"><em>fixes</em></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<a class="headerlink" href="#variables-that-generate-values-to-output" title="Permalink to this headline">¶</a></h3>
|
|
<p>Every <a class="reference internal" href="variable.html"><em>variables</em></a> defined in an input script generates
|
|
either a global scalar value or a per-atom vector (only atom-style
|
|
variables) when it is accessed. The formulas used to define equal-
|
|
and atom-style variables can contain references to the thermodynamic
|
|
keywords and to global and per-atom data generated by computes, fixes,
|
|
and other variables. The values generated by variables can be output
|
|
using the other commands described in this section.</p>
|
|
</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<a class="headerlink" href="#summary-table-of-output-options-and-data-flow-between-commands" title="Permalink to this headline">¶</a></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="30%" />
|
|
<col width="30%" />
|
|
<col width="1%" />
|
|
</colgroup>
|
|
<tbody valign="top">
|
|
<tr class="row-odd"><td>Command</td>
|
|
<td>Input</td>
|
|
<td>Output</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-even"><td><a class="reference internal" href="thermo_style.html"><em>thermo_style custom</em></a></td>
|
|
<td>global scalars</td>
|
|
<td>screen, log file</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-odd"><td><a class="reference internal" href="dump.html"><em>dump custom</em></a></td>
|
|
<td>per-atom vectors</td>
|
|
<td>dump file</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-even"><td><a class="reference internal" href="dump.html"><em>dump local</em></a></td>
|
|
<td>local vectors</td>
|
|
<td>dump file</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-odd"><td><a class="reference internal" href="fix_print.html"><em>fix print</em></a></td>
|
|
<td>global scalar from variable</td>
|
|
<td>screen, file</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-even"><td><a class="reference internal" href="print.html"><em>print</em></a></td>
|
|
<td>global scalar from variable</td>
|
|
<td>screen</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-odd"><td><a class="reference internal" href="compute.html"><em>computes</em></a></td>
|
|
<td>N/A</td>
|
|
<td>global/per-atom/local scalar/vector/array</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-even"><td><a class="reference internal" href="fix.html"><em>fixes</em></a></td>
|
|
<td>N/A</td>
|
|
<td>global/per-atom/local scalar/vector/array</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-odd"><td><a class="reference internal" href="variable.html"><em>variables</em></a></td>
|
|
<td>global scalars, per-atom vectors</td>
|
|
<td>global scalar, per-atom vector</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-even"><td><a class="reference internal" href="compute_reduce.html"><em>compute reduce</em></a></td>
|
|
<td>per-atom/local vectors</td>
|
|
<td>global scalar/vector</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-odd"><td><a class="reference internal" href="compute_slice.html"><em>compute slice</em></a></td>
|
|
<td>global vectors/arrays</td>
|
|
<td>global vector/array</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-even"><td><a class="reference internal" href="compute_property_atom.html"><em>compute property/atom</em></a></td>
|
|
<td>per-atom vectors</td>
|
|
<td>per-atom vector/array</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-odd"><td><a class="reference internal" href="compute_property_local.html"><em>compute property/local</em></a></td>
|
|
<td>local vectors</td>
|
|
<td>local vector/array</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-even"><td><a class="reference internal" href="fix_vector.html"><em>fix vector</em></a></td>
|
|
<td>global scalars</td>
|
|
<td>global vector</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-odd"><td><a class="reference internal" href="fix_ave_atom.html"><em>fix ave/atom</em></a></td>
|
|
<td>per-atom vectors</td>
|
|
<td>per-atom vector/array</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-even"><td><a class="reference internal" href="fix_ave_time.html"><em>fix ave/time</em></a></td>
|
|
<td>global scalars/vectors</td>
|
|
<td>global scalar/vector/array, file</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-odd"><td><a class="reference internal" href="fix_ave_spatial.html"><em>fix ave/spatial</em></a></td>
|
|
<td>per-atom vectors</td>
|
|
<td>global array, file</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-even"><td><a class="reference internal" href="fix_ave_histo.html"><em>fix ave/histo</em></a></td>
|
|
<td>global/per-atom/local scalars and vectors</td>
|
|
<td>global array, file</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-odd"><td><a class="reference internal" href="fix_ave_correlate.html"><em>fix ave/correlate</em></a></td>
|
|
<td>global scalars</td>
|
|
<td>global array, file</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-even"><td><a class="reference internal" href="fix_store_state.html"><em>fix store/state</em></a></td>
|
|
<td>per-atom vectors</td>
|
|
<td>per-atom vector/array</td>
|
|
<td> </td>
|
|
</tr>
|
|
<tr class="row-odd"><td> </td>
|
|
<td> </td>
|
|
<td> </td>
|
|
<td> </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<a class="headerlink" href="#thermostatting-barostatting-and-computing-temperature" title="Permalink to this headline">¶</a></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"><em>compute commands</em></a> calculate temperature, and the <a class="reference internal" href="compute_pressure.html"><em>compute pressure</em></a> command calculates pressure.</p>
|
|
<ul class="simple">
|
|
<li><a class="reference internal" href="compute_temp.html"><em>compute temp</em></a></li>
|
|
<li><a class="reference internal" href="compute_temp_sphere.html"><em>compute temp/sphere</em></a></li>
|
|
<li><a class="reference internal" href="compute_temp_asphere.html"><em>compute temp/asphere</em></a></li>
|
|
<li><a class="reference internal" href="compute_temp_com.html"><em>compute temp/com</em></a></li>
|
|
<li><a class="reference internal" href="compute_temp_deform.html"><em>compute temp/deform</em></a></li>
|
|
<li><a class="reference internal" href="compute_temp_partial.html"><em>compute temp/partial</em></a></li>
|
|
<li><a class="reference internal" href="compute_temp_profile.html"><em>compute temp/profile</em></a></li>
|
|
<li><a class="reference internal" href="compute_temp_ramp.html"><em>compute temp/ramp</em></a></li>
|
|
<li><a class="reference internal" href="compute_temp_region.html"><em>compute temp/region</em></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"><em>Compute temp/sphere</em></a> and <a class="reference internal" href="compute_temp_asphere.html"><em>compute temp/asphere</em></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"><em>fixes</em></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"><em>fix nvt</em></a></li>
|
|
<li><a class="reference internal" href="fix_nvt_sphere.html"><em>fix nvt/sphere</em></a></li>
|
|
<li><a class="reference internal" href="fix_nvt_asphere.html"><em>fix nvt/asphere</em></a></li>
|
|
<li><a class="reference internal" href="fix_nvt_sllod.html"><em>fix nvt/sllod</em></a></li>
|
|
<li><a class="reference internal" href="fix_temp_berendsen.html"><em>fix temp/berendsen</em></a></li>
|
|
<li><a class="reference internal" href="fix_temp_csvr.html"><em>fix temp/csvr</em></a></li>
|
|
<li><a class="reference internal" href="fix_langevin.html"><em>fix langevin</em></a></li>
|
|
<li><a class="reference internal" href="fix_temp_rescale.html"><em>fix temp/rescale</em></a></li>
|
|
<li><a class="reference internal" href="pair_dpd.html"><em>pair_style dpd/tstat</em></a></li>
|
|
</ul>
|
|
<p><a class="reference internal" href="fix_nh.html"><em>Fix nvt</em></a> only thermostats the translational velocity of
|
|
particles. <a class="reference internal" href="fix_nvt_sllod.html"><em>Fix nvt/sllod</em></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>NEMD simulations</span></a> section of this page for further details. <a class="reference internal" href="fix_nvt_sphere.html"><em>Fix nvt/sphere</em></a> and <a class="reference internal" href="fix_nvt_asphere.html"><em>fix nvt/asphere</em></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"><em>fix langevin</em></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"><em>fix_modify</em></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"><em>compute temp/partial</em></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"><em>compute temp/profile</em></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"><em>fix nve</em></a></li>
|
|
<li><a class="reference internal" href="fix_nve_sphere.html"><em>fix nve/sphere</em></a></li>
|
|
<li><a class="reference internal" href="fix_nve_asphere.html"><em>fix nve/asphere</em></a></li>
|
|
</ul>
|
|
<p>Barostatting in LAMMPS is also performed by <a class="reference internal" href="fix.html"><em>fixes</em></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"><em>fix npt</em></a></li>
|
|
<li><a class="reference internal" href="fix_npt_sphere.html"><em>fix npt/sphere</em></a></li>
|
|
<li><a class="reference internal" href="fix_npt_asphere.html"><em>fix npt/asphere</em></a></li>
|
|
<li><a class="reference internal" href="fix_nh.html"><em>fix nph</em></a></li>
|
|
<li><a class="reference internal" href="fix_press_berendsen.html"><em>fix press/berendsen</em></a></li>
|
|
</ul>
|
|
<p>The <a class="reference internal" href="fix_nh.html"><em>fix npt</em></a> commands include a Nose-Hoover thermostat
|
|
and barostat. <a class="reference internal" href="fix_nh.html"><em>Fix nph</em></a> is just a Nose/Hoover barostat;
|
|
it does no thermostatting. Both <a class="reference internal" href="fix_nh.html"><em>fix nph</em></a> and <a class="reference internal" href="fix_press_berendsen.html"><em>fix press/bernendsen</em></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"><em>fix npt</em></a> and <a class="reference internal" href="fix_nh.html"><em>fix nph</em></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"><em>Fix npt/sphere</em></a> and <a class="reference internal" href="fix_npt_asphere.html"><em>fix npt/asphere</em></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"><em>compute pressure</em></a> compute to calculate a current
|
|
pressure. By default, this compute is created with a simple <a class="reference internal" href="compute_temp.html"><em>compute temp</em></a> (see the last argument of the <a class="reference internal" href="compute_pressure.html"><em>compute pressure</em></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"><em>fix_modify</em></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"><em>fix npt</em></a> and <a class="reference internal" href="fix_nh.html"><em>fix nph</em></a>) perform time integration.
|
|
<a class="reference internal" href="fix_press_berendsen.html"><em>Fix press/berendsen</em></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"><em>thermo_style</em></a> command, often includes temperature
|
|
and pressure values. As explained on the doc page for the
|
|
<a class="reference internal" href="thermo_style.html"><em>thermo_style</em></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"><em>thermo_style custom</em></a> command. Or
|
|
you can use the <a class="reference internal" href="thermo_modify.html"><em>thermo_modify</em></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<a class="headerlink" href="#walls" title="Permalink to this headline">¶</a></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"><em>lattice</em></a> and <a class="reference internal" href="create_atoms.html"><em>create_atoms</em></a> commands,
|
|
or read in via the <a class="reference internal" href="read_data.html"><em>read_data</em></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"><em>fix nve</em></a> or <a class="reference internal" href="fix_nh.html"><em>fix nvt</em></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"><em>fix aveforce</em></a> - set force on particles to average value, so they move together</li>
|
|
<li><a class="reference internal" href="fix_setforce.html"><em>fix setforce</em></a> - set force on particles to a value, e.g. 0.0</li>
|
|
<li><a class="reference internal" href="fix_freeze.html"><em>fix freeze</em></a> - freeze particles for use as granular walls</li>
|
|
<li><a class="reference internal" href="fix_nve_noforce.html"><em>fix nve/noforce</em></a> - advect particles by their velocity, but without force</li>
|
|
<li><a class="reference internal" href="fix_move.html"><em>fix move</em></a> - prescribe motion of particles by a linear velocity, oscillation, rotation, variable</li>
|
|
</ul>
|
|
<p>The <a class="reference internal" href="fix_move.html"><em>fix move</em></a> command offers the most generality, since
|
|
the motion of individual particles can be specified with
|
|
<a class="reference internal" href="variable.html"><em>variable</em></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"><em>neigh_modify exclude</em></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"><em>bond</em></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"><em>Fix wall/gran</em></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"><em>fix wall/reflect</em></a> - reflective flat walls</li>
|
|
<li><a class="reference internal" href="fix_wall.html"><em>fix wall/lj93</em></a> - flat walls, with Lennard-Jones 9/3 potential</li>
|
|
<li><a class="reference internal" href="fix_wall.html"><em>fix wall/lj126</em></a> - flat walls, with Lennard-Jones 12/6 potential</li>
|
|
<li><a class="reference internal" href="fix_wall.html"><em>fix wall/colloid</em></a> - flat walls, with <a class="reference internal" href="pair_colloid.html"><em>pair_style colloid</em></a> potential</li>
|
|
<li><a class="reference internal" href="fix_wall.html"><em>fix wall/harmonic</em></a> - flat walls, with repulsive harmonic spring potential</li>
|
|
<li><a class="reference internal" href="fix_wall_region.html"><em>fix wall/region</em></a> - use region surface as wall</li>
|
|
<li><a class="reference internal" href="fix_wall_gran.html"><em>fix wall/gran</em></a> - flat or curved walls with <a class="reference internal" href="pair_gran.html"><em>pair_style granular</em></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"><em>fix wall/region</em></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"><em>Regions</em></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"><em>Regions</em></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"><em>fix wall/gran</em></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<a class="headerlink" href="#elastic-constants" title="Permalink to this headline">¶</a></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"><em>change_box</em></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"><em>this section</em></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>(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<a class="headerlink" href="#library-interface-to-lammps" title="Permalink to this headline">¶</a></h2>
|
|
<p>As described in <a class="reference internal" href="Section_start.html#start-5"><span>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>coupled manner</span></a> with other
|
|
codes, or driven through a <a class="reference internal" href="Section_python.html"><em>Python interface</em></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-python"><div class="highlight"><pre>void lammps_open(int, char **, MPI_Comm, void **)
|
|
void lammps_close(void *)
|
|
int lammps_version(void *)
|
|
void lammps_file(void *, char *)
|
|
char *lammps_command(void *, char *)
|
|
</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>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-python"><div class="highlight"><pre>void *lammps_extract_global(void *, char *)
|
|
void *lammps_extract_atom(void *, char *)
|
|
void *lammps_extract_compute(void *, char *, int, int)
|
|
void *lammps_extract_fix(void *, char *, int, int, int, int)
|
|
void *lammps_extract_variable(void *, char *, char *)
|
|
int lammps_set_variable(void *, char *, char *)
|
|
int lammps_get_natoms(void *)
|
|
void lammps_get_coords(void *, double *)
|
|
void lammps_put_coords(void *, double *)
|
|
</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"><em>Python interface</em></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<a class="headerlink" href="#calculating-thermal-conductivity" title="Permalink to this headline">¶</a></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>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>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>Ikeshoji and Hafskjold</span></a> for
|
|
details of this idea. Note that thermostatting fixes such as <a class="reference internal" href="fix_nh.html"><em>fix nvt</em></a>, <a class="reference internal" href="fix_langevin.html"><em>fix langevin</em></a>, and <a class="reference internal" href="fix_temp_rescale.html"><em>fix temp/rescale</em></a> store the cumulative energy they
|
|
add/subtract.</p>
|
|
<p>Alternatively, as a second method, the <a class="reference internal" href="fix_heat.html"><em>fix heat</em></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"><em>fix ave/spatial</em></a> and <a class="reference internal" href="compute_ke_atom.html"><em>compute ke/atom</em></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"><em>fix thermal/conductivity</em></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"><em>fix ave/spatial</em></a> and <a class="reference internal" href="compute_ke_atom.html"><em>compute ke/atom</em></a> commands. The fix tallies the
|
|
cumulative energy transfer that it performs. See the <a class="reference internal" href="fix_thermal_conductivity.html"><em>fix thermal/conductivity</em></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"><em>compute heat/flux</em></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"><em>fix ave/correlate</em></a> command to calculate the needed
|
|
auto-correlation. See the doc page for the <a class="reference internal" href="compute_heat_flux.html"><em>compute heat/flux</em></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<a class="headerlink" href="#calculating-viscosity" title="Permalink to this headline">¶</a></h2>
|
|
<p>The shear viscosity eta of a fluid can be measured in at least 4 ways
|
|
using various options in LAMMPS. See the examples/VISCOSITY 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-20"><span>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"><em>fix deform</em></a>
|
|
command, and using the <a class="reference internal" href="fix_nvt_sllod.html"><em>fix nvt/sllod</em></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"><em>fix ave/spatial</em></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"><em>compute pressure</em></a>
|
|
command, can also be monitored, which is the J term in the equation
|
|
above. See <a class="reference internal" href="#howto-13"><span>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"><em>fix viscosity</em></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"><em>fix ave/spatial</em></a> command.
|
|
The fix tallies the cummulative momentum transfer that it performs.
|
|
See the <a class="reference internal" href="fix_viscosity.html"><em>fix viscosity</em></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 steady-state
|
|
equilibrated simulation which is in contrast to the two preceding
|
|
non-equilibrium methods, where momentum flows continuously through the
|
|
simulation box.</p>
|
|
<p>Here is an example input script that calculates the viscosity of
|
|
liquid Ar via the GK formalism:</p>
|
|
<div class="highlight-python"><div class="highlight"><pre><span class="c"># Sample LAMMPS input script for viscosity of liquid Ar</span>
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>units real
|
|
variable T equal 86.4956
|
|
variable V equal vol
|
|
variable dt equal 4.0
|
|
variable p equal 400 # correlation length
|
|
variable s equal 5 # sample interval
|
|
variable d equal $p*$s # dump interval
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre><span class="c"># convert from LAMMPS real units to SI</span>
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>variable kB equal 1.3806504e-23 # [J/K/** Boltzmann
|
|
variable atm2Pa equal 101325.0
|
|
variable A2m equal 1.0e-10
|
|
variable fs2s equal 1.0e-15
|
|
variable convert equal ${atm2Pa}*${atm2Pa}*${fs2s}*${A2m}*${A2m}*${A2m}
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre><span class="c"># setup problem</span>
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>dimension 3
|
|
boundary p p p
|
|
lattice fcc 5.376 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
|
|
region box block 0 4 0 4 0 4
|
|
create_box 1 box
|
|
create_atoms 1 box
|
|
mass 1 39.948
|
|
pair_style lj/cut 13.0
|
|
pair_coeff * * 0.2381 3.405
|
|
timestep ${dt}
|
|
thermo $d
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre><span class="c"># equilibration and thermalization</span>
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>velocity all create $T 102486 mom yes rot yes dist gaussian
|
|
fix NVT all nvt temp $T $T 10 drag 0.2
|
|
run 8000
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre><span class="c"># viscosity calculation, switch to NVE if desired</span>
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre><span class="c">#unfix NVT</span>
|
|
<span class="c">#fix NVE all nve</span>
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>reset_timestep 0
|
|
variable pxy equal pxy
|
|
variable pxz equal pxz
|
|
variable pyz equal pyz
|
|
fix SS all ave/correlate $s $p $d &
|
|
v_pxy v_pxz v_pyz type auto file S0St.dat ave running
|
|
variable scale equal ${convert}/(${kB}*$T)*$V*$s*${dt}
|
|
variable v11 equal trap(f_SS[3])*${scale}
|
|
variable v22 equal trap(f_SS[4])*${scale}
|
|
variable v33 equal trap(f_SS[5])*${scale}
|
|
thermo_style custom step temp press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33
|
|
run 100000
|
|
variable v equal (v_v11+v_v22+v_v33)/3.0
|
|
variable ndens equal count(all)/vol
|
|
print "average viscosity: $v [Pa.s/** @ $T K, ${ndens} /A^3"
|
|
</pre></div>
|
|
</div>
|
|
<hr class="docutils" />
|
|
</div>
|
|
<div class="section" id="calculating-a-diffusion-coefficient">
|
|
<span id="howto-22"></span><h2>6.22. Calculating a diffusion coefficient<a class="headerlink" href="#calculating-a-diffusion-coefficient" title="Permalink to this headline">¶</a></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"><em>compute msd</em></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"><em>fix vector</em></a> command, and a line fit to the vector to
|
|
compute its slope via the <a class="reference internal" href="variable.html"><em>variable slope</em></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"><em>compute vacf</em></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"><em>fix vector</em></a> command,
|
|
and time integrated via the <a class="reference internal" href="variable.html"><em>variable trap</em></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<a class="headerlink" href="#using-chunks-to-calculate-system-properties" title="Permalink to this headline">¶</a></h2>
|
|
<p>In LAMMS, “chunks” are collections of atoms, as defined by the
|
|
<a class="reference internal" href="compute_chunk_atom.html"><em>compute chunk/atom</em></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"><em>fix ave/chunk</em></a></li>
|
|
<li>any of the <a class="reference internal" href="compute.html"><em>compute */chunk</em></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:<a class="headerlink" href="#compute-chunk-atom-command" title="Permalink to this headline">¶</a></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"><em>compute cluster/atom</em></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"><em>compute pe/atom</em></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"><em>compute centro/atom</em></a> or <a class="reference internal" href="compute_coord_atom.html"><em>compute coord/atom</em></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"><em>computes</em></a>, <a class="reference internal" href="fix.html"><em>fixes</em></a>, and
|
|
<a class="reference internal" href="variable.html"><em>variables</em></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"><em>Section_modify</em></a> of the documentation for how to
|
|
do this. You can also define a <a class="reference internal" href="variable.html"><em>per-atom variable</em></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:<a class="headerlink" href="#fix-ave-chunk-command" title="Permalink to this headline">¶</a></h3>
|
|
<p>This fix takes the ID of a <a class="reference internal" href="compute_chunk_atom.html"><em>compute chunk/atom</em></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"><em>computes</em></a>, <a class="reference internal" href="fix.html"><em>fixes</em></a>, and <a class="reference internal" href="variable.html"><em>per-atom variables</em></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="#id70"><span class="problematic" id="id71">*</span></a>/chunk commands:<a class="headerlink" href="#compute-chunk-commands" title="Permalink to this headline">¶</a></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"><em>compute com/chunk</em></a></li>
|
|
<li><a class="reference internal" href="compute_gyration_chunk.html"><em>compute gyration/chunk</em></a></li>
|
|
<li><a class="reference internal" href="compute_inertia_chunk.html"><em>compute inertia/chunk</em></a></li>
|
|
<li><a class="reference internal" href="compute_msd_chunk.html"><em>compute msd/chunk</em></a></li>
|
|
<li><a class="reference internal" href="compute_property_chunk.html"><em>compute property/chunk</em></a></li>
|
|
<li><a class="reference internal" href="compute_temp_chunk.html"><em>compute temp/chunk</em></a></li>
|
|
<li><a class="reference internal" href="compute_vcm_chunk.html"><em>compute torque/chunk</em></a></li>
|
|
<li><a class="reference internal" href="compute_vcm_chunk.html"><em>compute vcm/chunk</em></a></li>
|
|
</ul>
|
|
<p>They each take the ID of a <a class="reference internal" href="compute_chunk_atom.html"><em>compute chunk/atom</em></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"><em>compute property/chunk</em></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"><em>fix ave/chunk command</em></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"><em>fix ave/time</em></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"><em>fix ave/histo</em></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"><em>equal-style variables</em></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<a class="headerlink" href="#example-calculations-with-chunks" title="Permalink to this headline">¶</a></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-python"><div class="highlight"><pre>compute cc1 all chunk/atom bin/2d x 0.0 0.1 y lower 0.01 units reduced
|
|
fix 1 all ave/chunk 100 10 1000 cc1 vx vy file tmp.out
|
|
</pre></div>
|
|
</div>
|
|
<p>(2) Temperature in each spatial bin, after subtracting a flow
|
|
velocity:</p>
|
|
<div class="highlight-python"><div class="highlight"><pre>compute cc1 all chunk/atom bin/2d x 0.0 0.1 y lower 0.1 units reduced
|
|
compute vbias all temp/profile 1 0 0 y 10
|
|
fix 1 all ave/chunk 100 10 1000 cc1 temp bias vbias file tmp.out
|
|
</pre></div>
|
|
</div>
|
|
<ol class="arabic simple" start="3">
|
|
<li>Center of mass of each molecule:</li>
|
|
</ol>
|
|
<div class="highlight-python"><div class="highlight"><pre>compute cc1 all chunk/atom molecule
|
|
compute myChunk all com/chunk cc1
|
|
fix 1 all ave/time 100 1 100 c_myChunk file tmp.out mode vector
|
|
</pre></div>
|
|
</div>
|
|
<ol class="arabic simple" start="4">
|
|
<li>Total force on each molecule and ave/max across all molecules:</li>
|
|
</ol>
|
|
<div class="highlight-python"><div class="highlight"><pre>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**2**)
|
|
variable xmax equal max(f_1**2**)
|
|
thermo 1000
|
|
thermo_style custom step temp v_xave v_xmax
|
|
</pre></div>
|
|
</div>
|
|
<ol class="arabic simple" start="5">
|
|
<li>Histogram of cluster sizes:</li>
|
|
</ol>
|
|
<div class="highlight-python"><div class="highlight"><pre>compute cluster all cluster/atom 1.0
|
|
compute cc1 all chunk/atom c_cluster compress yes
|
|
compute size all property/chunk cc1 count
|
|
fix 1 all ave/histo 100 1 100 0 20 20 c_size mode vector ave running beyond ignore file tmp.histo
|
|
</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"><em>kspace_style pppm/disp</em></a> command<a class="headerlink" href="#setting-parameters-for-the-kspace-style-pppm-disp-command" title="Permalink to this headline">¶</a></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<a class="headerlink" href="#polarizable-models" title="Permalink to this headline">¶</a></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"><em>QEQ</em></a>
|
|
package,</li>
|
|
<li>the adiabatic core-shell method, implemented in the
|
|
<a class="reference internal" href="#howto-26"><span>CORESHELL</span></a> package,</li>
|
|
<li>the thermalized Drude dipole method, implemented in the
|
|
<a class="reference internal" href="#howto-27"><span>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"><em>fix qeq</em></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"><em>fix qeq/point</em></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"><em>fix qeq/dynamic</em></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<a class="headerlink" href="#adiabatic-core-shell-model" title="Permalink to this headline">¶</a></h2>
|
|
<p>The adiabatic core-shell model by <a class="reference internal" href="pair_cs.html#mitchellfinchham"><span>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"><em>atom_style</em></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-python"><div class="highlight"><pre>432 atoms # core and shell atoms
|
|
216 bonds # number of core/shell springs
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>4 atom types # 2 cores and 2 shells for Na and Cl
|
|
2 bond types
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>0.0 24.09597 xlo xhi
|
|
0.0 24.09597 ylo yhi
|
|
0.0 24.09597 zlo zhi
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre><span class="n">Masses</span> <span class="c"># core/shell mass ratio = 0.1</span>
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>1 20.690784 # Na core
|
|
2 31.90500 # Cl core
|
|
3 2.298976 # Na shell
|
|
4 3.54500 # Cl shell
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre><span class="n">Atoms</span>
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>1 1 2 1.5005 0.00000000 0.00000000 0.00000000 # core of core/shell pair 1
|
|
2 1 4 -2.5005 0.00000000 0.00000000 0.00000000 # shell of core/shell pair 1
|
|
3 2 1 1.5056 4.01599500 4.01599500 4.01599500 # core of core/shell pair 2
|
|
4 2 3 -0.5056 4.01599500 4.01599500 4.01599500 # shell of core/shell pair 2
|
|
(...)
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre><span class="n">Bonds</span> <span class="c"># Bond topology for spring forces</span>
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>1 2 1 2 # spring for core/shell pair 1
|
|
2 2 3 4 # spring for core/shell pair 2
|
|
(...)
|
|
</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"><em>special_bonds</em></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"><em>special_bonds</em></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"><em>special_bonds</em></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"><em>this doc page</em></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"><em>kspace_style</em></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-python"><div class="highlight"><pre>pair_style born/coul/long/cs 20.0 20.0
|
|
pair_coeff * * 0.0 1.000 0.00 0.00 0.00
|
|
pair_coeff 3 3 487.0 0.23768 0.00 1.05 0.50 #Na-Na
|
|
pair_coeff 3 4 145134.0 0.23768 0.00 6.99 8.70 #Na-Cl
|
|
pair_coeff 4 4 405774.0 0.23768 0.00 72.40 145.40 #Cl-Cl
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>bond_style harmonic
|
|
bond_coeff 1 63.014 0.0
|
|
bond_coeff 2 25.724 0.0
|
|
</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"><em>compute temp/cs</em></a> command can be used, in conjunction with
|
|
any of the thermostat fixes, such as <a class="reference internal" href="fix_nh.html"><em>fix nvt</em></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"><em>compute temp/cs</em></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"><em>group *type*</em></a> command.
|
|
Note that to perform thermostatting using this definition of
|
|
temperature, the <a class="reference internal" href="fix_modify.html"><em>fix modify temp</em></a> command should be
|
|
used to assign the compute to the thermostat fix. Likewise the
|
|
<a class="reference internal" href="thermo_modify.html"><em>thermo_modify temp</em></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-python"><div class="highlight"><pre>group cores type 1 2
|
|
group shells type 3 4
|
|
compute CSequ all temp/cs cores shells
|
|
fix thermoberendsen all temp/berendsen 1427 1427 0.4 # thermostat for the true physical system
|
|
fix thermostatequ all nve # integrator as needed for the berendsen thermostat
|
|
fix_modify thermoberendsen temp CSequ
|
|
thermo_modify temp CSequ # output of center-of-mass derived temperature
|
|
</pre></div>
|
|
</div>
|
|
<p>If <a class="reference internal" href="compute_temp_cs.html"><em>compute temp/cs</em></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"><em>fix momentum</em></a> command in combination with <a class="reference internal" href="compute_temp_cs.html"><em>compute temp/cs</em></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"><em>velocity create</em></a> command and assigning the <a class="reference internal" href="compute_temp_cs.html"><em>compute temp/cs</em></a> command to the <em>temp</em> keyword of the
|
|
<a class="reference internal" href="velocity.html"><em>velocity</em></a> commmand, e.g.</p>
|
|
<div class="highlight-python"><div class="highlight"><pre>velocity all create 1427 134 bias yes temp CSequ
|
|
velocity all scale 1427 temp CSequ
|
|
</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"><em>compute temp/cs</em></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"><em>compute chunk/atom</em></a> and <a class="reference internal" href="compute_temp_chunk.html"><em>compute temp/chunk</em></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"><em>variable</em></a> command. Or they can
|
|
be time/averaged and output using the <a class="reference internal" href="fix_ave_time.html"><em>fix ave/time</em></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"><em>fix property/atom</em></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"><em>read_data</em></a> command. This field can then be
|
|
accessed by the <a class="reference internal" href="compute_property_atom.html"><em>compute property/atom</em></a>
|
|
command, to use as input to the <a class="reference internal" href="compute_chunk_atom.html"><em>compute chunk/atom</em></a> command to define the core/shell
|
|
pairs as chunks.</p>
|
|
<p>For example,</p>
|
|
<div class="highlight-python"><div class="highlight"><pre>fix csinfo all property/atom i_CSID # property/atom command
|
|
read_data NaCl_CS_x0.1_prop.data fix csinfo NULL CS-Info # atom property added in the data-file
|
|
compute prop all property/atom i_CSID
|
|
compute cs_chunk all chunk/atom c_prop
|
|
compute cstherm all temp/chunk cs_chunk temp internal com yes cdof 3.0 # note the chosen degrees of freedom for the core/shell pairs
|
|
fix ave_chunk all ave/time 10 1 10 c_cstherm file chunk.dump mode vector
|
|
</pre></div>
|
|
</div>
|
|
<p>The additional section in the date file would be formatted like this:</p>
|
|
<div class="highlight-python"><div class="highlight"><pre><span class="n">CS</span><span class="o">-</span><span class="n">Info</span> <span class="c"># header of additional section</span>
|
|
</pre></div>
|
|
</div>
|
|
<div class="highlight-python"><div class="highlight"><pre>1 1 # column 1 = atom ID, column 2 = core/shell ID
|
|
2 1
|
|
3 2
|
|
4 2
|
|
5 3
|
|
6 3
|
|
7 4
|
|
8 4
|
|
(...)
|
|
</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<a class="headerlink" href="#drude-induced-dipoles" title="Permalink to this headline">¶</a></h2>
|
|
<p>The thermalized Drude model, similarly to the <a class="reference internal" href="#howto-26"><span>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="tutorial_drude.html#lamoureux"><span>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"><em>available here</em></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"><em>fix drude</em></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"><em>special_bonds</em></a> command. If using <a class="reference internal" href="fix_shake.html"><em>fix shake</em></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"><em>fix langevin/drude</em></a>
|
|
for a Langevin thermostat, or <a class="reference internal" href="fix_drude_transform.html"><em>fix drude/transform/*</em></a> for a Nose-Hoover
|
|
thermostat. The former requires use of the command <a class="reference internal" href="comm_modify.html"><em>comm_modify vel yes</em></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"><em>compute temp/drude</em></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"><em>pair style thole</em></a> in <a class="reference internal" href="pair_hybrid.html"><em>pair_style hybrid/overlay</em></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> |