Why does a bulk crystal get tilted in NPT MD simulations

Dear friends,
I am carrying molecular dynamics of a Lennard-Jones
fcc solid, in 3D. While equilibrating the system in the NPT ensemble,
I often find that the crystal gets tilted (or rotated) for some reason. Can
anyone provide an explanation for this.


What is your fix npt command? And your

system boundary conditions? Are you running

with an orthogonal or triclinic box?

And what do you mean by “gets tilted” and

how do you know that is happening? Visual,

numbers in a file, etc?


Dear Steve,
I use periodic boundary conditions and
an orthogonal simulation box. The same tilt of the crystal
also happens in an NVT simulation and also for
cubic simulation boxes. I also do not know whether
to actually call it a tilt because it is more as if a wave is traversing
across the crystal. I do not understand the origin
of this.


If the box size changes suddenly it can induce

a shock-wave in a solid. “Suddenly” is a subjective

term - it depends on how far your current pressure

is from the target, how rapidly your fix npt params

allow the volume to change, etc, etc.


I believe this is an issue with your initial setup (geometry/interactions).
If you are running with fix nvt your box geometry cannot change.
What is that you considered to be a "tilt" if the box remains unchanged?
Have you tried minimizing the initial structure before running dynamics?
Bets are this strange behavior is 99.9% your fault. Looking to prove the
opposite? Then one or two clear pics and a small reproducible input script
"may" prompt someone to give you further free assistance.


I have attached a snapshot of the configuration with the (100) orientation

of the fcc crystal along the z-direction. This was from a NVT simulation but

also for NPT simulation I sometimes get similar configurations i.e. as if a wave is
moving through the crystal. However, this does not happen all the time and

I get a crystal with no wavy structure. The snapshot corresponds to an inverse-256th

power potential (i.e. almost hard-sphere like) at the coexistence temperature.


config.pdf (250 KB)

I can explain what I see. I just don’t know how realistic this is for your particular case.
This wave pattern appears sometimes in hexagonal crystals under high pressure/temperature. You can probably find simulations
where people try to convert hexagonal carbon into diamond and find this kind of behavior. The deal is that the intra-layer
interactions tend to be much stronger than inter-layer ones and combined with the use of PBCs, this leads to a collective atomic motion where the sheets oscillate. Think about the case of one atomic layer under thermal noise and BCs (a drum sounds familiar?)
I just don’t know ( don’t care either, sorry) how much of a real thing or of and artifact your results are.
Hope this comment helps you in the way to finding your final answer,

Here is a link: http://www.nature.com/articles/srep05930 that could help you further.

"The snapshot corresponds to an inverse-256th
power potential (i.e. almost hard-sphere like) at the coexistence temperature."

You said you were running LJ. That is not LJ. You are going to have to
use an awfully small timestep. Even then, the model may exhibit
pathological behavior that is just a consequence of the extreme
collisional stiffness. More generally, depending on how a crystal is
initialized, it can exhibit all kinds of interesting elastic waves.
The best way to quickly remove these is with a Langevin thermostat.
With Nose-Hoover thermostat, they can sometimes persist for a long

Dear all,
thanks for all the helpful suggestions. It has greatly
expanded my knowledge.