# Free surface in one direction

I am simulating thin film alumina having nano thickness in Z direction. I want to create free surface in Z direction. I am using fixed boundary conditions and fix reflect wall in z direction for creating a free surface. X and Y are periodic boundary conditions.

Is this right way for creating a free surface. I get a melting point of alumina far above the required value using this method. Is there any other mehtod

I dont want to simulate bulk alumina but nanoscale thin film of alumina having nanothickness in z direction. For this I use periodic boundary conditon in x and y direction but fixed boundary condition in z direction. But using fixed Boundary, atoms get lost so to avoid them I use fix wall /reflect.

Hi Nikhil.

By fixing wall and reflecting atoms that tend to go out, you are artificially constraining the atoms to be in the simulation region and in effect you are not having a “free boundary”. Try the following in order to not lose atoms:

1. Choose realistic alumina density and create atomic positions.

2. Carry out minimization of the system using the “minimize” command.

3. Start with a very low temperature of the system (say 1 K) and raise the temperature gradually to desired temperature.

If the above does not work, try

1. Using the fix dt/reset command too.

If none of the above work, send a sample to the mail-list and perhaps someone may find time to go through it. Try to understand the implementation of BCs, etc in MD so that you can judge if the constraints you impose correspond to the conditions you are simulating.

Best’

Manoj

Thanks for your comments.

But in fix dt/reset commnad one has to input maximum distance for an atom to move in one timestep. How to decide that provided I have 1 nm thin film and I dont want atoms to get lost while using fixed boundary condition.

The thumb rule is that the fastest atom in the simulation should not move more than 1/20 times the inter-atomic distance. The physics/numerics behind this is that the inter-atomic potential should be properly sampled. Check out how the interatomic potential varies between two atoms and what would happen if your step size was say 1/2 the interatomic distance. As was suggested in this mail-list by someone else, spend time to understand the MD algorithm.

Manoj

Thanks for your comments.

But in fix dt/reset commnad one has to input maximum distance for an atom to
move in one timestep. How to decide that provided I have 1 nm thin film and
I dont want atoms to get lost while using fixed boundary condition.

since you don't post your input deck, people can only speculate what
is going on. most of them will assume that you know what you are doing
and that your initial input is reasonable. however, it seems to me,
that your input has a serious problem. if you have a solid with proper
box dimensions and fixed boundaries, you should not lose atoms.
rather than trying to band-aid your simulation with walls and messing
with time steps, you need to find the origin of what is going on. i am
willing to bet a substantial amount of money, that it is something
rather trivial (and possibly embarrassing).

but if you want competent and accurate help, you need to post a
complete input deck.

axel.

Hello Sir,

Here is my input script

units metal
atom_style full
boundary p p f

# Create atoms

lattice custom 4.76 a1 1 0.0 0.0 &
a2 -0.5 0.866 0.0 &
a3 0.0 0.0 2.7296 &
basis 0.306 0.0 0.25 basis 0.0 0.306 0.25 basis 0.694 0.694 0.25 basis …

region box prism 0.0 4.76 0.0 4.122 0 12.993 -2.38 0.0 0.0 units box
create_box 2 box
create_atoms 1 basis 1 2 basis …
mass 1 26.981539
mass 2 15.9994

#Define Potential
pair_style vashishta
pair_coeff * * AlO.vashishta Al O

timestep 0.0005
thermo 10
thermo_style custom step temp press vol

velocity all create 0.1 102486 mom yes rot yes dist gaussian

fix 3 all npt temp 0.1 0.1 1.0 x 0.0 0.0 1.0 y 0.0 0.0 1.0 couple xy
run 50000
unfix 3

variable tmp index 0.001 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
variable newtemp index 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
variable temperature equal v_tmp100
variable newtemperature equal v_newtemp
100
label looptemp
fix 4 all npt temp {temperature} {newtemperature} 1.0 x 0.0 0.0 100.0 y 0.0 0.0 100.0 couple xy
run 100000
next tmp
next newtemp

jump aluminamelt.in looptemp

Hello Sir,

Here is my input script

this is incomplete and thus i cannot run it. however, it looks like
you are placing atoms right on the box boundary.
with a fixed box, that is *asking* for losing atoms, as simple thermal
vibrations can move an atom outside the box.
applying a wall, isn't the solution. why not use shrinkwrap boundaries?

