Inconsistency in VASP Calculations for the Same Element from API

[Minju]

Hello everyone,

I am writing this post because I discovered something strange while performing calculations using VASP.

Let me summarize in two lines:
The data structure received through API in primitive and conventional formats is strange.
When processing these two with BoltzTraP2, a post-processing tool, there are serious value errors. A solution is needed.

Here’s a brief workflow:

1. Load structure data from MP_API using Materials Project
2. Convert the data using get_primitive_structure() function to generate a Primitive structure
3. Separately, using to_conventional() function on the API structure data generates a Conventional structure
4. At this point, the Final energy (= Total energy; has no physical meaning, but should have values within error range for the same material) from calculations of primitive and conventional structures is confirmed to be the same
5. However, I found that the Fermi energy is different between the two structures
6. As a result, the product of resistivity and mean free path obtained during interpolation using BoltzTraP2 comes out very differently. (Since the same material should have the same properties, the values should be the same)

Let me explain in more detail.
First, it seems there is a critical issue with get_primitive.
Previously, when other users and I created topics on similar subjects, the structure obtained through ‘get_primitive’ is quite different from the structure available on the Materials Project website.
The data from the website is actually more similar to the values obtained through ‘to_conventional’.
The structure obtained through ‘get_primitive’ even damages the symmetry of the cell structure.
(e.g., in HCP structure, α, β, γ are not 90, 90, 120, but rather 90, 90, 119.98, etc.)
As a result, calculation errors due to symmetry frequently occurred in both Quantum ESPRESSO and VASP calculations. This was also a common problem in BoltzTraP2, which is used for post-processing.

Anyway, having different Fermi energies in the same structure is indeed problematic. This is because such small differences in Fermi level can cause very large errors in the final results.

get_primitive of Co1F3_mp-559435

• static
  
• NSCF
  

INCAR (Static)

ALGO = Fast
EDIFF = 0.0002
ENCUT = 520
ISMEAR = -5
ISPIN = 2
LASPH = True
LCHARG = True
LDAU = True
LDAUJ = 0 0
LDAUL = 2 0
LDAUPRINT = 1
LDAUTYPE = 2
LDAUU = 3.32 0
LMAXMIX = 4
LORBIT = 11
LREAL = False
LWAVE = False
MAGMOM = 13.189 30.081
NELM = 300
NSW = 0
PREC = Accurate
SIGMA = 0.05
AMIX = 0.2
BMIX = 0.0001
KPAR = 4
NCORE = 12

INCAR(Non-self consistency field)

ALGO = Fast
EDIFF = 0.0002
ENCUT = 520
ICHARG = 11
ISMEAR = -5
ISPIN = 2
KSPACING = 0.05
LASPH = True
LCHARG = False
LDAU = True
LDAUJ = 0 0
LDAUL = 2 0
LDAUPRINT = 1
LDAUTYPE = 2
LDAUU = 3.32 0
LMAXMIX = 4
LORBIT = 11
LREAL = False
LWAVE = False
NELM = 100
NSW = 0
PREC = Accurate
SIGMA = 0.05
KPAR = 4
NCORE = 12

• POSCAR

Co1 F3
1.0
2.6472384033333332 1.5283830066666670 2.1613651899999997
0.0000023733333334 -3.0567679833333337 2.1613661900000003
2.6472397666666665 -1.5283839666666665 -2.1613641800000001
Co F
1 3
direct
0.9999965233333332 0.9999965233333331 0.0000034766666666 Co3+
0.5000006766666666 0.0001207133333332 0.0001185800000000 F-
0.0001207133333334 0.9998814199999999 0.4999993233333332 F-
0.9998810866666666 0.5000006766666668 0.9998789533333333 F-

to_conventional of Co1F3_mp-559435

• INCAR: static

ALGO = Fast
EDIFF = 0.0002
ENCUT = 520
ISMEAR = -5
ISPIN = 2
LASPH = True
LCHARG = True
LDAU = True
LDAUJ = 0 0
LDAUL = 2 0
LDAUPRINT = 1
LDAUTYPE = 2
LDAUU = 3.32 0
LMAXMIX = 4
LORBIT = 11
LREAL = False
LWAVE = False
MAGMOM = 4*0.6
NELM = 300
NSW = 0
PREC = Accurate
SIGMA = 0.05
AMIX = 0.2
BMIX = 0.0001
KPAR = 4
NCORE = 4

