I have been trying to understand the algorithms updating frozen / partially frozen atoms, virtual sites, and shell particles. These are handled differently in different integrator schemes, which (in the best case) is confusing and prone to the introduction of errors, and likely already yields slight differences between the integrator schemes. This issue has several goals: - [ ] Establish and document the protocol the update of these particles should follow - [ ] Verify whether the current implementation is correct - [ ] Understand whether there is potential for simplification or unification across the different integration algorithms. The three types of particles have characteristics that require special treatment during updates: * **Frozen particles** are immobile in some or all dimensions. Their position should stay invariant, their velocity should be zero. The inverse mass and inverse mass per dimension of completely frozen atoms is set to zero (infinitely heavy). The inverse mass per dimension of partially-frozen atoms is set to zero, inverse mass of the atom is 1/(mass of atom type). Frozen particles experience forces as any other particle, i.e. forces on them are expected to be non-zero. * **Virtual site** positions are recomputed from the position of other particles (after their propagation). Their velocity is recalculated from the difference of their position between steps. Forces on virtual sites are redistributed to their constructing particles. The forces are expected to be zero at the end of do_force. The inverse mass of virtual sites is set to zero. * **Shell particles** represent polarization. Their position is determined by minimization of the potential energy every step (during force calculation). *I am not familiar with the shell code. Looking over it, I think to see that (but I'd be happy if someone could confirm or correct me!):* * *Forces are non-zero* * *Inverse masses are zero* * *Velocities are never used, but printed to trajectory (or is there a check somewhere to exclude them?)* * *Positions are updated in do_force, printed to trajectory (?)* Treatment in current `do_md`: 1. Before start of the main loop if new runs (input not read from checkpoint), the velocities of frozen dimensions, virtual sites and shell particles are set to 0. 1. During removal of COM motion, some velocity might be reintroduced (except for frozen particles / dimensions which are explicitly exempt) 1. *Beginning of simulation loop* 1. Force calculation * If shell & vsite, vsite positions & velocities are updated. * If shell particles, shell particle positions are updated. 1. If md-vv, do half-step velocity update * md-vv update sets velocities of freeze dimensions, virtual sites and shell particles explicitly to zero. 1. If md-vv, call compute_globals, might include COM motion removal adding some velocity to zero-ed non-frozen velocities. 1. Write trajectory & checkpoint * For frozen atoms / dimensions, we expect the velocity to be zero * For v-sites, we expect the velocity to be equal to the velocity calculated at the last step (zero if first step) 1. Update positions (full step) and velocities (half-step for md-vv, full step for other algorithms) * md-vv, bd and sd explicitly check particle type, set velocities of freeze dimensions, virtual sites and shell particles to zero. * md-vv, bd and sd explicitly check particle type, set the positions of freeze dimensions, virtual sites and shell particles to their previous values (`xp[i][d] = x[i][d]`). * md has three update paths: * Simple SIMD * No Nose-Hoover, no Parrinello-Rahman, no acceleration, single temperature scaling value (one T-group or non-T-coupling step), no partially frozen atoms (SIMD version uses invMass instead of invMassPerDim to reduce cache pressure). * v-site velocities get set to zero before update * Does not check particle type - all atoms treated equally * Simple * No Nose-Hoover, no off-diagonal Parrinello-Rahman, no acceleration * v-site velocities get set to zero before update * Does not check particle type - all atoms treated equally * General * Only used if Nose-Hoover, off-diagonal Parrinello-Rahman, or acceleration groups are used * Does not zero v-site velocities * Does not check particle type - all atoms treated equally * As md does leap-frog integration, the basic operations for the simple and general code paths is (`i`: atom, `d`: dimension) ``` v[i][d] += dt*f[i][d]*invMassPerDim[i][d] x[i][d] += dt*v[i][d] ``` while SIMD does not differentiate by dimension ``` v[i] += dt*f[i]*invMass[i] x[i] += dt*v[i] ``` * Frozen atoms will not change velocity thanks to the inverse mass / inverse mass per dimension being zero and SIMD being disallowed for partially frozen atoms. Non-zero velocities are unlikely to be introduced anywhere, as COM motion removal is explicitly excluding frozen dimensions → positions should not change. * Virtual sites will not change velocity since their forces and inverse masses are zero. In the simple and simple SIMD implementation, v-site velocities are set to zero before the update. In the general implementation, the v-site velocities are not zero-ed, and will propagate the positions. * The shell particle velocities are not zero-ed in either leap-frog implementation. The inverse masses are zero (?), so the velocities will not change, but velocity introduced via COM motion removal might slightly move the particles. 1. Constraining: This might be a problem when using lincs or shake (which uses invMass) and partially frozen atoms. The constraining might displace the atoms along the frozen dimensions. Fully frozen atoms don't suffer that problem due to the zero inverse mass. 1. Finish update: If we have partially frozen atoms and constraints, we avoid copying back the updated positions (`xp`) to the position vector in the state object. In all other cases, we copy back all positions. *What about the velocities? Velocities are updated in lincs / shake as well, wouldn't we expect to have non-zero velocities in the frozen dimensions of partially frozen, constrained atoms?* 1. Virtual site recalculation: This calculates the virtual sites from the new positions, and the velocities from the displacement between the old and the new virtual site position. This only works if the positions were not changed during the update. *If we moved the virtual sites in the general leap-frog implementation, the velocities might be off.* 1. Compute globals / remove COM motion, might subtract COM motion from vsite / shell velocities. I see the following problems in the current implementation: * ~~If we use the general implementation (Nose-Hoover / off-diagonal Parrinello-Rahman, acceleration), the virtual site velocities are not zeroed (and they are not expected to be zero except in the first step!). This will not result in wrong virtual site positions (since they are recalculated from non-virtual atoms), but in wrong v-site velocities (since the algorithm expects unchanged positions since the last call).~~ → #3866, !979 * ~~md-vv will always print zero / almost zero velocities for the v-sites (due to the velocity half-step before trajectory printing), while the other algorithms will print the v-site velocities calculated at the last step. Md-vv should probably not zero the v-site velocities in the first half step.~~ → #3866, !979 * ~~Virtual sites and shell particles get changed by COM motion removal - this is likely to be a very small contribution in healthy simulations, but is it correct? If it is correct, md-vv needs to be investigated exactly, as it is zero-ing the velocities and removing COM at different places compared to the other algorithms. If it is not correct, the COM removal should expand its particle check to include v-sites and shell particles.~~ → #3866, !979, !980 * ~~Partially-frozen atoms might have non-zero velocities due to constraining, unless I've missed something.~~ → No, the frozen-dimension velocities are zero-ed in the constraining implementation. The frozen-dimension positions, as mentioned above, are not copied back, so this is all as expected. However, the initial constraining might move partially frozen atom *positions*. * ~~Shell particle velocities will accumulate the COM motion removal velocity changes, but be unused otherwise.~~ → !980 I think that otherwise, the frozen / partially frozen atoms and shell particles implementation looks ok. Since the leap-frog update does not care about shell particles, it seems that it wouldn't be necessary to check for it in md-vv / sd / bd either, somewhat unifying the code paths.