doc/src/fix_shake.txt

   1 "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
   2
   3 :link(lws,http://lammps.sandia.gov)
   4 :link(ld,Manual.html)
   5 :link(lc,Section_commands.html#comm)
   6
   7 :line
   8
   9 fix shake command :h3
  10 fix rattle command :h3
  11
  12 [Syntax:]
  13
  14 fix ID group-ID style tol iter N constraint values ... keyword value ... :pre
  15
  16 ID, group-ID are documented in "fix"_fix.html command :ulb,l
  17 style = shake or rattle = style name of this fix command :l
  18 tol = accuracy tolerance of SHAKE solution :l
  19 iter = max # of iterations in each SHAKE solution :l
  20 N = print SHAKE statistics every this many timesteps (0 = never) :l
  21 one or more constraint/value pairs are appended :l
  22 constraint = {b} or {a} or {t} or {m} :l
  23   {b} values = one or more bond types
  24   {a} values = one or more angle types
  25   {t} values = one or more atom types
  26   {m} value = one or more mass values :pre
  27 zero or more keyword/value pairs may be appended :l
  28 keyword = {mol} :l
  29   {mol} value = template-ID
  30     template-ID = ID of molecule template specified in a separate "molecule"_molecule.html command :pre
  31 :ule
  32
  33 [Examples:]
  34
  35 fix 1 sub shake 0.0001 20 10 b 4 19 a 3 5 2
  36 fix 1 sub shake 0.0001 20 10 t 5 6 m 1.0 a 31
  37 fix 1 sub shake 0.0001 20 10 t 5 6 m 1.0 a 31 mol myMol
  38 fix 1 sub rattle 0.0001 20 10 t 5 6 m 1.0 a 31
  39 fix 1 sub rattle 0.0001 20 10 t 5 6 m 1.0 a 31 mol myMol :pre
  40
  41 [Description:]
  42
  43 Apply bond and angle constraints to specified bonds and angles in the
  44 simulation by either the SHAKE or RATTLE algorithms.  This typically
  45 enables a longer timestep.
  46
  47 [SHAKE vs RATTLE:]
  48
  49 The SHAKE algorithm was invented for schemes such as standard Verlet
  50 timesteppnig, where only the coordinates are integrated and the
  51 velocities are approximated as finite differences to the trajectories
  52 ("Ryckaert et al. (1977)"_#Ryckaert).  If the velocities are
  53 integrated explicitly, as with velocity Verlet which is what LAMMPS
  54 uses as an integration method, a second set of constraining forces is
  55 required in order to eliminate velocity components along the bonds
  56 ("Andersen (1983)"_#Andersen).
  57
  58 In order to formulate individual constraints for SHAKE and RATTLE,
  59 focus on a single molecule whose bonds are constrained.  Let Ri and Vi
  60 be the position and velocity of atom {i} at time {n}, for
  61 {i}=1,...,{N}, where {N} is the number of sites of our reference
  62 molecule. The distance vector between sites {i} and {j} is given by
  63
  64 :c,image(Eqs/fix_rattle_rij.jpg)
  65
  66 The constraints can then be formulated as
  67
  68 :c,image(Eqs/fix_rattle_constraints.jpg)
  69
  70 The SHAKE algorithm satisfies the first condition, i.e. the sites at
  71 time {n+1} will have the desired separations Dij immediately after the
  72 coordinates are integrated.  If we also enforce the second condition,
  73 the velocity components along the bonds will vanish.  RATTLE satisfies
  74 both conditions.  As implemented in LAMMPS, fix rattle uses fix shake
  75 for satisfying the coordinate constraints. Therefore the settings and
  76 optional keywords are the same for both fixes, and all the information
  77 below about SHAKE is also relevant for RATTLE.
  78
  79 [SHAKE:]
  80
  81 Each timestep the specified bonds and angles are reset to their
  82 equilibrium lengths and angular values via the SHAKE algorithm
  83 ("Ryckaert et al. (1977)"_#Ryckaert).  This is done by applying an
  84 additional constraint force so that the new positions preserve the
  85 desired atom separations.  The equations for the additional force are
  86 solved via an iterative method that typically converges to an accurate
  87 solution in a few iterations.  The desired tolerance (e.g. 1.0e-4 = 1
  88 part in 10000) and maximum # of iterations are specified as arguments.
  89 Setting the N argument will print statistics to the screen and log
  90 file about regarding the lengths of bonds and angles that are being
  91 constrained.  Small delta values mean SHAKE is doing a good job.
  92
  93 In LAMMPS, only small clusters of atoms can be constrained.  This is
  94 so the constraint calculation for a cluster can be performed by a
  95 single processor, to enable good parallel performance.  A cluster is
  96 defined as a central atom connected to others in the cluster by
  97 constrained bonds.  LAMMPS allows for the following kinds of clusters
  98 to be constrained: one central atom bonded to 1 or 2 or 3 atoms, or
  99 one central atom bonded to 2 others and the angle between the 3 atoms
 100 also constrained.  This means water molecules or CH2 or CH3 groups may
 101 be constrained, but not all the C-C backbone bonds of a long polymer
 102 chain.
 103
 104 The {b} constraint lists bond types that will be constrained.  The {t}
 105 constraint lists atom types.  All bonds connected to an atom of the
 106 specified type will be constrained.  The {m} constraint lists atom
 107 masses.  All bonds connected to atoms of the specified masses will be
 108 constrained (within a fudge factor of MASSDELTA specified in
 109 fix_shake.cpp).  The {a} constraint lists angle types.  If both bonds
 110 in the angle are constrained then the angle will also be constrained
 111 if its type is in the list.
 112
 113 For all constraints, a particular bond is only constrained if both
 114 atoms in the bond are in the group specified with the SHAKE fix.
 115
 116 The degrees-of-freedom removed by SHAKE bonds and angles are accounted
 117 for in temperature and pressure computations.  Similarly, the SHAKE
 118 contribution to the pressure of the system (virial) is also accounted
 119 for.
 120
 121 NOTE: This command works by using the current forces on atoms to
 122 caculate an additional constraint force which when added will leave
 123 the atoms in positions that satisfy the SHAKE constraints (e.g. bond
 124 length) after the next time integration step.  If you define fixes
 125 (e.g. "fix efield"_fix_efield.html) that add additional force to the
 126 atoms after fix shake operates, then this fix will not take them into
 127 account and the time integration will typically not satisfy the SHAKE
 128 constraints.  The solution for this is to make sure that fix shake is
 129 defined in your input script after any other fixes which add or change
 130 forces (to atoms that fix shake operates on).
 131
 132 :line
 133
 134 The {mol} keyword should be used when other commands, such as "fix
 135 deposit"_fix_deposit.html or "fix pour"_fix_pour.html, add molecules
 136 on-the-fly during a simulation, and you wish to contrain the new
 137 molecules via SHAKE.  You specify a {template-ID} previously defined
 138 using the "molecule"_molecule.html command, which reads a file that
 139 defines the molecule.  You must use the same {template-ID} that the
 140 command adding molecules uses.  The coordinates, atom types, special
 141 bond restrictions, and SHAKE info can be specified in the molecule
 142 file.  See the "molecule"_molecule.html command for details.  The only
 143 settings required to be in this file (by this command) are the SHAKE
 144 info of atoms in the molecule.
 145
 146 :line
 147
 148 Styles with a suffix are functionally the same as the corresponding
 149 style without the suffix.  They have been optimized to run faster,
 150 depending on your available hardware, as discussed in
 151 "Section 5"_Section_accelerate.html of the manual.  The
 152 accelerated styles take the same arguments and should produce the same
 153 results, except for round-off and precision issues.
 154
 155 These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
 156 USER-OMP and OPT packages, respectively.  They are only enabled if
 157 LAMMPS was built with those packages.  See the "Making
 158 LAMMPS"_Section_start.html#start_3 section for more info.
 159
 160 You can specify the accelerated styles explicitly in your input script
 161 by including their suffix, or you can use the "-suffix command-line
 162 switch"_Section_start.html#start_7 when you invoke LAMMPS, or you can
 163 use the "suffix"_suffix.html command in your input script.
 164
 165 See "Section 5"_Section_accelerate.html of the manual for
 166 more instructions on how to use the accelerated styles effectively.
 167
 168 :line
 169
 170 [RATTLE:]
 171
 172 The velocity constraints lead to a linear system of equations which
 173 can be solved analytically.  The implementation of the algorithm in
 174 LAMMPS closely follows ("Andersen (1983)"_#Andersen).
 175
 176 NOTE: The fix rattle command modifies forces and velocities and thus
 177 should be defined after all other integration fixes in your input
 178 script.  If you define other fixes that modify velocities or forces
 179 after fix rattle operates, then fix rattle will not take them into
 180 account and the overall time integration will typically not satisfy
 181 the RATTLE constraints.  You can check whether the constraints work
 182 correctly by setting the value of RATTLE_DEBUG in src/fix_rattle.cpp
 183 to 1 and recompiling LAMMPS.
 184
 185 :line
 186
 187 [Restart, fix_modify, output, run start/stop, minimize info:]
 188
 189 No information about these fixes is written to "binary restart
 190 files"_restart.html.  None of the "fix_modify"_fix_modify.html options
 191 are relevant to these fixes.  No global or per-atom quantities are
 192 stored by these fixes for access by various "output
 193 commands"_Section_howto.html#howto_15.  No parameter of these fixes
 194 can be used with the {start/stop} keywords of the "run"_run.html
 195 command.  These fixes are not invoked during "energy
 196 minimization"_minimize.html.
 197
 198 [Restrictions:]
 199
 200 These fixes are part of the RIGID package.  They are only enabled if
 201 LAMMPS was built with that package.  See the "Making
 202 LAMMPS"_Section_start.html#start_3 section for more info.
 203
 204 For computational efficiency, there can only be one shake or rattle
 205 fix defined in a simulation.
 206
 207 If you use a tolerance that is too large or a max-iteration count that
 208 is too small, the constraints will not be enforced very strongly,
 209 which can lead to poor energy conservation.  You can test for this in
 210 your system by running a constant NVE simulation with a particular set
 211 of SHAKE parameters and monitoring the energy versus time.
 212
 213 SHAKE or RATTLE should not be used to contrain an angle at 180 degrees
 214 (e.g. linear CO2 molecule).  This causes numeric difficulties.
 215
 216 [Related commands:] none
 217
 218 [Default:] none
 219
 220 :line
 221
 222 :link(Ryckaert)
 223 [(Ryckaert)] J.-P. Ryckaert, G. Ciccotti and H. J. C. Berendsen,
 224 J of Comp Phys, 23, 327-341 (1977).
 225
 226 :link(Andersen)
 227 [(Andersen)] H. Andersen, J of Comp Phys, 52, 24-34 (1983).