• INCAR: NSCF

ALGO = Fast
EDIFF = 0.0002
ENCUT = 520
ICHARG = 11
ISMEAR = -5
ISPIN = 2
KSPACING = 0.05
LASPH = True
LCHARG = False
LDAU = True
LDAUJ = 0 0
LDAUL = 2 0
LDAUPRINT = 1
LDAUTYPE = 2
LDAUU = 3.32 0
LMAXMIX = 4
LORBIT = 11
LREAL = False
LWAVE = False
NELM = 100
NSW = 0
PREC = Accurate
SIGMA = 0.05
KPAR = 4
NCORE = 4

• POSCAR
  Co1 F3

1.0
-2.6472393272282879 -1.5283843381846158 2.1613651866667190
-2.6472393272282879 1.5283843381846154 -2.1613651866667194
0.0000000000000000 -3.0567686763692312 -2.1613651866667190
Co F
1 3
direct
0.9999999999999999 0.0000000000000001 0.0000000000000001 Co3+
0.0001198116666665 0.5000000000000001 0.0001198116666667 F-
0.9998801883333335 0.9998801883333335 0.4999999999999999 F-
0.5000000000000000 0.0001198116666666 0.9998801883333335 F-

[tsmathis]

I think this is another instance of differences between the tolerances used in MP for symmetry determination and the defaults in pymatgen.

Linking a similar post from a few days ago:

wrong space group in mp-1271793

The default tolerances used for symmetry determination are slightly different between the default in, e.g., pymatgen (0.01), and the software we use to publish data in the Materials Project (0.1): from mp_api.client import MPRester with MPRester() as mpr: summary_doc = mpr.materials.summary.search(material_ids=["mp-1271793"])[0] print(summary_doc.symmetry.symbol) >>> 'R-3m' print(summary_doc.get_space_group_info()) >>> ('C2/m', 12) print(summary_doc.get_space_group_info(symprec=0.1)) >>> (…

@Minju, does this get you going in the right direction?

[Minju]

Hello.
Thank you for your interest in my problem.

I forgot to include information about the Fermi energy value of the conventional structure,
but while I obtained a Fermi energy value of -0.2370 eV from the structure converted directly to primitive via API, the structure converted from API > conventional to primitive has a value of -0.5512 eV.

According to the similar post you mentioned earlier, symmetry issues can occur due to tolerance errors, but I’m not sure whether it’s reasonable for structural differences caused by such errors to change the Fermi energy to such a large extent.

Due to this degree of Fermi energy difference, my final calculation result (resistivity * mean free path value obtained using BoltzTraP2) ultimately shows approximately a 100-fold difference.

[Aaron_Kaplan]

The Fermi level is measurable only with respect to a reference state (usually the vacuum energy), and thus can change with the pseudopotentials you use, the magnetic configuration, etc… The Fermi level is not a physical observable without similar referencing (e.g., this paper)

You probably want to look at energy differences instead, like the bandgap, occupied bandwidth (Fermi level minus the bottom of the valence band), or effective mass

Another thing to note: your two static calculations use different initial magnetic configurations (MAGMOM) - are you sure that both calculations are in the same final magnetic state?

[Minju]

Hello Aaron.

I understand that the Fermi level has a relative value, and I know that it does not have a physical value like Total or Final energy.
As you mentioned, if the pseudopotentials are different, such differences could occur.
However, the POTCAR used in my calculations is the same.

And here, I’ll attach the SCF and NSCF output of both bulk to primitive and bulk to conventional to primitive structure.

• Bulk (API get structure) → Primitive
  SCF
  

  
  
  NSCF
  
  

• Bulk (API get structure) → Conventional → Primitive
  SCF
  

  
  
  NSCF
  
  

As you can see, both structures have the same (nearly identical) final magnetization.

Hmm… I’ll also attach the Free energy data as well.
As I mentioned in our previous discussion that the Free energy is also the same, let me show you that data once more for confirmation.

	B → P	B → C → P
SCF	-17.33810578	-17.78246006
NSCF	-17.34184624	-17.78248674

Based on my knowledge, I don’t understand why two structures with the same total energy and magnetization would have different Fermi energies.
If there is something I should know or if I missed something, I would greatly appreciate it if you could let me know.

Additionally, I also calculated the band gap that you mentioned.
i) B → P, SCF & NSCF / Fermi energy: -0.236951 eV

