git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13182 f3b2605a-c512-4ea7-a41b-209d697bcdaa

doc/fix_shake.html

@@ -15,11 +15,11 @@
 </H3>
 <P><B>Syntax:</B>
 </P>
-<PRE>fix ID group-ID shake tol iter N constraint values ... keyword value ... 
+<PRE>fix ID group-ID <I>shake/rattle</I> tol iter N constraint values ... keyword value ... 
 </PRE>
 <UL><LI>ID, group-ID are documented in <A HREF = "fix.html">fix</A> command 
 
-<LI>shake = style name of this fix command 
+<LI>shake or rattle = style name of this fix command 
 
 <LI>tol = accuracy tolerance of SHAKE solution 
 
@@ -51,11 +51,32 @@
 fix 1 sub shake 0.0001 20 10 t 5 6 m 1.0 a 31
 fix 1 sub shake 0.0001 20 10 t 5 6 m 1.0 a 31 mol myMol 
 </PRE>
+<PRE>fix 1 sub rattle 0.0001 20 10 t 5 6 m 1.0 a 31
+fix 1 sub rattle 0.0001 20 10 t 5 6 m 1.0 a 31 mol myMol 
+</PRE>
 <P><B>Description:</B>
 </P>
 <P>Apply bond and angle constraints to specified bonds and angles in the
 simulation.  This typically enables a longer timestep.
 </P>
+<P><B>SHAKE vs. RATTLE:</B>
+</P>
+<P>The SHAKE algorithm was invented for schemes such as Verlet, where only the coordinates 
+are integrated and the velocities are approximated as finite differences to the trajectories (<A HREF = "#Ryckaert">Ryckaert et al. (1977)</A>). 
+If the velocities are integrated explicitly, for example with velocity Verlet, 
+a second set of constraining forces is required in order to eliminate velocity components along the bonds (<A HREF = "#Andersen">Andersen (1983)</A>). 
+</P>
+<P>In order to formulate the individual constraints for SHAKE and RATTLE, let us focus on a single molecule whose 
+bonds are to be constrained. Let image(Eqs/fix_rattle_ri.jpg) and image(Eqs/fix_rattle_vi.jpg) be the position and velocity of atom <I>i</I> at time <I>n</I>, for <I>i</I>=1,...,<I>N</I>, where <I>N</I> is the number of sites of our reference molecule. The distance vector between sites <I>i</I> and <I>j</I> is given by image(Eqs/fix_rattle_rij.jpg). The constraints can then be formulated as
+</P>
+<CENTER><IMG SRC = "Eqs/fix_rattle_constraints.jpg">
+</CENTER>
+<P>The SHAKE algorithm satisfies the first condition, i.e. the sites at time <I>n+1</I> will have the desired separations image(Eqs/fix_rattle_dij.jpg) immediately after the coordinates are integrated. 
+If we also enforce the second condition, the velocity components along the bonds will vanish. 
+RATTLE satisfies both conditions and fix rattle is implemented in a way that it uses fix shake for dealing with the coordinate constraints. Therefore all the information about SHAKE below is also relevant for RATTLE.
+</P>
+<P><B>SHAKE:</B>
+</P>
 <P>Each timestep the specified bonds and angles are reset to their
 equilibrium lengths and angular values via the well-known SHAKE
 algorithm.  This is done by applying an additional constraint force so
@@ -144,6 +165,24 @@ more instructions on how to use the accelerated styles effectively.
 </P>
 <HR>
 
+<P><B>RATTLE:</B>
+</P>
+<P>The velocity constraints lead to a linear system of equations which
+can be solved analytically. The implementation of the algorithm in LAMMPS
+closely follows <A HREF = "#Andersen">Andersen (1983)</A>. 
+</P>
+<P>IMPORTANT NOTE: This command modifies forces and velocities and 
+it has to come after all other integration fixes. If you define other
+fixes that modify velocities or forces after fix rattle operates, 
+then this fix will not take them into account and the time integration
+will typically not satisfy the RATTLE constraints. You can check whether
+the constraints work correctly by setting the value of RATTLE_DEBUG in fix_rattle.cpp
+to 1 and recompiling LAMMPS.
+</P>
+<P>IMPORTANT NOTE: RATTLE does not support cuda at the moment.
+</P>
+<HR>
+
 <P><B>Restart, fix_modify, output, run start/stop, minimize info:</B>
 </P>
 <P>No information about this fix is written to <A HREF = "restart.html">binary restart
@@ -173,4 +212,16 @@ of SHAKE parameters and monitoring the energy versus time.
 </P>
 <P><B>Default:</B> none
 </P>
+<HR>
+
+<A NAME = "Ryckaert"></A>
+
+<P><B>(Ryckaert)</B> J.-P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, 
+Journal of Computational Physics, 23, 327–341 (1977).
+</P>
+<A NAME = "Andersen"></A>
+
+<P><B>(Andersen)</B> H. Andersen,
+Journal of Computational Physics, 52, 24-34 (1983).
+</P>
 </HTML>

