# How to simulate phase separation in a discrete system using LAMMPS?

I would like to simulate phase separation in LAMMPS and need advice what kind of two particle systems and interactions will yield phase separation.

For example, a meta stable phase separation is achieved with a continuous phase flow model governed by Cahn–Hilliard equation, where $$c$$ is the concentration, $$\mu$$ is the chemical potential, and $$\gamma$$ is a constant: \begin{align} \frac{\partial c}{\partial t} &= D\nabla^2\mu& \mu &= \left(c^3-c-\gamma\nabla^2 c\right)\\ \end{align}

I tried to use a discrete two atom system in LAMMPS. The first type of the particles has an inter-LJ interaction and the second type has an inter-soft-potential interaction. Between the distinct two types there is also a soft potential:

Initial configuration (interleaved lattice of the two particles):

Steady state configuration at $$T=0.1$$ (reduced LJ units):

Configuration file:

units lj

timestep 0.001

dimension 2
boundary p p p
atom_style atomic
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no

variable temp equal 0.1

lattice hex 0.5
region simbox block 0 50 0 25 -0.1 0.1
create_box 2 simbox

lattice hex 0.5
create_atoms 1 region simbox

lattice hex 0.5 origin 0.5 0.5 0
create_atoms 2 region simbox

mass 1 1
mass 2 1.0

velocity all create ${temp} 1234567 dist gaussian pair_style hybrid lj/cut 2.5 soft 0.5 pair_coeff 1 1 lj/cut 2.0 0.8 4.0 pair_coeff 1 2 soft 0.5 pair_coeff 2 2 soft 0.5 fix 1 all nvt temp $${temp}$${temp} 0.1 dump 1 all atom 500 out/dump_temp_${temp}.lammpstrj
write_data data.all
thermo 500
run 500000


1. How do I make my system to look more similar to the continuous model suggested by Cahn and Hiliard? I want to get a more prominent phase separation like in the continuous model. I am not sure how to calibrate my system/interactions.
2. Is there a different simple system I could use to simulate phase separation?
3. If you think Matlab could be more suitable to simulate phase separation, please elaborate.
• According to IUPAC (which I consider to be the authority on terminological issues) 'phase separation' is defined as: "The process by which a single solid or liquid phase separates into two or more new phases." Again, based on looking at your system, you have hard spheres sticking to each other in a simulation space larger than the system, comprising a single phase in vacuum. – Bdrs Apr 18 '18 at 13:51
• @Bdrs I ended up using 2 particle system. But still looking for a more natural system that will experience phase separation. – 0x90 Apr 18 '18 at 13:53
• Great! Try water/1-octanol. That's a popular one, and you should be able to find plenty of empirical data about it. – Bdrs Apr 18 '18 at 13:55
• What is it with that first picture... it seriously messes with my head, no clue why though. Does anyone else have this feeling? – Martin - マーチン Nov 27 '18 at 16:04
• @Martin-マーチン what do you mean? – 0x90 Nov 27 '18 at 16:06

The Cahn-Hilliard model, for which you gave a snapshot, represents the two fluid phases on a much larger length scale than the atomic scale. It's a continuum model. So, whatever particle-based model you use, you would need an enormous system to get close to similar-looking pictures. The interfaces will always look more irregular, in a practical particle-based simulation. Don't expect to see nice circular droplets.

Secondly, even that Cahn-Hilliard simulation is not yet at equilibrium. The droplets are in the process of coalescing, and the final state will be a single large droplet of one phase in the other (or possibly the two phases will form slabs, with planar interfaces between them, it depends somewhat on the boundary conditions and dimensions of the surrounding box). The driving force is the reduction of the total interfacial area. This whole process takes time. Something similar should happen in the particle-based simulations, but because the length scales are so different, the time scales will also be different.

