2015-07-30 22:38:28 +08:00
|
|
|
<HTML>
|
|
|
|
|
<script type="text/javascript"
|
|
|
|
|
src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
|
|
|
|
|
</script>
|
|
|
|
|
<script type="text/x-mathjax-config">
|
|
|
|
|
MathJax.Hub.Config({ TeX: { equationNumbers: {autoNumber: "AMS"} } });
|
|
|
|
|
</script>
|
|
|
|
|
|
|
|
|
|
<H3>Tutorial for Thermalized Drude oscillators in LAMMPS
|
|
|
|
|
</H3>
|
|
|
|
|
<P>This tutorial explains how to use Drude oscillators in LAMMPS to
|
|
|
|
|
simulate polarizable systems using the USER-DRUDE package. As an
|
|
|
|
|
illustration, the input files for a simulation of 250 phenol molecules
|
|
|
|
|
are documented. First of all, LAMMPS has to be compiled with the
|
|
|
|
|
USER-DRUDE package activated. Then, the data file and input scripts
|
|
|
|
|
have to be modified to include the Drude dipoles and how to handle
|
|
|
|
|
them.
|
|
|
|
|
</P>
|
|
|
|
|
<HR>
|
|
|
|
|
|
|
|
|
|
<P><B>Overview of Drude induced dipoles</B>
|
|
|
|
|
</P>
|
|
|
|
|
<P>Polarizable atoms acquire an induced electric dipole moment under the
|
2015-07-25 01:57:41 +08:00
|
|
|
action of an external electric field, for example the electric field
|
|
|
|
|
created by the surrounding particles. Drude oscillators represent
|
|
|
|
|
these dipoles by two fixed charges: the core (DC) and the Drude
|
|
|
|
|
particle (DP) bound by a harmonic potential. The Drude particle can be
|
|
|
|
|
thought of as the electron cloud whose center can be displaced from
|
2015-07-30 22:38:28 +08:00
|
|
|
the position of the the corresponding nucleus.
|
|
|
|
|
</P>
|
|
|
|
|
<P>The sum of the masses of a core-Drude pair should be the mass of the
|
|
|
|
|
initial (unsplit) atom, \(m_C + m_D = m\). The sum of their charges
|
|
|
|
|
should be the charge of the initial (unsplit) atom, \(q_C + q_D = q\).
|
2015-07-25 01:57:41 +08:00
|
|
|
A harmonic potential between the core and Drude partners should be
|
2015-07-30 22:38:28 +08:00
|
|
|
present, with force constant \(k_D\) and an equilibrium distance of
|
|
|
|
|
zero. The (half-)stiffness of the <A HREF = "bond_harmonic.html">harmonic bond</A>
|
|
|
|
|
\(K_D = k_D/2\) and the Drude charge \(q_D\) are related to the atom
|
|
|
|
|
polarizability \(\alpha\) by
|
|
|
|
|
</P>
|
|
|
|
|
<P>\begin{equation} K_D = \frac 1 2\, \frac {q_D^2} \alpha
|
|
|
|
|
\end{equation}
|
|
|
|
|
</P>
|
|
|
|
|
<P>Ideally, the mass of the Drude particle should be small, and the
|
2015-07-25 01:57:41 +08:00
|
|
|
stiffness of the harmonic bond should be large, so that the Drude
|
|
|
|
|
particle remains close ot the core. The values of Drude mass, Drude
|
|
|
|
|
charge, and force constant can be chosen following different
|
|
|
|
|
strategies, as in the following examples of polarizable force
|
2015-07-30 22:38:28 +08:00
|
|
|
fields.
|
|
|
|
|
</P>
|
|
|
|
|
<UL><LI><A HREF = "#Lamoureux">Lamoureux and Roux</A> suggest adopting a global
|
|
|
|
|
half-stiffness, \(K_D\) = 500 kcal/(mol Å<sup>2</sup>) —
|
|
|
|
|
which corresponds to a force constant \(k_D\) = 4184 kJ/(mol
|
|
|
|
|
Å<sup>2</sup>) — for all types of core-Drude bond, a
|
|
|
|
|
global mass \(m_D\) = 0.4 g/mol (or u) for all types of Drude
|
2015-07-25 01:57:41 +08:00
|
|
|
particle, and to calculate the Drude charges for individual atom types
|
|
|
|
|
from the atom polarizabilities using equation (1). This choice is
|
2015-07-30 22:38:28 +08:00
|
|
|
followed in the polarizable CHARMM force field.
|
|
|
|
|
|
|
|
|
|
<LI><A HREF = "#Schroeder">Schroeder and Steinhauser</A> suggest adopting a global
|
|
|
|
|
charge \(q_D\) = -1.0e and a global mass \(m_D\) = 0.1 g/mol (or u)
|
2015-07-25 01:57:41 +08:00
|
|
|
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
|
2015-07-30 22:38:28 +08:00
|
|
|
hydrogen atoms are treated as non-polarizable.
|
|
|
|
|
</UL>
|
|
|
|
|
<P>The motion of of the Drude particles can be calculated by minimizing
|
2015-07-25 01:57:41 +08:00
|
|
|
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
|
2015-07-30 22:38:28 +08:00
|
|
|
can be remediated by the "cold Drude" technique (<A HREF = "#Lamoureux">Lamoureux and
|
|
|
|
|
Roux</A>).
|
|
|
|
|
</P>
|
|
|
|
|
<P>Two closely related models are used to represent polarization through
|
|
|
|
|
"charges on a spring": the core-shell model and the Drude
|
2015-07-25 01:57:41 +08:00
|
|
|
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
|
2015-07-30 22:38:28 +08:00
|
|
|
important features:
|
|
|
|
|
</P>
|
|
|
|
|
<OL><LI>The possibility to thermostat the additional degrees of freedom
|
|
|
|
|
associated with the induced dipoles at very low temperature, in terms
|
2015-07-25 01:57:41 +08:00
|
|
|
of the reduced coordinates of the Drude particles with respect to
|
|
|
|
|
their cores. This makes the trajectory close to that of relaxed
|
2015-07-30 22:38:28 +08:00
|
|
|
induced dipoles.
|
|
|
|
|
|
|
|
|
|
<LI>The Drude dipoles on covalently bonded atoms interact too strongly
|
2015-07-25 01:57:41 +08:00
|
|
|
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
|
2015-07-30 22:38:28 +08:00
|
|
|
can be done by <A HREF = "#Thole">Thole</A> functions.
|
|
|
|
|
</OL>
|
|
|
|
|
<HR>
|
|
|
|
|
|
|
|
|
|
<P><B>Preparation of the data file</B>
|
|
|
|
|
</P>
|
|
|
|
|
<P>The data file is similar to a standard LAMMPS data file for
|
|
|
|
|
<I>atom_style full</I>. The DPs and the <I>harmonic bonds</I> connecting them
|
|
|
|
|
to their DC should appear in the data file as normal atoms and bonds.
|
|
|
|
|
</P>
|
|
|
|
|
<P>You can use the <I>polarizer</I> tool (Python script distributed with the
|
2015-07-25 01:57:41 +08:00
|
|
|
USER-DRUDE package) to convert a non-polarizable data file (here
|
2015-07-30 22:38:28 +08:00
|
|
|
<I>data.102494.lmp</I>) to a polarizable data file (<I>data-p.lmp</I>)
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>polarizer -q -f phenol.dff data.102494.lmp data-p.lmp
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>This will automatically insert the new atoms and bonds.
|
2015-07-25 01:57:41 +08:00
|
|
|
The masses and charges of DCs and DPs are computed
|
2015-07-30 22:38:28 +08:00
|
|
|
from <I>phenol.dff</I>, as well as the DC-DP bond constants. The file
|
|
|
|
|
<I>phenol.dff</I> contains the polarizabilities of the atom types
|
|
|
|
|
and the mass of the Drude particles, for instance:
|
|
|
|
|
</P>
|
|
|
|
|
<PRE># units: kJ/mol, A, deg
|
2015-07-25 01:57:41 +08:00
|
|
|
# kforce is in the form k/2 r_D^2
|
|
|
|
|
# type m_D/u q_D/e k_D alpha/A3 thole
|
|
|
|
|
OH 0.4 -1.0 4184.0 0.63 0.67
|
|
|
|
|
CA 0.4 -1.0 4184.0 1.36 2.51
|
2015-07-30 22:38:28 +08:00
|
|
|
CAI 0.4 -1.0 4184.0 1.09 2.51
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>The hydrogen atoms are absent from this file, so they will be treated
|
2015-07-25 01:57:41 +08:00
|
|
|
as non-polarizable atoms. In the non-polarizable data file
|
2015-07-30 22:38:28 +08:00
|
|
|
<I>data.102494.lmp</I>, atom names corresponding to the atom type numbers
|
|
|
|
|
have to be specified as comments at the end of lines of the <I>Masses</I>
|
2015-07-25 01:57:41 +08:00
|
|
|
section. You probably need to edit it to add these names. It should
|
2015-07-30 22:38:28 +08:00
|
|
|
look like
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>Masses
|
|
|
|
|
</PRE>
|
|
|
|
|
<PRE>1 12.011 # CAI
|
2015-07-25 01:57:41 +08:00
|
|
|
2 12.011 # CA
|
|
|
|
|
3 15.999 # OH
|
|
|
|
|
4 1.008 # HA
|
2015-07-30 22:38:28 +08:00
|
|
|
5 1.008 # HO
|
|
|
|
|
</PRE>
|
|
|
|
|
<HR>
|
|
|
|
|
|
|
|
|
|
<P><B>Basic input file</B>
|
|
|
|
|
</P>
|
|
|
|
|
<P>The atom style should be set to (or derive from) <I>full</I>, so that you
|
|
|
|
|
can define atomic charges and molecular bonds, angles, dihedrals...
|
|
|
|
|
</P>
|
|
|
|
|
<P>The <I>polarizer</I> tool also outputs certain lines related to the input
|
2015-07-25 01:57:41 +08:00
|
|
|
script (the use of these lines will be explained below). In order for
|
|
|
|
|
LAMMPS to recognize that you are using Drude oscillators, you should
|
2015-07-30 22:38:28 +08:00
|
|
|
use the fix <I>drude</I>. The command is
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>fix DRUDE all drude C C C N N D D D
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>The N, C, D following the <I>drude</I> keyword have the following meaning:
|
2015-07-27 22:54:15 +08:00
|
|
|
There is one tag for each atom type. This tag is C for DCs, D for DPs
|
|
|
|
|
and N for non-polarizable atoms. Here the atom types 1 to 3 (C and O
|
|
|
|
|
atoms) are DC, atom types 4 and 5 (H atoms) are non-polarizable and
|
2015-07-30 22:38:28 +08:00
|
|
|
the atom types 6 to 8 are the newly created DPs.
|
|
|
|
|
</P>
|
|
|
|
|
<P>By recognizing the fix <I>drude</I>, LAMMPS will find and store matching
|
2015-07-25 01:57:41 +08:00
|
|
|
DC-DP pairs and will treat DP as equivalent to their DC in the
|
2015-07-30 22:38:28 +08:00
|
|
|
<I>special bonds</I> relations. It may be necessary to extend the space
|
2015-07-25 01:57:41 +08:00
|
|
|
for storing such special relations. In this case extra space should
|
2015-07-30 22:38:28 +08:00
|
|
|
be reserved by using the <I>extra</I> keyword of the <I>special_bonds</I>
|
2015-07-25 01:57:41 +08:00
|
|
|
command. With our phenol, there is 1 more special neighbor for which
|
|
|
|
|
space is required. Otherwise LAMMPS crashes and gives the required
|
2015-07-30 22:38:28 +08:00
|
|
|
value.
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>special_bonds lj/coul 0.0 0.0 0.5 extra 1
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>Let us assume we want to run a simple NVT simulation at 300 K. Note
|
2015-07-25 01:57:41 +08:00
|
|
|
that Drude oscillators need to be thermalized at a low temperature in
|
|
|
|
|
order to approximate a self-consistent field (SCF), therefore it is not
|
|
|
|
|
possible to simulate an NVE ensemble with this package. Since dipoles
|
2015-07-30 22:38:28 +08:00
|
|
|
are approximated by a charged DC-DP pair, the <I>pair_style</I> must
|
|
|
|
|
include Coulomb interactions, for instance <I>lj/cut/coul/long</I> with
|
|
|
|
|
<I>kspace_style pppm</I>. For example, with a cutoff of 10. and a precision
|
|
|
|
|
1.e-4:
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>pair_style lj/cut/coul/long 10.0
|
|
|
|
|
kspace_style pppm 1.0e-4
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>As compared to the non-polarizable input file, <I>pair_coeff</I> lines need
|
2015-07-25 01:57:41 +08:00
|
|
|
to be added for the DPs. Since the DPs have no Lennard-Jones
|
2015-07-30 22:38:28 +08:00
|
|
|
interactions, their <I>epsilon</I> is 0. so the only <I>pair_coeff</I> line
|
|
|
|
|
that needs to be added is
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>pair_coeff * 6* 0.0 0.0 # All-DPs
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>Now for the thermalization, the simplest choice is to use the <A HREF = "fix_langevin_drude.html">fix
|
|
|
|
|
langevin/drude</A>.
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>fix LANG all langevin/drude 300. 100 12435 1. 20 13977
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>This applies a Langevin thermostat at temperature 300. to the centers
|
2015-07-25 01:57:41 +08:00
|
|
|
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
|
2015-07-30 22:38:28 +08:00
|
|
|
atoms need to be in this fix's group. LAMMPS will thermostate the DPs
|
2015-07-25 01:57:41 +08:00
|
|
|
together with their DC. For this, ghost atoms need to know their
|
2015-07-30 22:38:28 +08:00
|
|
|
velocities. Thus you need to add the following command:
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>comm_modify vel yes
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>In order to avoid that the center of mass of the whole system
|
2015-07-25 01:57:41 +08:00
|
|
|
drifts due to the random forces of the Langevin thermostat on DCs, you
|
2015-07-30 22:38:28 +08:00
|
|
|
can add the <I>zero yes</I> option at the end of the fix line.
|
|
|
|
|
</P>
|
|
|
|
|
<P>If the fix <I>shake</I> is used to constrain the C-H bonds, it should be
|
|
|
|
|
invoked after the fix <I>langevin/drude</I> for more accuracy.
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>fix SHAKE ATOMS shake 0.0001 20 0 t 4 5
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>IMPORTANT NOTE: The group of the fix <I>shake</I> must not include the DPs.
|
|
|
|
|
If the group <I>ATOMS</I> is defined by non-DPs atom types, you could use
|
|
|
|
|
</P>
|
|
|
|
|
<P>Since the fix <I>langevin/drude</I> does not perform time integration (just
|
2015-07-25 01:57:41 +08:00
|
|
|
modification of forces but no position/velocity updates), the fix
|
2015-07-30 22:38:28 +08:00
|
|
|
<I>nve</I> should be used in conjunction.
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>fix NVE all nve
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>Finally, do not forget to update the atom type elements if you use
|
|
|
|
|
them in a <I>dump_modify ... element ...</I> command, by adding the element
|
|
|
|
|
type of the DPs. Here for instance
|
|
|
|
|
</P>
|
|
|
|
|
<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>
|
|
|
|
|
<P>The input file should now be ready for use!
|
|
|
|
|
</P>
|
|
|
|
|
<P>You will notice that the global temperature <I>thermo_temp</I> computed by
|
2015-07-25 01:57:41 +08:00
|
|
|
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
|
2015-07-30 22:38:28 +08:00
|
|
|
to their DCs, you should use the <A HREF = "compute_temp_drude.html">compute
|
|
|
|
|
temp_drude</A>
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>compute TDRUDE all temp/drude
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>And then output the correct temperatures of the Drude oscillators
|
|
|
|
|
using <I>thermo_style custom</I> with respectively <I>c_TDRUDE[1]</I> and
|
|
|
|
|
<I>c_TDRUDE[2]</I>. These should be close to 300.0 and 1.0 on average.
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>thermo_style custom step temp c_TDRUDE[1] c_TDRUDE[2]
|
|
|
|
|
</PRE>
|
|
|
|
|
<HR>
|
|
|
|
|
|
|
|
|
|
<P><B>Thole screening</B>
|
|
|
|
|
</P>
|
|
|
|
|
<P>Dipolar interactions represented by point charges on springs may not
|
2015-07-25 01:57:41 +08:00
|
|
|
be stable, for example if the atomic polarizability is too high for
|
|
|
|
|
instance, a DP can escape from its DC and be captured by another DC,
|
|
|
|
|
which makes the force and energy diverge and the simulation
|
|
|
|
|
crash. Even without reaching this extreme case, the correlation
|
|
|
|
|
between nearby dipoles on the same molecule may be exagerated. Often,
|
|
|
|
|
special bond relations prevent bonded neighboring atoms to see the
|
2015-07-30 22:38:28 +08:00
|
|
|
charge of each other's DP, so that the problem does not always appear.
|
2015-07-25 01:57:41 +08:00
|
|
|
It is possible to use screened dipole dipole interactions by using the
|
2015-07-30 22:38:28 +08:00
|
|
|
<A HREF = "pair_thole.html"><I>pair_style thole</I></A>. This is implemented as a
|
2015-07-25 01:57:41 +08:00
|
|
|
correction to the Coulomb pair_styles, which dampens at short distance
|
|
|
|
|
the interactions between the charges representing the induced dipoles.
|
2015-07-30 22:38:28 +08:00
|
|
|
It is to be used as <I>hybrid/overlay</I> with any standard <I>coul</I> pair
|
|
|
|
|
style. In our example, we would use
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>pair_style hybrid/overlay lj/cut/coul/long 10.0 thole 2.6 10.0
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>This tells LAMMPS that we are using two pair_styles. The first one is
|
|
|
|
|
as above (<I>lj/cut/coul/long 10.0</I>). The second one is a <I>thole</I>
|
|
|
|
|
pair_style with default screening factor 2.6 (<A HREF = "#Noskov">Noskov</A>) and
|
|
|
|
|
cutoff 10.0.
|
|
|
|
|
</P>
|
|
|
|
|
<P>Since <I>hybrid/overlay</I> does not support mixing rules, the interaction
|
|
|
|
|
coefficients of all the pairs of atom types with i < j should be
|
|
|
|
|
explicitly defined. The output of the <I>polarizer</I> script can be used
|
|
|
|
|
to complete the <I>pair_coeff</I> section of the input file. In our
|
|
|
|
|
example, this will look like:
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>pair_coeff 1 1 lj/cut/coul/long 0.0700 3.550
|
2015-07-25 01:57:41 +08:00
|
|
|
pair_coeff 1 2 lj/cut/coul/long 0.0700 3.550
|
|
|
|
|
pair_coeff 1 3 lj/cut/coul/long 0.1091 3.310
|
|
|
|
|
pair_coeff 1 4 lj/cut/coul/long 0.0458 2.985
|
|
|
|
|
pair_coeff 2 2 lj/cut/coul/long 0.0700 3.550
|
|
|
|
|
pair_coeff 2 3 lj/cut/coul/long 0.1091 3.310
|
|
|
|
|
pair_coeff 2 4 lj/cut/coul/long 0.0458 2.985
|
|
|
|
|
pair_coeff 3 3 lj/cut/coul/long 0.1700 3.070
|
|
|
|
|
pair_coeff 3 4 lj/cut/coul/long 0.0714 2.745
|
|
|
|
|
pair_coeff 4 4 lj/cut/coul/long 0.0300 2.420
|
|
|
|
|
pair_coeff * 5 lj/cut/coul/long 0.0000 0.000
|
|
|
|
|
pair_coeff * 6* lj/cut/coul/long 0.0000 0.000
|
|
|
|
|
pair_coeff 1 1 thole 1.090 2.510
|
|
|
|
|
pair_coeff 1 2 thole 1.218 2.510
|
|
|
|
|
pair_coeff 1 3 thole 0.829 1.590
|
|
|
|
|
pair_coeff 1 6 thole 1.090 2.510
|
|
|
|
|
pair_coeff 1 7 thole 1.218 2.510
|
|
|
|
|
pair_coeff 1 8 thole 0.829 1.590
|
|
|
|
|
pair_coeff 2 2 thole 1.360 2.510
|
|
|
|
|
pair_coeff 2 3 thole 0.926 1.590
|
|
|
|
|
pair_coeff 2 6 thole 1.218 2.510
|
|
|
|
|
pair_coeff 2 7 thole 1.360 2.510
|
|
|
|
|
pair_coeff 2 8 thole 0.926 1.590
|
|
|
|
|
pair_coeff 3 3 thole 0.630 0.670
|
|
|
|
|
pair_coeff 3 6 thole 0.829 1.590
|
|
|
|
|
pair_coeff 3 7 thole 0.926 1.590
|
|
|
|
|
pair_coeff 3 8 thole 0.630 0.670
|
|
|
|
|
pair_coeff 6 6 thole 1.090 2.510
|
|
|
|
|
pair_coeff 6 7 thole 1.218 2.510
|
|
|
|
|
pair_coeff 6 8 thole 0.829 1.590
|
|
|
|
|
pair_coeff 7 7 thole 1.360 2.510
|
|
|
|
|
pair_coeff 7 8 thole 0.926 1.590
|
2015-07-30 22:38:28 +08:00
|
|
|
pair_coeff 8 8 thole 0.630 0.670
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>For the <I>thole</I> pair style the coefficients are
|
|
|
|
|
</P>
|
|
|
|
|
<OL><LI>the atom polarizability in units of cubic length
|
|
|
|
|
|
|
|
|
|
<LI>the screening factor of the Thole function (optional, default value
|
|
|
|
|
specified by the pair_style command)
|
|
|
|
|
|
|
|
|
|
<LI>the cutoff (optional, default value defined by the pair_style command)
|
|
|
|
|
|
|
|
|
|
</OL>
|
|
|
|
|
<P>The special neighbors have charge-charge and charge-dipole
|
|
|
|
|
interactions screened by the <I>coul</I> factors of the <I>special_bonds</I>
|
2015-07-27 22:54:15 +08:00
|
|
|
command (0.0, 0.0, and 0.5 in the example above). Without using the
|
2015-07-30 22:38:28 +08:00
|
|
|
pair_style <I>thole</I>, dipole-dipole interactions are screened by the
|
|
|
|
|
same factor. By using the pair_style <I>thole</I>, dipole-dipole
|
|
|
|
|
interactions are screened by Thole's function, whatever their special
|
2015-07-25 01:57:41 +08:00
|
|
|
relationship (except within each DC-DP pair of course). Consider for
|
2015-07-30 22:38:28 +08:00
|
|
|
example 1-2 neighbors: using the pair_style <I>thole</I>, their dipoles
|
|
|
|
|
will see each other (despite the <I>coul</I> factor being 0.) and the
|
|
|
|
|
interactions between these dipoles will be damped by Thole's function.
|
|
|
|
|
</P>
|
|
|
|
|
<HR>
|
|
|
|
|
|
|
|
|
|
<P><B>Thermostats and barostats</B>
|
|
|
|
|
</P>
|
|
|
|
|
<P>Using a Nose-Hoover barostat with the <I>langevin/drude</I> thermostat is
|
|
|
|
|
straightforward using fix <I>nph</I> instead of <I>nve</I>. For example:
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>fix NPH all nph iso 1. 1. 500
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>It is also possible to use a Nose-Hoover instead of a Langevin
|
|
|
|
|
thermostat. This requires to use <A HREF = "fix_drude_transform.html"><I>fix
|
|
|
|
|
drude/transform</I></A> just before and after the
|
|
|
|
|
time intergation fixes. The <I>fix drude/transform/direct</I> converts the
|
2015-07-25 01:57:41 +08:00
|
|
|
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
|
2015-07-30 22:38:28 +08:00
|
|
|
with respect to their DC. The <I>fix drude/transform/inverse</I> performs
|
2015-07-25 01:57:41 +08:00
|
|
|
the reverse transformation. For a NVT simulation, with the DCs and
|
2015-07-30 22:38:28 +08:00
|
|
|
atoms at 300 K and the DPs at 1 K relative to their DC one would use
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>fix DIRECT all drude/transform/direct
|
2015-07-25 01:57:41 +08:00
|
|
|
fix NVT1 ATOMS nvt temp 300. 300. 100
|
|
|
|
|
fix NVT2 DRUDES nvt temp 1. 1. 20
|
2015-07-30 22:38:28 +08:00
|
|
|
fix INVERSE all drude/transform/inverse
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>For our phenol example, the groups would be defined as
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>group ATOMS type 1 2 3 4 5 # DCs and non-polarizable atoms
|
2015-07-25 01:57:41 +08:00
|
|
|
group CORES type 1 2 3 # DCs
|
2015-07-30 22:38:28 +08:00
|
|
|
group DRUDES type 6 7 8 # DPs
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>Note that with the fixes <I>drude/transform</I>, it is not required to
|
|
|
|
|
specify <I>comm_modify vel yes</I> because the fixes do it anyway (several
|
2015-07-25 01:57:41 +08:00
|
|
|
times and for the forces also). To avoid the flying ice cube artifact
|
2015-07-30 22:38:28 +08:00
|
|
|
<A HREF = "#Lamoureux">(Lamoureux)</A>, where the atoms progressively freeze and the
|
|
|
|
|
center of mass of the whole system drifts faster and faster, the <I>fix
|
|
|
|
|
momentum</I> can be used. For instance:
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>fix MOMENTUM all momentum 100 linear 1 1 1
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>It is a bit more tricky to run a NPT simulation with Nose-Hoover
|
2015-07-25 01:57:41 +08:00
|
|
|
barostat and thermostat. First, the volume should be integrated only
|
2015-07-30 22:38:28 +08:00
|
|
|
once. So the fix for DCs and atoms should be <I>npt</I> while the fix for
|
|
|
|
|
DPs should be <I>nvt</I> (or vice versa). Second, the <I>fix npt</I> computes a
|
2015-07-25 01:57:41 +08:00
|
|
|
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
|
2015-07-30 22:38:28 +08:00
|
|
|
the <I>fix_modify</I> command for this. In the end, the block of
|
|
|
|
|
instructions for thermostating and barostating will look like
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>compute TATOMS ATOMS temp
|
2015-07-25 01:57:41 +08:00
|
|
|
fix DIRECT all drude/transform/direct
|
|
|
|
|
fix NPT ATOMS npt temp 300. 300. 100 iso 1. 1. 500
|
|
|
|
|
fix_modify NPT temp TATOMS press thermo_press
|
|
|
|
|
fix NVT DRUDES nvt temp 1. 1. 20
|
2015-07-30 22:38:28 +08:00
|
|
|
fix INVERSE all drude/transform/inverse
|
|
|
|
|
</PRE>
|
|
|
|
|
<HR>
|
|
|
|
|
|
|
|
|
|
<P><B>Rigid bodies</B>
|
|
|
|
|
</P>
|
|
|
|
|
<P>You may want to simulate molecules as rigid bodies (but polarizable).
|
|
|
|
|
Common cases are water models such as <A HREF = "#SWM4-NDP">SWM4-NDP</A>, which is a
|
2015-07-25 01:57:41 +08:00
|
|
|
kind of polarizable TIP4P water. The rigid bodies and the DPs should
|
|
|
|
|
be integrated separately, even with the Langevin thermostat. Let us
|
2015-07-30 22:38:28 +08:00
|
|
|
review the different thermostats and ensemble combinations.
|
|
|
|
|
</P>
|
|
|
|
|
<P>NVT ensemble using Langevin thermostat:
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>comm_modify vel yes
|
2015-07-25 01:57:41 +08:00
|
|
|
fix LANG all langevin/drude 300. 100 12435 1. 20 13977
|
|
|
|
|
fix RIGID ATOMS rigid/nve/small molecule
|
2015-07-30 22:38:28 +08:00
|
|
|
fix NVE DRUDES nve
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>NVT ensemble using Nose-Hoover thermostat:
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>fix DIRECT all drude/transform/direct
|
2015-07-25 01:57:41 +08:00
|
|
|
fix RIGID ATOMS rigid/nvt/small molecule temp 300. 300. 100
|
|
|
|
|
fix NVT DRUDES nvt temp 1. 1. 20
|
2015-07-30 22:38:28 +08:00
|
|
|
fix INVERSE all drude/transform/inverse
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>NPT ensemble with Langevin thermostat:
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>comm_modify vel yes
|
2015-07-25 01:57:41 +08:00
|
|
|
fix LANG all langevin/drude 300. 100 12435 1. 20 13977
|
|
|
|
|
fix RIGID ATOMS rigid/nph/small molecule iso 1. 1. 500
|
2015-07-30 22:38:28 +08:00
|
|
|
fix NVE DRUDES nve
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>NPT ensemble using Nose-Hoover thermostat:
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>compute TATOM ATOMS temp
|
2015-07-25 01:57:41 +08:00
|
|
|
fix DIRECT all drude/transform/direct
|
|
|
|
|
fix RIGID ATOMS rigid/npt/small molecule temp 300. 300. 100 iso 1. 1. 500
|
|
|
|
|
fix_modify RIGID temp TATOM press thermo_press
|
|
|
|
|
fix NVT DRUDES nvt temp 1. 1. 20
|
2015-07-30 22:38:28 +08:00
|
|
|
fix INVERSE all drude/transform/inverse
|
|
|
|
|
</PRE>
|
|
|
|
|
<HR>
|
|
|
|
|
|
|
|
|
|
<A NAME = "Lamoureux"></A>
|
|
|
|
|
|
|
|
|
|
<P><B>(Lamoureux)</B> Lamoureux and Roux, J Chem Phys, 119, 3025-3039 (2003)
|
|
|
|
|
</P>
|
|
|
|
|
<A NAME = "Schroeder"></A>
|
|
|
|
|
|
|
|
|
|
<P><B>(Schroeder)</B> Schröder and Steinhauser, J Chem Phys, 133,
|
|
|
|
|
154511 (2010).
|
|
|
|
|
</P>
|
|
|
|
|
<A NAME = "Jiang"></A>
|
|
|
|
|
|
|
|
|
|
<P><B>(Jiang)</B> Jiang, Hardy, Phillips, MacKerell, Schulten, and Roux,
|
|
|
|
|
J Phys Chem Lett, 2, 87-92 (2011).
|
|
|
|
|
</P>
|
|
|
|
|
<A NAME = "Thole"></A>
|
|
|
|
|
|
|
|
|
|
<P><B>(Thole)</B> Chem Phys, 59, 341 (1981).
|
|
|
|
|
</P>
|
|
|
|
|
<A NAME = "Noskov"></A>
|
|
|
|
|
|
|
|
|
|
<P><B>(Noskov)</B> Noskov, Lamoureux and Roux, J Phys Chem B, 109, 6705 (2005).
|
|
|
|
|
</P>
|
|
|
|
|
<A NAME = "SWM4-NDP"></A>
|
|
|
|
|
|
|
|
|
|
<P><B>(SWM4-NDP)</B> Lamoureux, Harder, Vorobyov, Roux, MacKerell, Chem Phys
|
|
|
|
|
Let, 418, 245-249 (2006)
|
|
|
|
|
</P>
|
|
|
|
|
</HTML>
|
