creating slabs

Dear Lammps users,

I wish to simulate a simple solid-liquid interface in which both the solid and liquid species are represented by simple Lennard-Jones interactions with epsilon parameters 10 and 1 respectively.

I wish to generate the solid slab with minimum free-energy and surround it with the liquid species and am wondering what the best way to do this in Lammps is? Ultimately I want to run my simulations in the NVE ensemble.

Many thanks, Jeff.

There are many ways. (I myself wrote a script to construct a data
file containing atoms arranged in a crystal structure.) Probably the
simplest way makes use of the "create_atoms" and "lattice" commands.
Check out:
http://lammps.sandia.gov/doc/create_atoms.html
http://lammps.sandia.gov/doc/lattice.html
but that's not why I posted a reply...

---- main post ----

I wanted to warn you that equilibrating the liquid around the slab at
constant pressure can be very tricky, depending on how you plan to
immobilize the slab. One way to do this is to define two "groups", a
mobile and an immobile group (consisting of the atoms in your slab).
With this approach, you only want to integrate the equations of motion
on one of the groups (the mobile group).
http://lammps.sandia.gov/doc/group.html
http://lammps.sandia.gov/doc/fix_nh.html
(Don't use isotropic pressure rescaling. Use "z" instead of "iso".
See example below.)

This approach uses "fix npt ... dilate". Unfortunately, this feature
of lammps is currently broken, but you can obtain a patch to fix it on
this thread:
http://sourceforge.net/mailarchive/message.php?msg_id=29497398
Up-to-date example files here:
http://sourceforge.net/mailarchive/message.php?msg_id=29497399

---- Details ----
Getting the pressure right may require redefining the way the
temperature is calculated to eliminate the degrees of freedom you have
removed. (I also found. Whether or not you need to use the "fix
rigid" command remains a point of contention.
http://lammps.sandia.gov/threads/msg23402.html
(Look at: comparison_high_pressure.png. The initial geometry of the
system in that example was pre-equilibrated at 10000 bar, and the
simulation was continued at 10000 bar. If you use fix rigid, it seems
to do a better job of maintaining the original system volume, although
it causes lammps to report the pressure incorrectly. The change in
free volume is a few percent, because about 75% of the atoms are
immobilized in this example.)

-------- Here is an excerpt of the files posted in that thread -------
------- (so you don't have to read the entire messy discussion) -------

# I'm not sure if the next two lines below are necessary, but they don't hurt:
velocity immobile zero angular
velocity immobile zero linear

# Again, the next line (fix ... rigid) is optional and is somewhat
controversial:

fix Ffreezestuff immobile rigid single force * off off off torque * off off off

# The developers (Steve and Trung) don't think it is necessary to include this
# line, but I use it because. See image attached to thread in discussion above.

neigh_modify exclude group immobile immobile

compute tempMobile mobile temp

compute pressMobile all pressure tempMobile

thermo_style custom step c_tempMobile c_pressMobile temp press vol

fix Fmovestuff mobile npt temp 300.0 300.0 100.0 z 10000.0 10000.0
1000.0 dilate mobile

#(The command above beginning in "fix" and ending in "mobile"
# should fit on one line. {Sometimes email breaks it into multiple lines.})
# "mobile" is the name of the group of atoms which can move.
# (You might name it something different.)
# The pressure in this example is 10000 bar

fix_modify Fmovestuff temp tempMobile

I hope this helps.

Andrew