which corresponds to a force constant <spanclass="math">\(k_D\)</span> = 4184 kJ/(mol
&Aring;<sup>2</sup>) &mdash; for all types of core-Drude bond, a
global mass <spanclass="math">\(m_D\)</span> = 0.4 g/mol (or u) for all types of Drude
particle, and to calculate the Drude charges for individual atom types
from the atom polarizabilities using equation (1). This choice is
followed in the polarizable CHARMM force field.</li>
<li><aclass="reference internal"href="#schroeder"><spanclass="std std-ref">Schroeder and Steinhauser</span></a> suggest adopting a global
charge <spanclass="math">\(q_D\)</span> = -1.0e and a global mass <spanclass="math">\(m_D\)</span> = 0.1 g/mol (or u)
for all Drude particles, and to calculate the force constant for each
type of core-Drude bond from equation (1). The timesteps used by these
authors are between 0.5 and 2 fs, with the degrees of freedom of the
Drude oscillators kept cold at 1 K. In both these force fields
hydrogen atoms are treated as non-polarizable.</li>
</ul>
<p>The motion of of the Drude particles can be calculated by minimizing
the energy of the induced dipoles at each timestep, by an interative,
self-consistent procedure. The Drude particles can be massless and
therefore do not contribute to the kinetic energy. However, the
relaxed method is computationall slow. An extended-lagrangian method
can be used to calculate the positions of the Drude particles, but
this requires them to have mass. It is important in this case to
decouple the degrees of freedom associated with the Drude oscillators
from those of the normal atoms. Thermalizing the Drude dipoles at
temperatures comparable to the rest of the simulation leads to several
problems (kinetic energy transfer, very short timestep, etc.), which
can be remediated by the “cold Drude” technique (<aclass="reference internal"href="#lamoureux"><spanclass="std std-ref">Lamoureux and Roux</span></a>).</p>
<p>Two closely related models are used to represent polarization through
“charges on a spring”: the core-shell model and the Drude
model. Although the basic idea is the same, the core-shell model is
normally used for ionic/crystalline materials, whereas the Drude model
is normally used for molecular systems and fluid states. In ionic
crystals the symmetry around each ion and the distance between them
are such that the core-shell model is sufficiently stable. But to be
applicable to molecular/covalent systems the Drude model includes two
important features:</p>
<olclass="arabic simple">
<li>The possibility to thermostat the additional degrees of freedom</li>
</ol>
<blockquote>
<div>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.</div></blockquote>
<olclass="arabic simple">
<li>The Drude dipoles on covalently bonded atoms interact too strongly
due to the short distances, so an atom may capture the Drude particle
(shell) of a neighbor, or the induced dipoles within the same molecule
may align too much. To avoid this, damping at short of the
interactions between the point charges composing the induced dipole
can be done by <aclass="reference internal"href="#thole"><spanclass="std std-ref">Thole</span></a> functions.</li>
</ol>
<hrclass="docutils"/>
<p><strong>Preparation of the data file</strong></p>
<p>The data file is similar to a standard LAMMPS data file for
<em>atom_style full</em>. The DPs and the <em>harmonic bonds</em> connecting them
to their DC should appear in the data file as normal atoms and bonds.</p>
<p>You can use the <em>polarizer</em> tool (Python script distributed with the
USER-DRUDE package) to convert a non-polarizable data file (here
<em>data.102494.lmp</em>) to a polarizable data file (<em>data-p.lmp</em>)</p>
<p>Now for the thermalization, the simplest choice is to use the <aclass="reference internal"href="fix_langevin_drude.html"><spanclass="doc">fix langevin/drude</span></a>.</p>
<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"><spanclass="doc">*fix drude/transform*</span></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>
<p>For our phenol example, the groups would be defined as</p>
<divclass="highlight-default"><divclass="highlight"><pre><span></span><spanclass="n">group</span><spanclass="n">ATOMS</span><spanclass="nb">type</span><spanclass="mi">1</span><spanclass="mi">2</span><spanclass="mi">3</span><spanclass="mi">4</span><spanclass="mi">5</span><spanclass="c1"># DCs and non-polarizable atoms</span>
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>.