geometry minimization of a highly skewed periodic system

LAMMPS LAMMPS Mailing List Mirror
1

Hello,

I am attempting to optimize the structure of a 3d crystal whose lattice parameters are (say) a=b=c=10.0 ang and angles of alpha=beta=90 and gamma=120.0 deg using fix minimize and fix box/relax commands. In such a case the tilt factor xy = bcos(gamma) is exactly (xhi-xlo)/2.
During the minimization the tilt factor sometimes exceeds the -(xhi-xlo)/2 < xy < (xhi -xlo)/2 limit and lammps give the expected “ERROR: Triclinic box skew is too large (domain.cpp:144)” error.

Currently, I am using multiple unitcells in order circumvent this error - for the xy tilt use multiple unitcells in the x direction to get a larger xhi value. I’d really not want to continue doing this.

More importantly, I am interested in geometry optimizations of systems that are more skewed that gamma = 120 deg. For example a crystal with a=b=c and alpha=beta=gamma=60 deg.

Is there a more elegant solution to this? Or am I missing something very trivial? Is there a way to bypass the tilt limitations?

Thank You,

Regards,
Ambar

2

First, see section 6.12 of the manual. The limitation on
the amount of skew is not a restriction, e.g. your 120 degree
question. You can always reformulate the shape of the
box into an identical system with a skew within the limits.

If the skew is changing during the minimization we could
add a fudge factor to allow titls to exceed the maximum, That
already occurs for fix npt (dynamics). I thought Aidan might
have added it to fix box/relax, or has thought about doing
it, so maybe he wants to comment on your min question.

Steve

3

Firstly, thank you for the reply.
I have already read through the section 6.12, but I guess I do not understand how a given unit cell can give different xhi,yhi,zhi,xy,xz,yz values.

Please let me know where my logic is incorrect.

  1. Say I have a crystal with a=b=c and alpha=beta=gamma = 60 deg
  2. I CHOOSE xlo,ylo and zlo to be zero
  3. From the formulas of section 6.12 lx = a
  4. xy (=b/2=a/2) and xz are defined from beta and gamma
  5. This give ly, yz and lz.
  6. I do understand xy can be defined equivalently as xy + n*(xhi-xlo)/2, where n = integer, but that does not ensure that xy will now like between +/- (xhi-xlo)/2

Am I missing something obvious here? Should I rotate the unitcell about an axis such that ‘a direction’ DOES NOT lie along the x direction anymore? I.E. choose such that lx ~= a.

Thank you,

Regards,
Ambar

4

6. I do understand xy can be defined equivalently as xy + n*(xhi-xlo)/2, where n = integer, but that >does not ensure that xy will now like between +/- (xhi-xlo)/2

What's an example where you can't do this? Perhaps if
the xdim is much larger than the ydim? In that case
you could simply flip the x,y dimensions in your crystal
definition.

Steve

5

Hello,

Thanks for the reply.

In a trclinic system with a ~=b, ~=c I understand that the cell should be oriented such that the largest lattice constant lies in the x direction while the smallest lies in the z direction. Given any unit cell one can always rotate it by n*90 deg about the x,y and z axis to find such an orientation.

My question was about a system in which a = b = c and alpha=beta=gamma=(say) 50 deg.

In this case rotation of the unit cell does not matter because all the lattice constants are equal and one of them has to be oriented about the x direction.

Assuming a= b =c = 1;

lx = a = 1;

xy = b*cos(gamma) = 0.64; >0.5

xz = c*cos(beta) = 0.64; > 0.5

etc.

What should be done in this case?

Also, I know that a highly skewed system will be inefficient,.but is there a hack that will prevent this error from occurring? Perhaps just commenting out the relevant lines in domain.cpp?

I have been thinking about ways to overcome this but I am all out of ideas…

Thanks,
Ambar

6

I'll talk to Aidan about this (he may be out of town thru the end
of the month). There should be a way to allow this pretty
easily. However, I believe what you are calling the unit
cell is only a matter of convention. If a = b = c, then a tilt
of 0.64 is exactly equivalent to a unit cell with a tilt of -0.36
is it not?

Steve

7

Hello,
Currently I have just commented out the relevant lines in domain.cpp.

8

Those changes are probably fine for both Qs, so long as you don't enable
other commands, like fix npt or fix box/relax that will do additional checking.

Was I correct that you can reformulate your system with tilt factors less
than 0.5?

Steve

9

Hi Steve,
I might be wrong, but I don’t think tilts of of 0.64 and -0.36 are equivalent (even for the simplest case of a=b=c=1, alpha = 50; beta = gamma = 90 deg)
I have attached a simple matlab script that converts a,b,c,alpha,beta,gamma to xhi,yhi,zhi,xlo,ylo,zlo,xy,xz,yz and vice versa using mode = 1 or 2.

Here is the output.

script_lammps_triclinic.m (3.18 KB)

10

I'm not seeing this geometrically. From section 6.12 of the manual:

To avoid extremely tilted boxes (which would be computationally
inefficient), no tilt factor can skew the box more than half the
distance of the parallel box length, which is the 1st dimension in the
tilt factor (x for xz). For example, if xlo = 2 and xhi = 12, then the
x box length is 10 and the xy tilt factor must be between -5 and 5.
Similarly, both xz and yz must be between -(xhi-xlo)/2 and
+(yhi-ylo)/2. Note that this is not a limitation, since if the maximum
tilt factor is 5 (as in this example), then configurations with tilt =
..., -15, -5, 5, 15, 25, ... are geometrically all equivalent.

