constant (incorrect) c/a ratio during energy minimization for hcp Ti




Hello everyone, I seem to be running into a problem while performing relaxations by energy minimization in LAMMPS for hcp Ti. I am sure this is not a bug or software related problem, it is most likely something wrong with my LAMMPS script and interpretation. For this structure I am expecting a c/a of about 1.58 (for the 6 atom unit cell), however in my efforts to yield this value in simulation, I get a ratio that is substantially higher and relatively constant throughout energy minimization. Being able to yield the proper lattice constants will be a great help for me moving forward on bigger projects, and I am greatly appreciative of any help you guys might be able to provide.

-Log file:

units metal
boundary p p p
atom_style atomic

lattice hcp 2.9508
region box block 0 4 0 4 0 4 create_box 1 box
create_atoms 1 box
mass 1 47.867

pair_style meam/spline
pair_coeff * * Ti.meam.spline Ti
neighbor 2.0 bin
neigh_modify delay 10 check yes

reset_timestep 0
fix 1 all box/relax iso 0 vmax 0.01
thermo 1
thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms
dump 1 all cfg 1 dump.emin_*.cfg mass type xs ys zs c_pe
dump_modify 1 element Ti
min_style cg
minimize 0 0 5000 10000

(WARNING: Resetting reneighboring criteria during minimization (…/min.cpp:173) Memory usage per processor = 4.62926 Mbytes)

Loop time of 0.2021 on 8 procs for 7 steps with 256 atoms

Minimization stats:
Stopping criterion = linesearch alpha is zero
Energy initial, next-to-last, final =
-1233.52693295 -1236.22330454 -1236.22330454
Force two-norm initial, final = 372.841 0.00897025
Force max component initial, final = 372.841 0.00897025
Final line search alpha, max atom move = 0.00012207 1.095e-06
Iterations, force evaluations = 7 53
Lx, Ly, Lz initial:
11.8032 20.443742 19.274545
Lx, Ly, Lz final:
11.637455 20.1566664 19.003885


When I visualize the atoms (atomeye), it is clear that the output ‘Ly’ runs along the ‘a’ lattice vector direction (of a 6 atom HCP unit cell). And with a box that is 4 unit cells long in each direction, I assume that ‘Lx’ is 8 ‘a’ lattice vectors long in this case (since it takes two lattice vectors to traverse across the floor of an HCP unit cell). Also with the visualization, I find that ‘Lz’ reflects the c lattice vector direction, being 4 ‘c’ vectors long in this case. By this convention I get a c/a of about 1.8856 throughout the simulation.
Now of course with the LAMMPS defined 4 atom hcp unit cell, ‘Lx’ represents the a direction instead of ‘Ly’, and by this convention I get a c/a ratio of about 1.63 throughout the simulation (the LAMMPS defined ratio for hcp).
Aside from this trouble of converting from one unit cell orientation to the other, I guess my original question boils down to this: with the above input script, and overall abilities of energy minimization in LAMMPS, should I expect the c/a ratio to decrease from the default 1.63 value to the expected value of about 1.58? I have tried different energy tolerances and although the iteration count changes, I seem to still get a rather constant c/a for either unit cell orientation.
Also, in regards to the unit cell orientations, which value (1.8856, 1.63) should I take as the true c/a ratio when comparing to the expected value of 1.58 (for the 6 atom unit cell)? I am aware I could use a custom lattice to layout the exact unit cell I would expect, but I am hoping this exercise will both straighten out my confusion using the default hcp lattice, and also illustrate the use of energy minimization in relaxations.
Any advice in regards to either of these questions would be greatly appreciated and very helpful to me moving forward. Thank you very much to anybody taking time to read this.