triclinic box skew is too large for a minimized trigonal crystal

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

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

Unless you use fix box/relax, a minimization doesn't change anything
about the box, lx or xy.


Thanks Steve for the answer. There is indeed a line “fix 3 all box/relax aniso 0.0” before the first minimize in the in.elastic script of Aidan. And lx does change, … but not xy (which is normal in the usual use of Aidan’s example for fcc crystals, but a priori not for trigonal crystals).

Sorry for the noise. I now understand that I should either add xy 0.0 to the fix line or use “fix 3 all box/relax tri 0.0”…