I was examining a system with a problematic equilibration (more details given at the end, but probably unnecessary to answer the question), and came upon what is presumably a pretty basic conceptual question that I unfortunately didn’t find an answer to on my own.

I am using LAMMPS to simulate a NPT ensemble for an ionic crystal with some small organic molecules, using a force field made up of two different pairwise potentials (lj and Buckingham), and some molecular terms (bonds, angles, tetrahedrals). When plotting the pairwise energy, the molecular energy, and the total potential energy (pair+molec+long range interactions), I noticed a trend that I don’t understand well:  graph of the different energies

Here is what I observe:

The system reaches the target temperature (300K) almost instantly at the scale of the graph, and then the geometry take a (very) long time to relax. While doing so, the pairwise energy goes down, but the molecular energy compensates it pretty much exactly to keep the potential energy constant.

I do not understand this well. In particular, why is the molecular energy going up? Upon a quick examination of the system, it seems like the bond lengths are simply vibrating but not drifting higher or lower. What is the source of this energy change? And why does it have to compensate the change in pairwise energy? I don’t see a good reason why the potential energy would be conserved in a NPT. On a side note, kinetic energy and total system energy are also conserved throughout  

Here is what I would have qualitatively expected to observe:

  • Pairwise energy does go down as the relaxation/equilibrations happens: this would be because the atoms are adjusting and getting into their “optimal average positions” energy wise under the action of the interactions, minimizing this pairwise “binding energy” (similar to how a dft simulation looks for an energy minimum. Of course, in a statistical ensemble, there is no such thing, but I would still reasonably expect the pairwise energy to decrease as the equilibration proceeds)
  • Molecular energy remains pretty much the same, granted that no significant change is going on in the average molecule geometry
  • Potential energy decreases with the pairwise one
  • Kinetic energy remains more or less the same once the target temperature is reached, implying a decrease of the total energy of the system (which is alright since it’s not NVE, we are not expecting energy to be a conserved quantity)

Basically, I just don’t understand why the molecular energy is increasing, and why it’s precisely compensating the pairwise energy. Could anyone point me in the right direction? Many thanks!

  In case it helps, my system is $\ce{MAPbI3}$, an ionic crystal of $\ce{Pb2+}$, $\ce{I-}$, and methyllamonium cation. The dynamics of the cation (notably rotations) are relevant to the overall relaxation to a ~ cubic system. The force field used is taken from the literature and has been used by several authors. I am observing a slow relaxation of the crystal lattice over ~1.5e6 timesteps of 0.5fs.  

EDIT: I am using a harmonic model for the bond which has the form $E_{bond}=K(r-r_0)^2$, with $r_0$ the parametrized bond length, and $r$ the actual bond length. If $r_0 \approx \langle r\rangle$ (which should be the case since it is the point of parameterizing the bond after all), then we have:

$$\langle E_{bond} \rangle \approx K\langle(r-\langle r \rangle)^2\rangle \approx K\textrm{var}(r),\tag{1}$$

and the increasing bond energy can be explained by an increasing variance in bond length during equilibration, even if the average bond length remains the same. I still do not understand how this is related to the decrease in pairwise energy and why the potential energy ends up being conserved

  • $\begingroup$ Do you observe any significant structural difference before and after the slow relaxation? $\endgroup$
    – wzkchem5
    Commented Sep 3, 2021 at 15:25
  • $\begingroup$ yes indeed. The starting point is a minimized geometry at 0K, so it undergoes serious changes when thermalized a 300K. I could include the evolution of the lattice constants if it helps. Those changes happen on the long timescale mentioned, along with the energy changes in the graph $\endgroup$ Commented Sep 3, 2021 at 15:46
  • $\begingroup$ While you will see any comments on your post, you will need to ping other users to notify them of your comment. For example, I'll ping @wzkchem5 $\endgroup$
    – Tyberius
    Commented Dec 1, 2021 at 19:34
  • $\begingroup$ @Tyberius I actually saw the OP's reply, but I had nothing else to comment, so I didn't respond $\endgroup$
    – wzkchem5
    Commented Dec 1, 2021 at 20:14
  • $\begingroup$ @wzkchem5 No problem, just going through the unanswered queue and checking if some comments exchanges got lost. $\endgroup$
    – Tyberius
    Commented Dec 1, 2021 at 20:15


You must log in to answer this question.

Browse other questions tagged .