Is it reasonable that the energy increases after the NPT relaxation?


The relevant codes are as follows
boundry p p s
fix 1 all npt temp 300 300 0.02 x 0.0 0.0 0.2 y 0.0 0.0 0.2

You are the researcher here and it is your project, so the question would rather have to be: what do supposed the answer should be and what is your explanation for that.

This does not provide much useful information, since what you observe is something that comes from the entire system and depends on the entire setup and particularly what happens to your system during the simulation for example, what happens to the temperature or whether the system expands or shrinks.

Hmm interesting question (even if not a LAMMPS one)… at a first moment, I think it is possible in general. I would say it depends on what was the initial thermodynamic state (I mean the one before the one reached as you relax the system) and what is the shape of your potential.

For example, if you have a system of one particle type that interacts uniquely via the pairwise non-bonded potential shown below, you would have r_min be the equilibrium position at 0K. As you increase the temperature progressively (at constant pressure), your equilibrium relative distance between the pairs should shift in the direction of increase (shown in the plot). This is because of the asymmetry of the potential (you can reason why it happens by thinking carefully about Newtonian mechanics). Naturally, a larger equilibrium position will mean that the new equilibrium volume of the system is overall larger, which fits what you were observing (I mean, I am assuming your system is expanding).

So for example if you were initially in a thermodynamic state where the temperature was 300K and the pressure was such that was keeping the system at a pairwise distance smaller than r_eq at 300K, you would then observe an increase of Ep. Since the Ek is the same (temperature is the same in the initial and relaxed state), your overall internal energy would increase.

A similar reasoning could hold if the assymetry of your well was in the other direction (which is also possible).

PS: I am assuming your energy is not increasing indefinately and your simulation doesnt explode at some point.

Also, this is just an idea of interpretation of the result - as Axel said you need to make a collective assessment of the evolution of different properties and reason it in the context of your system specifically (potential form and etc).

It is generally reasonable (and is likely the case if the simulation starts from a structure at a local minimum of its potential energy surface, e.g. when a minimize is done before the MD starts). At finite temperatures the particles don’t stay at their local minimums in most time, so the total potential energy would be higher.

For your specific scenario, however, I’m more concerned about the fact that your potential energy is not converging. That means the system is not reaching the thermal equilibrium in that timescale (several ns). There could be several reasons for this:

  1. the system equilibrates very slowly due to its nature (e.g. having slow degrees of freedom);
  2. the system is going through some continuous transition and does not equilibrate at all;
  3. there is something wrong with your setup.

Generally speaking, it’s your task to find out the exact cause and determine if it’s acceptable for your purpose.