also, what is the point of that crazy loop? fix npt can do a linear
temperature ramp out of the box.

it also has been pointed out to your already, that this kind of ramp
is not a good way to determine a melting point.
have you never heard about hysteresis and nucleation?

axel.

Thanks Sir for your comments. Using shrink wrap boundaries , drops the density near to zero after some time steps. The above method helped me to find melting of bulk alumina using periodic boundary condition in all directions. I just want to go one step ahead with thin alumina film.

Complete script is here:

units metal
atom_style full
boundary p p f

# Create atoms

lattice custom 4.76 a1 1 0.0 0.0 &
a2 -0.5 0.866 0.0 &
a3 0.0 0.0 2.7296 &
basis 0.306 0.0 0.25 basis 0.0 0.306 0.25 basis 0.694 0.694 0.25 basis 0.694 0.0 0.75 basis 0.0 0.694 0.75 &
basis 0.306 0.306 0.75 basis 0.0 0.0 0.352 basis 0.0 0.0 0.648 basis 0.0 0.0 0.852 basis 0.0 0.0 0.148 &
basis 0.639 0.667 0.917 basis 0.333 0.973 0.917 basis 0.027 0.361 0.917 basis 0.027 0.667 0.417 basis 0.333 0.361 0.417 &
basis 0.639 0.973 0.417 basis 0.333 0.667 0.019 basis 0.333 0.667 0.315 basis 0.333 0.667 0.519 basis 0.333 0.667 0.815 &
basis 0.973 0.333 0.583 basis 0.667 0.639 0.583 basis 0.361 0.027 0.583 basis 0.361 0.333 0.083 basis 0.667 0.027 0.083 &
basis 0.973 0.639 0.083 basis 0.667 0.333 0.685 basis 0.667 0.333 0.981 basis 0.667 0.333 0.185 basis 0.667 0.333 0.481

region box prism 0.0 4.76 0.0 4.122 0 12.993 -2.38 0.0 0.0 units box
create_box 2 box
create_atoms 1 box basis 1 2 basis 2 2 basis 3 2 basis 4 2 basis 5 2 basis 6 2 basis 7 1 basis 8 1 basis 9 1 basis 10 1 &
basis 11 2 basis 12 2 basis 13 2 basis 14 2 basis 15 2 basis 16 2 basis 17 1 basis 18 1 basis 19 1 basis 20 1 &
basis 21 2 basis 22 2 basis 23 2 basis 24 2 basis 25 2 basis 26 2 basis 27 1 basis 28 1 basis 29 1 basis 30 1

mass 1 26.981539
mass 2 15.9994

#Define Potential
pair_style vashishta
pair_coeff * * AlO.vashishta Al O

timestep 0.0005
thermo 1
thermo_style custom step temp press vol etotal ke pe

velocity all create 0.1 102486 mom yes rot yes dist gaussian

fix 3 all npt temp 0.1 0.1 1.0 x 0.0 0.0 1.0 y 0.0 0.0 1.0 couple xy
run 50000
unfix 3

variable tmp index 0.001 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
variable newtemp index 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
variable temperature equal v_tmp100
variable newtemperature equal v_newtemp
100
label looptemp
fix 4 all npt temp {temperature} {newtemperature} 1.0 x 0.0 0.0 100.0 y 0.0 0.0 100.0 couple xy
run 100000
next tmp
next newtemp

jump aluminamelt.in looptemp

AlO.vashishta (1.53 KB)

Thanks Sir for your comments. Using shrink wrap boundaries , drops the
density near to zero after some time steps. The above method helped me to

have you visualized the trajectory? shrinking density would mean that
atoms "fly off the surface" and thus you are creating an unstable
surface.

find melting of bulk alumina using periodic boundary condition in all
directions. I just want to go one step ahead with thin alumina film.

it doesn't make sense that you need such a complex loop for that. a
single ramp should do as well.
since you are heating up rather rapidly, you shouldn't use a
nose-hoover thermostat in the first place, as those are not meant to
transfer large amounts of kinetic energy quickly.

overall, what you are saying doesn't make too much sense. if all you
care is just to melt your material, you can just initialize and keep a
sufficiently high temperature and be done with it. if you care about
the transition, you would have to be *much* more careful and worry
about heating rates, hysteresis and so on.

axel.