# Different results with cubic box and cylindrical box

Dear all,

I am investigating collision cascades in Fe. To identify interstitials and vacancies I use the Voronoi package. So far I was using a cubic box (eg 100100100 lattice param) with periodic boundary conditions. In order to avoid overlap and influence of shock wave on the development of the cascade, I use a Berendsen thermostat (T = 300K) two layer thick around the cubic box to extract the energy of the box. The results I obtain are as expected. Everything goes smooth.

Few days ago I had the idea to use a cylindrical box instead of a cubic one. Since I know (more or less) the direction of the cascade, I thought I could create a cylinder along the direction of the cascade in order to avoid large cubic boxes in which many atoms are not playing any role (in particular far from the center of the cascade). Of course, the cylinder must be large enough to avoid border effects.

To do so, I did the following:

• As before, I create a cubic box (a large one) with Fe atoms.

• Inside this cubic box, I define a sufficiently large cylinder of appropriate radius and height that will contain the cascade. And I rotate it. All this with the region command.

• Then, I create a cylinder a bit larger (+ 2 layers in radius and height).

• The thermostat is the result of subtracting the inner cylinder from the outer one. It is an empty cylinder with a thickness of 2 atomic layers.

• Then, I delete the rest of the atoms around the largest cylinder.

• Then, the rest as for the cubic box. Equilibration, MD steps, thermostat, etc…

All these steps are summarized below in LAMMPS language:

# Angle of rotation

variable theta equal -0.882228287749

# Definition of inner cylinder

region cylindin cylinder z 59.4231067727 73.9774529773 29.6879405819 26.5231503881 96.230594021 units lattice rotate v_theta 59.4231067727 73.9774529773 61.3768722045 -0.727716646765 0.258175294724 0.0

group cylindingr region cylindin

# Definition of outer cylinder

region cylindout cylinder z 59.4231067727 73.9774529773 31.6879405819 24.5231503881 98.230594021 units lattice rotate v_theta 59.4231067727 73.9774529773 61.3768722045 -0.727716646765 0.258175294724 0.0

group cylindcrystal region cylindout

# Define group of atoms around the cylinder, to be deleted

group waste subtract crystal cylindcrystal

# Delete atoms around the cylinder

delete_atoms group waste

# Define the cylindrical thermostat

group thermocylind subtract cylindcrystal cylindingr

The results I obtain with this cylindrical box differ a lot from what I obtain with the cubic box (which are correct). For instance, the number of defects does not reach a constant value after a long time (eg 20 ps), as it should and as it does in the cubic box. Using VMD I could see that the defects always reach the border of the cylinder, even though the cylinder is very large. In the cubic box this does not occur. Defects stay confined in a certain region and do not reach borders.
I tried large cylinders, thinking it could be due to some border effects. But it is not that. I observe the same results even for large cylinders.

Any idea why I obtain very different results when I use a cylindrical box ?
Any help will be appreciated.

With best regards,
Christophe

What happens if you use a cylinder with an active (inner)

region larger than your original box? What happens

if you don’t use a periodic box, but use a 2-layer box (like

the cylinder) where there is an inner active region

and an outer (frozen?) region? Periodic boundaries

are not the same as non-periodic.

Steve

Hi Steve,

What happens if you use a cylinder with an active (inner)
region larger than your original box?

Instead of a cylinder, I created a block with region command inside the
original box. In principle simpler. When it arrives at the first run of
Voronoi: segmentation fault.

compute v1 all voronoi/atom occupation
compute r0 all reduce sum c_v1[1]
compute r1 all reduce sum c_v1[2]

thermo_style custom c_r0 c_r1
run 0

Hi Steve,

What happens if you use a cylinder with an active (inner)
region larger than your original box?

Instead of a cylinder, I created a block with region command inside the
original box. In principle simpler. When it arrives at the first run of
Voronoi: segmentation fault.

which version of LAMMPS exactly is this with?
a bunch of possible segfault reasons, particularly for systems like
you describe were addressed in some of the recent patches.

axel.

>
> Hi Steve,
>
> Thank you for your email.
>
>>
>> What happens if you use a cylinder with an active (inner)
>> region larger than your original box?
>
>
>
> Instead of a cylinder, I created a block with region command inside the
> original box. In principle simpler. When it arrives at the first run of
> Voronoi: segmentation fault.

