[lammps-users] Help with interdiffusion simulation

Hi everyone,

My ultimate goal is to model interdiffusion at an interface between
titanium and aluminum.

To start, I am simulating the mixing of atoms between two aluminum
blocks that are positioned side by side. Boundary conditions are
periodic in y and z, and shrink-wrapped in x. (I have no experience
with MD simulations, please advise if these assumptions are correct for
my application.)

My input file is as follows:

units metal
atom_style atomic

boundary s p p

lattice fcc 4.1262

region total_domain block -5.5 5 0 2 0 2
create_box 2 total_domain

region al_block block -5.5 -0.5 0 2 0 2
create_atoms 1 region al_block

#lattice bcc 2.870 # iron
#lattce hcp 2.950 # titanium
region ti_block block 0 5 0 2 0 2
create_atoms 2 region ti_block

pair_style eam/alloy
pair_coeff * * ../../potentials/Al_mm.eam.fs Al Al
#pair_coeff * * ../../potentials/AlFe_mm.eam.fs Al Al

velocity all create 1700.0 314159

neighbor 1.0 bin
neigh_modify every 1 delay 5 check yes

fix 1 all nvt 1700.0 1700.0 1.0

timestep 0.001
thermo_style custom step pe press temp
thermo 5000

dump atom_pos all atom 5000 block_interface.dump
dump_modify atom_pos image yes scale no

run 20000

Here is my problem:

When I run the simulation using the potentials in Al_mm.eam.fs, there
are no issues. However when I use the potentials from AlFe_mm.eam.fs,
my simulation blows up and my atoms are lost.

Any insights as to what is happening would be greatly appreciated.

Thanks,
John

Hi, Everyone:

I have two questions about compute_group_group.cpp file.

  1. If I understand correctly, LAMMPS has two types of atoms on each processor: “real” atom and “ghost” atoms. Does the following code build a neighbor list for all of the “real” atoms on each processor?

neighbor->build_one(list->index);

  1. I really have a hard time to understand the following code.

if (j < nall) factor_coul = factor_lj = 1.0;
else {
factor_coul = special_coul[j/nall];
factor_lj = special_lj[j/nall];
j %= nall;

Does “j” represent the global index of an atom? Why is “j” related to nall, which is the total number of “real” and “ghost” atoms on a processor? I searched the archive. It was mentioned these codes were used to turn of LJ and ELE interactions between bonded and connected atoms. Can anyone provide any hints to help me out?

Thanks a lot for your help in advance,
Leo

hi leo,

2009/12/3 Leo Ho <[email protected]>:

Hi, Everyone:

I have two questions about compute_group_group.cpp file.

1) If I understand correctly, LAMMPS has two types of atoms on each
processor: "real" atom and "ghost" atoms. Does the following code build a
neighbor list for all of the "real" atoms on each processor?

neighbor->build_one(list->index);

this will add a _request_ to build the neighbor lists in
a given style, if needed. the neighbor list build is done
at a given point in the simulation and it is done in parallel
and the code will create the lists as they are requested,
optionally just copying information or extracting sublists
from already created lists.

2) I really have a hard time to understand the following code.

