doc/src/fix_langevin_drude.txt

   1 <script type="text/javascript"
   2   src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
   3 </script>
   4 <script type="text/x-mathjax-config">
   5   MathJax.Hub.Config({ TeX: { equationNumbers: {autoNumber: "AMS"} } });
   6 </script>
   7
   8 "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
   9
  10 :link(lws,http://lammps.sandia.gov)
  11 :link(ld,Manual.html)
  12 :link(lc,Section_commands.html#comm)
  13
  14 :line
  15
  16 fix langevin/drude command :h3
  17
  18 [Syntax:]
  19
  20 fix ID group-ID langevin/drude Tcom damp_com seed_com Tdrude damp_drude seed_drude keyword values ... :pre
  21
  22 ID, group-ID are documented in "fix"_fix.html command :ulb,l
  23 langevin/drude = style name of this fix command :l
  24 Tcom = desired temperature of the centers of mass (temperature units) :l
  25 damp_com = damping parameter for the thermostat on centers of mass (time units) :l
  26 seed_com = random number seed to use for white noise of the thermostat on centers of mass (positive integer) :l
  27 Tdrude = desired temperature of the Drude oscillators (temperature units) :l
  28 damp_drude = damping parameter for the thermostat on Drude oscillators (time units) :l
  29 seed_drude = random number seed to use for white noise of the thermostat on Drude oscillators (positive integer) :l
  30 zero or more keyword/value pairs may be appended :l
  31 keyword = {zero} :l
  32   {zero} value = {no} or {yes}
  33     {no} = do not set total random force on centers of mass to zero
  34     {yes} = set total random force on centers of mass to zero :pre
  35 :ule
  36
  37 [Examples:]
  38
  39 fix 3 all langevin/drude 300.0 100.0 19377 1.0 20.0 83451
  40 fix 1 all langevin/drude 298.15 100.0 19377 5.0 10.0 83451 zero yes :pre
  41
  42 [Description:]
  43
  44 Apply two Langevin thermostats as described in "(Jiang)"_#Jiang for
  45 thermalizing the reduced degrees of freedom of Drude oscillators.
  46 This link describes how to use the "thermalized Drude oscillator
  47 model"_tutorial_drude.html in LAMMPS and polarizable models in LAMMPS
  48 are discussed in "this Section"_Section_howto.html#howto_25.
  49
  50 Drude oscillators are a way to simulate polarizables atoms, by
  51 splitting them into a core and a Drude particle bound by a harmonic
  52 bond.  The thermalization works by transforming the particles degrees
  53 of freedom by these equations.  In these equations upper case denotes
  54 atomic or center of mass values and lower case denotes Drude particle
  55 or dipole values. Primes denote the transformed (reduced) values,
  56 while bare letters denote the original values.
  57
  58 Velocities:
  59 \begin\{equation\} V' = \frac \{M\, V + m\, v\} \{M'\} \end\{equation\}
  60 \begin\{equation\} v' = v - V \end\{equation\}
  61 Masses:
  62 \begin\{equation\} M' = M + m \end\{equation\}
  63 \begin\{equation\} m' = \frac \{M\, m \} \{M'\} \end\{equation\}
  64 The Langevin forces are computed as
  65 \begin\{equation\} F' = - \frac \{M'\} \{\mathtt\{damp\_com\}\}\, V' + F_r' \end\{equation\}
  66 \begin\{equation\} f' = - \frac \{m'\} \{\mathtt\{damp\_drude\}\}\, v' + f_r' \end\{equation\}
  67 \(F_r'\) is a random force proportional to
  68 \(\sqrt \{ \frac \{2\, k_B \mathtt\{Tcom\}\, m'\}
  69                  \{\mathrm dt\, \mathtt\{damp\_com\} \}
  70         \} \). :b
  71 \(f_r'\) is a random force proportional to
  72 \(\sqrt \{ \frac \{2\, k_B \mathtt\{Tdrude\}\, m'\}
  73                  \{\mathrm dt\, \mathtt\{damp\_drude\} \}
  74         \} \). :b
  75 Then the real forces acting on the particles are computed from the inverse
  76 transform:
  77 \begin\{equation\} F = \frac M \{M'\}\, F' - f' \end\{equation\}
  78 \begin\{equation\} f = \frac m \{M'\}\, F' + f' \end\{equation\}
  79
  80 This fix also thermostates non-polarizable atoms in the group at
  81 temperature {Tcom}, as if they had a massless Drude partner.  The
  82 Drude particles themselves need not be in the group. The center of
  83 mass and the dipole are thermostated iff the core atom is in the
  84 group.
  85
  86 Note that the thermostat effect of this fix is applied to only the
  87 translational degrees of freedom of the particles, which is an
  88 important consideration if finite-size particles, which have
  89 rotational degrees of freedom, are being thermostated. The
  90 translational degrees of freedom can also have a bias velocity removed
  91 from them before thermostating takes place; see the description below.
  92
  93 NOTE: Like the "fix langevin"_fix_langevin.html command, this fix does
  94 NOT perform time integration. It only modifies forces to effect
  95 thermostating. Thus you must use a separate time integration fix, like
  96 "fix nve"_fix_nve.html or "fix nph"_fix_nh.html to actually update the
  97 velocities and positions of atoms using the modified forces.
  98 Likewise, this fix should not normally be used on atoms that also have
  99 their temperature controlled by another fix - e.g. by "fix
 100 nvt"_fix_nh.html or "fix temp/rescale"_fix_temp_rescale.html commands.
 101
 102 See "this howto section"_Section_howto.html#howto_16 of the manual for
 103 a discussion of different ways to compute temperature and perform
 104 thermostating.
 105
 106 :line
 107
 108 This fix requires each atom know whether it is a Drude particle or
 109 not.  You must therefore use the "fix drude"_fix_drude.html command to
 110 specify the Drude status of each atom type.
 111
 112 NOTE: only the Drude core atoms need to be in the group specified for
 113 this fix. A Drude electron will be transformed together with its cores
 114 even if it is not itself in the group.  It is safe to include Drude
 115 electrons or non-polarizable atoms in the group. The non-polarizable
 116 atoms will simply be thermostatted as if they had a massless Drude
 117 partner (electron).
 118
 119 NOTE: Ghost atoms need to know their velocity for this fix to act
 120 correctly.  You must use the "comm_modify"_comm_modify.html command to
 121 enable this, e.g.
 122
 123 comm_modify vel yes :pre
 124
 125 :line
 126
 127 {Tcom} is the target temperature of the centers of mass, which would
 128 be used to thermostate the non-polarizable atoms.  {Tdrude} is the
 129 (normally low) target temperature of the core-Drude particle pairs
 130 (dipoles).  {Tcom} and {Tdrude} can be specified as an equal-style
 131 "variable"_variable.html.  If the value is a variable, it should be
 132 specified as v_name, where name is the variable name. In this case,
 133 the variable will be evaluated each timestep, and its value used to
 134 determine the target temperature.
 135
 136 Equal-style variables can specify formulas with various mathematical
 137 functions, and include "thermo_style"_thermo_style.html command
 138 keywords for the simulation box parameters and timestep and elapsed
 139 time.  Thus it is easy to specify a time-dependent temperature.
 140
 141 Like other fixes that perform thermostating, this fix can be used with
 142 "compute commands"_compute.html that remove a "bias" from the atom
 143 velocities.  E.g. removing the center-of-mass velocity from a group of
 144 atoms.  This is not done by default, but only if the
 145 "fix_modify"_fix_modify.html command is used to assign a temperature
 146 compute to this fix that includes such a bias term.  See the doc pages
 147 for individual "compute commands"_compute.html to determine which ones
 148 include a bias.  In this case, the thermostat works in the following
 149 manner: bias is removed from each atom, thermostating is performed on
 150 the remaining thermal degrees of freedom, and the bias is added back
 151 in.  NOTE: this feature has not been tested.
 152
 153 Note: The temperature thermostating the core-Drude particle pairs
 154 should be chosen low enough, so as to mimic as closely as possible the
 155 self-consistent minimization. It must however be high enough, so that
 156 the dipoles can follow the local electric field exerted by the
 157 neighbouring atoms. The optimal value probably depends on the
 158 temperature of the centers of mass and on the mass of the Drude
 159 particles.
 160
 161 {damp_com} is the characteristic time for reaching thermal equilibrium
 162 of the centers of mass.  For example, a value of 100.0 means to relax
 163 the temperature of the centers of mass in a timespan of (roughly) 100
 164 time units (tau or fmsec or psec - see the "units"_units.html
 165 command).  {damp_drude} is the characteristic time for reaching
 166 thermal equilibrium of the dipoles. It is typically a few timesteps.
 167
 168 The number {seed_com} and {seed_drude} are positive integers. They set
 169 the seeds of the Marsaglia random number generators used for
 170 generating the random forces on centers of mass and on the
 171 dipoles. Each processor uses the input seed to generate its own unique
 172 seed and its own stream of random numbers.  Thus the dynamics of the
 173 system will not be identical on two runs on different numbers of
 174 processors.
 175
 176 The keyword {zero} can be used to eliminate drift due to the
 177 thermostat on centers of mass. Because the random forces on different
 178 centers of mass are independent, they do not sum exactly to zero.  As
 179 a result, this fix applies a small random force to the entire system,
 180 and the momentum of the total center of mass of the system undergoes a
 181 slow random walk.  If the keyword {zero} is set to {yes}, the total
 182 random force on the centers of mass is set exactly to zero by
 183 subtracting off an equal part of it from each center of mass in the
 184 group. As a result, the total center of mass of a system with zero
 185 initial momentum will not drift over time.
 186
 187 The actual temperatures of cores and Drude particles, in
 188 center-of-mass and relative coordinates, respectively, can be
 189 calculated using the "compute temp/drude"_compute_temp_drude.html
 190 command.
 191
 192 :line
 193
 194 Usage example for rigid bodies in the NPT ensemble:
 195
 196 comm_modify vel yes
 197 fix TEMP all langevin/drude 300. 100. 1256 1. 20. 13977 zero yes
 198 fix NPH ATOMS rigid/nph/small molecule iso 1. 1. 500.
 199 fix NVE DRUDES nve
 200 compute TDRUDE all temp/drude
 201 thermo_style custom step cpu etotal ke pe ebond ecoul elong press vol temp c_TDRUDE\[1\] c_TDRUDE\[2\] :pre
 202
 203 Comments:
 204
 205 Drude particles should not be in the rigid group, otherwise the Drude
 206 oscillators will be frozen and the system will lose its
 207 polarizability. :ulb,l
 208
 209 {zero yes} avoids a drift of the center of mass of
 210 the system, but is a bit slower. :l
 211
 212 Use two different random seeds to avoid unphysical correlations. :l
 213
 214 Temperature is controlled by the fix {langevin/drude}, so the
 215 time-integration fixes do not thermostate.  Don't forget to
 216 time-integrate both cores and Drude particles. :l
 217
 218 Pressure is time-integrated only once by using {nve} for Drude
 219 particles and {nph} for atoms/cores (or vice versa). Do not use {nph}
 220 for both. :l
 221
 222 The temperatures of cores and Drude particles are calculated by
 223 "compute temp/drude"_compute_temp_drude.html :l
 224
 225 Contrary to the alternative thermostating using Nose-Hoover thermostat
 226 fix {npt} and "fix drude/transform"_fix_drude_transform.html, the
 227 {fix_modify} command is not required here, because the fix {nph}
 228 computes the global pressure even if its group is {ATOMS}. This is
 229 what we want. If we thermostated {ATOMS} using {npt}, the pressure
 230 should be the global one, but the temperature should be only that of
 231 the cores. That's why the command {fix_modify} should be called in
 232 that case. :l
 233 :ule
 234
 235
 236 :line
 237
 238 [Restart, fix_modify, output, run start/stop, minimize info:]
 239
 240 No information about this fix is written to "binary restart
 241 files"_restart.html.  Because the state of the random number generator
 242 is not saved in restart files, this means you cannot do "exact"
 243 restarts with this fix, where the simulation continues on the same as
 244 if no restart had taken place.  However, in a statistical sense, a
 245 restarted simulation should produce the same behavior.
 246
 247 The "fix_modify"_fix_modify.html {temp} option is supported by this
 248 fix.  You can use it to assign a temperature "compute"_compute.html
 249 you have defined to this fix which will be used in its thermostating
 250 procedure, as described above. For consistency, the group used by the
 251 compute should include the group of this fix and the Drude particles.
 252
 253 This fix is not invoked during "energy minimization"_minimize.html.
 254
 255 [Restrictions:] none
 256
 257 [Related commands:]
 258
 259 "fix langevin"_fix_langevin.html,
 260 "fix drude"_fix_drude.html,
 261 "fix drude/transform"_fix_drude_transform.html,
 262 "compute temp/drude"_compute_temp_drude.html,
 263 "pair_style thole"_pair_thole.html
 264
 265 [Default:]
 266
 267 The option defaults are zero = no.
 268
 269 :line
 270
 271 :link(Jiang)
 272 [(Jiang)] Jiang, Hardy, Phillips, MacKerell, Schulten, and Roux, J
 273 Phys Chem Lett, 2, 87-92 (2011).