4
$\begingroup$

I'm trying to follow simple example of determining the lattice constant (https://wiki.fysik.dtu.dk/ase/tutorials/eos/eos.html), but using QE instead of EMT. Here is ASE code I use:

import numpy as np
from ase.build import bulk
from ase.calculators.espresso import Espresso, EspressoProfile
from ase.io.trajectory import Trajectory
from ase.io.espresso import write_espresso_in

from ase.units import kJ

import os

pseudo_dir = '/mnt/d/PhD/DFT/SSSP_1.3.0_PBE_efficiency'
profile = EspressoProfile(
    command='mpirun -n 2 /home/mateja/qe-7.3.1/bin/pw.x',
    pseudo_dir=pseudo_dir
) 
input_data = {
    'prefix': 'Ag_lattice',
    'pseudo_dir': pseudo_dir,
    'outdir': '/mnt/d/PhD/DFT/LatticeConstant',
    'calculation': 'scf',
    'restart_mode': 'from_scratch',
    'tprnfor': True,
    'etot_conv_thr': 1e-5,
    'forc_conv_thr': 1e-4,
    'ecutwfc': 60,
    'ecutrho': 480,
    'input_dft': 'rpbe',
    'vdw_corr': 'dft-d3',
    'occupations': 'smearing',
    'degauss': 0.01,
    'smearing': 'gaussian',
    'conv_thr': 1e-8,
    'mixing_mode': 'local-TF',
    'mixing_beta': 0.35,
    'diagonalization': 'david',
    'ion_dynamics': 'bfgs',
    'bfgs_ndim': 6,
    'startingwfc': 'random',
}

pseudopotentials = {
    'Ag': 'Ag_ONCV_PBE-1.0.oncvpsp.upf'  # Ensure this file is in the pseudo_dir
}

def run_espresso(atoms, prefix, input_pars):

    calc = Espresso(profile=profile, pseudopotentials=pseudopotentials, tstress=True, tprnfor=True,input_data=input_pars, kpts=(1, 1, 1))
    atoms.calc = calc
    # Write the input file
    input_file = f'{prefix}.in'
    output_file = f'{prefix}.out'
    write_espresso_in(input_file,atoms, input_pars, pseudopotentials=pseudopotentials, format='espresso-in')

    return atoms.get_potential_energy()


a = 4.0  
ag = bulk('Ag', 'fcc', a=a)
traj = Trajectory('Ag.traj', 'w', atoms=ag)
for x in np.linspace(0.95, 1.05, 5):
    run_espresso(ag, 'Ag_lattice', input_data)
    traj.write(ag)
    
    
from ase.eos import EquationOfState
from ase.io import read
from ase.units import kJ

configs = read('Ag.traj@0:5')

volumes = [ag.get_volume() for ag in configs]
energies = [ag.get_potential_energy() for ag in configs]
eos = EquationOfState(volumes, energies)
v0, e0, B = eos.fit()
print(B / kJ * 1.0e24, 'GPa')
eos.plot('Ag-eos.png',show=False)

Why are my results so off?

enter image description here

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ EMT is an interatomic potential (see here). QE does plane-wave DFT, which is expected to be more accurate than an interatomic potential because of explicit treatment of the electrons in a quantum mechanical way. So, they should give you different-ish answers. You could check what the experimental lattice constant of fcc-Ag is and compare with what you get from EMT and QE (via the EOS analysis). $\endgroup$
    – CW Tan
    Commented Jul 3 at 13:10
  • 3
    $\begingroup$ That being said, there are various things you could change within DFT that will affect the outcome. E.g. why are you using a D3 correction ('vdw_corr': 'dft-d3') for bulk silver? Your single k-point kpts=(1, 1, 1) might not be sufficiently converged. $\endgroup$
    – CW Tan
    Commented Jul 3 at 13:13
  • 1
    $\begingroup$ I am not going to have the time to write a full answer right now, but D3 must be included for bulk silver if they will use D3 later for a surface calculation for example. You can't get a lattice without D3 then use it with D3 later. $\endgroup$
    – Tristan Maxson
    Commented Jul 4 at 7:08
  • 1
    $\begingroup$ Also, (1, 1, 1) kpt grid is definitely incorrect. I would expect at least a (12, 12, 12) for a primitive cell of Ag to be needed. $\endgroup$
    – Tristan Maxson
    Commented Jul 4 at 7:08
  • $\begingroup$ Thank you for your answer, bigger kpts grid definitely gives better results $\endgroup$
    – Mateja Jovanovic
    Commented Jul 4 at 9:48

0

Reset to default

You must log in to answer this question.