<p>Now for the thermalization, the simplest choice is to use the <aclass="reference internal"href="fix_langevin_drude.html"><em>fix langevin/drude</em></a>.</p>
<divclass="highlight-python"><divclass="highlight"><pre>fix LANG all langevin/drude 300. 100 12435 1. 20 13977
</pre></div>
</div>
<p>This applies a Langevin thermostat at temperature 300. to the centers
of mass of the DC-DP pairs, with relaxation time 100 and with random
seed 12345. This fix applies also a Langevin thermostat at temperature
1. to the relative motion of the DPs around their DCs, with relaxation
time 20 and random seed 13977. Only the DCs and non-polarizable
atoms need to be in this fix’s group. LAMMPS will thermostate the DPs
together with their DC. For this, ghost atoms need to know their
velocities. Thus you need to add the following command:</p>
<divclass="highlight-python"><divclass="highlight"><pre>comm_modify vel yes
</pre></div>
</div>
<p>In order to avoid that the center of mass of the whole system
drifts due to the random forces of the Langevin thermostat on DCs, you
can add the <em>zero yes</em> option at the end of the fix line.</p>
<p>If the fix <em>shake</em> is used to constrain the C-H bonds, it should be
invoked after the fix <em>langevin/drude</em> for more accuracy.</p>
<divclass="highlight-python"><divclass="highlight"><pre>fix SHAKE ATOMS shake 0.0001 20 0 t 4 5
</pre></div>
</div>
<divclass="admonition note">
<pclass="first admonition-title">Note</p>
<pclass="last">The group of the fix <em>shake</em> must not include the DPs. If the
group <em>ATOMS</em> is defined by non-DPs atom types, you could use</p>
</div>
<p>Since the fix <em>langevin/drude</em> does not perform time integration (just
modification of forces but no position/velocity updates), the fix
<em>nve</em> should be used in conjunction.</p>
<divclass="highlight-python"><divclass="highlight"><pre>fix NVE all nve
</pre></div>
</div>
<p>Finally, do not forget to update the atom type elements if you use
them in a <em>dump_modify ... element ...</em> command, by adding the element
type of the DPs. Here for instance</p>
<divclass="highlight-python"><divclass="highlight"><pre>dump DUMP all custom 10 dump.lammpstrj id mol type element x y z ix iy iz
dump_modify DUMP element C C O H H D D D
</pre></div>
</div>
<p>The input file should now be ready for use!</p>
<p>You will notice that the global temperature <em>thermo_temp</em> computed by
LAMMPS is not 300. K as wanted. This is because LAMMPS treats DPs as
standard atoms in his default compute. If you want to output the
temperatures of the DC-DP pair centers of mass and of the DPs relative
to their DCs, you should use the <aclass="reference internal"href="compute_temp_drude.html"><em>compute temp_drude</em></a></p>
<divclass="highlight-python"><divclass="highlight"><pre>compute TDRUDE all temp/drude
</pre></div>
</div>
<p>And then output the correct temperatures of the Drude oscillators
using <em>thermo_style custom</em> with respectively <em>c_TDRUDE[1]</em> and
<em>c_TDRUDE[2]</em>. These should be close to 300.0 and 1.0 on average.</p>
<p>For the <em>thole</em> pair style the coefficients are</p>
<olclass="arabic simple">
<li>the atom polarizability in units of cubic length</li>
<li>the screening factor of the Thole function (optional, default value
specified by the pair_style command)</li>
</ol>
<ulclass="simple">
<li>the cutoff (optional, default value defined by the pair_style command)</li>
</ul>
<p>The special neighbors have charge-charge and charge-dipole
interactions screened by the <em>coul</em> factors of the <em>special_bonds</em>
command (0.0, 0.0, and 0.5 in the example above). Without using the
pair_style <em>thole</em>, dipole-dipole interactions are screened by the
same factor. By using the pair_style <em>thole</em>, dipole-dipole
interactions are screened by Thole’s function, whatever their special
relationship (except within each DC-DP pair of course). Consider for
example 1-2 neighbors: using the pair_style <em>thole</em>, their dipoles
will see each other (despite the <em>coul</em> factor being 0.) and the
interactions between these dipoles will be damped by Thole’s function.</p>
<hrclass="docutils"/>
<p><strong>Thermostats and barostats</strong></p>
<p>Using a Nose-Hoover barostat with the <em>langevin/drude</em> thermostat is
straightforward using fix <em>nph</em> instead of <em>nve</em>. For example:</p>
<divclass="highlight-python"><divclass="highlight"><pre>fix NPH all nph iso 1. 1. 500
</pre></div>
</div>
<p>It is also possible to use a Nose-Hoover instead of a Langevin
thermostat. This requires to use <aclass="reference internal"href="fix_drude_transform.html"><em>*fix drude/transform*</em></a> just before and after the
time intergation fixes. The <em>fix drude/transform/direct</em> converts the
atomic masses, positions, velocities and forces into a reduced
representation, where the DCs transform into the centers of mass of
the DC-DP pairs and the DPs transform into their relative position
with respect to their DC. The <em>fix drude/transform/inverse</em> performs
the reverse transformation. For a NVT simulation, with the DCs and
atoms at 300 K and the DPs at 1 K relative to their DC one would use</p>
<divclass="highlight-python"><divclass="highlight"><pre>fix DIRECT all drude/transform/direct
fix NVT1 ATOMS nvt temp 300. 300. 100
fix NVT2 DRUDES nvt temp 1. 1. 20
fix INVERSE all drude/transform/inverse
</pre></div>
</div>
<p>For our phenol example, the groups would be defined as</p>
<divclass="highlight-python"><divclass="highlight"><pre>group ATOMS type 1 2 3 4 5 # DCs and non-polarizable atoms
group CORES type 1 2 3 # DCs
group DRUDES type 6 7 8 # DPs
</pre></div>
</div>
<p>Note that with the fixes <em>drude/transform</em>, it is not required to
specify <em>comm_modify vel yes</em> because the fixes do it anyway (several
times and for the forces also). To avoid the flying ice cube artifact
<aclass="reference internal"href="#lamoureux"><span>(Lamoureux)</span></a>, where the atoms progressively freeze and the
center of mass of the whole system drifts faster and faster, the <em>fix
momentum</em> can be used. For instance:</p>
<divclass="highlight-python"><divclass="highlight"><pre>fix MOMENTUM all momentum 100 linear 1 1 1
</pre></div>
</div>
<p>It is a bit more tricky to run a NPT simulation with Nose-Hoover
barostat and thermostat. First, the volume should be integrated only
once. So the fix for DCs and atoms should be <em>npt</em> while the fix for
DPs should be <em>nvt</em> (or vice versa). Second, the <em>fix npt</em> computes a
global pressure and thus a global temperature whatever the fix group.
We do want the pressure to correspond to the whole system, but we want
the temperature to correspond to the fix group only. We must then use
the <em>fix_modify</em> command for this. In the end, the block of
instructions for thermostating and barostating will look like</p>
Built with <ahref="http://sphinx-doc.org/">Sphinx</a> using a <ahref="https://github.com/snide/sphinx_rtd_theme">theme</a> provided by <ahref="https://readthedocs.org">Read the Docs</a>.