I am trying to generate a phase diagram of pressure versus temperature in a NPT ensemble using molecular dynamics method. I am using the Lennard Jones potential with 108 particles and Berendsen thermostat and barostat to fix the parameter.

In order to observe phase transition I plot the radial distribution function (RDF) and also density versus temperature graph for various fixed temperature. However it is hard to differentiate liquid and gaseous phase from RDF.

Another peculiar thing is that there should be a discontinuity in density as the transition temperature is reached, which I observe only for solid liquid transition but not for the liquid vapor transition. The change in density should be huge for liquid vapor transition. Any idea on what possibly goes wrong?

Periodic boundary
cut off radius = 2.5
Truncated and shifted lennard jones potential
number of particle = 108
number of steps = 2000 timestep = 0.004
Berendsen barostat of scale coeefficient of 0.1 is used to scale the size of simulation box

The density versus temperature graph for pressure of 1.0MPa:

The solid liquid transition can be seen but not the liquid vapor transition.

  • $\begingroup$ Your temperature might be above (or close) to the critical temperature of the triple point... $\endgroup$
    – Fabian
    Dec 27, 2012 at 21:48
  • $\begingroup$ I gradually increase my temperature from 80K to 160K, so it will surely show a distinct phase transition $\endgroup$
    – Gosere
    Dec 28, 2012 at 1:30
  • 1
    $\begingroup$ if you would give all the parameters (the parameters of the Lennard Jones, which pressure, which temperature) in the question then maybe somebody can help. $\endgroup$
    – Fabian
    Dec 28, 2012 at 7:22
  • $\begingroup$ In addition, specify your cell and boundary conditions, potential truncation method, what is the average cell size corresponding to you preassure. Why don't you also just post your "density versus temperature graph"? Apperantly it's your own programm, am I right? $\endgroup$
    – Yrogirg
    Dec 28, 2012 at 10:20

2 Answers 2


First, some might be interested in your reduced pressure. I guess you are working argon parameters of $\varepsilon = 119.8 \text K$ and $\sigma = 3.405 \overset{\circ}{\text A}$, thus your pressure of 1.0 MPa (10 atmospheres) in a reduced form:

$$p^* = \frac{p \sigma^3}{\varepsilon} \approx 0.024 $$

Anyway, let's see where are you. From your plot I've deduced several poins of yours:

  ρ     T
 0.9    .7
 0.7    .75
 0.6    .9
 0.5   1.05
 0.4   1.1
 0.3   1.15
 0.2   1.25

and marked them on the density-temperature phase diagram of a Lennard-Jonesium:

enter image description here

I took the figure from the article I've already mentioned in your aother question, the article is by JP Hansen, L Verlet "Phase Transitions of the Lennard-Jones System." As you see, generally you start from a liquid+solid, then as you lower the density (or increase the temperature) solid melts into liquid, then somehow you reach the region of vapor-liquid coexistence where liquid gradually evaporates until you end up with mostly vapor phase.

Most of your process lies in a region of vapor-liquid coexistence, some part in the region of solid and liquid coexistence. The switch between the two occurs near $\rho =0.8$ and indeed your plot shows a peculiarity at these point. You should better resolve this region you may be near the triple point as Fabian suggests.

To sum up, $\rho(T)$ should be continuous, it is $\frac{d \rho}{d T}$ that should not. You plot seems to exhibit such a behavior, but I cannot tell for sure. You have to accurately measure the transition between the liquid and liquid+vapor to find the discontinuity of $\frac{d \rho}{d T}$. The problem is it can be hard to hit into liquid, since liquid for LJ (and argon) occupies only a small gap of $T$ and $\rho$. Why have you put $p$ to be 10 atm, maybe you wanted to have 1 atm?

You cannot rely on an RDF in the region of vapor-liquid coexistence, it has hardly any sense there --- basically you have bubbles of liquid floating in a vapor. Though due too the small size of your system it's hard to talk about bubbles.

You might want to study the system clearly exhibiting vapor-liquid interface, that is when the one part of your cell has vapor, the other one has liquid. Thus you could actually see how the amount of liquid diminishes due to evaporating. Given $z$ is a vertical axis, perpendicular to the surface of the liquid you can study $\rho(z)$ instead of RDF.

Last, but not least, obviously, I could not plot your points accurately on the phase diagram, be cautious not to rely on these too much, that's just an estimation. And the phase diagram itself is intended to reflect the true LJ, not the shifted and truncated one.

  • $\begingroup$ Is it suitable then to fix the density and varying the temperature to observe the transition? I am thinking of using compressin factor or potential energy as an indicator of phase transition $\endgroup$
    – Gosere
    Dec 31, 2012 at 13:03
  • $\begingroup$ Also is it possible to reproduce the above Verlet's result by using the shifted and truncated LJ potential (with potential and pressure correction terms)? $\endgroup$
    – Gosere
    Dec 31, 2012 at 13:30
  • $\begingroup$ Is the above diagram dependent on initial configuration of particles? I notice that the pressure I get at a particular temperature and density is different for simple cubic structure and face centered cubic (FCC) structure. $\endgroup$
    – Gosere
    Jan 3, 2013 at 2:41
  • $\begingroup$ @Sheng Sorry, I'm already moving beyond my expertise, I've never been deep in molecular simulation of fluids, just some basic things and I have never studied phase transitions, that seems to be quite and area. You should really study the state of the art on the subject and ask someone who knows better. Is that your diploma? If yes try to consult with your science advisor. But things are not so bad --- I'd say you have to just double your knowledge in the area (just double, no more) to be able to do it. Do you have all the books you need? $\endgroup$
    – Yrogirg
    Jan 3, 2013 at 17:03

For those who may read this later: you cannot find the liquid and vapor phases independently in this way. The reason is that for densities in between the liquid and vapor densities, the system will phase segregate into droplets within the simulation, and so your total density will include contributions from both phases. Consequently you can change the pressure or density continuously just by changing the occupancy of each phase. If you really want to see the separate phases with low particle numbers, Gibbs Ensemble Monte Carlo is probably the best bet. This method sets up two subsystems and allows them to exchange volume and particles with each other. Total free energy will be minimized when each subsystem is in a single phase, one liquid and the other vapor. It is easy to read off the densities from there.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.