(r1,r2)into r,(p,p2)into p,etc.for simplicity: p(t+8t)=p(t)+tf(t)+xg(t) r(t+6t)=r(t)8tp(t+t)/m choosing A such that:0=x(r(t +6t)) (10a) p(t+t)=p(t+号6t)+6tf(t+6t)+μg(t+6t) choosing u such that:0=(r(t+6t),p(t +6t)) (10b) Step(10a)may be implemented by defining unconstrained variables (t+26t)=p(t)+26tf(t),(t+6t)=r(t)+6tpt+2t)/m then solving the nonlinear equation for A X(t+6t)=x((t+t)+λ6tg(t)/m)=0 and substituting back p(t+t)=p(t+t)+xg(t);r(t+6t)=(t+6t)+6tAg(t)/m Step(10b)may be handled by defining p(t+t)=p(t+号t)+号tf(t+6t) solving the equation for the second Lagrange multiplier u X(t+6t)=(r(t+6t),(t+t)+g(t+6t)）=0 (which is actually linear,since x(r,p)=-g(r).p/m)and substituting back p(t+6t)=(6t)+ug(t+t) In pseudo-code this scheme may be written do step 1,nstep pp (dt/2)*f r =r dt*p/m lambda_g shake(r) pp lambda g r=r dt*lambda_g/m f force(r) p =p (dt/2)*f mu_g rattle(r,p) p =p mug enddo The routine called shake here calculates the constraint forces Ag;necessary to ensure that the end-of-step positions ri satisfy Eq.(9a).For a system of many constraints,this calculation is usually performed in an iterative fashion,so as to satisfy each constraint in turn until convergence;the original SHAKE algorithm was framed in this way.These constraint forces are incorporated into both the end-of-step positions and the mid-step mo- menta.The routine called rattle calculates a new set of constraint forces ug;to ensure that the end-of-step momenta satisfy Eq.(9b).This also may be carried out iteratively.(r1, r2) into r, (p1 , p2 ) into p, etc. for simplicity: p(t + 1 2 δt) = p(t) + 1 2 δtf(t) + λg(t) r(t + δt) = r(t) + δtp(t + 1 2 δt)/m choosing λ such that: 0 = χ(r(t + δt)) (10a) p(t + δt) = p(t + 1 2 δt) + 1 2 δtf(t + δt) + µg(t + δt) choosing µ such that: 0 = χ˙(r(t + δt), p(t + δt)) (10b) Step (10a) may be implemented by defining unconstrained variables p¯(t + 1 2 δt) = p(t) + 1 2 δtf(t), r¯(t + δt) = r(t) + δtp¯(t + 1 2 δt)/m then solving the nonlinear equation for λ χ(t + δt) = χ ￾ r¯(t + δt) + λδtg(t)/m
= 0 and substituting back p(t + 1 2 δt) = p¯(t + 1 2 δt) + λg(t), r(t + δt) = r¯(t + δt) + δtλg(t)/m Step (10b) may be handled by defining p¯(t + δt) = p(t + 1 2 δt) + 1 2 δtf(t + δt) solving the equation for the second Lagrange multiplier µ χ˙(t + δt) = χ˙ ￾ r(t + δt), p¯(t + δt) + µg(t + δt)
= 0 (which is actually linear, since χ˙(r, p) = −g(r) · p/m) and substituting back p(t + δt) = p¯(δt) + µg(t + δt) In pseudo-code this scheme may be written do step = 1, nstep p = p + (dt/2)*f r = r + dt*p/m lambda_g = shake(r) p = p + lambda_g r = r + dt*lambda_g/m f = force(r) p = p + (dt/2)*f mu_g = rattle(r,p) p = p + mu_g enddo The routine called shake here calculates the constraint forces λgi necessary to ensure that the end-of-step positions ri satisfy Eq. (9a). For a system of many constraints, this calculation is usually performed in an iterative fashion, so as to satisfy each constraint in turn until convergence; the original SHAKE algorithm was framed in this way. These constraint forces are incorporated into both the end-of-step positions and the mid-step mo￾menta. The routine called rattle calculates a new set of constraint forces µgi to ensure that the end-of-step momenta satisfy Eq. (9b). This also may be carried out iteratively. 9