Skip to content

Commit

Permalink
Made package installable via pip (pymatgen.io.validation) and added p…
Browse files Browse the repository at this point in the history
…otcar validation using potcar_summary_stats from PMG, temporarily disabling mypy on pre-commit due to clash of module / package name in PMG
  • Loading branch information
Aaron Kaplan authored and Aaron Kaplan committed Oct 27, 2023
1 parent b1447cf commit 2b654d3
Show file tree
Hide file tree
Showing 9 changed files with 694 additions and 480 deletions.
15 changes: 5 additions & 10 deletions .pre-commit-config.yaml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
exclude: ^(docs|.*test_files|cmd_line|dev_scripts)

default_language_version:
python: python3.8
python: python3.11

ci:
autoupdate_schedule: monthly
Expand All @@ -20,11 +20,6 @@ repos:
- id: end-of-file-fixer
- id: trailing-whitespace

- repo: https://github.com/PyCQA/flake8
rev: 6.1.0
hooks:
- id: flake8

- repo: https://github.com/asottile/pyupgrade
rev: v3.10.1
hooks:
Expand All @@ -42,7 +37,7 @@ repos:
- --expand-star-imports
- --ignore-init-module-imports

- repo: https://github.com/pre-commit/mirrors-mypy
rev: v1.5.1
hooks:
- id: mypy
# - repo: https://github.com/pre-commit/mirrors-mypy
# rev: v1.5.1
# hooks:
# - id: mypy
4 changes: 4 additions & 0 deletions pymatgen/io/validation/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,3 +3,7 @@
simulations and calculations. That is, it checks calculation details against a reference
to ensure that data is compatible with some standard.
"""

from __future__ import annotations

from .validation import ValidationDoc
83 changes: 42 additions & 41 deletions pymatgen/io/validation/check_common_errors.py
Original file line number Diff line number Diff line change
@@ -1,16 +1,17 @@
import numpy as np


def _check_common_errors(
reasons,
warnings,
task_doc,
calcs_reversed,
parameters,
incar,
structure,
valid_max_allowed_scf_gradient,
ionic_steps,
num_ionic_steps_to_avg_drift_over
reasons,
warnings,
task_doc,
calcs_reversed,
parameters,
incar,
structure,
valid_max_allowed_scf_gradient,
ionic_steps,
num_ionic_steps_to_avg_drift_over,
):
# Check for cases where both GGA and METAGGA are set. This should *not* be allowed, as it can erroneously change
# the outputted energy significantly. See https://github.com/materialsproject/atomate2/issues/453#issuecomment-1699605867
Expand All @@ -22,7 +23,6 @@ def _check_common_errors(
"for more information."
)


# check if structure electronically converged
final_esteps = ionic_steps[-1]["electronic_steps"] if incar.get("ALGO", "Normal").lower() != "chi" else 0
# In a response function run there is no ionic steps, there is no scf step
Expand All @@ -31,45 +31,48 @@ def _check_common_errors(
to_check = {"e_wo_entrp", "e_fr_energy", "e_0_energy"}
while set(final_esteps[i]) == to_check:
i += 1
is_converged = (i + 1 != parameters.get("NELM", 60))
is_converged = i + 1 != parameters.get("NELM", 60)
elif len(final_esteps) < parameters.get("NELM", 60):
is_converged = True
else:
is_converged = False

if not is_converged:
reasons.append(f"CONVERGENCE --> Did not achieve electronic convergence in the final ionic step. NELM should be increased.")


reasons.append(
"CONVERGENCE --> Did not achieve electronic convergence in the final ionic step. NELM should be increased."
)

# Check if drift force is too large
try:
all_drift_forces = calcs_reversed[0]["output"]["outcar"]["drift"]
if len(all_drift_forces) < num_ionic_steps_to_avg_drift_over:
drift_forces_to_avg_over = all_drift_forces
else:
drift_forces_to_avg_over = all_drift_forces[::-1][: num_ionic_steps_to_avg_drift_over]
drift_forces_to_avg_over = all_drift_forces[::-1][:num_ionic_steps_to_avg_drift_over]

drift_mags_to_avg_over = [np.linalg.norm(drift_forces) for drift_forces in drift_forces_to_avg_over]
cur_avg_drift_mag = np.average(drift_mags_to_avg_over)

valid_max_drift = 0.05
if cur_avg_drift_mag > valid_max_drift:
reasons.append(f"CONVERGENCE --> Excessive drift of {round(cur_avg_drift_mag,4)} eV/A is greater than allowed "\
f"value of {valid_max_drift} eV/A.")
except Exception as e:
reasons.append(
f"CONVERGENCE --> Excessive drift of {round(cur_avg_drift_mag,4)} eV/A is greater than allowed "
f"value of {valid_max_drift} eV/A."
)
except Exception:
warnings.append("Drift forces not contained in calcs_reversed! Can not check for excessive drift.")



# Check for excessively positive final energies (which usually indicates a bad structure)
valid_max_energy_per_atom = 50
cur_final_energy_per_atom = task_doc.output.energy_per_atom
if cur_final_energy_per_atom > valid_max_energy_per_atom:
reasons.append(f"LARGE POSITIVE FINAL ENERGY --> Final energy is {round(cur_final_energy_per_atom,4)} eV/atom, which is "\
f"greater than the maximum allowed value of {valid_max_energy_per_atom} eV/atom.")

reasons.append(
f"LARGE POSITIVE FINAL ENERGY --> Final energy is {round(cur_final_energy_per_atom,4)} eV/atom, which is "
f"greater than the maximum allowed value of {valid_max_energy_per_atom} eV/atom."
)

# Check for excessively large final magnetic moments
### TODO: make this also work for elements Gd and Eu, which have magmoms >5 in at least one of their pure structures
# TODO: make this also work for elements Gd and Eu, which have magmoms >5 in at least one of their pure structures
valid_max_magmoms = {"Gd": 10, "Eu": 10}
cur_magmoms = [abs(mag["tot"]) for mag in calcs_reversed[0]["output"]["outcar"]["magnetization"]]
bad_site_magmom_msgs = []
Expand All @@ -87,36 +90,34 @@ def _check_common_errors(
f"at least one {cur_site_ele} site with magmom greater than {cur_site_max_allowed_magmom}"
)
if len(bad_site_magmom_msgs) > 0:
msg = f"MAGNETISM --> Final structure contains sites with magnetic moments that are very likely erroneous. This includes: "
msg = msg + "; ".join(set(bad_site_magmom_msgs)) + "."
msg = "MAGNETISM --> Final structure contains sites with magnetic moments that are very likely erroneous. This includes: "
msg = msg + "; ".join(set(bad_site_magmom_msgs)) + "."
reasons.append(msg)


# Check for a SCF gradient that is too large (usually indicates unstable calculations)
# NOTE: do NOT use `e_0_energy`, as there is a bug in the vasprun.xml when printing that variable
# (see https://www.vasp.at/forum/viewtopic.php?t=16942 for more details).
skip = abs(parameters.get("NELMDL", -5)) - 1
energies = [
d["e_fr_energy"]
for d in ionic_steps[-1]["electronic_steps"]
]
energies = [d["e_fr_energy"] for d in ionic_steps[-1]["electronic_steps"]]
if len(energies) > skip:
cur_max_gradient = np.max(np.gradient(energies)[skip:])
cur_max_gradient_per_atom = cur_max_gradient / structure.num_sites
if cur_max_gradient_per_atom > valid_max_allowed_scf_gradient:
warnings.append(f"STABILITY --> The max SCF gradient is {round(cur_max_gradient_per_atom,4)} eV/atom, which is larger than the typical max expected value of {valid_max_allowed_scf_gradient} eV/atom. "\
f"This sometimes indicates an unstable calculation.")
warnings.append(
f"STABILITY --> The max SCF gradient is {round(cur_max_gradient_per_atom,4)} eV/atom, which is larger than the typical max expected value of {valid_max_allowed_scf_gradient} eV/atom. "
f"This sometimes indicates an unstable calculation."
)
else:
warnings.append(
"Not enough electronic steps to compute valid gradient"
" and compare with max SCF gradient tolerance."
"Not enough electronic steps to compute valid gradient" " and compare with max SCF gradient tolerance."
)


# Check for Am and Po elements. These currently do not have proper elemental entries
# and will not get treated properly by the thermo builder.
chemsys = task_doc.chemsys
if ("Am" in chemsys) or ("Po" in chemsys):
reasons.append("COMPOSITION --> Your structure contains the elements Am and/or Po, which are not currently being accepted.")
reasons.append(
"COMPOSITION --> Your structure contains the elements Am and/or Po, which are not currently being accepted."
)

return reasons
return reasons
26 changes: 13 additions & 13 deletions pymatgen/io/validation/check_for_excess_empty_space.py
Original file line number Diff line number Diff line change
@@ -1,30 +1,30 @@
from pymatgen.analysis.local_env import VoronoiNN
import numpy as np


def check_for_excess_empty_space(structure):
##### Check 1: find large gaps along one of the defined lattice vectors
# Check 1: find large gaps along one of the defined lattice vectors
lattice_vec_lengths = structure.lattice.lengths
fcoords = np.array(structure.frac_coords)

lattice_vec_1_frac_coords = list(np.sort(fcoords[:,0]))
lattice_vec_1_frac_coords = list(np.sort(fcoords[:, 0]))
lattice_vec_1_frac_coords.append(lattice_vec_1_frac_coords[0] + 1)
lattice_vec_1_diffs = np.diff(lattice_vec_1_frac_coords)
max_lattice_vec_1_empty_dist = max(lattice_vec_1_diffs) * lattice_vec_lengths[0]

lattice_vec_2_frac_coords = list(np.sort(fcoords[:,1]))
lattice_vec_2_frac_coords = list(np.sort(fcoords[:, 1]))
lattice_vec_2_frac_coords.append(lattice_vec_2_frac_coords[0] + 1)
lattice_vec_2_diffs = np.diff(lattice_vec_2_frac_coords)
max_lattice_vec_2_empty_dist = max(lattice_vec_2_diffs) * lattice_vec_lengths[1]

lattice_vec_3_frac_coords = list(np.sort(fcoords[:,2]))
lattice_vec_3_frac_coords = list(np.sort(fcoords[:, 2]))
lattice_vec_3_frac_coords.append(lattice_vec_3_frac_coords[0] + 1)
lattice_vec_3_diffs = np.diff(lattice_vec_3_frac_coords)
max_lattice_vec_3_empty_dist = max(lattice_vec_3_diffs) * lattice_vec_lengths[2]

max_empty_distance = max(max_lattice_vec_1_empty_dist, max_lattice_vec_2_empty_dist, max_lattice_vec_3_empty_dist)


##### Check 2: get max voronoi polyhedra volume in structure
max_empty_distance = max(max_lattice_vec_1_empty_dist, max_lattice_vec_2_empty_dist, max_lattice_vec_3_empty_dist)

# Check 2: get max voronoi polyhedra volume in structure
def get_max_voronoi_polyhedra_volume(structure):
max_voronoi_polyhedra_vol = 0
vnn = VoronoiNN().get_all_voronoi_polyhedra(structure)
Expand All @@ -39,20 +39,20 @@ def get_max_voronoi_polyhedra_volume(structure):
try:
max_voronoi_polyhedra_vol = get_max_voronoi_polyhedra_volume(structure)
except Exception as e:

if "No Voronoi neighbors found for site - try increasing cutoff".lower() in str(e).lower():
try:
structure.make_supercell(2) # to circumvent weird issue with voronoi class, though this decreases performance significantly.
structure.make_supercell(
2
) # to circumvent weird issue with voronoi class, though this decreases performance significantly.
max_voronoi_polyhedra_vol = get_max_voronoi_polyhedra_volume(structure)
except Exception as e:
except Exception:
pass

if "infinite vertex in the Voronoi construction".lower() in str(e).lower():
print(f"{str(e)} As a result, this structure is marked as having excess empty space.")
max_voronoi_polyhedra_vol = np.inf


if (max_voronoi_polyhedra_vol > 25) or (max_voronoi_polyhedra_vol > 5 and max_empty_distance > 7.5):
return True
else:
return False
return False
Loading

0 comments on commit 2b654d3

Please sign in to comment.