I was trying to model liquid-liquid interface between two immiscible liquids. In order to model it

I decided on using the repulsive part of the LJ potential by changing the cutoff radius to ~2^(1/6) for interactions

between the two liquids. For all other interactions with wall and between the same liquid molecules I have

used a cutoff radius of 2.2.

I am not sure what I am doing wrong but the interface seems to breakdown on visualizing it and the two types of liquid

end up mixing.

Would appreciate any help/suggestions. Have attached the code below.

Thank you.

BR/Joseph

# 3-d LJ flow simulation

units lj

dimension 3

boundary p s p

atom_style atomic

neighbor 0.3 bin

neigh_modify delay 1

# create geometry

region simulation_box block 0.0 24.136 0.0 54.82 0.0 7.22 units box

create_box 4 simulation_box

region fluid_1 block INF 11.0 2.39 52.5 INF INF units box

region fluid_2 block 12.0 INF 2.39 52.5 INF INF units box

region btm_wall block INF INF INF 0.95 INF INF units box

region top_wall block INF INF 53.86 INF INF INF units box

lattice fcc 0.81 orient x 1 1 -2 orient y 1 1 1 orient z 1 -1 0

create_atoms 1 region fluid_1 units box

lattice fcc 0.81 orient x 1 1 -2 orient y 1 1 1 orient z 1 -1 0

create_atoms 2 region fluid_2 units box

lattice fcc 3.24 orient x -1 1 0 orient y 1 1 1 orient z 1 1 -2

create_atoms 3 region btm_wall units box

create_atoms 4 region top_wall units box

mass 1 1.0

mass 2 1.0

mass 3 1.0

mass 4 1.0

# LJ potentials

pair_style lj/cut 2.2

pair_coeff 1 1 1.0 1.0

pair_coeff 2 2 1.0 1.0

pair_coeff 3 3 1.0 1.0

pair_coeff 4 4 1.0 1.0

pair_coeff 1 3*4 0.6 0.75

pair_coeff 2 3*4 0.6 0.75

pair_coeff 1 2 0.6 0.75 1.16

pair_modify shift yes

# define groups

group lower region btm_wall

group upper region top_wall

group boundary union lower upper

group flow subtract all boundary

set group flow type 1

set group lower type 2

set group upper type 3

# initial velocities

compute mobile flow temp

velocity flow create 1.10 482748 temp mobile units box

thermo_modify temp mobile

# thermostating

compute mobile_thermo all temp/partial 0 0 1

fix 3 all langevin 1.10 1.10 1.0 699483

fix 1 flow nve

fix 2 lower nve

fix_modify 3 temp mobile_thermo

thermo_modify temp mobile_thermo

# Couette flow

velocity lower set 0.0 0.0 0.0 units box

velocity upper set 0.0 0.0 0.0 units box

#fix 8 upper move wiggle 10.0 0.0 0.0 40.0 units box #<<<<< CHANGE AMPLITUDE AND TIMEPERIOD

fix 4 boundary setforce 0.0 0.0 0.0

# Run

timestep 0.002

thermo 1000

dump 10 all atom 100 dump.flow

dump 12 all xyz 100 dump_xyz.flow

run 6000