which version of LAMMPS exactly is this with?
a bunch of possible segfault reasons, particularly for systems like
you describe were addressed in some of the recent patches.

I am using version 16 Feb 16.

Christophe

>
> Hi Steve,
>
> Thank you for your email.
>
>>
>> What happens if you use a cylinder with an active (inner)
>> region larger than your original box?
>
>
>
> Instead of a cylinder, I created a block with region command inside the
> original box. In principle simpler. When it arrives at the first run of
> Voronoi: segmentation fault.

which version of LAMMPS exactly is this with?
a bunch of possible segfault reasons, particularly for systems like
you describe were addressed in some of the recent patches.

I am using version 16 Feb 16.

that is too old for the bugfixes i was referring to. please try the
very latest, 14 Jun 2016 (as of now).

axel.

Hi Steve,

What happens if you use a cylinder with an active (inner)
region larger than your original box?

Instead of a cylinder, I created a block with region command inside the
original box. In principle simpler. When it arrives at the first run of
Voronoi: segmentation fault.

which version of LAMMPS exactly is this with?
a bunch of possible segfault reasons, particularly for systems like
you describe were addressed in some of the recent patches.

I am using version 16 Feb 16.

that is too old for the bugfixes i was referring to. please try the
very latest, 14 Jun 2016 (as of now).

Ok, thanks for the advice. I will try and let you know if it solves the problem.

Christophe

Hi all,

As Axel suggested, I installed version 14 May 2016.
The segmentation fault that appeared at the first run of the compute voronoi occupation is now gone.
However, I still have the problem of some defects that appear at the faces of the block, which is not expected.

In summary, what I do is the following:
I define a large box with PBC (each face). Then, inside this box, I define a block region with a rotation. See fig1 in attachment. Then, I delete all the atoms (grey points) around this block (blue).
My goal I study the formation of vacancies and interstitials inside this block during a collision cascade, reducing this way the number of atoms.

The result of the cascade inside the block is shown in fig2. As you can see, some vacancies (white balls) and interstitials (red balls: dumbbells) appear on the faces of the block, which is not expected. This occurs even for large size of the block. Thus, this must be an artefact since defects should appear only in a limited zone in the center (more or less) of the block.
As Steve commented few emails ago, my system is no longer periodic, which could affect the results. Thus, if defects appear on the faces of the block, it must be due to a combination of non-periodic conditions + Voronoi package results.
If I run the same cascade in a block with periodic conditions, there is no defects on the faces of the box.
My question thus is: once I have defined a block with a rotation inside the original box, is it possible to make its faces with periodic conditions ?
Or is there another way–other than region block rotate command–to create a block box with a rotation and with periodic BC from the beginning ?

Many thanks in advance and saludos !
Christophe

Hi all,

As Axel suggested, I installed version 14 May 2016.
The segmentation fault that appeared at the first run of the compute voronoi
occupation is now gone.
However, I still have the problem of some defects that appear at the faces
of the block, which is not expected.

if you have open surfaces, reconstruction is quite expected. this also
has a significant effect on phonons in your system. in short, you are
doing a very different simulation and there is no reason to assume,
that it should produce the same results as a bulk system. that line of
thinking may work to some degree on a macroscopic level, but it
definitely has no place in atomistic simulations.

In summary, what I do is the following:
I define a large box with PBC (each face). Then, inside this box, I define a
block region with a rotation. See fig1 in attachment. Then, I delete all the
atoms (grey points) around this block (blue).
My goal I study the formation of vacancies and interstitials inside this
block during a collision cascade, reducing this way the number of atoms.
The result of the cascade inside the block is shown in fig2. As you can see,
some vacancies (white balls) and interstitials (red balls: dumbbells) appear
on the faces of the block, which is not expected. This occurs even for large
size of the block. Thus, this must be an artefact since defects should
appear only in a limited zone in the center (more or less) of the block.
As Steve commented few emails ago, my system is no longer periodic, which
could affect the results. Thus, if defects appear on the faces of the block,
it must be due to a combination of non-periodic conditions + Voronoi package
results.

voronoi analysis doesn't affect the trajectory of atoms.

If I run the same cascade in a block with periodic conditions, there is no
defects on the faces of the box.
My question thus is: once I have defined a block with a rotation inside the
original box, is it possible to make its faces with periodic conditions ?

