PourbaixEntry from ComputedEntry - what energy corrections to use?


I am trying to calculate the electrochemical stability map for a material that is not on materials project. I used VASP to calculate the energy of the material, and from this I was able to make a ComputedEntry and then convert it into a PourbaixEntry easily by using:

entry_ce = ce.ComputedEntry(comp_dict, E)
entry_pb = pb.PourbaixEntry(entry_ce)

However, I know that the PourbaixEntry needs the energy E to be the “formation energ(y) with respect to hydrogen and oxygen gas” so I don’t think directly using my DFT energy will work for this situation. I know that when I create the ComputedEntry, there are methods for adding energy corrections. Could anyone provide some guidance on how I could “correct” my DFT energy from VASP such that it would be compatible with the PourbaixEntry? Thank you!

Hi @mdresser , I can provide some guidance on this. In order to be compatible with existing MP calculations, you should take your ComputedEntry and process it with MaterialsProjectAqueoustCompatibility (see code here). With default arguments, this will first apply energy adjustments to your solid DFT energy (see MP docs for details), then it will apply corrections to O2, H2, and other elements that are liquid or gas at room temperature to ensure consistency with the experimental formation energy of water.

If you create a list that contains your ComputedEntry plus the MP derived ComputedEntry for O2 and H2O, you can just do:

from pymatgen.entries.compatibility import MaterialsProjectAqueousCompatibility

entries = [<list of entries including H2O and O2 entries from the API>]
aqcompat = MaterialsProjectAqueousCompatibility()
entries = aqcompat.process_entries(entries)

After that you should be able to create your PourbaixEntry from the corrected ComputedEntry and proceed with your analysis. Note that depending on your material, the above may not result in any change to your DFT energy. You can always inspect the energy_adjustments attribute of any entry to see what corrections have been applied.

Hope that helps!

Hi @rkingsbury thank you for your reply! This is very helpful. I just wanted to clarify two things.

When making a ComputedEntry for O2 and H2O, you mention they should be from the API. It looks like mp-12597 is the most stable O2 structure. Thus, I perform a query and extract the final_energy associated with mp-12597. Is this the energy I should then use when I make the ComputedEntry? Or do you think I should instead extract something like the final_energy_per_atom, or formation_energy_per_atom, etc.? I am referencing this source for the different energies.

Using the final_energy, I get an error when running the last line of the code you recommended:

import pymatgen.entries.computed_entries as ce
from pymatgen.entries.compatibility import MaterialsProjectAqueousCompatibility

entry_ce = ce.ComputedEntry(comp_dict, E)
entry_O2 = ce.ComputedEntry('O8',-39.58364375)
entry_H2O = ce.ComputedEntry('H8O4',-59.56501858)

entries = [entry_ce, entry_O2, entry_H2O]

aqcompat = MaterialsProjectAqueousCompatibility()
entries = aqcompat.process_entries(entries)

The error is: KeyError: 'potcar_symbols'
Do you have any idea why this error is occurring? From looking at the source code and my full error message, the problem occurs at line 166 here. But I’m really unsure how I could fix the problem. Thank you!

Hi @mdresser , yes the whole infrastructure for adjusting energies is rather intricate and tricky to understand.

For any entries you will get from MP, I suggest you just directly download ComputedEntry via, e.g. MPRester.get_entries_in_chemsys() or one of the similar methods, that way you don’t have to build the entries yourself.

In general when you make a ComputedEntry you just need to make sure that the energy you provide corresponds to the composition. So for example if the O2 energy is -5 eV/atom, you could make a ComputedEntry with composition O2 and -10 eV, or composition O4 and -20 eV, etc. Just make sure they’re consistent. When you create a ComputedEntry, you want to use the DFT energy, not any corrected energy (corrections are added later).

The potcar_hash error is a bit subtle. When we process ComputedEntry using MaterialsProjectCompatibility2020 to adjust solid energies, that code checks the POTCARS that were used in the calculation to make sure they are consistent with the ones we used to parameterize our energy corrections. To do that, it looks for a ComputedEntry.data["potcar_symbols"] key, which is not populated when you make an entry by hand.

As long as you have run your calculation using MP parameters (i.e., by using MPRelaxSet to set up the calculation), the easiest fix would be to download some existing ComputedEntry from MP using the API that have the same composition or the same elements as your material, and copy the respective parts of the .data dict into your entry manually. See here for more details. Does that make sense?

Hi @rkingsbury thank you for your response. I think I understand the approach.

I am trying to relax the structure using MPRelaxSet as you mentioned, and I am running into an error. The code I am using is:

import pymatgen.io.vasp.sets as vsp
import pymatgen.core.structure as ps
import pymatgen.cli.pmg_potcar as pot
import os

cwd = os.getcwd()
poscar = cwd + '/POSCAR'
struc = ps.Structure.from_file(poscar)
inputs = vsp.MPRelaxSet(struc).write_input(cwd,potcar_spec=True)
pot.gen_potcar(cwd, 'POTCAR.spec')

I get the following error:
OSError: You do not have the right POTCAR with functional PBE and label Na_pv in your VASP_PSP_DIR

Do you know how I could fix this issue? I had this error before using the gen_potcar method and then started using that because I thought it would fix the issue. Thank you for any advice!

Hi @mdresser, when MPRelaxSet writes input files it checks the POTCAR files on your computer to make sure they match the ones that we use in MP calculations. If they don’t match, you get an error such as this one.

I would double check that you have downloaded the correct POTCARs (maybe re-download them) and make sure that you have set the VASP_PSP_DIR environment variable. See these lines.

Hi @rkingsbury, thank you for all of your help so far. I was able to download the correct POTCARs and perform the DFT relaxation with MPRelaxSet. Now I have modified my code, and I still get an error (although it is now a slightly different error). My code is:

import pymatgen.entries.computed_entries as ce
from pymatgen import MPRester
from pymatgen.entries.compatibility import MaterialsProjectAqueousCompatibility

entry_ce = ce.ComputedEntry(comp_dict, E)
entry_ce.parameters["potcar_symbols"] = ['PBE Na_pv','PBE Fe_pv', 'PBE Ni_pv','PBE O']
with MPRester(api_key) as m:
    entry_O2 = m.get_entries('mp-12957') 
    entry_H2O = m.get_entries('mp-697111')

# define list of entries
entries = [entry_ce, entry_O2, entry_H2O]

# Compute the 'AqueousCompatibility'
aqcompat = MaterialsProjectAqueousCompatibility()
entries = aqcompat.process_entries(entries)

The error occurs from the last line of code and says: AttributeError: 'list' object has no attribute 'parameters'

Do you have any recommendations for fixing this? Thank you very much!

Hmmm, make sure that entry_CE, entry_O2, and entry_H2O are actually ComputedEntry and not single-item lists. The error arises because something is calling ComputedEntry.parameters but the ComputedEntry it’s expecting is actually a list.

I suspect m.get_entries might always return a list, so I’d start there.

Thank you, this was indeed the issue!