5
$\begingroup$

Using the below code I am trying to calculate the formation energy per atom and bandgap given the cif file. the values that are present in the material project. but I got these values which are wrong for formation energy per atom. can anyone point out where I am doing wrong.or the correct code to calculate it

Formation Energy per Atom: -138.67884553

Bandgap : -4.9528

import os
import pandas as pd
from ase.io import read
from ase.calculators.vasp import Vasp
from ase.optimize import BFGS
import xml.etree.ElementTree as ET
import time

def run_vasp_calculation(file_path):
    start_time = time.time()
    base_name = os.path.basename(file_path).replace('.cif', '')
    dir_path = f"/{base_name}"
    if not os.path.exists(dir_path):
        os.makedirs(dir_path)

    structure = read(file_path)
    calculator = Vasp(xc='PBE',
                      encut=400,
                      kpts=(3, 3, 3),
                      ibrion=2,
                      nsw=50,
                      ismear=0,
                      sigma=0.05,
                      lreal='Auto',
                      ispin=2,
                      ediff=1e-5,
                      algo='Fast',
                      command="mpirun -np 1 vasp_gpu",
                      directory=dir_path)
    structure.set_calculator(calculator)

    try:
        opt = BFGS(structure)
        opt.run(fmax=0.01)
        energy = structure.get_total_energy()
        num_atoms = structure.get_number_of_atoms()
        calculation_time = time.time() - start_time
    except Exception as e:
        print(f"Error during optimization for {file_path}: {e}")
        return None, None, None

    return energy, num_atoms, calculation_time

def extract_bandgap(directory):
    try:
        tree = ET.parse(os.path.join(directory, 'vasprun.xml'))
        root = tree.getroot()
        bandgap = root.find(".//bandgap").text
        return float(bandgap)
    except Exception as e:
        print(f"Error reading bandgap from {directory}: {e}")
        return None

def process_cifs(cif_directory):
    cif_files = [os.path.join(cif_directory, f) for f in os.listdir(cif_directory) if f.endswith('.cif')]
    data = []
    
    for cif_file in cif_files:
        print(f"Processing {cif_file}")
        start_time = time.time()
        compound_energy, num_atoms, calc_time = run_vasp_calculation(cif_file)
        if compound_energy is not None:
            formation_energy_per_atom = compound_energy / num_atoms
            bandgap = extract_bandgap(f"/data/{os.path.basename(cif_file).replace('.cif', '')}")
            data.append([cif_file, compound_energy, formation_energy_per_atom, bandgap, calc_time])
            print(f'Formation Energy per Atom for {cif_file}: {formation_energy_per_atom} eV/atom, Calculation Time: {calc_time} seconds')
        else:
            print("Failed to calculate compound energy.")
        overall_time = time.time() - start_time
        print(f'Total Time for processing {cif_file}: {overall_time} seconds')

    df = pd.DataFrame(data, columns=['CIF', 'Energy', 'Formation Energy per Atom', 'Bandgap', 'Calculation Time'])
    df.to_csv('summary.csv', index=False)

# Example usage
process_cifs('cifs')

cif I am using

# generated using pymatgen
data_Nd3Al11
_symmetry_space_group_name_H-M   Immm
_cell_length_a   4.36307790
_cell_length_b   9.97432892
_cell_length_c   12.95781714
_cell_angle_alpha   90.00000000
_cell_angle_beta   90.00000000
_cell_angle_gamma   90.00000000
_symmetry_Int_Tables_number   71
_chemical_formula_structural   Nd3Al11
_chemical_formula_sum   'Nd6 Al22'
_cell_volume   563.90831666
_cell_formula_units_Z   2
loop_
 _symmetry_equiv_pos_site_id
 _symmetry_equiv_pos_as_xyz
  1  'x, y, z'
  2  '-x, -y, -z'
  3  '-x, -y, z'
  4  'x, y, -z'
  5  'x, -y, -z'
  6  '-x, y, z'
  7  '-x, y, -z'
  8  'x, -y, z'
  9  'x+1/2, y+1/2, z+1/2'
  10  '-x+1/2, -y+1/2, -z+1/2'
  11  '-x+1/2, -y+1/2, z+1/2'
  12  'x+1/2, y+1/2, -z+1/2'
  13  'x+1/2, -y+1/2, -z+1/2'
  14  '-x+1/2, y+1/2, z+1/2'
  15  '-x+1/2, y+1/2, -z+1/2'
  16  'x+1/2, -y+1/2, z+1/2'
loop_
 _atom_site_type_symbol
 _atom_site_label
 _atom_site_symmetry_multiplicity
 _atom_site_fract_x
 _atom_site_fract_y
 _atom_site_fract_z
 _atom_site_occupancy
  Nd  Nd0  4  0.00000000  0.00000000  0.31720041  1
  Nd  Nd1  2  0.00000000  0.00000000  0.00000000  1
  Al  Al2  8  0.00000000  0.27452085  0.13593955  1
  Al  Al3  8  0.00000000  0.36810718  0.33399153  1
  Al  Al4  4  0.00000000  0.21533047  0.50000000  1
  Al  Al5  2  0.00000000  0.50000000  0.00000000  1

Modified code

def load_elemental_energies(save_file='elemental_energies.csv'):
    return pd.read_csv(save_file).set_index('Element')['Energy per Atom'].to_dict()

