Body Rotation in a Molecular Dynamics Simulation

I have been trying to simulate particles trapped in a parabolic potential in two dimensions. Due to the circular symmetry of the system, in the ground state at zero temperature, the particles are seen to arrange themselves in concentric rings. As I increase the temperature gradually, the expected result is to observe fluctuations about the mean position. However, due to the circular symmetry in the system, the particles are seen to undergo a body rotation. Can anyone please suggest a way to circumvent this self-rotation? It is imperative that the total energy remains constant during the simulation.

• It sounds as if angular momentum conservation about the symmetry axis of the paraboloid potential is being violated (rather than energy conservation). And if you have a constant temperature thermostat switched on, you should not expect energy conservation anyway. I doubt that anyone will be able to guess the reason, without knowing more details. Are you using a package, to perform the MD, or your own program? How is the thermostat implemented (it may simply be a source of external overall torque). – user197851 Jul 12 '18 at 20:32
• I have written my own program in FORTRAN. The thermostat implemented is the simple velocity re-scaling thermostat. – Mami Jul 14 '18 at 2:51

It's hard to guess, as I said in my comment, but I'm offering this answer anyway. Most molecular dynamics thermostats do not respect angular momentum conservation. This is because most MD simulations use cubic periodic boundary conditions, in which angular momentum is not conserved (the boundary conditions break spherical symmetry). Therefore, I reckon that your thermostat imparts some angular momentum about the symmetry axis, which causes the system to rotate. Either that or your method for giving the particles an initial set of velocities results in nonzero angular momentum.

If that is the cause of the problem, it may be possible to apply a thermostat which avoids it. There are thermostats which are Galilean invariant, and preserve linear momentum. Some of them also preserve angular momentum. For example the Lowe-Andersen thermostat thermalizes randomly selected pairs of atoms $ij$: $\mathbf{v}_i \rightarrow \mathbf{v}_i+\Delta\mathbf{v}$ and $\mathbf{v}_j \rightarrow \mathbf{v}_j-\Delta\mathbf{v}$, choosing $\Delta\mathbf{v}$ to be parallel to the vector $\mathbf{r}_i-\mathbf{r}_j$. The magnitude $|\Delta\mathbf{v}|$ is chosen from the appropriate probability distribution. Then the change in angular momentum is $m(\mathbf{r}_i-\mathbf{r}_j)\times\Delta\mathbf{v}=0$. For simplicity I'm assuming equal-mass atoms.

Otherwise, my advice would simply be to choose initial velocities randomly, apply a global correction to cancel out the total angular momentum, and allow the system to equilibrate without adjusting the velocities, or at worst, simply scale them all up and down by a common factor to achieve the desired temperature.

[Edit following comment]

To remove the initial angular momentum: here is one scheme, but it is easy to devise other versions of it. I assume all masses are equal, $m$, and that the motion is in the $xy$-plane. For each particle $\mathbf{r}_i=(x_i,y_i)$ calculate the distance from the origin $r_i=|\mathbf{r}_i|=\sqrt{x_i^2+y_i^2}$, and construct the unit tangential vector $\mathbf{t}_i=(-y_i,x_i)/r_i$ (see diagram). The idea is to adjust the velocity of every particle by $\Delta v \mathbf{t}_i$, in the tangential direction. The magnitude of the adjustment, $\Delta v$, is to be determined. By construction, the change in angular momentum about the origin will be $m\Delta v r_i$ for each particle, $m\Delta v \sum_i r_i$ in total, and we just need to choose $\Delta v$ to make this cancel the initial angular momentum.

The angular momentum ($z$-component) is given by the vector cross product formula, which we can easily rewrite in terms of the tangential vectors: $$L_z = \left(m \sum_i \mathbf{r}_i \times \mathbf{v}_i\right)_z = m\sum_i ( x_i v_{i,y} - y_i v_{i,x} ) = m\sum_i r_i ( v_{i,x} t_{i,x} + v_{i,y}t_{i,y}) = m\sum_i r_i ( \mathbf{v}_i\cdot\mathbf{t}_i )$$ From this we can calculate the necessary $\Delta v$ to cancel it out $$m\Delta v \sum_i r_i = -m\sum_i r_i ( \mathbf{v}_i\cdot\mathbf{t}_i ) \qquad\Rightarrow\qquad \Delta v = -\frac{\sum_i r_i ( \mathbf{v}_i\cdot\mathbf{t}_i )}{\sum_i r_i}$$ and then you simply replace each particle velocity by $$\mathbf{v}_i' = \mathbf{v}_i + \Delta v \mathbf{t}_i$$ Do this once, at the start, after you have chosen random initial velocities. Double check that it worked, by printing out the new total angular momentum: it should be zero. Thereafter, provided your thermostat simply entails multiplying all velocities by the same common factor, the angular momentum should be conserved at zero. The only numerical problem that I can see will arise if you have placed a particle exactly at the origin (because then $\mathbf{t}_i$ is undefined). You should avoid doing that (or exclude such a particle from the calculation).

• Thank you for your answer. I have initialised the positions and velocity randomly and have used a velocity re-scaling thermostat the equilibrate the system. Could you please elaborate on the "global correction" that I could use to eliminate the angular momentum? Shouldn't the total angular momentum be conserved as the system is rotationally invariant? – Mami Jul 14 '18 at 2:49
• @Mami I have edited my answer accordingly, hope this helps. – user197851 Jul 14 '18 at 9:17