1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458
|
from datetime import datetime
from typing import Dict, List, Union, Optional
import numpy as np
from pydantic import ConfigDict, Field, ImportString
from pymatgen.core.structure import Structure
from pymatgen.io.vasp.inputs import Kpoints
from pymatgen.io.vasp.sets import VaspInputSet
from emmet.core.settings import EmmetSettings
from emmet.core.base import EmmetBaseModel
from emmet.core.mpid import MPID
from emmet.core.utils import DocEnum
from emmet.core.tasks import TaskDoc
from emmet.core.vasp.calc_types.enums import CalcType, TaskType
from emmet.core.vasp.task_valid import TaskDocument
SETTINGS = EmmetSettings()
class DeprecationMessage(DocEnum):
MANUAL = "M", "Manual deprecation"
SYMMETRY = (
"S001",
"Could not determine crystalline space group, needed for input set check.",
)
KPTS = "C001", "Too few KPoints"
KSPACING = "C002", "KSpacing not high enough"
ENCUT = "C002", "ENCUT too low"
FORCES = "C003", "Forces too large"
MAG = "C004", "At least one site magnetization is too large"
POTCAR = (
"C005",
"At least one POTCAR used does not agree with the pymatgen input set",
)
CONVERGENCE = "E001", "Calculation did not converge"
MAX_SCF = "E002", "Max SCF gradient too large"
LDAU = "I001", "LDAU Parameters don't match the inputset"
SET = ("I002", "Cannot validate due to missing or problematic input set")
UNKNOWN = "U001", "Cannot validate due to unknown calc type"
class ValidationDoc(EmmetBaseModel):
"""
Validation document for a VASP calculation
"""
task_id: MPID = Field(..., description="The task_id for this validation document")
valid: bool = Field(False, description="Whether this task is valid or not")
last_updated: datetime = Field(
description="Last updated date for this document",
default_factory=datetime.utcnow,
)
reasons: Optional[List[Union[DeprecationMessage, str]]] = Field(
None, description="List of deprecation tags detailing why this task isn't valid"
)
warnings: List[str] = Field(
[], description="List of potential warnings about this calculation"
)
data: Dict = Field(
description="Dictioary of data used to perform validation."
" Useful for post-mortem analysis"
)
model_config = ConfigDict(extra="allow")
nelements: Optional[int] = Field(None, description="Number of elements.")
symmetry_number: Optional[int] = Field(
None,
title="Space Group Number",
description="The spacegroup number for the lattice.",
)
@classmethod
def from_task_doc(
cls,
task_doc: Union[TaskDoc, TaskDocument],
kpts_tolerance: float = SETTINGS.VASP_KPTS_TOLERANCE,
kspacing_tolerance: float = SETTINGS.VASP_KSPACING_TOLERANCE,
input_sets: Dict[str, ImportString] = SETTINGS.VASP_DEFAULT_INPUT_SETS,
LDAU_fields: List[str] = SETTINGS.VASP_CHECKED_LDAU_FIELDS,
max_allowed_scf_gradient: float = SETTINGS.VASP_MAX_SCF_GRADIENT,
max_magmoms: Dict[str, float] = SETTINGS.VASP_MAX_MAGMOM,
potcar_stats: Optional[Dict[CalcType, Dict[str, str]]] = None,
) -> "ValidationDoc":
"""
Determines if a calculation is valid based on expected input parameters from a pymatgen inputset
Args:
task_doc: the task document to process
kpts_tolerance: the tolerance to allow kpts to lag behind the input set settings
kspacing_tolerance: the tolerance to allow kspacing to lag behind the input set settings
input_sets: a dictionary of task_types -> pymatgen input set for validation
pseudo_dir: directory of pseudopotential directory to ensure correct hashes
LDAU_fields: LDAU fields to check for consistency
max_allowed_scf_gradient: maximum uphill gradient allowed for SCF steps after the
initial equillibriation period
potcar_stats: Dictionary of potcar stat data. Mapping is calculation type -> potcar symbol -> hash value.
"""
nelements = task_doc.nelements or None
symmetry_number = task_doc.symmetry.number if task_doc.symmetry else None
bandgap = task_doc.output.bandgap
calc_type = task_doc.calc_type
task_type = task_doc.task_type
run_type = task_doc.run_type
inputs = task_doc.orig_inputs
chemsys = task_doc.chemsys
calcs_reversed = [
calc if not hasattr(calc, "model_dump") else calc.model_dump()
for calc in task_doc.calcs_reversed
]
if calcs_reversed[0].get("input", {}).get("structure", None):
structure = calcs_reversed[0]["input"]["structure"]
else:
structure = task_doc.input.structure or task_doc.output.structure
if isinstance(structure, dict):
structure = Structure.from_dict(structure)
reasons = []
data = {} # type: ignore
warnings: List[str] = []
if str(calc_type) in input_sets:
try:
valid_input_set = _get_input_set(
run_type, task_type, calc_type, structure, input_sets, bandgap
)
except (TypeError, KeyError, ValueError):
reasons.append(DeprecationMessage.SET)
valid_input_set = None
try:
# Sometimes spglib can't determine space group with the default
# `symprec` and `angle_tolerance`. In these cases,
# `Structure.get_space_group_info()` fails
valid_input_set.structure.get_space_group_info()
except Exception:
reasons.append(DeprecationMessage.SYMMETRY)
valid_input_set = None
if valid_input_set:
# Checking POTCAR summary_stats if a directory is supplied
if potcar_stats:
if _potcar_stats_check(task_doc, potcar_stats):
if task_type in [
TaskType.NSCF_Line,
TaskType.NSCF_Uniform,
TaskType.DFPT_Dielectric,
TaskType.Dielectric,
]:
warnings.append(DeprecationMessage.POTCAR.__doc__) # type: ignore
else:
reasons.append(DeprecationMessage.POTCAR)
# Checking K-Points
# Calculations that use KSPACING will not have a .kpoints attr
if task_type != TaskType.NSCF_Line:
# Not validating k-point data for line-mode calculations as constructing
# the k-path is too costly for the builder and the uniform input set is used.
if valid_input_set.kpoints is not None:
if _kpoint_check(
valid_input_set,
inputs,
calcs_reversed,
data,
kpts_tolerance,
):
reasons.append(DeprecationMessage.KPTS)
else:
# warnings
_kspacing_warnings(
valid_input_set, inputs, data, warnings, kspacing_tolerance
)
# warn, but don't invalidate if wrong ISMEAR
valid_ismear = valid_input_set.incar.get("ISMEAR", 1)
incar = inputs.get("incar", {})
curr_ismear = incar.get("ISMEAR", 1)
if curr_ismear != valid_ismear:
warnings.append(
f"Inappropriate smearing settings. Set to {curr_ismear},"
f" but should be {valid_ismear}"
)
# Checking ENCUT
encut = incar.get("ENCUT")
valid_encut = valid_input_set.incar["ENCUT"]
data["encut_ratio"] = float(encut) / valid_encut # type: ignore
if data["encut_ratio"] < 1:
reasons.append(DeprecationMessage.ENCUT)
# U-value checks
if _u_value_checks(task_doc, valid_input_set, warnings):
reasons.append(DeprecationMessage.LDAU)
# Check the max upwards SCF step
if _scf_upward_check(
calcs_reversed, inputs, data, max_allowed_scf_gradient, warnings
):
reasons.append(DeprecationMessage.MAX_SCF)
# Check for Am and Po elements. These currently do not have proper elemental entries
# and will not get treated properly by the thermo builder.
if ("Am" in chemsys) or ("Po" in chemsys):
reasons.append(DeprecationMessage.MANUAL)
# Check for magmom anomalies for specific elements
if _magmom_check(calcs_reversed, structure, max_magmoms=max_magmoms):
reasons.append(DeprecationMessage.MAG)
else:
if "Unrecognized" in str(calc_type):
reasons.append(DeprecationMessage.UNKNOWN)
else:
reasons.append(DeprecationMessage.SET)
doc = ValidationDoc(
task_id=task_doc.task_id,
calc_type=calc_type,
run_type=task_doc.run_type,
valid=len(reasons) == 0,
reasons=reasons,
data=data,
warnings=warnings,
nelements=nelements,
symmetry_number=symmetry_number,
)
return doc
def _get_input_set(run_type, task_type, calc_type, structure, input_sets, bandgap):
# Ensure inputsets get proper additional input values
if "SCAN" in run_type.value:
valid_input_set: VaspInputSet = input_sets[str(calc_type)](structure, bandgap=bandgap) # type: ignore
elif task_type == TaskType.NSCF_Uniform or task_type == TaskType.NSCF_Line:
# Constructing the k-path for line-mode calculations is too costly, so
# the uniform input set is used instead and k-points are not checked.
valid_input_set = input_sets[str(calc_type)](structure, mode="uniform")
elif task_type == TaskType.NMR_Electric_Field_Gradient:
valid_input_set = input_sets[str(calc_type)](structure, mode="efg")
else:
valid_input_set = input_sets[str(calc_type)](structure)
return valid_input_set
def _scf_upward_check(calcs_reversed, inputs, data, max_allowed_scf_gradient, warnings):
skip = abs(inputs.get("incar", {}).get("NLEMDL", -5)) - 1
energies = [
d["e_fr_energy"]
for d in calcs_reversed[0]["output"]["ionic_steps"][-1]["electronic_steps"]
]
if len(energies) > skip:
max_gradient = np.max(np.gradient(energies)[skip:])
data["max_gradient"] = max_gradient
if max_gradient > max_allowed_scf_gradient:
return True
else:
warnings.append(
"Not enough electronic steps to compute valid gradient"
" and compare with max SCF gradient tolerance"
)
return False
def _u_value_checks(task_doc, valid_input_set, warnings):
# NOTE: Reverting to old method of just using input.hubbards which is wrong in many instances
input_hubbards = {} if task_doc.input.hubbards is None else task_doc.input.hubbards
if valid_input_set.incar.get("LDAU", False) or len(input_hubbards) > 0:
# Assemble required input_set LDAU params into dictionary
input_set_hubbards = dict(
zip(
valid_input_set.poscar.site_symbols,
valid_input_set.incar.get("LDAUU", []),
)
)
all_elements = list(set(input_set_hubbards.keys()) | set(input_hubbards.keys()))
diff_ldau_params = {
el: (input_set_hubbards.get(el, 0), input_hubbards.get(el, 0))
for el in all_elements
if not np.allclose(input_set_hubbards.get(el, 0), input_hubbards.get(el, 0))
}
if len(diff_ldau_params) > 0:
warnings.extend(
[
f"U-value for {el} should be {good} but was {bad}"
for el, (good, bad) in diff_ldau_params.items()
]
)
return True
return False
def _kpoint_check(input_set, inputs, calcs_reversed, data, kpts_tolerance):
"""
Checks to make sure the total number of kpoints is correct
"""
valid_num_kpts = input_set.kpoints.num_kpts or np.prod(input_set.kpoints.kpts[0])
if calcs_reversed:
input_dict = calcs_reversed[0].get("input", {})
if not input_dict:
input_dict = inputs
else:
input_dict = inputs
kpoints = input_dict.get("kpoints", {})
if isinstance(kpoints, Kpoints):
kpoints = kpoints.as_dict()
elif kpoints is None:
kpoints = {}
num_kpts = kpoints.get("nkpoints", 0) or np.prod(kpoints.get("kpoints", [1, 1, 1]))
data["kpts_ratio"] = num_kpts / valid_num_kpts
return data["kpts_ratio"] < kpts_tolerance
def _kspacing_warnings(input_set, inputs, data, warnings, kspacing_tolerance):
"""
Issues warnings based on KSPACING values
"""
valid_kspacing = input_set.incar.get("KSPACING", 0)
if kspacing := inputs.get("incar", {}).get("KSPACING"):
data["kspacing_delta"] = kspacing - valid_kspacing
# larger KSPACING means fewer k-points
if data["kspacing_delta"] > kspacing_tolerance:
warnings.append(
f"KSPACING is greater than input set: {data['kspacing_delta']}"
f" lower than {kspacing_tolerance} "
)
elif data["kspacing_delta"] < kspacing_tolerance:
warnings.append(
f"KSPACING is lower than input set: {data['kspacing_delta']}"
f" lower than {kspacing_tolerance} "
)
def _potcar_stats_check(task_doc, potcar_stats: dict):
"""
Checks to make sure the POTCAR summary stats is equal to the correct
value from the pymatgen input set.
"""
data_tol = 1.0e-6
try:
potcar_details = task_doc.calcs_reversed[0].model_dump()["input"]["potcar_spec"]
except KeyError:
# Assume it is an old calculation without potcar_spec data and treat it as passing POTCAR hash check
return False
use_legacy_hash_check = False
if any(entry.get("summary_stats", None) is None for entry in potcar_details):
# potcar_spec doesn't include summary_stats kwarg needed to check potcars
# fall back to header hash checking
use_legacy_hash_check = True
all_match = True
for entry in potcar_details:
if not entry["titel"]:
all_match = False
break
symbol = entry["titel"].split(" ")[1]
ref_summ_stats = potcar_stats[str(task_doc.calc_type)].get(symbol, None)
if not ref_summ_stats:
# Symbol differs from reference set - deprecate
all_match = False
break
if use_legacy_hash_check:
all_match = any(
all(
entry[key] == ref_stat[key]
for key in (
"hash",
"titel",
)
)
for ref_stat in ref_summ_stats
)
else:
all_match = False
for ref_stat in ref_summ_stats:
key_match = all(
set(ref_stat["keywords"][key])
== set(entry["summary_stats"]["keywords"][key])
for key in ["header", "data"]
)
data_match = False
if key_match:
data_match = all(
abs(
ref_stat["stats"][key][stat]
- entry["summary_stats"]["stats"][key][stat]
)
< data_tol
for stat in ["MEAN", "ABSMEAN", "VAR", "MIN", "MAX"]
for key in ["header", "data"]
)
all_match = key_match and data_match
if all_match:
# Found at least one match to reference POTCAR summary stats,
# that suffices for the check
break
if not all_match:
break
return not all_match
def _magmom_check(
calcs_reversed: list, structure: Structure, max_magmoms: dict[str, float]
):
"""
Checks for maximum magnetization values for specific elements.
Returns True if the maximum absolute value outlined below is exceded for the associated element.
"""
if (outcar := calcs_reversed[0]["output"]["outcar"]) and (
mag_info := outcar.get("magnetization", [])
):
return any(
abs(mag_info[isite].get("tot", 0.0))
> abs(max_magmoms.get(site.label, np.inf))
for isite, site in enumerate(structure)
)
return False
def _get_unsorted_symbol_set(structure: Structure):
"""
Have to build structure_symbol set manually to ensure
we get the right order since pymatgen sorts its symbol_set list.
"""
return list(
{
str(sp): 1 for site in structure for sp, v in site.species.items() if v != 0
}.keys()
)
|