if (j < nall) factor_coul = factor_lj = 1.0;
else {
factor_coul = special_coul[j/nall];
factor_lj = special_lj[j/nall];
j %= nall;

Does "j" represent the global index of an atom? Why is "j" related to nall,

no. j is the index in the list of all atoms that are present on the processor
(i.e. local and ghost). the relation to nall is given by the fact that
nall is added
to j in order to indicate that for this atom non-bonded interactions have to
be scaled by a given factor that is to be taken from the special_lj array.
this is for the so-called 1-2, 1-3, and 1-4 exclusions.

so if the i-j entry in the pairlist is a 1-2 bonded interaction, then nall
is added to j while creating the neighborlists; if it is a 1-3 interaction,
then 2*nall is added and if it is a 1-4 interaction, the 3*nall is added.

hope that helps,
    axel.

If your system is all Al, why are you trying to use an alloy file
for Al and Fe? I would contact the author of that potential
file (or see the paper), and check if what you are trying to
do is valid.

Steve

Hi, Axel,

Thank you for your prompt reply! Your answer definite helps me a lot.

It seems that the “i” and “j” in this file mean the unique local index of atoms (local and ghost atoms). One question I am not sure yet is: Does the “i” mean the unique local index of local atoms while the “j” means the unique local index of ghost atoms?
or put it in another way, Does the “ilist” array ONLY contain the index of each local atom while “jlist” array ONLY contains the index of each ghost atom? I really want to know the exact information in “ilist” and “jlist” after the “build_one” subroutine is called.

Thanks,
Leo

2009/12/3 Leo Ho <[email protected]>:

Hi, Axel,

Thank you for your prompt reply! Your answer definite helps me a lot.

It seems that the "i" and "j" in this file mean the unique local index of
atoms (local and ghost atoms). One question I am not sure yet is: Does the

atoms are indexed by their index in the individual domain.
the _global_ index is called "tag".

"i" mean the unique local index of local atoms while the "j" means the
unique local index of ghost atoms?

no. index covers local and ghost atoms.

local atoms have the index 0 to nlocal-1, ghost atoms nlocal to nall-1.

or put it in another way, Does the "ilist" array ONLY contain the index of
each local atom while "jlist" array ONLY contains the index of each ghost

ilist contains the _neighbors_ of "i". that is usually a subset of nall.

atom? I really want to know the exact information in "ilist" and "jlist"
after the "build_one" subroutine is called.

those are "extended verlet lists", i.e. neighbor lists that include
neighbors up
to the cutoff plus a half the "skin" distance. please read the LAMMPS papers
and presentations, that explain this in great detail and much better than i ever
could do it.

cheers,
    axel.

Hi, Alex:

Thank you for your reply! I guess I didn’t make my point very clear. The reason I asked this question is that I want to make sure this subroutine doesn’t have any duplicate counting of interaction forces.

There are two loops in the interact() member function of this file. The outer loop picks atom “i”, then the inner loop search the neighbor atoms (atom “j”) of each atom “i”. My question is "Does EACH atom “i” in the outer loop belongs to local atoms while EACH atom “j” in the inner loop belongs to ghost atoms? If the answer is “yes”, then the current program may result in duplicate counting of interaction forces. If the answer is “no”, that means we may run into both local atoms and ghost atoms when we traverse the outer loop.

Do you have any comments on this question? I will also appreciate any reference specifically talking about how neighbor lists are built in LAMMPS.

Cheers,
Leo

2009/12/3 Leo Ho <[email protected]>:

Hi, Alex:

Thank you for your reply! I guess I didn't make my point very clear. The
reason I asked this question is that I want to make sure this subroutine
doesn't have any duplicate counting of interaction forces.

from where do you get the idea that something would be counted twice?

There are two loops in the interact() member function of this file. The
outer loop picks atom "i", then the inner loop search the neighbor atoms
(atom "j") of each atom "i". My question is "Does EACH atom "i" in the outer
loop belongs to local atoms while EACH atom "j" in the inner loop belongs to
ghost atoms? If the answer is "yes", then the current program may result in

no. how do you get this idea? this doesn't make sense.

duplicate counting of interaction forces. If the answer is "no", that means
we may run into both local atoms and ghost atoms when we traverse the outer
loop.

there is no double counting.

Do you have any comments on this question? I will also appreciate any
reference specifically talking about how neighbor lists are built in LAMMPS.

i already mentioned the papers on LAMMPS. they have references, too!
what else do you want? sorry, but nobody has the time to spoonfeed this
information to you. the compute code that you are looking at is _very_
obvious, iff you understand the principle way how LAMMPS is structured.
this principle is in the LAMMPS papers. please read them and be enlightened.

axel.

Hi, Alex:

Thanks for your reply. I will read more papers and code about LAMMPS before I jump into any conclusion.

Cheers,
Leo

Steve: Thank you for your reply. I realized my mistake; I should be
using pair_style eam/fs rather than eam/alloy.

All: I have some general questions about how to set up conditions for my
simulation.

Instead of a pure Al system, I now have a block of Al sitting next to a
block of Fe. As before, the boundary is periodic in y and z, and
shrink-wrapped in x. (Interface between Al and Fe is on the y-z plane.)
I have added an equilibration step prior to production run where I scale
the velocities to 300K and zero the linear and angular momentum of the
system.

1) How do I deal with the fact that the lattice size is different
between Al and Fe? Due to the ratio of the lattice lengths, I end up
with fractions of a lattice at the periodic boundaries. Do I just pick
a simulation box size large enough such that the inconsistencies at the
boundary becomes negligible?

2) Do I need to impose any constraints on my boundaries in the
x-direction?

3) What distance should separate the two blocks?

Any thoughts would be appreciated.

Thanks,
John

1) How do I deal with the fact that the lattice size is different
between Al and Fe? Due to the ratio of the lattice lengths, I end up
with fractions of a lattice at the periodic boundaries. Do I just pick
a simulation box size large enough such that the inconsistencies at the
boundary becomes negligible?
2) Do I need to impose any constraints on my boundaries in the
x-direction?
3) What distance should separate the two blocks?

All of these are questions about your model you need to
answer to solve the physics you are interested in. There
is no one right answer. I suggest you find a paper doing
something similar and start there.

Steve