Adopt the value of the 12th band, which is ultimately occupied among the band energies recorded at the last K-point
-0.236951 - (-1.9024) = 1.665449
-0.214707 - (-1.8836) = 1.668893

ii) B → C → P, SCF & NSCF / Fermi energy: -0.897697 eV

Adopt the value of the 13th band, which is ultimately occupied, from the band energies recorded at the last K-point
-0.897697 - (-1.3030) = 0.405303
-0.551184 - (-1.2489) = 0.697716
At the Fermi level, there is a significant difference between the two structures, and when converting from Conventional to Primitive, there is also a difference between SCF and NSCF.

Thank you for taking the time out of your busy schedule.

[Aaron_Kaplan]

I’m not sure about the way you’re reading off the gap, why not just look at the DOS (DOSCAR or vasprun.xml)?

Two things you could check:

1. For an insulator, the Fermi level could be placed anywhere from the VBM to just below the CBM. There are different conventions for where to place the Fermi level, see the conventions available in VASP.
   It’s possible that VASP is placing the Fermi level somewhere between the VBM and CBM with no other significant difference to the electronic structure.
   (By contrast, the chemical potential at 0K should be at the midgap, i.e., the negative of half the sum of ionization potential and electron affinity)

2. Do the initial structures in both cases match?

from pymatgen.analysis.structure_matcher import StructureMatcher
from pymatgen.core import Structure

matcher = StructureMatcher()
s_bp = Structure.from_file("POSCAR_bp") # the API to primitive
s_bcp = Structure.from_file("POSCAR_bcp") # the API to conventional to primitive

for s in (s_bp, s_bcp):
  for k in s.site_properties:
    s.remove_site_property(k)
matcher.fit(s_bp, s_bcp)

If that prints True, then the two structures are very likely symmetrically equivalent. I tried this just now with mp-559435 and got True

[Minju]

/home/0/.local/lib/python3.11/site-packages/pymatgen/core/structure.py:3087: EncodingWarning: We strongly encourage explicit encoding, and we would use UTF-8 by default as per PEP 686
with zopen(filename, mode=“rt”, errors=“replace”) as file:
/home/0/.local/lib/python3.11/site-packages/pymatgen/core/structure.py:3087: EncodingWarning: We strongly encourage explicit encoding, and we would use UTF-8 by default as per PEP 686
with zopen(filename, mode=“rt”, errors=“replace”) as file:

기존 당신이 저에게 주신 코드는 encoding 에러가 발생했습니다.
이에 저는 아래와 같이 코드를 수정하였습니다.

import warnings

warnings.filterwarnings("ignore", message=r".*explicit `encoding`.*")

from pymatgen.analysis.structure_matcher import StructureMatcher
from pymatgen.core import Structure

matcher = StructureMatcher()
s_bp = Structure.from_file("POSCAR_bp")
s_bcp = Structure.from_file("POSCAR_bcp")

for s in (s_bp, s_bcp):
    for k in s.site_properties:
        s.remove_site_property(k)

print(matcher.fit(s_bp, s_bcp))

위 코드를 실행해서 제가 얻은 값은 True 입니다.

POSCAR_bp

Co1 F3
1.0
2.6472384033333332 1.5283830066666670 2.1613651899999997
0.0000023733333334 -3.0567679833333337 2.1613661900000003
2.6472397666666665 -1.5283839666666665 -2.1613641800000001
Co F
1 3
direct
0.9999965233333332 0.9999965233333331 0.0000034766666666 Co3+
0.5000006766666666 0.0001207133333332 0.0001185800000000 F-
0.0001207133333334 0.9998814199999999 0.4999993233333332 F-
0.9998810866666666 0.5000006766666668 0.9998789533333333 F-

POSCAR_bcp

Co1 F3
1.0
-2.6472393272282879 -1.5283843381846158 2.1613651866667190
-2.6472393272282879 1.5283843381846154 -2.1613651866667194
0.0000000000000000 -3.0567686763692312 -2.1613651866667190
Co F
1 3
direct
0.9999999999999999 0.0000000000000001 0.0000000000000001 Co3+
0.0001198116666665 0.5000000000000001 0.0001198116666667 F-
0.9998801883333335 0.9998801883333335 0.4999999999999999 F-
0.5000000000000000 0.0001198116666666 0.9998801883333335 F-

