I am playing a bit with few-layer twisted graphene, and a few more interlayer potentials than I normally use (and a much newer LAMMPS as well).
For twisted bilayers of a large size (twist angle 1/3 degree, 59402 x2 carbon atoms) for various potentials (e.g., REBO+ILP, LCBOP+KC, …) the graphene sheets themselves become unstable (at zero temperature, i.e., a straight minimization). This is clearly due to the normals tilting to the xy-plane, giving a long range interaction in the graphene layers, which then leads to more bending of the layers, and more interaction, thus running away…
In past simulations I have often fixed the z-coordinate, or worked for smaller systems, but there are some interesting questions here! I could of course suppress this by using the KC/z interaction, but I was wondering if anyone has a deeper understanding of such instabilities or is a aware of a discussion in the literature.
Any input appreciated.
This is more a question for a discussion with other researchers doing research using interlayer potentials. You may try to contact them through contact information provided in the LAMMPS manual or the corresponding source code or the corresponding publications. While LAMMPS does provide the potentials, any question involving those becomes a LAMMPS topic only if it is about the specifics of the implementation or inconsistencies of the implementation and the documentation or any of the referenced publications.
MatSci has the “ScienceTalk” category for more open discussions.