I want to simulate a two-phase system, which is wall-fluid1-fluid2-wall (basically a Couette flow with two phase) and all follow the Lennard Jones potential. What is the best way to set up the system? for one fluid I used “region” and “create_atoms” with no hassle (walls have the same density as the fluid) but having different phases with different densities and arrange them on a lattice makes the set-up complicated as the lattice length sizes are different so I get much higher/lower densities at a phase than the density that I expect. I hope I was clear in explaining the issue. Would you help me with this I really appreciate it.

I want to simulate a two-phase system, which is wall-fluid1-fluid2-wall
(basically a Couette flow with two phase) and all follow the Lennard Jones
potential. What is the best way to set up the system? for one fluid I used
"region" and "create_atoms" with no hassle (walls have the same density as
the fluid) but having different phases with different densities and arrange
them on a lattice makes the set-up complicated as the lattice length sizes
are different so I get much higher/lower densities at a phase than the
density that I expect. I hope I was clear in explaining the issue. Would you
help me with this I really appreciate it.

since you want to simulate fluids, you don't need to worry
too much about perfect match of the lattice(s), you will
have to equilibrate your system *anyway*, and the resulting
densities should be a function or the proper choice of parameters
and not the initial placement. it still helps to get close to the
final equilibrated state at the setup phase.

that being said, i don't see where there is a problem with
the lattice command. you can just use a suitable lattice
command for the first system and use region to select the
placement of those atoms and then use a second lattice
command to define the second set of atoms. you need to
make sure that all boxes, regions are either defined consistently
using the first lattice command or that you use "units box"
or else it can be a bit confusing.

from then one, you need to carefully equilibrate.
liquid/liquid system are a bit tricky. you don't
want to disturb them too much, and the setup
is going to be pretty disturbing. so careful
monitoring of the system and a quench or two
might help.

Thank you Axel for your thorough reply. I have one more question, it’s actually a confusion that I have. I tried to test a very simple set-up to understand what is going on so I just build a system which is consisted of two phase (phase1 density = 0.755 and phase 2 density = 0.81) and used just one lattice for each phase as follows:

lattice fcc 0.755
region phase1 block 0 1 0.001 1 0 1 units lattice
region phase2 block 0 1 1.001 2 0 1 units lattice
region total union 2 phase1 phase2
create_box 2 total

lattice fcc 0.755
create_atoms 1 region phase1

lattice fcc 0.81
create_atoms 2 region phase2

By looking at the initial set-up (without turning on any interaction between particles) I got the following coordinate for them:

So the density for type 1 is 4/(1.7433)^3 = 0.755 but for type 2 I get 9/(1.7433)^3 = 1.698 which is way off from what I expect ( 0.81). I sort of can understand the way lammps sets particles and these results are inevitable since type 2 lattice sides are smaller so basically 5 more particles which are associated with adjacent lattices that are included in the region contribute to these large number of particles in phase 2…so would you please tell me in which part my problem lies? By the way I’m going to do NVT simulation which I assume the volume of my system should remain constant and so the number of particles created initially divided by the volume of the region should give me a density close to what I expect (0.755 for phase 1 and 0.81 for phase 2) and two fluids are immiscible.

Thank you Axel for your thorough reply. I have one more question, it's
actually a confusion that I have. I tried to test a very simple set-up to
understand what is going on so I just build a system which is consisted of
two phase (phase1 density = 0.755 and phase 2 density = 0.81) and used just
one lattice for each phase as follows:

lattice fcc 0.755
region phase1 block 0 1 0.001 1 0 1 units
lattice
region phase2 block 0 1 1.001 2 0 1 units
lattice
region total union 2 phase1 phase2
create_box 2 total

lattice fcc 0.755
create_atoms 1 region phase1

lattice fcc 0.81
create_atoms 2 region phase2

By looking at the initial set-up (without turning on any interaction between
particles) I got the following coordinate for them:

So the density for type 1 is 4/(1.7433)^3 = 0.755 but for type 2 I get
9/(1.7433)^3 = 1.698 which is way off from what I expect ( 0.81). I sort of
can understand the way lammps sets particles and these results are
inevitable since type 2 lattice sides are smaller so basically 5 more
particles which are associated with adjacent lattices that are included in
the region contribute to these large number of particles in phase 2...so
would you please tell me in which part my problem lies? By the way I'm going

yes. you define the total volume and the volume of the second phase
in terms of the lattice constant of the first system. since the two are
not commensurate, there has to be a mismatch. it would probably be
better to compute volumes in terms of box units.

with so few atoms, you will have a hard time to get both
densities represented properly.

to do NVT simulation which I assume the volume of my system should remain
constant and so the number of particles created initially divided by the
volume of the region should give me a density close to what I expect (0.755
for phase 1 and 0.81 for phase 2) and two fluids are immiscible.