2015-07-30 22:53:28 +08:00
<!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 > < no title> — LAMMPS 15 May 2015 version 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 15 May 2015 version documentation" href = "index.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 >
< 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" > < a class = "reference internal" href = "Section_howto.html" > 6. How-to discussions< / a > < / 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 > < no title> < / 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 >
< div role = "main" class = "document" itemscope = "itemscope" itemtype = "http://schema.org/Article" >
< div itemprop = "articleBody" >
< p > < strong > Overview of Drude induced dipoles< / strong > < / 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:53: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:53:28 +08:00
present, with force constant (k_D) and an equilibrium distance of
zero. The (half-)stiffness of the < a class = "reference internal" href = "bond_harmonic.html" > < em > harmonic bond< / em > < / 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:53:28 +08:00
fields.< / p >
< ul class = "simple" >
< li > < a class = "reference internal" href = "#lamoureux" > < span > Lamoureux and Roux< / span > < / a > suggest adopting a global
half-stiffness, (K_D) = 500 kcal/(mol & Aring;< sup> 2< /sup> ) & mdash;
which corresponds to a force constant (k_D) = 4184 kJ/(mol
& Aring;< sup> 2< /sup> ) & mdash; 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:53:28 +08:00
followed in the polarizable CHARMM force field.< / li >
< li > < a class = "reference internal" href = "#schroeder" > < span > Schroeder and Steinhauser< / span > < / 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:53:28 +08:00
hydrogen atoms are treated as non-polarizable.< / li >
< / 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:53:28 +08:00
can be remediated by the “ cold Drude” technique (< a class = "reference internal" href = "#lamoureux" > < span > 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
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:53:28 +08:00
important features:< / p >
< ol class = "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
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:53:28 +08:00
induced dipoles.< / div > < / blockquote >
< ol class = "arabic simple" >
< 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:53:28 +08:00
can be done by < a class = "reference internal" href = "#thole" > < span > Thole< / span > < / a > functions.< / li >
< / ol >
< hr class = "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
2015-07-25 01:57:41 +08:00
USER-DRUDE package) to convert a non-polarizable data file (here
2015-07-30 22:53:28 +08:00
< em > data.102494.lmp< / em > ) to a polarizable data file (< em > data-p.lmp< / em > )< / p >
< div class = "highlight-python" > < div class = "highlight" > < pre > polarizer -q -f phenol.dff data.102494.lmp data-p.lmp
< / pre > < / div >
< / div >
< 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:53:28 +08:00
from < em > phenol.dff< / em > , as well as the DC-DP bond constants. The file
< em > phenol.dff< / em > contains the polarizabilities of the atom types
and the mass of the Drude particles, for instance:< / p >
< div class = "highlight-python" > < div class = "highlight" > < 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:53:28 +08:00
CAI 0.4 -1.0 4184.0 1.09 2.51
< / pre > < / div >
< / div >
< 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:53:28 +08:00
< em > data.102494.lmp< / em > , atom names corresponding to the atom type numbers
have to be specified as comments at the end of lines of the < em > Masses< / em >
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:53:28 +08:00
look like< / p >
< div class = "highlight-python" > < div class = "highlight" > < pre > < span class = "n" > Masses< / span >
< / pre > < / div >
< / div >
< div class = "highlight-python" > < div class = "highlight" > < 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:53:28 +08:00
5 1.008 # HO
< / pre > < / div >
< / div >
< hr class = "docutils" / >
< p > < strong > Basic input file< / strong > < / p >
< p > The atom style should be set to (or derive from) < em > full< / em > , so that you
can define atomic charges and molecular bonds, angles, dihedrals...< / p >
< p > The < em > polarizer< / em > 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:53:28 +08:00
use the fix < em > drude< / em > . The command is< / p >
< div class = "highlight-python" > < div class = "highlight" > < pre > fix DRUDE all drude C C C N N D D D
< / pre > < / div >
< / div >
< p > The N, C, D following the < em > drude< / em > 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:53:28 +08:00
the atom types 6 to 8 are the newly created DPs.< / p >
< p > By recognizing the fix < em > drude< / em > , 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:53:28 +08:00
< em > special bonds< / em > 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:53:28 +08:00
be reserved by using the < em > extra< / em > keyword of the < em > special_bonds< / em >
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:53:28 +08:00
value.< / p >
< div class = "highlight-python" > < div class = "highlight" > < pre > special_bonds lj/coul 0.0 0.0 0.5 extra 1
< / pre > < / div >
< / div >
< 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:53:28 +08:00
are approximated by a charged DC-DP pair, the < em > pair_style< / em > must
include Coulomb interactions, for instance < em > lj/cut/coul/long< / em > with
< em > kspace_style pppm< / em > . For example, with a cutoff of 10. and a precision
1.e-4:< / p >
< div class = "highlight-python" > < div class = "highlight" > < pre > pair_style lj/cut/coul/long 10.0
kspace_style pppm 1.0e-4
< / pre > < / div >
< / div >
< p > As compared to the non-polarizable input file, < em > pair_coeff< / em > 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:53:28 +08:00
interactions, their < em > epsilon< / em > is 0. so the only < em > pair_coeff< / em > line
that needs to be added is< / p >
< div class = "highlight-python" > < div class = "highlight" > < pre > pair_coeff * 6* 0.0 0.0 # All-DPs
< / pre > < / div >
< / div >
< p > Now for the thermalization, the simplest choice is to use the < a class = "reference internal" href = "fix_langevin_drude.html" > < em > fix langevin/drude< / em > < / a > .< / p >
< div class = "highlight-python" > < div class = "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
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:53: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:53:28 +08:00
velocities. Thus you need to add the following command:< / p >
< div class = "highlight-python" > < div class = "highlight" > < pre > comm_modify vel yes
< / pre > < / div >
< / div >
< 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:53:28 +08:00
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 >
< div class = "highlight-python" > < div class = "highlight" > < pre > fix SHAKE ATOMS shake 0.0001 20 0 t 4 5
< / pre > < / div >
< / div >
< div class = "admonition warning" >
< p class = "first admonition-title" > Warning< / p >
< p class = "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
2015-07-25 01:57:41 +08:00
modification of forces but no position/velocity updates), the fix
2015-07-30 22:53:28 +08:00
< em > nve< / em > should be used in conjunction.< / p >
< div class = "highlight-python" > < div class = "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 >
< div class = "highlight-python" > < div class = "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
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:53:28 +08:00
to their DCs, you should use the < a class = "reference internal" href = "compute_temp_drude.html" > < em > compute temp_drude< / em > < / a > < / p >
< div class = "highlight-python" > < div class = "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 >
< div class = "highlight-python" > < div class = "highlight" > < pre > thermo_style custom step temp c_TDRUDE[1] c_TDRUDE[2]
< / pre > < / div >
< / div >
< hr class = "docutils" / >
< p > < strong > Thole screening< / strong > < / 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:53: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:53:28 +08:00
< a class = "reference internal" href = "pair_thole.html" > < em > *pair_style thole*< / em > < / 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:53:28 +08:00
It is to be used as < em > hybrid/overlay< / em > with any standard < em > coul< / em > pair
style. In our example, we would use< / p >
< div class = "highlight-python" > < div class = "highlight" > < pre > pair_style hybrid/overlay lj/cut/coul/long 10.0 thole 2.6 10.0
< / pre > < / div >
< / div >
< p > This tells LAMMPS that we are using two pair_styles. The first one is
as above (< em > lj/cut/coul/long 10.0< / em > ). The second one is a < em > thole< / em >
pair_style with default screening factor 2.6 (< a class = "reference internal" href = "#noskov" > < span > Noskov< / span > < / a > ) and
cutoff 10.0.< / p >
< p > Since < em > hybrid/overlay< / em > 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 < em > polarizer< / em > script can be used
to complete the < em > pair_coeff< / em > section of the input file. In our
example, this will look like:< / p >
< div class = "highlight-python" > < div class = "highlight" > < 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:53:28 +08:00
pair_coeff 8 8 thole 0.630 0.670
< / pre > < / div >
< / div >
< p > For the < em > thole< / em > pair style the coefficients are< / p >
< ol class = "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 >
< ul class = "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 >
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:53:28 +08:00
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
2015-07-25 01:57:41 +08:00
relationship (except within each DC-DP pair of course). Consider for
2015-07-30 22:53:28 +08:00
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 >
< hr class = "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 >
< div class = "highlight-python" > < div class = "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 < a class = "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
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:53:28 +08:00
with respect to their DC. The < em > fix drude/transform/inverse< / em > performs
2015-07-25 01:57:41 +08:00
the reverse transformation. For a NVT simulation, with the DCs and
2015-07-30 22:53:28 +08:00
atoms at 300 K and the DPs at 1 K relative to their DC one would use< / p >
< div class = "highlight-python" > < div class = "highlight" > < 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:53:28 +08:00
fix INVERSE all drude/transform/inverse
< / pre > < / div >
< / div >
< p > For our phenol example, the groups would be defined as< / p >
< div class = "highlight-python" > < div class = "highlight" > < 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:53:28 +08:00
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
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:53:28 +08:00
< a class = "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 >
< div class = "highlight-python" > < div class = "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
2015-07-25 01:57:41 +08:00
barostat and thermostat. First, the volume should be integrated only
2015-07-30 22:53:28 +08:00
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
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:53:28 +08:00
the < em > fix_modify< / em > command for this. In the end, the block of
instructions for thermostating and barostating will look like< / p >
< div class = "highlight-python" > < div class = "highlight" > < 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:53:28 +08:00
fix INVERSE all drude/transform/inverse
< / pre > < / div >
< / div >
< hr class = "docutils" / >
< p > < strong > Rigid bodies< / strong > < / p >
< p > You may want to simulate molecules as rigid bodies (but polarizable).
Common cases are water models such as < a class = "reference internal" href = "#swm4-ndp" > < span > SWM4-NDP< / span > < / 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:53:28 +08:00
review the different thermostats and ensemble combinations.< / p >
< p > NVT ensemble using Langevin thermostat:< / p >
< div class = "highlight-python" > < div class = "highlight" > < 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:53:28 +08:00
fix NVE DRUDES nve
< / pre > < / div >
< / div >
< p > NVT ensemble using Nose-Hoover thermostat:< / p >
< div class = "highlight-python" > < div class = "highlight" > < 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:53:28 +08:00
fix INVERSE all drude/transform/inverse
< / pre > < / div >
< / div >
< p > NPT ensemble with Langevin thermostat:< / p >
< div class = "highlight-python" > < div class = "highlight" > < 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:53:28 +08:00
fix NVE DRUDES nve
< / pre > < / div >
< / div >
< p > NPT ensemble using Nose-Hoover thermostat:< / p >
< div class = "highlight-python" > < div class = "highlight" > < 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:53:28 +08:00
fix INVERSE all drude/transform/inverse
< / pre > < / div >
< / div >
< hr class = "docutils" / >
< p id = "lamoureux" > < strong > (Lamoureux)< / strong > Lamoureux and Roux, J Chem Phys, 119, 3025-3039 (2003)< / p >
< p id = "schroeder" > < strong > (Schroeder)< / strong > Schr& ouml;der and Steinhauser, J Chem Phys, 133,
154511 (2010).< / p >
< dl class = "docutils" id = "jiang" >
< dt > < strong > (Jiang)< / strong > Jiang, Hardy, Phillips, MacKerell, Schulten, and Roux,< / dt >
< dd > J Phys Chem Lett, 2, 87-92 (2011).< / dd >
< / dl >
< p id = "thole" > < strong > (Thole)< / strong > Chem Phys, 59, 341 (1981).< / p >
< p id = "noskov" > < strong > (Noskov)< / strong > Noskov, Lamoureux and Roux, J Phys Chem B, 109, 6705 (2005).< / p >
< p id = "swm4-ndp" > < strong > (SWM4-NDP)< / strong > Lamoureux, Harder, Vorobyov, Roux, MacKerell, Chem Phys
Let, 418, 245-249 (2006)< / p >
< / div >
< / div >
< footer >
< hr / >
< div role = "contentinfo" >
< p >
© Copyright .
< / 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:'15 May 2015 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 >
