Polycrystalline Surfaces

Hi there,

I have been trying to simulate polycrystalline surfaces of copper. I’m exposing a copper surface to methane molecules and other small hydrocarbons including atomic carbon. To do so, I use the lattice command to create a first lattice orientation, then I create copper atoms in the region I define. This initially create the 211 surface. Then I do the same, but create the 001 surface of Copper. When Creating the region, I let the two overlap, then I use the delete atoms command to create the polycrystalline surface. This simulation takes a long time to run as it is, and sometimes, if I make the runtime too long, it will terminate telling me that some atoms have been lost on a proc. I was wondering if you have any suggestions that could help me improve the speed of the simulations and help me resolve the issue of lost atoms. My Code is pasted below (this one creates a polycrystalline surface with 211 and 001 in the presence of methane:

units metal

atom_style full

dimension 3

boundary p p p

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

region box block 0 20 0 20 0 16 #Define a box for the super-cell

create_box 3 box bond/types 1 angle/types 1 extra/bond/per/atom 4 extra/angle/per/atom 4

region in block 0 10 0 10 0 9.33 side in

create_atoms 1 region in

molecule CH4 methane3.mol

mass 1 63.54 #Define mass of Cu atom

group one region in

lattice fcc 3.6149 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1

region in2 block 0 20 9 20 0 9.33 side in

create_atoms 1 region in2

mass 2 12.01 #Define mass of C atom

group two region in2

pair_style comb3 polar_off #Make use of comb3 potential

pair_coeff * * ffield.comb3 Cu C H #Define comb3 for Cu and C and H

neighbor 0.5 bin

neigh_modify every 1 delay 1 check yes

mass 2 12.01 #Define mass of C and H atoms

mass 3 1.01

delete_atoms overlap 0.3 all all

lattice diamond 3.567 # Lattice constant of C

region out block 0 10 0 10 0 9.33 side out

create_atoms 1 random 10 1 out mol CH4 1

bond_style harmonic

bond_coeff * 5.0 1.0

compute eng all pe/atom

compute eatoms all reduce sum c_eng

timestep 0.00025

min_style cg

minimize 1e-4 1e-6 100 1000

fix 1 all npt temp 1 300 0.025 iso 0 0 0.25

thermo_style custom step temp etotal pe evdwl ecoul c_eatoms press vol lx ly lz xz

thermo_modify norm yes

thermo 1000

Hi there,

I have been trying to simulate polycrystalline surfaces of copper. I’m
exposing a copper surface to methane molecules and other small hydrocarbons
including atomic carbon. To do so, I use the lattice command to create a
first lattice orientation, then I create copper atoms in the region I
define. This initially create the 211 surface. Then I do the same, but
create the 001 surface of Copper. When Creating the region, I let the two
overlap, then I use the delete atoms command to create the polycrystalline
surface. This simulation takes a long time to run as it is, and sometimes,
if I make the runtime too long, it will terminate telling me that some atoms
have been lost on a proc. I was wondering if you have any suggestions that
could help me improve the speed of the simulations and help me resolve the
issue of lost atoms. My Code is pasted below (this one creates a
polycrystalline surface with 211 and 001 in the presence of methane:

some comments:

- COMB/COMB3 are computationally demanding potentials, see
http://lammps.sandia.gov/bench.html#potentials
- using a cutoff of 0.3 seems very small for removing overlapping
atoms. since you do a minimization before running MD, any tension
resulting from close neighbors should be mostly relaxed, but there may
be an optimal value somewhere, and my guess would be that it is larger
than that.
- you are using bonds/angles with a manybody potential. that is a
very, *very* bad idea. bonds are implicit in these potentials. using
bonds triggers excluding "special" neighbors, which in turn will
result in bogus forces. exclusions can only work for pairwise additive
potentials.

axel.

Hi there,

If I don’t use any bonds/angles options it gives me errors. How would I be able to create molecules in my simulation box without using those options.
PS. Sorry for previous email.

Best Regards

Hi there,

I have been trying to simulate polycrystalline surfaces of copper. I’m
exposing a copper surface to methane molecules and other small hydrocarbons
including atomic carbon. To do so, I use the lattice command to create a
first lattice orientation, then I create copper atoms in the region I
define. This initially create the 211 surface. Then I do the same, but
create the 001 surface of Copper. When Creating the region, I let the two
overlap, then I use the delete atoms command to create the polycrystalline
surface. This simulation takes a long time to run as it is, and sometimes,
if I make the runtime too long, it will terminate telling me that some atoms
have been lost on a proc. I was wondering if you have any suggestions that
could help me improve the speed of the simulations and help me resolve the
issue of lost atoms. My Code is pasted below (this one creates a
polycrystalline surface with 211 and 001 in the presence of methane:

some comments:

- COMB/COMB3 are computationally demanding potentials, see
http://lammps.sandia.gov/bench.html#potentials
- using a cutoff of 0.3 seems very small for removing overlapping
atoms. since you do a minimization before running MD, any tension
resulting from close neighbors should be mostly relaxed, but there may
be an optimal value somewhere, and my guess would be that it is larger
than that.
- you are using bonds/angles with a manybody potential. that is a
very, *very* bad idea. bonds are implicit in these potentials. using
bonds triggers excluding "special" neighbors, which in turn will
result in bogus forces. exclusions can only work for pairwise additive
potentials.

axel.

Hi there,

If I don’t use any bonds/angles options it gives me errors. How would I be able to create molecules in my simulation box without using those options.

you are doing things backwards.

before even trying to set up a simulation, you have to *understand*
what kinds of interactions your force field of choice can compute, how
they are implemented, and - most importantly - whether the
parameterization you use is suitable for the system you are aiming to
simulated. the implementation then determines, how you have to specify
your input and what information has to be included or not. as i
already mentioned in my previous response, mixing manybody potentials
in LAMMPS with additional explicit bonds will lead to bogus results.
that is for two reasons: 1) explicit bonds and angles will trigger
removal of "special" pairs from the neighbor lists and since manybody
potentials are not pairwise additive this will create inconsistent and
incorrect force computations as certain n-tuple interactions will no
longer be computed , 2) bonded (or more general short range)
interactions are determined implicitly and usually in a fashion that
bonds can be broken and atoms can change the number of bonded
interactions (and the bonds then change their interaction strength).
this added level of detail and the resulting computational effort is
one of the reasons, why these potentials are so much more time
consuming.

you obviously have not done this due diligence, as then you would not
ask your questions. i thus strongly suggest you make up for that
first. there is no point in continuing with COMB3, if the
parameterization is not suitable. in case it is suitable, the
necessary changes are obvious (remove any explicit bonded interactions
from the topology), and then the errors should go away. it is still
highly advisable to do some tests and validate your input. COMB/COMB3
is not a trivial potential and must not be used in a "black box" way.

axel.