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