더 나아가 band gap을 아래의 코드를 활용해 계산하였습니다.
I referenced your other answer.

from pymatgen.io.vasp import Vasprun

vasprun = Vasprun("bcp_vasprun.xml") # or the path to your vasprun.xml
gap = vasprun.get_band_structure(efermi="smart").get_band_gap()
print("band gap from vasprun:", gap)

For B > C > P structure
band gap from vasprun: {‘direct’: False, ‘transition’: ‘(0.500,0.500,0.500)-(0.000,0.000,0.000)’, ‘energy’: 0.20450000000000002}

For B > P structure
band gap from vasprun: {‘energy’: 0.0, ‘direct’: False, ‘transition’: None}

Now, regarding the question you asked,

1. Since the material I’m examining is a metal, using the default value LEGACY might be inaccurate. (Though I don’t fully understand, is it correct that the default sets the Fermi level below VBM and CBM like insulators?)

So I modified it to EFERMI = MIDGAP and ran the codes again, but
the error keeps occurring as shown in the image.

2. The initial structures are identical in both cases.

3. The Band Gap is as shown above.
   According to those results, does this mean that Co1F3, which should be a metal, is behaving like a semiconductor (with a band gap of ~0.2045 eV) in the b>c>p structure?

How can it be divided into conductor/semiconductor when the structures are identical and the Total energy (though it has no physical meaning) is the same?

I tried to keep it as short as possible, but it ended up being long again.
I apologize.

[Aaron_Kaplan]

You can ignore the encoding error, it’s not a meaningful warning (also please translate text beforehand)

EFERMI = LEGACY should place the fermi level at the lowest energy which occupies N electrons, so the VBM in a gapped material

The only other possibility I see is that the calculation is underconverged (EDIFF is too large) and thus sensitive to the parallelization settings. Can you rerun the calculation with identical INCARs (you will have to change the MAGMOM to reflect difference in sites in conventional and MP cells)?

[Minju]

Oh, I’d written those texts in English, but It seems my device automatically translated it again, which caused the problem.
My apologize

Now, I use some codes for making INCAR and POSCAR etc., this code makes little bit difference between bcp (bulk > conventional > primitive) and bp (bulk > primitive) structure.

When I make bp structure, it recalls MAGMOM from its site product (I understood this is after relaxation magnetization).
However, in bcp structureit recalls its default value from pymatgen, MPStatic.

In conclusion, it was determined that MAGMOM causes the largest difference.
It’s quite strange, but despite the total charge and total magnetization being similar, the difference in MAGMOM causes a difference in Fermi energy.
(+ I will attach the photo below, but with that option, the EFERMI = MIDGAP option doesn’t work. Um… I’m not sure what the problem is, but seeing that IERR=5 was output, it seems to recognize that the INCAR is written incorrectly. For now, I’ve set this problem aside.)

Could there be a reason that would cause such changes?

• Different MAGMOM with Co1F3_bcp_MAG_Default INCAR
  
  

Co1F3_bp: Bulk > Primitive; MAGMOM from Site property
Co1F3_bcp_Mag_Siteproduct: Bulk > Conventional > Primitive; MAGMOM from Site property
Co1F3_bcp_Mag_Default: Bulk > Conventional > Primitive; MAGMOM from MPStaticSet

[Aaron_Kaplan]

No worries! EFERMI is an INCAR option only in newer versions of VASP, maybe just double check that you’re using a supported version?

Could there be a reason that would cause such changes?

The electronic structure depends on the magnetic moments, but beyond that, it could be a VASP-related quirk. You might get a clearer answer from posting on VASP’s forum and referencing this discussion.

[Minju]

That’s indeed the case.
When I first thought about it, I believed this cause was due to symmetry differences resulting from structural differences between bcp and bp, but now it seems like it’s a VASP issue.
Actually, there are many cases where post-processing programs like BoltzTraP2 don’t work due to sym errors because of such structural differences, so I inquired about this, but fortunately it’s not that kind of problem.

Thank you very much for your help.
Best regards,
GH Maeng