different energies in "energy final" and the last line of the log file

Konstantin_Khromov · August 5, 2016, 12:57pm

Dear lammps users,

I'm trying to calculate vacancy formation energy on an edge dislocation
core in FCC Al under uniaxial tensile and compressive stresses.

To this end I first relax atomic positions for a dislocation without a
vacancy, setting uniaxial stress to +-50 MPa along a Burgers vector.

The log file lists potential energy, stresses, simulation volume size -
one line per ten minimization steps. After minimization is finished
"energy final" is printed as well.

In all the lammps simulations I've made so far the two values: potential
energy at the last line of the corresponding log list and the "energy
final" values coincided, or at most differed by very insignificant
amount less that 0.01 eV. This is not the case in my current dislocation
under pressure modelling.

Here is the table:

Py "last line" Energy "energy final" Energy

-50 MPa -152265.06 -152261.660793
0 MPa -152265.12 -152265.11734
+50 MPa -152265.03 -152268.70936

My first idea looking at this table was: taking into account the facts that
1. for all pressures "last line" Energies are very close,
2 for zero pressure "last line" and "energy final" Energies are very close,
3. for +- 50 MPa "energy final" energies are shifted by about the same
amount but in different directions,
one may suppose that "last line" energy does not include (Pressure) x
(Volume change) term and "energy final" does.

But unfortunately such supposition seems not to be true. I made some
further experiments:
1. The same +- 50 MPa for ideal FCC Al crystal without dislocations with
the same simulation volume size. In this case "last line" and "energy
final" are very close.
2. Deliberately setting somewhat wrong lattice constant thus introducing
extra pressure, "last line" and "energy final" are the same again.

So, can somebody give me a clue, what might this discrepancy between
"last line" and "energy final" mean?

Konstantin

============================= -50 MPa input====================
# Input file for Vacancy Formation Energy

#clear
#units metal
#dimension 3
#boundary p p p
#atom_style atomic

read_restart tmp.restart.230
# ------------------------ FORCE FIELDS -----------------------
  pair_style eam/alloy
  pair_coeff * * ../../Al99.eam.alloy Al
#---------------------------Settings----------------------------
  compute csym all centro/atom fcc
  compute eng all pe/atom
  compute eatoms all reduce sum c_eng
#fix 1 all box/relax x 0.0 y 0.0 z 0.0 vmax 0.001
fix 1 all box/relax x 0.0 y -500.0 z 0.0 vmax 0.001

reset_timestep 0

thermo 10
thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms

dump 1 all custom 1 dump.relax.2.* id type x y z c_csym c_eng fx fy fz
#dump 1 all custom 400 dump.relax.2.* x y z

#min_style sd
#min_style hftn
min_style cg
minimize 1e-25 1e-25 500000 500000

=============================== -50 MPa log ==========================
LAMMPS (14 May 2016)
# Input file for Vacancy Formation Energy

#clear
#units metal
#dimension 3
#boundary p p p
#atom_style atomic

read_restart tmp.restart.230
   orthogonal box = (0 0 0) to (74.4033 91.641 112.237)
   1 by 2 by 2 MPI processor grid
   45360 atoms
# ------------------------ FORCE FIELDS -----------------------
  pair_style eam/alloy
  pair_coeff * * ../../Al99.eam.alloy Al
#---------------------------Settings----------------------------
  compute csym all centro/atom fcc
  compute eng all pe/atom
  compute eatoms all reduce sum c_eng
#fix 1 all box/relax x 0.0 y 0.0 z 0.0 vmax 0.001
fix 1 all box/relax x 0.0 y -500.0 z 0.0 vmax 0.001

reset_timestep 0

thermo 10
thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms

dump 1 all custom 1 dump.relax.2.* id type x y z c_csym c_eng fx fy fz
#dump 1 all custom 400 dump.relax.2.* x y z