Given all that, your first choice of models may be suitable. However, you need to choose a temperature which is high enough for both the phases to be liquid. You haven't done this! The triple point temperature of 2D Lennard-Jones is $$k_BT/\epsilon\approx 0.4$$. At your temperature of $$0.1$$, the Lennard-Jones atoms will be trying to crystallize; in practice they will get "stuck" in irregular semi-crystalline semi-amorphous clumps. Nonetheless, if you raise the temperature above $$0.4$$, you might see liquid-liquid coexistence for this system. I can't be sure, and you may need to experiment.

I have some suggestions for a simpler approach, if all you want to do is see liquid-liquid phase separation using LAMMPS. Firstly, it is OK for both the fluids to be identical. For example, calling the two species A and B, the A-A and B-B interactions can both be Lennard-Jones. If you choose the A-B interaction to be different, and unfavourable compared to Lennard-Jones (for example, the completely repulsive Weeks-Chandler-Andersen version of Lennard-Jones), you should see phase separation. The advantage of this approach is that the Lennard-Jones phase diagram is well known, making it fairly easy to choose a density and temperature at which both the phases will be liquid.

Even simpler, you can choose A and B to interact through a coarse-grained potential such as the dissipative particle dynamics DPD potential $$v(r) = \begin{cases} \frac{1}{2}a(1-r)^2 & r<1 \\ 0 & r>1 \end{cases}$$ where the cutoff distance is taken as the unit of length. LAMMPS has a DPD option, and you can find a perspective article on DPD by Español and Warren J Chem Phys, 146, 150901 (2017) also available as a preprint. There are also tutorials on DPD readily available on the internet. DPD has the advantage of being very fast; also it is sort-of intermediate between an atomic model and a continuum one, since the particles in DPD really represent regions of fluid, not individual molecules.

Once again, the A-A and B-B interactions can be made identical, for example $$a=25$$, with a typical temperature $$k_BT=1$$ and a density $$\rho=3$$ in 3D; you would need to experiment for suitable parameters in 2D, but the precise values may not be crucial. The A-B interaction could be the same, but with a larger value of $$a$$, for instance $$a=50$$. This approach has already been used to model fluid-fluid separation, for instance in papers by Novik and Coveney Phys Rev E, 54, 5134 (1996) also available as a preprint, and Phys Rev E, 61, 435 (2000) also available as preprint. Those papers are a bit involved, and there are many complicated technical issues with the DPD method. It is not essential to go into those in depth. If you are not essentially interested in the hydrodynamics, it is quite satisfactory to use the DPD potential in a standard molecular dynamics simulation with a thermostat of your choice.

• I simulated the data from $T=0.01$ to $T=0.1$ only bellow $T=0.5$ a structure that resembles phase separation emerges. For temperature higher than $T=0.4$ both particles have a continuous phase. – 0x90 Nov 13 '18 at 0:08
• I don't quite understand your temperatures. You say that you saw phase separation below 0.5 but you did not see it above 0.4? Can you please clarify? More generally, it is hard to predict what you will see for a model which has not been studied before. I have suggested, in my answer, using the DPD potential, which has been studied before. Suitable parameters are suggested in the papers I cited, and you could also try the values I suggested. I suspect that there is a wide range of DPD parameters that will work. This may give you a better chance of seeing what you want to see. – LonelyProf Nov 13 '18 at 7:50
• it's a typo. There is a $T_c$. It's around 0.5 – 0x90 Nov 13 '18 at 13:14
• OK. The only other suggestion I have for that system is to make the 1-2 interaction more repulsive, instead of being just the same as the soft potential. For example WCA Lennard-Jones, which is fairly easy to do in LAMMPS: lammps.sandia.gov/threads/msg02715.html This would discourage the mixing more strongly. Maybe push $T_c$ up a bit? Perhaps worth trying at least. – LonelyProf Nov 13 '18 at 13:30
• Short answer is "no". You would need a way of predicting the thermodynamic free energy of each of the phases, as a function of temperature, pressure, composition, knowing just the interaction potentials. This is not a trivial task to do, and is a whole field of science in itself. Doing particle-based simulations (as you are doing) is actually one of the standard ways of finding out which phase is stable. – LonelyProf Nov 13 '18 at 14:12