Phonon Dispersion of SrTiO3, EuTiO3 using Quantum Espresso (DFT-U)

Question

the last few days I've tried to calculate the phonon-dispersion of SrTiO3 (STO) and EuTiO3 (ETO) using the PHonon package of QE. SCF, DOS, PDOS and Bandsctructure calculation worked well on these materials, I've used a Hubbard-correction for both. My scf.in file for STO looks as follows:

STO:

&CONTROL
    calculation = 'scf'
    outdir = './out/'
    prefix      = 'STO_gamma'
    pseudo_dir  = './nc-sr-04_pbe_stringent_upf/'
/

&SYSTEM
  ! a = 3.94513 # old value before vc-relax.x
  A = 3.93795
  degauss =   0.005
  ecutwfc =   80
  lda_plus_u=.true.
  lda_plus_u_kind = 0
  hubbard_u(1) = 0
  hubbard_u(2) = 5.6
  hubbard_u(3) = 4
  ibrav = 0
  nbnd = 33
  nat = 5
  nosym = .false.
  nspin = 2
  ntyp = 3
  occupations = 'smearing'
  starting_magnetization(1) =   1.0d-01
  starting_magnetization(2) =   1.0d-01
  starting_magnetization(3) =   4.16d-01
/

&ELECTRONS
    conv_thr =  1.d-10
    mixing_beta = 0.7
/

ATOMIC_SPECIES
Sr    87.62000  Sr.upf
Ti   47.86700  Ti.upf
O   15.99940  O.upf

CELL_PARAMETERS {alat}
   0.998180509   0.000000000   0.000000000
   0.000000000   0.998180509   0.000000000
   0.000000000   0.000000000   0.998180509

ATOMIC_POSITIONS {crystal}
Sr           0.0000000000       0.0000000000       0.0000000000
Ti           0.5000000000       0.5000000000       0.5000000000
O            0.5000000000       0.0000000000       0.5000000000
O            0.5000000000       0.5000000000       0.0000000000
O            0.0000000000       0.5000000000       0.5000000000

K_POINTS {automatic}
16 16 16

The pesudopotentials are the norm conserving potentials from pseudo-dojo.org. The bandstructure shows an energy gap of about 3.2 eV with these settings.

The ph.x input looks like:

phonon calculation at Gamma point.
&inputph
  outdir = './out/'
  prefix      = 'STO_gamma'
  tr2_ph = 1.0d-3
  nmix_ph = 10
  epsil = .false.
  amass(1) = 87.62000
  amass(2) = 47.86700
  amass(3)  = 15.99940
  fildyn = 'STO_phonon.dyn'
/
0 0 0

The output generated by dynmat.x after applying the acustic sum rule = crystal looks like:

# mode   [cm-1]    [THz]      IR
    1  -6235.76 -186.9435    0.0000
    2  -6235.76 -186.9435    0.0000
    3  -6235.76 -186.9435    0.0000
    4  -4352.85 -130.4953    0.0000
    5  -4352.85 -130.4953    0.0000
    6  -4352.85 -130.4953    0.0000
    7  -4105.78 -123.0883    0.0000
    8  -4105.78 -123.0883    0.0000
    9  -4105.78 -123.0883    0.0000
   10  -2138.73  -64.1174    0.0000
   11  -2138.73  -64.1174    0.0000
   12  -2138.73  -64.1174    0.0000
   13      0.00    0.0000    0.0000
   14      0.00    0.0000    0.0000
   15      0.00    0.0000    0.0000

This calculation alone took over 3 hours only for the gamma point. The input file for my scf calculation of ETO looks like

ETO:

&CONTROL
  calculation = 'scf'
  outdir = './out/'
  prefix = 'EuTiO3'
  pseudo_dir = './pseudo/'
/

&SYSTEM
  a = 3.9611945300
  degauss = 0.01
  ecutrho =   600
  ecutwfc =   75
  lda_plus_u=.true.
  lda_plus_u_kind = 0
  hubbard_u(1) = 9
  hubbard_u(2) = 0
  hubbard_u(3) = 0
  ibrav = 1
  nat = 5
  nosym = .false.
  nspin = 2
  ntyp = 3
  occupations = 'smearing'
  smearing = 'mv'
  starting_magnetization(1) =   4.1176470588d-01
  starting_magnetization(2) =   1.0000000000d-01
  starting_magnetization(3) =   4.1666666667d-01
/

&ELECTRONS
  conv_thr =   1.0000000000d-06
  electron_maxstep = 1000
  mixing_beta =   4.0000000000d-01
/

ATOMIC_SPECIES
Eu     151.964 Eu.GGA-PBE-paw-v1.0.UPF
O      15.9994 O.pbe-n-kjpaw_psl.0.1.UPF
Ti     47.867 ti_pbe_v1.4.uspp.F.UPF

ATOMIC_POSITIONS {crystal}
Eu           0.0000000000       0.0000000000       0.0000000000 
Ti           0.5000000000       0.5000000000       0.5000000000 
O            0.5000000000       0.5000000000       0.0000000000 
O            0.5000000000       0.0000000000       0.5000000000 
O            0.0000000000       0.5000000000       0.5000000000 

K_POINTS {automatic}
6 6 6 0 0 0

Bandstructure with band gap of about 1.05 eV: The ph.x input is:

!phonons of EuTiO3 at Gamma
&inputph
  prefix = 'EuTiO3'
  outdir = './out/'
  tr2_ph = 1.0d-14
  amass(1) = 151.964
  amass(2) = 15.9994
  amass(3) =  47.867
  ! epsil = .true.
  ! lraman .true.
  fildyn='EuTiO3.dynG'
 /
0.000000000000000   0.000000000000000   0.000000000000000

Which got me the following phonon frequencies (again asr = crystal):

# mode   [cm-1]    [THz]      IR
    1   -636.90  -19.0939    0.0000
    2   -636.90  -19.0939    0.0000
    3   -636.90  -19.0939    0.0000
    4   -547.18  -16.4041    0.0000
    5   -547.18  -16.4041    0.0000
    6   -547.18  -16.4041    0.0000
    7   -401.41  -12.0340    0.0000
    8   -401.41  -12.0340    0.0000
    9   -401.41  -12.0340    0.0000
   10     -0.00   -0.0000    0.0000
   11     -0.00   -0.0000    0.0000
   12      0.00    0.0000    0.0000
   13     64.61    1.9371    0.0000
   14     64.61    1.9371    0.0000
   15     64.61    1.9371    0.0000

I have a question about these results:

Why is it not possible to set epsil = .true., lraman .true. ?

I hope you are still around and thank you for reading :)

Have a nice day,

Simon

Answer 1

Follow up: Using a fixed occupation it was possible to get the dielectric constant since then QE won't treat the system as a metal anymore. Since STO is not stable at 0K, ph.x only showed imaginary frequencies (see https://arxiv.org/pdf/1506.02382.pdf). As for ETO it might be similar but I am not sure yet.