[lammps-users] Cascade Simulation : help needed

Hi,
I was trying to do cascade simulation for Fe with EAM potential using
LAMMPS. But each time the simulation is getting aborted due to an error
related with loss of atoms. I have used periodic boundary conditins as
well as shrink wrapped b.c. Both times it ended with same error. I am
using NVE ensemble, Nose-Hoover thermostat. Can anybody plz help me in
this regard ? Can you post an example INPUT file for such kind of
simulation?
with regards.
Prithwish Nandi

Hi,
I was trying to do cascade simulation for Fe with EAM potential using
LAMMPS. But each time the simulation is getting aborted due to an error
related with loss of atoms. I have used periodic boundary conditins as

lost atoms are invariably a sign of an error in the input
(geometries, choice of parameters, etc.). this pops up frequently.
please check the documentation and the mailing list archives for
the possible reasons and solutions.

well as shrink wrapped b.c. Both times it ended with same error. I am
using NVE ensemble, Nose-Hoover thermostat. Can anybody plz help me in

??? if you use a nose-hoover, you don't have an NVE ensemble.

this regard ? Can you post an example INPUT file for such kind of
simulation?

how about you post _your_ input and people comment
about possible problems with it?

cheers,
   axel.

Hi.

Periodic b.cs on a cascade simulation does not make sense because the
cascade gets replicated in the direction the p.b.cs are applied in and
energetic atoms participating in the cascade would start interacting
across the boundaries, unless you use large sized targets, such that
the cascades do not reach the boundaries.

Shrink wrapped boundary conditions do not take into consideration the
background effect of atoms on the cascading atoms that get to the
shrink wrapped boundaries, again, unless you use large sized targets
such that the cascades do not reach the boundaries.

Even if you take a large sized target, the right thing to do for
cascade simulations is to use slab boundary conditions (which is not a
"boundary" choice in LAMMPS as far as I know). It can however be
implemented using the input commands LAMMPS provides you with by
dividing the simulation volume into regions. In slab boundary
conditions you take a large enough size of sample such that the
cascade does not reach the boundaries. Beyond the boundaries you
create a boundary region wherein the atoms are scaled to the target
bulk temperature. For details see Kai Nordlunds notes (lecture2.ps,
page 9, to be specific) at
http://beam.acclab.helsinki.fi/~knordlun/atomistiset/lecture2.ps
You might find the whole series of notes at the above link useful.

Best Regards,
Manoj

2008/11/3 Prithwish <[email protected]>:

In principle, I agree with Manoj. However, there is a simpler way that has been used with some degree of success in the cascade literature: select a domain large enough to contain the cascade with some space to separate it from the boundaries (this is tough in general… especially since the extent of a cascade can vary widely), with periodic boundary conditions (so you don’t need to worry about relaxation). Designate some thin boundary region to be the ‘thermostat’ - apply something simple in there like velocity scaling and keep the rest of the domain NVE integration. This will have the same net effect as what Manoj describes, but won’t require any detailed coding. The drawback is that you really need a fairly large domain, but you should be using a pretty big domain regardless.

Also - avoid system-wide thermostats (i.e. Nose-Hoover over the entire domain). They simply will not give you the correct physics for the sort of system you are studying.

Below is what I have used (for a cascade in 3C-SiC, but principles should still apply pretty well).

Hope this helps,

Dave

units metal # sets units to ‘metal’ units - Atomic units w/ ps timescale
dimension 3
boundary p p p # periodic BCs on x,y plane faces, z dimension fixed volume

atom_style atomic # default attribute style
neighbor 1.0 bin # sets maximum neighbor search radius to cutoff+value, using bin-sort algorithm
neigh_modify delay 5 check yes # checks if neighbor list should be rebuilt every 5 steps

Create geometry using internal stuff

lattice diamond 4.3765
region box block 0 50 0 50 0 50
create_box 2 box # initialize box
create_atoms 2 box basis 1 2 basis 2 2 basis 3 2 basis 4 2 basis 5 1 basis 6 1 basis 7 1 basis 8 1 # create atoms
region rPKA sphere 25 24.5 49.5 .1

mass 1 12.0 # type 1 = C
mass 2 28.1 # type 2 = Si

Define potential