Is your argument below also true for gamma, which is essentially the
xy tilt factor in the above paragraph?

If the box length is 1, then isn't every config of the system
with tilt factors that vary by 1, equivalent?

E.g. tilt factors of ...,-1.36,-0.36,0.64,1.64, ... are all the same
system, so long as you choose an appropriate unit cell and
basis atoms. By this I mean that you can define a unit cell,
say with tilt factor of 0.64, and there is a different unit cell
with tilt factor -0.36 (different atoms, different positions within
the new unit cell), that produces the exact same atoms and
periodicity of the infinite lattice as does the first unit cell.

There is an additional coupling constraint of yz to xz if
I recall, but I believe that doesn't prevent all 3 tilt factors to
be remapped between -0.5 and 0.5 while preserving an equivalent
system.

What am I missing?

Steve

11

Hi,
I am not sure I understand how different tilts give the same unit cell.
Perhaps a simpler case to analyze will be a 2-d system, a=b=1 and gamma = 50 deg
For this system

12

I am not sure I understand how different tilts give the same unit cell.

They don't. But they do give identical simulations with different unit
cells. There are an infinite number of equivalent unit cells. The -0.5 to 0.5
restriction is saying: pick the unit cell that meets that restriction, b/c it
will be more efficient for LAMMPS. It may not be more conventient for
the user, so Aidan and I will think about relaxing it, when he gets back
around Aug 1.

One way to think about this is as follows. Imagine you replicate
your init cell to form a big simulation box with the same shape.
With corner pts A,B,C,D (in 2d). Then replicate that a few times
in x for periodic BC. So you get somthing like:

H --- A ---- B --- F

      > > >

J --- D ---- C ---- G

The angle ADC does not have to be orghongal if your unit cell is tilted.

You can run the simulation with ABCD as your simulatoin box. But you
could also run with the (more tilted) BFCD as your box, or with HACD
as your box. If you replicate either of those boxes you tile all of
space with the exact same set of atoms, and will run an identical
simulatoin. If AB is of length 1, then those 3 simulation boxes
differ in tilt factor by +/- 1. There are an infinite # of such tilted
boxes, all equivalent, each made of differntly tilted unit cells.
LAMMPS requires you to choose the one that is least tilted.

Thta is what is meant by this info on the doc page:

To avoid extremely tilted boxes (which would be computationally
inefficient), no tilt factor can skew the box more than half the
distance of the parallel box length, which is the 1st dimension in the
tilt factor (x for xz). For example, if xlo = 2 and xhi = 12, then the
x box length is 10 and the xy tilt factor must be between -5 and 5.
Similarly, both xz and yz must be between -(xhi-xlo)/2 and
+(yhi-ylo)/2. Note that this is not a limitation, since if the maximum
tilt factor is 5 (as in this example), then configurations with tilt =
..., -15, -5, 5, 15, 25, ... are geometrically all equivalent.

Steve

13

Hi Steve,
That makes perfect sense.
The reason why I do not want to do this is because I have DFT data for the experimentally reported unit cells.

I have a related question -
If the smallest periodic unit i.e. the primitive unitcell is triclinic (with or without a high tilt), then any other unitcell will have to be larger than the primitive unit cell. Is this correct?
Thank you for the replies, this does clear up my original question.

Finally, it seems that commenting out the relevant lines in domain.cpp is enough for removing the 0.5 restriction. When I looked through fix_box_relax, I did not find an error_check for the tilt. Does this seem about right?

Thank you,

Regards,
Ambarish

14

Comments below.

Steve

Hi Steve,
That makes perfect sense.
The reason why I do not want to do this is because I have DFT data for the
experimentally reported unit cells.

I have a related question -
If the smallest periodic unit i.e. the primitive unitcell is triclinic (with
or without a high tilt), then _any_ other unitcell will _have_ to be larger
than the primitive unit cell. Is this correct?
Thank you for the replies, this does clear up my original question.

There are many unit cells that have the same volume. In my previous
email, all 3 simultaion boxes (and corresponding unit cells) have the
same volume. There are also, of course, other unit cells with larger
volume. (e.,g. double the primitive unit cell)

Finally, it seems that commenting out the relevant lines in domain.cpp is
enough for removing the 0.5 restriction. When I looked through
fix_box_relax, I did not find an error_check for the tilt. Does this seem
about right?

I think that's correct. Aidan may not yet have implemented the restriction
in fix box/relax the same as in fix npt. The conceptual Q is if/how
we want to allow users to specify highly skilled tilt factors, even for static
simluations. It means someone (not you), will abuse it to do something
nutty.

15

Perfect.
Thank you so much.

Regards,
Ambarish

16

final note: fix_box_relax.cpp does not perform box flips, and so there is no check on the amount of skew. This could result in a highly skewed box, but in practice this does not seem to occur.

Aidan

17

Just posted a 9 Aug patch that allow arbitrary tilt values,
either at setup or during dynamics, if you override
the default restrictions.

Users are now free to shoot themselves in the foot by
using extremely skewed simulation boxes.

Steve