doc/fix_shake.txt

@@ -11,10 +11,10 @@ fix shake/cuda command :h3
 
 [Syntax:]
 
-fix ID group-ID shake tol iter N constraint values ... keyword value ... :pre
+fix ID group-ID {shake/rattle} tol iter N constraint values ... keyword value ... :pre
 
 ID, group-ID are documented in "fix"_fix.html command :ulb,l
-shake = style name of this fix command :l
+shake or rattle = style name of this fix command :l
 tol = accuracy tolerance of SHAKE solution :l
 iter = max # of iterations in each SHAKE solution :l
 N = print SHAKE statistics every this many timesteps (0 = never) :l
@@ -35,12 +35,34 @@ keyword = {mol} :l
 fix 1 sub shake 0.0001 20 10 b 4 19 a 3 5 2
 fix 1 sub shake 0.0001 20 10 t 5 6 m 1.0 a 31
 fix 1 sub shake 0.0001 20 10 t 5 6 m 1.0 a 31 mol myMol :pre
+fix 1 sub rattle 0.0001 20 10 t 5 6 m 1.0 a 31
+fix 1 sub rattle 0.0001 20 10 t 5 6 m 1.0 a 31 mol myMol :pre
 
 [Description:]
 
 Apply bond and angle constraints to specified bonds and angles in the
 simulation.  This typically enables a longer timestep.
 
+
+[SHAKE vs. RATTLE:]
+
+The SHAKE algorithm was invented for schemes such as Verlet, where only the coordinates 
+are integrated and the velocities are approximated as finite differences to the trajectories ("Ryckaert et al. (1977)"_#Ryckaert). 
+If the velocities are integrated explicitly, for example with velocity Verlet, 
+a second set of constraining forces is required in order to eliminate velocity components along the bonds ("Andersen (1983)"_#Andersen). 
+
+In order to formulate the individual constraints for SHAKE and RATTLE, let us focus on a single molecule whose 
+bonds are to be constrained. Let image(Eqs/fix_rattle_ri.jpg) and image(Eqs/fix_rattle_vi.jpg) be the position and velocity of atom {i} at time {n}, for {i}=1,...,{N}, where {N} is the number of sites of our reference molecule. The distance vector between sites {i} and {j} is given by image(Eqs/fix_rattle_rij.jpg). The constraints can then be formulated as
+
+:c,image(Eqs/fix_rattle_constraints.jpg)
+
+The SHAKE algorithm satisfies the first condition, i.e. the sites at time {n+1} will have the desired separations image(Eqs/fix_rattle_dij.jpg) immediately after the coordinates are integrated. 
+If we also enforce the second condition, the velocity components along the bonds will vanish. 
+RATTLE satisfies both conditions and fix rattle is implemented in a way that it uses fix shake for dealing with the coordinate constraints. Therefore all the information about SHAKE below is also relevant for RATTLE.
+
+[SHAKE:]
+
+
 Each timestep the specified bonds and angles are reset to their
 equilibrium lengths and angular values via the well-known SHAKE
 algorithm.  This is done by applying an additional constraint force so
@@ -92,6 +114,7 @@ satisfy the SHAKE constraints.  The solution for this is to make sure
 that fix shake is defined in your input script after any other fixes
 which add or change forces (to atoms that fix shake operates on).
 
+
 :line
 
 The {mol} keyword should be used when other commands, such as "fix
@@ -127,6 +150,27 @@ use the "suffix"_suffix.html command in your input script.
 See "Section_accelerate"_Section_accelerate.html of the manual for
 more instructions on how to use the accelerated styles effectively.
 
+:line
+
+
+[RATTLE:]
+
+
+The velocity constraints lead to a linear system of equations which
+can be solved analytically. The implementation of the algorithm in LAMMPS
+closely follows "Andersen (1983)"_#Andersen. 
+	
+IMPORTANT NOTE: This command modifies forces and velocities and 
+it has to come after all other integration fixes. If you define other
+fixes that modify velocities or forces after fix rattle operates, 
+then this fix will not take them into account and the time integration
+will typically not satisfy the RATTLE constraints. You can check whether
+the constraints work correctly by setting the value of RATTLE_DEBUG in fix_rattle.cpp
+to 1 and recompiling LAMMPS.
+
+IMPORTANT NOTE: RATTLE does not support cuda at the moment.
+
+
 :line
 
 [Restart, fix_modify, output, run start/stop, minimize info:]
@@ -157,3 +201,15 @@ of SHAKE parameters and monitoring the energy versus time.
 [Related commands:] none
 
 [Default:] none
+
+:line 
+ 
+:link(Ryckaert)
+[(Ryckaert)] J.-P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, 
+Journal of Computational Physics, 23, 327–341 (1977).
+ 
+:link(Andersen)
+[(Andersen)] H. Andersen,
+Journal of Computational Physics, 52, 24-34 (1983).
+
+