Dear LAMMPS users and developpers,

It seems to me that minimize does not update the tilt factors after a minimization and that it can cause problems:

If I start from a trigonal crystal (e.g. quartz) with lattice parameters and coordinates of the basis atoms taken from an online database (e.g. mincryst) it seems to me I must use a xy tilt that is precisely -0.5*(xhi - xlo).

Then, if I try to use minimize with “boundary p p p” in order to get the equilibrium positions of the ions in the unit cell corresponding to the potential I use (e.g. BKS), I get at the end of the minimization coordinates and lx that have changed by less than one percent (and xy that has not changed at all). Consequently, if I do a write-restart than a read-restart as in the in.elastic script of Aidan, I get “ERROR: triclinic box skew is too large”.

Is that normal that minimize does not attempt to change xy or am I misunderstanding the problem?

(I note that I have read several times the help for region, create_box, section 6.12 of the manual, and various threads on the users list found by searching for “lammps-users triclinic box skew is too large”, but I must admit, I have not always understood everything).

Please find hereafter my modified init.mod to be used with Aidan’s in.elastic, potential.mod (modified for BKS potential) and displace.mod.

I note that I have also tried to read coordinates from a file, but that gives the same problem.

units metal

atom_style charge

dimension 3

boundary p p p

# lattice parameters and basis from Mincryst card 7096

variable a equal 4.923

variable c equal 5.409

variable a1x equal $a/c
variable a2x equal -{a1x}/2

variable a2y equal ${a1x}*.8660254040

lattice custom c origin 0.0 0.0 0.0 &
a1 {a1x} 0 0 &

a2 {a2x} {a2y} 0 &

a3 0 0 1 &

basis 0.468 0.0 0.0 &

basis 0.414 0.264 0.121 &

basis 0.0 0.468 0.6667 &

basis 0.532 0.532 0.3333 &

basis 0.736 0.15 0.7877 &

basis 0.85 0.586 0.4543 &

basis 0.264 0.414 0.5456 &

basis 0.15 0.736 0.8789 &

basis 0.586 0.85 0.2123

variable n equal 1

variable mydx equal n/(1-{a2x}/{a1x})
variable mydxo2 equal -{mydx}/2

region simbox prism 0 ${mydx} 0 $n 0 n {mydxo2} 0. 0. units lattice

create_box 2 simbox

create_atoms 1 box &

basis 1 1 &

basis 2 2 &

basis 3 1 &

basis 4 1 &

basis 5 2 &

basis 6 2 &

basis 7 2 &

basis 8 2 &

basis 9 2

mass 1 28.09

mass 2 16.0

group siliconatoms type 1

group oxygenatoms type 2

set group siliconatoms charge 2.4

set group oxygenatoms charge -1.2

# metal units, elastic constants in GPa

#units metal

variable cfac equal 1.0e-4

variable cunits string GPa

# real units, elastic constants in GPa

#units real

#variable cfac equal 1.01325e-4

#variable cunits string GPa

# Define minimization parameters

variable etol equal 0.0

variable ftol equal 1.0e-5

variable maxiter equal 600

variable maxeval equal 2000

variable dmax equal 1.0e-2