How to do a geometrical optimization of hexagonal Iron

We will consider the ferromagnetic state

1. Introduction

Iron in the hcp structure has a hexagonal lattice with two parameters, a and c:

a1 = (a, 0, 0),       a2 = (−a/2 ,a√3/2 , 0)    a3 = (0, 0, c)

and two atoms in the unit cell, at positions      d1 = (0, 0, 0),      d2 = (0, a√3/2, c/2)

The non cubic lattice can be described in many equivalent ways:

• ibrav=4,       celldm(1)= a,    celldm(3)= c/a, as in sample file fe_hcp.scf.in

• ibrav=4,   a= a,   c= c,    both in Å, not a.u.

• ibrav=0, celldm(1)= a, lattice vectors in card
CELL PARAMETERS (alat)
1.000000 0.000000 0.000000
-0.500000 0.866025 0.000000
0.000000 0.000000 1.633000

• Using space group=194, and Wyckoff positions 2c (see sample file fe hcp.wy.in)
....

Useful hint: use xcrysden --pwi fe_hcp.scf.in to verify the structure.

fe_hcp.scf.in

&control
    prefix='fe'
 /
 &system
    ibrav= 4, celldm(1)=5.0, celldm(3)=1.633, nat=2, ntyp= 1,
    ecutwfc = 25.0, ecutrho = 200.0,
    occupations='smearing', smearing='mv', degauss=0.02,
    nbnd=12,
    nspin=2, 
    starting_magnetization(1)=0.4
 /
 &electrons
    mixing_beta = 0.5 
 /
ATOMIC_SPECIES
 Fe  1.  Fe.pbe-nd-rrkjus.UPF
ATOMIC_POSITIONS alat
 Fe 0.0 0.00       0.0 
 Fe 0.0 0.5773503  0.8165
K_POINTS automatic
8 8 6 0 0 0

fe_hcp.wy.in

&control
    prefix='fe'
 /
 &system
    ! Specify lattice as space group number + Wyckoff positions !
    ! On input, nat = number of Wyckoff positions, not of atoms !
    space_group=194, celldm(1)=5.0, celldm(3)=1.633, nat=1, ntyp= 1,
    ecutwfc = 25.0, ecutrho = 200.0,
    occupations='smearing', smearing='mv', degauss=0.02,
    nbnd=12,
    nspin=2, 
    starting_magnetization(1)=0.4
 /
 &electrons
    mixing_beta = 0.5 
 /
ATOMIC_SPECIES
 Fe  1.  Fe.pbe-nd-rrkjus.UPF
# crystal_sg: use Wyckoff positions with needed parameters (none in this case)
ATOMIC_POSITIONS crystal_sg
 Fe 2c
K_POINTS automatic
8 8 6 0 0 0

2. Equation of state for non cubic structures

For a few selected values of the unit cell volume,

1. Perform a variable-cell calculation (sample file fe_hcp.vcr.in)

• set calculation=’vc-relax’
• add two namelists, &ions ... / and &cell ... /, after &electrons
• in &cell ... /, set cell dofree=’shape’ (keeps volume fixed)

2. Inspect output of fe_hcp.vcr.out

• several scf steps are performed, forces (zero by symmetry) and stresses computed
• the energy decrease and the stress becomes more and more isotropic as the minimum is approached
• a final scf step is performed with plane waves computed for the final cell
• the final cell is printed after the last CELL PARAMETERS card

3. Calculation

We need 2 files

fe_hcp.vcr.in

 &control
    calculation='vc-relax'
    outdir='./OUT_HCP'
    pseudo_dir = './',
    prefix='fe'
 /
 &system
    ibrav= 4, celldm(1)=5.0, celldm(3)=1.633, nat=2, ntyp= 1,
    ecutwfc = 25.0, ecutrho = 200.0,
    occupations='smearing', smearing='mv', degauss=0.02,
    nbnd=12,
    nspin=2, 
    starting_magnetization(1)=0.4
 /
 &electrons
    mixing_beta = 0.5 
 /
 &ions
 /
 &cell
    cell_dofree='shape'
 /
ATOMIC_SPECIES
 Fe  1.  Fe.pbe-nd-rrkjus.UPF
ATOMIC_POSITIONS alat
 Fe 0.0 0.00       0.0 
 Fe 0.0 0.5773503  0.8165
K_POINTS automatic
8 8 6 0 0 0

Fe.pbe-nd-rrkjus.UPF

You can download the file using the following command

$ wget https://www.quantum-espresso.org/upf_files/Fe.pbe-nd-rrkjus.UPF --no-check-certificate

algerien1970@linux-wipm:~/abinitio/QE-tutorials/Fe-hex> ls
fe_hcp.vcr.in  Fe.pbe-nd-rrkjus.UPF

algerien1970@linux-wipm:~/abinitio/QE-tutorials/Fe-hex> pw.x < fe_hcp.vcr.in |tee fe_hcp.vcr.out

.
.
End of BFGS Geometry Optimization

     Final enthalpy           =    -111.5771652771 Ry

     File ./OUT_HCP/fe.bfgs deleted, as requested
Begin final coordinates
     new unit-cell volume =    176.77744 a.u.^3 (    26.19571 Ang^3 )
     density =      0.12678 g/cm^3

CELL_PARAMETERS (alat=  5.00000000)
   0.980620254   0.000000000  -0.000000000
  -0.490310127   0.849242051  -0.000000000
  -0.000000000  -0.000000000   1.698182913

ATOMIC_POSITIONS (alat)
Fe            0.0000000000       -0.0000000000        0.0000000000
Fe           -0.0000000000        0.5661613978        0.8490914563
End final coordinates