if you make that periodic again, it will be like your original system,
only smaller.

Or is there another way--other than region block rotate command--to create a
block box with a rotation and with periodic BC from the beginning ?

i think that your overall approach is very misguided at a fundamental
level. you are spending a lot of effort trying to save some CPU time
and are disregarding the physics (and even more so the facts of
spatial geometry) to a degree that is quite scary. your simulations
clearly tell you, that what you are trying to do, i.e. cutting away
atoms, has a significant impact and thus is not applicable.

have you made an estimate of how much simulation time you would be
able save compared to how much time you have already spent on this?

axel.

Hi all,

As Axel suggested, I installed version 14 May 2016.
The segmentation fault that appeared at the first run of the compute voronoi
occupation is now gone.
However, I still have the problem of some defects that appear at the faces
of the block, which is not expected.

if you have open surfaces, reconstruction is quite expected. this also
has a significant effect on phonons in your system.

Good to know.

in short, you are
doing a very different simulation and there is no reason to assume,
that it should produce the same results as a bulk system. that line of
thinking may work to some degree on a macroscopic level, but it
definitely has no place in atomistic simulations.

In summary, what I do is the following:
I define a large box with PBC (each face). Then, inside this box, I define a
block region with a rotation. See fig1 in attachment. Then, I delete all the
atoms (grey points) around this block (blue).
My goal I study the formation of vacancies and interstitials inside this
block during a collision cascade, reducing this way the number of atoms.
The result of the cascade inside the block is shown in fig2. As you can see,
some vacancies (white balls) and interstitials (red balls: dumbbells) appear
on the faces of the block, which is not expected. This occurs even for large
size of the block. Thus, this must be an artefact since defects should
appear only in a limited zone in the center (more or less) of the block.
As Steve commented few emails ago, my system is no longer periodic, which
could affect the results. Thus, if defects appear on the faces of the block,
it must be due to a combination of non-periodic conditions + Voronoi package
results.

voronoi analysis doesn’t affect the trajectory of atoms.

If I run the same cascade in a block with periodic conditions, there is no
defects on the faces of the box.
My question thus is: once I have defined a block with a rotation inside the
original box, is it possible to make its faces with periodic conditions ?

if you make that periodic again, it will be like your original system,
only smaller.

And how can this be done in LAMMPS ? Can I make the block periodic again ?

Or is there another way–other than region block rotate command–to create a
block box with a rotation and with periodic BC from the beginning ?

i think that your overall approach is very misguided at a fundamental
level. you are spending a lot of effort trying to save some CPU time
and are disregarding the physics (and even more so the facts of
spatial geometry) to a degree that is quite scary.

I am not disregarding the physics. Besides this idea, I have tested others and each time I have compared the method to full MD simulations (number of Frenkel pairs vs time, number of defects in clusters, size distribution of clusters…). I hope to be able to publish something by the end of the year.

clearly tell you, that what you are trying to do, i.e. cutting away
atoms, has a significant impact and thus is not applicable.

Yes, I do agree with you. It was an idea I had, I tried it and I found out that it has, as you say, a significant impact on the result. Not everything always works. I am just trying to figure out why it does not work and if I can improve it.

have you made an estimate of how much simulation time you would be
able save compared to how much time you have already spent on this?

I have dedicated something like 10 days. Not big deal. The case I showed is very small and of course, the CPU time I could save in this case is insignificant. My real goal is to simulate high-energy PKA (200 keV - 1 MeV). In general, such simulations require thousands of hours on a supercomputer. If we could reduce the CPU time, don’t you think it would worth it ?

Christophe

I’m not following all the details of this thread.
However, one thing to note. Compute voronoi
should not affect the results at all. It is simply
a diagnostic. It should also work fine at
surfaces, although you obviously may get large
Voronoi volumes if the volume extends from a free
surface into open space.

Periodic vs non-periodic can obviously give
different results. If you put a cylinder in a periodic
box (don’t see a good reason to do that), you should
make sure the surface of the cylinder is far enough
from the box wall to avoid interactions with another
image of the cylinder. Other than that there is nothing
special about a cylinder vs block of material
from the LAMMPS perspective.

Steve

if you make that periodic again, it will be like your original system,
only smaller.

