I am simulating an h-BN bilayer using the Tersoff potential for intralayer interactions and LJ/cut for interlayer interactions. However, after 2 ps of simulation, I observe an unintended lateral translation of the layers. is there any way to prevent this motion and ensure the bilayer remains stable in its initial configuration?
LAMMPS (29 Aug 2024)
the data and lammps script are attached for your reference
Thank you for your time and valuable insights. I look forward to your kind suggestions.
You are using a very aggressive cutoff on the plain LJ interaction with 4.5 Angstrom. That will create a lot of noise, which in turn can induce a drift since there can be significant residual forces in the atoms at the cutoff that may not cancel. At the same time, your timestep is uncharacteristically small (you should definitely be able to use 1fs possibly even 2fs with good error conservation). By having to use 100x more timesteps you also generate 100x more noise due to the short cutoff.
Also, you are using both fix momentum and fix recenter and fix momentum with a frequency of 1 which will guarantee that your dynamics will be tainted. Since fix recenter only adjusts the frame of reference it does not help at all. It is meant for cases where you would induce a motion and then want to follow a specific feature that is moving. Fix momentum cannot correct problems with the force field. Its purpose is to remove slowly accumulating moment, e.g. over 1000 or more MD steps.
The fact that you do get a drift is a sign that your simulation or force field settings are not very good. Since you remove the initial momentum when using the velocity command, there is no reason that your system should accumulate a drift as quickly as you report. Suppressing this drift is not going to fix that problem, only hide it. You still have a simulation that produces data you cannot trust, even if it looks like you want it. Such a simulation is not really science but more of a computer animation (and a costly one).
That said, the way to enforce that groups of atoms remain close to their initial position are either fix spring or fix spring/self. But same as fix momentum, these affect the dynamics of your simulation, so you have to first improve your input that you need only very little residual drift, if any, that can then be addressed with very minimal interference, i.e. a small force constant of the spring.
Thank you so much, Dear @akohlmey, for your prompt response. I have implemented all the suggested changes in my LAMMPS script, but the issue persists. My ultimate goal is to create multiple layers to calculate thermal conductivity. I would greatly appreciate any further guidance you can provide. I want to reproduce the result of this paper (Phonon Thermal Transport across Multilayer Graphene/Hexagonal Boron Nitride van der Waals Heterostructures)
Given how fast you are posting, all that you have done is that you have just blindly changed your input according to my comments, but you have not made any attempt to understand where the drift is coming from. I cannot tell you from remote (and I don’t care enough about the kind of research you are doing, my focus is elsewhere). That makes it your job as a researcher to figure out and understand your simulation. For example, it could be that you have stacked the layers incorrectly and the layer-to-layer preferred geometry (for the potential you are using) is different, i.e. shifted. Also, LAMMPS has specific potentials that can handle the interlayer interactions more accurately in the INTERLAYER package. Perhaps, those are needed. This require careful attention to detail, and careful testing, though.
At any rate, from here on this is no longer a LAMMPS issue but a research issue. I have discussed the various options and it is up to you to try and figure out which is the physically best way to model your system. That is not always something that can be sorted out in a few minutes.
