I am trying to use Lammps to calculate the torque on a twisted bilayer graphene system. My rotation center is at the corner of the triclinic cell (at the AA stacking, in case you are familiar with this terminology).
This compute calculates the three components of the torque vector for each chunk, due to the forces on the individual atoms in the chunk around the center-of-mass of the chunk. The calculation includes all effects due to atoms passing through periodic boundaries.
I’m mostly confused about the periodic boundary part. How can a center of mass be defined for a periodic system? Is it just using the center of mass of the supercell?
Is it possible to define another center than the center of mass for the torque in LAMMPS, like the actual rotation center of my system? I am getting a different torque for the bottom and top layer so I suspect the center of mass has something to do with it, so I’d like to make some checks where I control the choice of the reference for the torque.
Unless I am missing something completely and the torque the way it is calculated in lammps does not allow me to probe the physics I am interested in, i.e. the desire for twisted systems to rotate back towards zero twist alignment.
Looking through the code, “all effects due to atoms passing through periodic boundaries” means that particles’ unwrapped positions are used instead of wrapped positions. That is:
the COM is calculated by summing over all particles’ unwrapped positions
each particle’s torque is calculated as a cross-product of force with the particle’s unwrapped displacement-from-COM
Interestingly, particles are chunked according to their wrapped positions. I am not sure if this is consistent. But at least it’s precisely defined!
If the torque about a reference point is \tau_0, the torque \tau_{\Delta \mathbf{r}} about some point displaced \Delta \mathbf{r} from the reference point is the original torque minus the displacement crossed with the total force: \tau_{\Delta \mathbf{r}} = \tau_0 - \Delta \mathbf{r} \times \sum \mathbf{F}. This is left as a fact from Wikipedia an exercise for the reader.
It should be straightforward (if tedious) to define per-atom torque components using wrapped instead of unwrapped lever arms as per-atom variables, and then average those over chunks as v_ variable inputs to fix ave/chunk.
Thank you. Let me read through your answer in detail.
The reason I am wondering if Lammps allows me to calculate what I want is that the forces after relaxation by definition will be nearly zero (and by extension, the torque if it indeed uses the previous forces to calculate the torque).
But in my case, i’d like to have the forces acting on the atoms in such a way that a twisted system rotates back to zero twist angle. But due to the constraints imposed by the chosen periodic cell, it cannot rotate back. But maybe it is these constraints that counter the forces that I am interested in, thus leading to the global zero force on each atom reported by LAMMPS at the end of the minimization run?
When talking about the unwrapped positions, if I take the relaxed positions after minimization as input for a torque calculation where I only allow for ONE energy minimization step, the atoms won’t move across any periodic boundary, thus the input positions will simply count as both the wrapped and unwrapped positions, no?
Basing this on
Unwrapped means that if the atom has passed through a periodic boundary one or more times, the value is printed for what the coordinate would be if it had not been wrapped back into the periodic box. Note that these coordinates may thus be far outside the box size stored with the snapshot.
With my theory hat on – if any directions perpendicular to a rotation axis have non-zero projection into a periodic direction, I fail to see how rotation about that axis (through any point) can be a valid continuous deformation of the periodic system as a whole.
If you are instead rotating every image about its own pivot point you must be inducing tremendous discontinuities along the former periodic boundaries, which would dominate any genuine torques in your system.
Ok, so it does make sense that my resultant torques calculated using the Lammps command seem to be very small in my current calculations.
I guess I must think a bit deeper about what observable I can look at to quantify the tendency one twisted bilayer system has to rotate back into zero-twist alignment. Of course, it is obvious from the total energy curve vs twist angle showing a downwards slope toward zero twist angle (where i can take dE/dtheta to get the torque constant), but I was hoping to get the info in a different way too.
Thanks for your comments.
I checked, ixiy and iz are all zero, as I am starting from relaxed positions where the convergence criterion on the forces is already satisfied to begin with.
