Body Rotation in a Molecular Dynamics Simulation

Question

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.

Answer 1

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).