Water Molecules in Couette Flow

Dear all,

I am new to LAMMPS and have been stuck on this problem for some time. I was unable to find a solution in the LAMMPS manual, and I am hoping someone could help me. I am trying to perform a simulation in LAMMPS in which a collection of water molecules are bound by two plates, where the plates begin moving in opposite directions half way through the simulation. I aim to obtain the density and velocity field from this interaction. An initial coordinate file has been created for the water molecules using TopoTools in VMD. The water molecules follow the TIP3P water model style and interact with the walls and each other is LJ-like.

Once I begin my simulation, I allow my water molecules to move around freely to reach an equilibrium point. They easily fill the entire simulation box (bounded by the two plates and periodic in the x and z direction); however, after a while the water molecules seem to condense together, usually on one of the plates. This is a problem because the water is no longer in contact with both plates so when the two plates begin moving, there is no shear on the water molecules.

I am wondering if I am using the correct pair_coeff for a water-carbon interaction (where my wall is made of carbon), the correct thermostating technique, or if there is another obvious problem with my script.

The input script is shown below for the system.

Basic Set-up

units real
dimension 3
boundary p s p

atom_style full
neighbor 0.3 bin
neigh_modify delay 1 one 30000 page 300000000

region simulation_box block -10.5 10.3 -18.0 18.0 -10.3 10.3 units box
create_box 4 simulation_box bond/types 1 angle/types 1 extra/bond/per/atom 2 extra/angle/per/atom 1

region btm_wall block INF INF INF -15.2 INF INF units box
region top_wall block INF INF 14.9 INF INF INF units box

lattice fcc 2.50 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

read_data SolvantBox_10x15x10.data add append

pair_style hybrid lj/cut/coul/cut 13.0 lj/cut 8.0
bond_style harmonic
angle_style harmonic

mass 1 1.0080000 #HT
mass 2 15.999400 #OT
mass 3 12.010700 #Btm_Wall-Carbon
mass 4 12.010700 #Top_Wall-Carbon

pair_coeff 1 1 lj/cut/coul/cut 0.00000 0.0000 #HT-HT
pair_coeff 1 2 lj/cut/coul/cut 0.00000 0.0000 #HT-OT
pair_coeff 1 3 lj/cut 0.50000 1.0000 #HT-Wall
pair_coeff 1 4 lj/cut 0.50000 1.0000 #HT-Wall
pair_coeff 2 2 lj/cut/coul/cut 0.10200 3.1880 #OT-OT
pair_coeff 2 3 lj/cut 0.74950 1.1800 #OT-Wall
pair_coeff 2 4 lj/cut 0.74950 1.1800 #OT-Wall
pair_coeff 3 3 lj/cut 0.00000 0.0000 #Wall-Wall
pair_coeff 3 4 lj/cut 0.00000 0.0000 #Wall-Wall
pair_coeff 4 4 lj/cut 0.00000 0.0000 #Wall-Wall

bond_coeff 1 450.0 0.9572 #HT-OT

angle_coeff 1 55.0 104.52 #HT-OT-HT

group water type 1 2
group lower type 3
group upper type 4
group boundary union lower upper

initial velocities

compute mobile water temp
velocity water create 310.0 12345 temp mobile units box
thermo_modify temp mobile

thermostating

compute mobile_thermo water temp/partial 0 0 1
fix lngv_2 water langevin 310.0 310.0 50 123456
fix 1 all nve
fix_modify lngv_2 temp mobile_thermo
thermo_modify temp mobile_thermo

equlibrating flow

velocity lower set 0 0.0 0.0 units box
velocity upper set 0.0 0.0 0.0 units box
fix 3 boundary setforce 0.0 0.0 0.0

timestep 0.01
thermo 500
dump 1 all atom 500 H2O_SolvantBox.dump
run 500000

Couette flow

velocity lower set -0.025 0.0 0.0 units box
velocity upper set 0.025 0.0 0.0 units box
fix still boundary setforce 0.0 0.0 0.0

compute chunky all chunk/atom bin/1d y 0.0 0.5
fix profile all ave/chunk 1 200000 1000000 chunky vx density/number density/mass file H2O_SolvantBox.profile

run 1000000

Any response would be greatly appreciated, and if you wish to receive my read_data file please let me know.

Kind regards,

Matthew Brown

Physics students at the University of Florida

Dear all,

I am new to LAMMPS and have been stuck on this problem for some time. I was
unable to find a solution in the LAMMPS manual, and I am hoping someone
could help me. I am trying to perform a simulation in LAMMPS in which a
collection of water molecules are bound by two plates, where the plates
begin moving in opposite directions half way through the simulation. I aim
to obtain the density and velocity field from this interaction. An initial
coordinate file has been created for the water molecules using TopoTools in
VMD. The water molecules follow the TIP3P water model style and interact
with the walls and each other is LJ-like.

Once I begin my simulation, I allow my water molecules to move around freely
to reach an equilibrium point. They easily fill the entire simulation box
(bounded by the two plates and periodic in the x and z direction); however,
after a while the water molecules seem to condense together, usually on one
of the plates. This is a problem because the water is no longer in contact
with both plates so when the two plates begin moving, there is no shear on
the water molecules.

I am wondering if I am using the correct pair_coeff for a water-carbon
interaction (where my wall is made of carbon), the correct thermostating
technique, or if there is another obvious problem with my script.

please note that what you are asking about are not really questions
about LAMMPS, but rather questions about how to set up and run MD
simulations correctly. thus it is something that you should primarily
discuss with your adviser and/or colleagues that already have
experience in running MD simulations.

as for pair_coeff values, those are determined by the force field you
have selected and thus something that you will have to verify by
yourself. this is not "the one and only correct" parameter, but it
depends on the overall approximation defined in your force field. the
rest are all very minor considerations.

what i would worry about is whether you have chosen a suitable initial
water density to start with. the behavior you describe suggests that
you didn't.

axel.