And how can this be done in LAMMPS ? Can I make the block periodic again ?

you have it backwards. recall your physics classes when you were
taught about frame of reference. then you should easily realize, that
converting your rotated block to a periodic system is nothing else but
and impact vectors.

now look up a little bit theory of crystallography and solid state
physics and also recall that you cannot rotate to arbitrary
orientations, but the resulting periodicity of that rotated lattice
has to be commensurate with the box dimensions. otherwise you would
create even worse artifacts in those areas where the periodic
continuation is supposed to occur.

> Or is there another way--other than region block rotate command--to
> create a
> block box with a rotation and with periodic BC from the beginning ?

i think that your overall approach is very misguided at a fundamental
level. you are spending a lot of effort trying to save some CPU time
and are disregarding the physics (and even more so the facts of
spatial geometry) to a degree that is quite scary.

I am not disregarding the physics. Besides this idea, I have tested others
and each time I have compared the method to full MD simulations (number of
Frenkel pairs vs time, number of defects in clusters, size distribution of
clusters...). I hope to be able to publish something by the end of the year.

clearly tell you, that what you are trying to do, i.e. cutting away
atoms, has a significant impact and thus is not applicable.

Yes, I do agree with you. It was an idea I had, I tried it and I found out
that it has, as you say, a significant impact on the result. Not everything
always works. I am just trying to figure out why it does not work and if I
can improve it.

it cannot work because you are neglecting rather fundamental issues of
physics, geometry and solid state physics.
neither periodic boundary conditions nor surface effects in atomistic
simulations are a new thing. there is a large body of research on
that, not to mention text book material.

now from a simple consideration of the physics of your system and from
what i recall of methodology of reducing the CPU time of simulations,
i see two possible approaches:

1) an isolated system with suitably chosen "stochastic boundaries"
that will avoid reconstruction and minimize artifacts from surfaces.
this is quite commonly done in slab calculations to have a finite size
slab, but mimic a very large body underneath.
2) an implementation of r-RESPA multi-timestepping, that recomputes
forces for the region(s) with fast modifications more frequent than
others.

for 2) we already have an implementation in LAMMPS, but that would
only work properly with pair-wise additive potentials. for your
material you are most likely using a manybody potential, thus
implementing a compatible multi-timestep scheme would be rather tricky
and require deep knowledge of LAMMPS and good programming skills.

1) will also have some artifacts, but it is much less than what you
are currently getting. what you need to do is: increase your region of
interest by 4-5 layers of atoms in each direction. then run a
dissipative thermostat with stronger than usual coupling on the 2-3
inner layers of those added layers and immobilize the outer 2 layers
at the equilibrated bulk lattice positions. immobilizing the outer
layers will avoid all reconstruction at the surface, impact of surface
tension, reflection of phonons and so on, and the dissipative
thermalization should mimic the thermal coupling with a bulk system
despite being next to immobile (i.e. 0K) atoms and should (help to)
squash incoming phonons.

axel.

if you make that periodic again, it will be like your original system,
only smaller.

And how can this be done in LAMMPS ? Can I make the block periodic again ?

you have it backwards. recall your physics classes when you were
taught about frame of reference. then you should easily realize, that
converting your rotated block to a periodic system is nothing else but
and impact vectors.

now look up a little bit theory of crystallography and solid state
physics and also recall that you cannot rotate to arbitrary
orientations, but the resulting periodicity of that rotated lattice
has to be commensurate with the box dimensions. otherwise you would
create even worse artifacts in those areas where the periodic
continuation is supposed to occur.

I confess that I always had some difficulties with crystallography but I see what you mean. I wrongly assumed that I could make periodic again any shape.

Or is there another way–other than region block rotate command–to
create a
block box with a rotation and with periodic BC from the beginning ?

i think that your overall approach is very misguided at a fundamental
level. you are spending a lot of effort trying to save some CPU time
and are disregarding the physics (and even more so the facts of
spatial geometry) to a degree that is quite scary.

I am not disregarding the physics. Besides this idea, I have tested others
and each time I have compared the method to full MD simulations (number of
Frenkel pairs vs time, number of defects in clusters, size distribution of
clusters…). I hope to be able to publish something by the end of the year.

clearly tell you, that what you are trying to do, i.e. cutting away
atoms, has a significant impact and thus is not applicable.