#min_style sd
#min_style hftn
min_style cg
minimize 1e-25 1e-25 500000 500000
WARNING: Resetting reneighboring criteria during minimization
(../min.cpp:168)
Neighbor list info ...
   2 neighbor list requests
   update every 1 steps, delay 0 steps, check yes
   max neighbors/atom: 2000, page size: 100000
   master list distance cutoff = 8.28721
   ghost atom cutoff = 8.28721
   binsize = 4.1436 -> bins = 18 23 28
Memory usage per processor = 10.3951 Mbytes
Step PotEng Lx Ly Lz Press Pxx Pyy Pzz eatoms
        0 -152200.55 74.403251 91.641039 112.23689
-11855.853 -9455.8598 -17988.508 -8123.1911 -152200.55
       10 -152260.06 74.1996 90.724629 112.03954
-419.568 842.70847 -3316.4051 1214.9927 -152260.06
       20 -152264.86 74.327298 90.338085 112.27821
-85.131029 85.731173 -426.77391 85.649645 -152264.86
       30 -152264.87 74.330271 90.339113 112.28318
-163.50448 5.5088942 -497.31138 1.2890504 -152264.87
       40 -152264.88 74.330168 90.338961 112.28332
-163.50159 6.9628234 -498.8159 1.3482988 -152264.88
       50 -152264.9 74.330706 90.338522 112.28306
-163.59221 2.0970711 -495.51682 2.6431209 -152264.9
.................

     7380 -152265.06 74.335704 90.331457 112.28355
-164.3806 1.0867458 -497.95564 3.7270898 -152265.06
     7390 -152265.06 74.335605 90.331148 112.28393
-163.30988 2.8165779 -495.00709 2.2608655 -152265.06
     7400 -152265.06 74.335894 90.33112 112.28356
-163.6151 0.1080608 -495.07686 4.1234828 -152265.06
     7410 -152265.06 74.33566 90.331197 112.28381
-163.56737 2.1707065 -495.4881 2.6152899 -152265.06
     7420 -152265.06 74.33664 90.330737 112.28303
-164.61913 -6.9254688 -493.31943 6.3875139 -152265.06
     7430 -152265.06 74.335668 90.331172 112.28381
-163.44074 2.2497751 -495.23289 2.660889 -152265.06
     7440 -152265.06 74.335551 90.331632 112.28385
-166.43118 0.98132387 -500.53198 0.25712749 -152265.06
     7446 -152265.06 74.335678 90.331182 112.2838
-163.5457 2.1119319 -495.33747 2.5884246 -152265.06
Loop time of 909.591 on 4 procs for 7446 steps with 45360 atoms

95.5% CPU use with 4 MPI tasks x no OpenMP threads

Minimization stats:
   Stopping criterion = linesearch alpha is zero
   Energy initial, next-to-last, final =
          -152200.54572 -152261.660793 -152261.660793
   Force two-norm initial, final = 10258.3 0.062062
   Force max component initial, final = 8353.34 0.0380085
   Final line search alpha, max atom move = 2.12807e-06 8.08847e-08
   Iterations, force evaluations = 7446 14974

athomps · August 5, 2016, 4:17pm

You are comparing potential energy to minimization energy. If you read the doc page for minimize, you will see that these are not the same thing. Specifically, when you use fix box/relax, you are are adding an extra energy term(s) given on the doc page for fix box relax. Finally, if you look at the LAMMPS example where fix box/relax is demonstrated, examples/in.min.box, you will learn how to print out the full minimization energy in thermo output. In general, when you encounter something that seems inexplicable, you should first check to see if it is explained in the manual. I will update the note in fix box/relax explaining this issue more explicitly, as you are not the first person to get confused by this.

"This fix computes a global scalar which can be accessed by various output commands. The scalar is given by the energy expression shown above. The energy values reported at the end of a minimization run under “Minimization stats” include this energy, and so differ from what LAMMPS normally reports as potential energy. This fix does not support the fix_modify energy yes keyword, because that would result in double-counting of the fix energy in the minimization energy. Instead, the fix energy can be explicitly added to the potential energy using a variable command e.g. “variable emin equal pe+f_1” or “variable emin equal pe+f_1/atoms”

Aidan