strange energy and structure for hcp Zr crystal

Dear all,

I got something strange when doing geometry optimization of crystal Zr of HCP structure using eam potential. I used LAMMPS (15 Jan 2011), and the potential file Cu-Zr_2.eam.fs is available on http://www.ctcms.nist.gov/~cbecker/Cu-Zr.html.

The control file in.eam.zr looks as follows:
units metal
boundary p p p
atom_style atomic
read_data geoin.lmp
pair_style eam/fs
pair_coeff * * /home/eugene/code/lammps/lammps-15Jan11/potentials/Cu-Zr_2.eam.fs Zr
group zr type 1
thermo 1
thermo_style custom step temp pe ke etotal vol
velocity all create 0.0 111121
dump 1 all atom 1 dump.lammpstrj
dump_modify 1 scale no sort id
minimize 1.0e-6 0.001 1000 10000

The input data file for geometry geoin.lmp looks as follows:
  Position data for * system

            2 atoms
            1 atom types
  0.0 3.23191850299478 xlo xhi
  0.0 2.79897057788360 ylo yhi
  0.0 5.14700000000000 zlo zhi
     -1.61587775 0.00000000 0.00000000 xy xz yz

   Atoms

          1 1 0.00005432 1.86598045 1.28675000
          2 1 1.61598648 0.93299021 3.86025000

  Masses

    1 91.2240

The " 1 91.2240" is the last line. Some relevant messages in log files are:

Step Temp PotEng KinEng TotEng Volume
        0 0 -12.935023 0 -12.935023 46.559993
        1 0 -12.935023 0 -12.935023 46.559993
Loop time of 0.000221973 on 1 procs for 1 steps with 2 atoms

Minimization stats:
   Stopping criterion = linesearch alpha is zero
   Energy initial, next-to-last, final =
         -12.9350227629 -12.9350227629 -12.9350227629
   Force two-norm initial, final = 0.000126032 0.000126032
   Force max component initial, final = 7.71812e-05 7.71812e-05
   Final line search alpha, max atom move = 0.5 3.85906e-05
   Iterations, force evaluations = 1 2

What's strange is the unit cell shape in the dump file. The initial (zero) snapshot looks:
ITEM: TIMESTEP

ITEM: NUMBER OF ATOMS
2
ITEM: BOX BOUNDS xy xz yz
-1.61588 3.23192 -1.61588
0 2.79897 0
0 5.147 0
ITEM: ATOMS id type x y z
1 1 5.432e-05 1.86598 1.28675
2 1 1.61599 0.93299 3.86025

Note xlo changes from 0 in the input geometry to negative in dump file. It's also negative in the next snapshot. What's more bizarre is if I expand the cell to 2x2x2 cells containing 16 atoms, which looks like:
  Position data for * system

           16 atoms
            1 atom types
  0.0 6.46383700598955 xlo xhi
  0.0 5.59794115576721 ylo yhi
  0.0 10.2940000000000 zlo zhi
     -3.23175551 0.00000000 0.00000000 xy xz yz

   Atoms

          1 1 0.00005432 1.86598045 1.28675000
          2 1 1.61598648 0.93299021 3.86025000
          3 1 3.23197282 0.25010270 1.28675000
          4 1 4.84790499 -0.68288754 3.86025000
          5 1 0.00005432 4.66495103 1.28675000
          6 1 1.61598648 3.73196079 3.86025000
          7 1 3.23197282 3.04907327 1.28675000
          8 1 4.84790499 2.11608304 3.86025000
          9 1 0.00005432 1.86598045 6.43375000
         10 1 1.61598648 0.93299021 9.00725000
         11 1 3.23197282 0.25010270 6.43375000
         12 1 4.84790499 -0.68288754 9.00725000
         13 1 0.00005432 4.66495103 6.43375000
         14 1 1.61598648 3.73196079 9.00725000
         15 1 3.23197282 3.04907327 6.43375000
         16 1 4.84790499 2.11608304 9.00725000

  Masses

    1 91.2240

then the log file showed that the potential in the first several steps were positive, looks like:
Step Temp PotEng KinEng TotEng Volume
        0 0 237.95627 0 237.95627 372.47994
        1 0 68.051335 0 68.051335 372.47994
        2 0 -9.4961439 0 -9.4961439 372.47994
        3 0 -48.813207 0 -48.813207 372.47994
       ...
       21 0 -100.17936 0 -100.17936 372.47994
Again, xlo in dump file is negative.

Could someone help? Thanks.

Eugene

What's strange is the unit cell shape in the dump file. The initial
(zero) snapshot looks:
ITEM: TIMESTEP

ITEM: NUMBER OF ATOMS
2
ITEM: BOX BOUNDS xy xz yz
-1.61588 3.23192 -1.61588
0 2.79897 0
0 5.147 0
ITEM: ATOMS id type x y z
1 1 5.432e-05 1.86598 1.28675
2 1 1.61599 0.93299 3.86025

Note xlo changes from 0 in the input geometry to negative in dump
file. It's also negative in the next snapshot. What's more bizarre is

this is due to the fact that the xlo in the dump file
and the xlo in the data file are not exactly the same
property. this is a bit of an unfortunate choice, but
the lo/hi values in the dump are "boundaries", i.e.
the maximum/minimum of _any_ position in that direction.

then the log file showed that the potential in the first several steps
were positive, looks like:
Step Temp PotEng KinEng TotEng Volume
0 0 237.95627 0 237.95627 372.47994
1 0 68.051335 0 68.051335 372.47994
2 0 -9.4961439 0 -9.4961439 372.47994
3 0 -48.813207 0 -48.813207 372.47994
...
21 0 -100.17936 0 -100.17936 372.47994
Again, xlo in dump file is negative.

again, that is ok.

axel.

To augment what Axel said, the BOX BOUNDS in
the dump file are described on the dump doc page and
in Section 4.12 of the manual. They are formatted
that way for the convenience of post-processing viz.

Steve

Thanks a lot, Axel and Steve. Now the only worry is the initial potential energy is too high for the large box, a possible indication of unoptimized geometry. I expected it was linearly scaled from the small box. I will look into the geometry.

Eugene

Dear Axel and Steve,

Using the information in Section 4.12, I managed to build orthogonal box for hcp, which has 2 times volume of the primitive one. For this new box the energy converges immediately. Thanks.

Eugene