[QE-users] Bulk modulus and Total energy MgH2

[RAQUEL YANES RODRIGUEZ]

Dear Quantum Espresso users,

I am trying to obtain the bulk modulus of MgH2 by fitting the total  
energy versus volume curve to the Murnaghan equation of state. MgH2  
has a tetragonal crystal system with P42/mnm symmetry. From the  
literature, it appears that B0 should be around 40-50 GPa. However, I  
am obtaining a value of 5 GPa. The input I am using is as follows:

------------------------------------------------------------------------------------------------

&CONTROL
                        title = 'ES-calculation' ,
                  calculation = 'relax' ,
                 restart_mode = 'from_scratch' ,
                       outdir = './' ,
                   pseudo_dir = './' ,
                   wf_collect = .FALSE. ,
                      disk_io = 'none' ,
  /
  &SYSTEM
                        ibrav = 6,
                            A = 4.00,
                            C = 2.672,
                         nat = 6 ,
                         ntyp = 2 ,
                      ecutwfc = 80 ,
                      ecutrho = 320 ,
                    input_dft = 'PBE' ,
                  occupations = 'smearing',
                     smearing = 'gaussian',
                      degauss = 0.02,
  /
  &ELECTRONS
                     conv_thr = 1.0d-8,
                  mixing_beta = 0.7,
  /
  &IONS
                 ion_dynamics = 'bfgs' ,
             trust_radius_max = 0.2 ,
             trust_radius_min = 1.0d-4 ,
  /
ATOMIC_SPECIES
     Mg   24.3050  Mg.pbe-spnl-kjpaw_psl.1.0.0.UPF
     H    1.00000  H.pbe-kjpaw_psl.1.0.0.UPF
ATOMIC_POSITIONS crystal
  H      0.19520000000000      0.80480000000000      0.50000000000000
  H      0.69520000000000      0.69520000000000     -0.00000000000000
  H      0.80480000000000      0.19520000000000      0.50000000000000
K_POINTS automatic
   3 3 5   1 1 1

------------------------------------------------------------------------------------------------

The output of this example ends like this:

!    total energy              =      -6.61308602 Ry
      Harris-Foulkes estimate   =      -6.61308602 Ry
      estimated scf accuracy    <          5.1E-12 Ry

      total all-electron energy =        -6.613060 Ry

      The total energy is the sum of the following terms:

      one-electron contribution =      -0.79116562 Ry
      hartree contribution      =       1.35759033 Ry
      xc contribution           =      -3.89869029 Ry
      ewald contribution        =      -3.27197209 Ry
      one-center paw contrib.   =      -0.00084285 Ry
      smearing contrib. (-TS)   =      -0.00800550 Ry

      convergence has been achieved in   5 iterations

      Forces acting on atoms (cartesian axes, Ry/au):

      atom    1 type  2   force =    -0.00001212    0.00001212    0.00000000
      atom    2 type  2   force =    -0.00000576   -0.00000576    0.00000000
      atom    3 type  2   force =     0.00001212   -0.00001212    0.00000000
      atom    4 type  2   force =     0.00000576    0.00000576    0.00000000
      atom    5 type  1   force =     0.00000000    0.00000000    0.00000000
      atom    6 type  1   force =     0.00000000    0.00000000    0.00000000

      Total force =     0.000027     Total SCF correction =     0.000001

      bfgs converged in  48 scf cycles and  47 bfgs steps
      (criteria: energy <  1.0E-04 Ry, force <  1.0E-03 Ry/Bohr)

      End of BFGS Geometry Optimization

      Final energy   =      -6.6130860246 Ry

------------------------------------------------------------------------------------------------

The obtained values of volume (where V=a^2*c) and total energy are:

Vol (ang^3)      Total energy (eV)
28.640500       -89.84334
33.836204	-89.93096
36.654496	-89.95390
39.625092	-89.96478
42.752000       -89.96531
46.039228	-89.95738
49.490784	-89.94272
56.902912	-89.89999
65.020448	-89.84715

------------------------------------------------------------------------------------------------

The fit to the Murnaghan equation seems to be correct, which makes me  
question whether the calculation of the energies might be wrong.  
According to literature, the total energy should be around -17.822 eV  
after optimization. However, my values are around -89 eV (-6.61 Ry). I  
cannot find where the problem lies. I have verified that ecutwfc,  
ecutrho, and the k-points are converged. I can’t think of what else to  
do. Could someone please verify the total energy of this structure or  
guide me on where the issue might be?

Thank you in advance!