Yes, I do agree with you. It was an idea I had, I tried it and I found out
that it has, as you say, a significant impact on the result. Not everything
always works. I am just trying to figure out why it does not work and if I
can improve it.

it cannot work because you are neglecting rather fundamental issues of
physics, geometry and solid state physics.
neither periodic boundary conditions nor surface effects in atomistic
simulations are a new thing. there is a large body of research on
that, not to mention text book material.

now from a simple consideration of the physics of your system and from
what i recall of methodology of reducing the CPU time of simulations,
i see two possible approaches:

1. an isolated system with suitably chosen “stochastic boundaries”
that will avoid reconstruction and minimize artifacts from surfaces.
this is quite commonly done in slab calculations to have a finite size
slab, but mimic a very large body underneath.
2. an implementation of r-RESPA multi-timestepping, that recomputes
forces for the region(s) with fast modifications more frequent than
others.

for 2) we already have an implementation in LAMMPS, but that would
only work properly with pair-wise additive potentials. for your
material you are most likely using a manybody potential, thus
implementing a compatible multi-timestep scheme would be rather tricky
and require deep knowledge of LAMMPS and good programming skills.

I am using pair_style eam/fs (Ackland04 or Mendelev07). I mainly work on Fe and W.

1. will also have some artifacts, but it is much less than what you
are currently getting. what you need to do is: increase your region of
interest by 4-5 layers of atoms in each direction. then run a
dissipative thermostat with stronger than usual coupling on the 2-3
inner layers of those added layers and immobilize the outer 2 layers
at the equilibrated bulk lattice positions. immobilizing the outer
layers will avoid all reconstruction at the surface, impact of surface
tension, reflection of phonons and so on, and the dissipative
thermalization should mimic the thermal coupling with a bulk system
despite being next to immobile (i.e. 0K) atoms and should (help to)
squash incoming phonons.

Before I tested this idea with the rotated block/cylinder, I was using something similar. The only difference is that the outer layer was not so elaborated. I was using a block with periodic conditions in each direction and between 5-10 extra-layers around the region of interest in each direction. From this extra-layers, two were the thermostat (Langevin) with a damping parameter set to 0.1.
I will try what you suggest.

Christophe

Hi Axel,

I tried what you suggested, ie, an inner thermostat (2 layers) strongly dissipative at room temp and an outer thermostat (2 layers) at 0 K.

fix thermostatout thermal_atoms_out temp/berendsen 0.0 0.0 0.1

fix thermostatin thermal_atoms_in temp/berendsen 300.0 300.0 0.1

You were right, it works pretty well. See picture in attachment. The two thermostats are represented with two different blue tones. Red balls = SIA and brown balls = Vacancies. With these thermostats there is almost no defects that form at the surface of the block (with non-periodic BC) due to surface reconstruction, as you expected. Or maybe once in a while, 1 or 2 sometimes, but much less than before.

Do you think the way I defined the thermostats is ok ? Or do you think it can be improved ?

With best regards,
Christophe

Hi Axel,

I tried what you suggested, ie, an inner thermostat (2 layers) strongly
dissipative at room temp and an outer thermostat (2 layers) at 0 K.

​fix temp/berendsen is *not* a dissipative thermostat. on the contrary. try
fix langevin or fix temp/csld of fix gld or fix gle.

fix thermostatout thermal_atoms_out temp/berendsen 0.0 0.0
0.1
fix thermostatin thermal_atoms_in temp/berendsen 300.0 300.0
0.1

​there is no need to have a thermostat for maintaining 0K. simply do not do
time integration for those atoms. or use a combination of setting
velocities and and forces to zero for those atoms, if you cannot for some
reason avoid time integration.

You were right, it works pretty well. See picture in attachment. The two
thermostats are represented with two different blue tones. Red balls = SIA
and brown balls = Vacancies. With these thermostats there is almost no
defects that form at the surface of the block (with non-periodic BC) due to
surface reconstruction, as you expected. Or maybe once in a while, 1 or 2
sometimes, but much less than before.

​there *cannot* be surface reconstructions or surface defects, if the
surface atoms cannot move. ​

Do you think the way I defined the thermostats is ok ? Or do you think it
can be improved ?

​it can be improved. see my comments above.

axel.​