pair_style zbl/tersoff .529 .00552635 1.0
pair_coeff * * SiC_Devanathan_JNM_1998_zbl.tersoff C Si

Groups

region rallatoms block INF INF INF INF INF INF
region rinterior block 2 48 2 48 2 INF
region rexterior block 2 48 2 48 2 INF side out

group interior region rinterior
group exterior region rexterior
group PKA region rPKA

Initialization

compute 1 all temp
compute 2 interior coord/atom 2.2
compute 3 interior ke/atom
#velocity all create 1200.0 34986 # initialize atom velocities to temperature, seed random num generator

fix 1 all nve
fix 2 exterior temp/rescale 1 1200.0 1200.0 0.5 1.0 region rexterior
fix 3 interior temp/rescale 1 1200.0 1200.0 0.5 1.0 region rinterior

timestep 0.001
thermo 100
thermo_style custom step temp pe etotal press
run 1000 # run the equilibriation phase

unfix 3 # remove the rescaling fix from the interior
run 1000 # continue to equilibriate

set PKA velocity to correspond to ~ 10 keV - 2620.549 Ang/ps

set group PKA vx 109.5108
set group PKA vy 301.1547
set group PKA vz -2600.8815

set timestep to smaller value for initial phase of collisions (.01 fs for .2 ps)

timestep 0.00001
thermo 100
dump 1 interior custom 100 test_init_col.dump x y z c_3 tag type c_2
dump 2 PKA custom 100 PKA_traj_init_col.dump x y z c_3 tag type c_2
run 20000 # run the collisional phase for .2 ps

run intermediate phase with intermediate timestep (.1 fs for 1 ps)

timestep 0.0001
undump 1
undump 2
dump 3 interior custom 100 test_inter_evolve.dump x y z c_3 tag type c_2
dump 4 PKA custom 100 PKA_traj_inter_evolve.dump x y z c_3 tag type c_2
run 10000

after initial phase, let evolve for remainder of time using .001 ps timestep (1 fs for 10 ps)

undump 3 # close the previous dumps
undump 4
timestep .001
dump 5 interior custom 100 test_final_evolve.dump x y z c_3 tag type c_2
dump 6 PKA custom 100 PKA_traj_final_evolve.dump x y z c_3 tag type c_2
run 10000

A cascade simulation implies impact by a hi energy ion. Such
simulations can require small timesteps or good potentials to
allow for close approach of 2 atoms. Else you blow atoms out
of the simulation box. See the fix dt/reset and fix nve/limit
commands for ideas.

Steve

Hi,
Thanks everybody. The discussion in this thread helped me to gain a very
very useful insight to model problems related with cascade simulation.
Anyway, I am going to implement the ideas in my next simulation. I shall
reply back very soon describing the status.

With regards.

Prithwish Nandi
Research Scholar
Theory & Computational Studies Section
Indira Gandhi Center for Atomic Research
Kalpakkam - 603 102

Hi,
At first, many many thanks for your response to my post. Your suggestion
is very useful. I reduced the timestep and the issue of "lost atoms"
problem is no more there. Other options whatever you mentioned in your
reply, I have not implemented yet.

I also specially want to thank Dr. Manoj Warrier and David Farrell for
their help.

I have another question. This is all about the analysis of the defect
structure i.e. I want to calculate the number of vacancies and
interstitials using the dump file. For that, I have written a fortran
code. The algorithm is as follows : I have formed a reference structure.
Then I comapare the reference structure with the structure obtained from
the simulation. I defined a cutoff radius ( 0.5*NN distance). If no atom
is there inside this sphere, I declare it as a vacancy and if more than
one atom is there, declaring the sphere as containing interstitials. But
this method is not satisfactory. I am also trying the algorithm based on
voronoi tessellation. But at moment this seems to be too difficult.

My question : Does LAMMPS give the occupancies of each wigner-seitz
cell? (i.e. voronoi tessellation based analysis )

Is there any other way to determine the defect microstructure ?

With highest regards.

Prithwish

You can look at compute coord/atom which calculates the coordination
number of each atom (out to some specified cutoff). LAMMPS stores
no reference structure to compute defects with respect to. You'd have
to write a fix yourself that stored this lattice-dependent structure and
then compute what you want.

Steve