def calculate_formation_energy(structure, total_energy, elemental_energies):
    element_counts = {elem: int(count) for elem, count in re.findall(r'([A-Z][a-z]*)(\d*)', structure.get_chemical_formula()) if count.isdigit()}
    total_atoms = sum(element_counts.values())
    formation_energy = total_energy
    print(f"Debug: Total Energy of the Compound: {total_energy}, Total Atoms: {total_atoms}")
    for elem, count in element_counts.items():
        if elem in elemental_energies:  # Ensure element key exists in reference data
            element_energy = count * elemental_energies[elem]
            formation_energy -= element_energy
            print(f"Debug: Element: {elem}, Count: {count}, Elemental Energy: {element_energy}")
    formation_energy_per_atom = formation_energy #/ total_atoms
    print(f"Debug: Formation Energy Per Atom: {formation_energy_per_atom}")
    return formation_energy_per_atom

elemental energy calculation

def calculate_elemental_energies(elemental_cif_directory, save_file='elemental_energies.csv'):
    if os.path.exists(save_file):
        return pd.read_csv(save_file).set_index('Element')['Energy per Atom'].to_dict()

    elemental_energies = {}
    elemental_files = [os.path.join(elemental_cif_directory, f) for f in os.listdir(elemental_cif_directory) if f.endswith('.cif')]
    
    for elemental_file in elemental_files:
        element_name = os.path.basename(elemental_file).split('.')[0]
        energy, num_atoms, _ = run_vasp_calculation(elemental_file, elemental_cif_directory)
        if energy is not None and num_atoms != 0:
            elemental_energies[element_name] = energy / num_atoms
    
    # Save the elemental energies to a CSV file
    pd.DataFrame(list(elemental_energies.items()), columns=['Element', 'Energy per Atom']).to_csv(save_file, index=False)
    return elemental_energies
````
$\endgroup$
2
  • $\begingroup$ What are the correct values that you are looking for? How large is the deviation? Adding these info to your post, and also a link to the MaterialsProject CIF file that you are working on would be helpful. $\endgroup$ Commented May 5 at 10:17
  • 1
    $\begingroup$ @AbdulMuhaymin-FreePalestine the correct value for formation energy is near to this value -0.379 and the PBE bandgap should be near to this value -0.009 eV. and I have updated the post with the cif I am using. $\endgroup$
    – harsh
    Commented May 5 at 10:34

1 Answer 1

5
$\begingroup$

To put it simply, your code has some issues. It doesn't calculate the formation energy per atom. Basically, what you did is:

# This is NOT formation energy
energy = structure.get_total_energy()
num_atoms = structure.get_number_of_atoms()
formation_energy_per_atom = energy / num_atoms

This is not the definition of formation energy per atom. You are just dividing the total energy by number of atoms. In order to calculate the formation energy, you need the energies of isolated atoms too.

So, first calculate the total energy of pure metallic $\mathrm{Al}$. Let's say the structure has $m$ atoms of $\mathrm{Al}$. Then find the total energy of aluminium per atom, $E_{\textrm{Al, per atom}}=E_{\textrm{Al, total}}/m$. Similarly, find the total energy of pure neodymium per atom, $E_{\textrm{Nd, per atom}}=E_{\textrm{Nd, total}}/n$ where the $\mathrm{Nd}$ structure has $n$ atoms. You can find the CIF of $\mathrm{Al}$ and $\mathrm{Nd}$ here and here, respectively.

Now you have the total energies of the elements. You now need the total energy of the compound structure (you have already calculated it in the above code). You need to know how many $\mathrm{Al}$ atoms and how many $\mathrm{Nd}$ atoms you have in your crystal. Let's say that the number is $a$ and $b$ for $\mathrm{Al}$ and $\mathrm{Nd}$, respectively.

Then you can do the calculation using the following formula: $$ E_{\textrm{formation per atom}} = \frac{E_{\textrm{total}} - aE_{\textrm{Al, per atom}} - bE_{\textrm{Nd, per atom}}}{a+b} $$ This should fix your formation energy. I couldn't figure out the issue with band gap though. Maybe check the vasprun.xml file directly or using tools like VASP-DOS_extractor. It seems something else is not correct since the band gap value should be positive for a material with gap. Negative band gap value means that the compound is metal and in that case you should not worry about the band gap.

Lastly, in Materials Project, the formation energies are calculated using the r2SCAN functional (metaGGA) which is very accurate but computationally expensive. However, here you used PBE functional which is at the level of GGA. So you should not expect the exact same value as given in Materials Project but they should be close.

$\endgroup$
6
  • $\begingroup$ rather than calculating the energy of the isolated atoms by doing dft which is computationally expensive as I am doing for 100's of crystals. can I find these values some where for the entire elements in the periodic table so that I can use them as a lookup and divide those values with the compound energy to calculate the formation energy per atom $\endgroup$
    – harsh
    Commented May 6 at 11:02
  • $\begingroup$ @harsh not really. The absolute value of the energy of system, be it an isolated atom or a crystal, depends on the code, functional, pseudopotential, etc. So there is no universal values for these. You must get the elemental energies using the same pseudopotential, same code, etc. Calculating isolated atom energies are expensive but notice that the method I described is not for isolated atom. They are for metallic crystal, and then you can divide the energy by the number of atoms to get the atomic energy. These calculations are very cheap computationally, only a few atoms and very small cell. $\endgroup$ Commented May 6 at 15:51
  • $\begingroup$ hey I have added the modified code in the post can you take a look and tell me where it correct or wrong. still my values are off by a very big margin. $\endgroup$
    – harsh
    Commented May 7 at 6:50
  • $\begingroup$ dont you think the formula must be something like $𝐸_{formation per atom}=(𝐸_{total}−𝑎𝐸_{Al, per atom}−𝑏𝐸_{Nd, per atom}) / a+b$ $\endgroup$
    – harsh
    Commented May 7 at 7:04
  • 1
    $\begingroup$ @harsh your question might get rolled back. Generally we don't allow questions to be significantly modified after getting answered. $\endgroup$ Commented May 7 at 11:37

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .