#!/usr/bin/env python3 """ Calculates equilibrium crystal structure and energy at zero temperature and pressure. Supports arbitrary crystal prototypes specified according to the AFLOW prototype label standard for portable models and standard simulator models. Computes `binding-energy-crystal `_ and `crystal-structure-npt `_ OpenKIM properties. The result is obtained by a symmetry-constrained minimization. If :data:`num_param_sets`>1, multiple relaxations from different initial guesses will be performed. If the relaxed structures are unique, multiple pairs of these properties will be written. Inputs read from ``stdin`` =========================== These are usually fed by a test's ``runner`` from a file rendered from a test generator (:ref:`doc.test_gen_format`) using ``test_template/pipeline.stdin.tpl.genie``). .. data:: species Space-separated list of chemical elements :type: list[str] .. data:: prototype_label AFLOW prototype label :type: str .. data:: parameter_names Space-separated list of parameter names :type: list[str] .. data:: modeltype Model type (only 'standard' supported at this time) :type: str .. data:: num_param_sets Number of parameter sets :type: int .. data:: param_values_superlist Newline-separated list of space-separated lists of parameter values :type: list[list[str]] .. data:: modelname KIM model name :type: str Overview of operation: ---------------------- 1) Inputs are read and processed 2) For each parameter set: a) Use :func:`util.aflow_util.AFLOW.build_atoms_from_prototype` to build the :class:`ase.atoms.Atoms` object b) Use :class:`ase.spacegroup.symmetrize.FixSymmetry` to constrain the atoms to the crystal symmetry c) Perform symmetry-constrained minimization using :class:`ase.optimize.BFGS` d) Save the structure to a file for later 3) Use :func:`util.aflow_util.AFLOW.compare_materials_dir` to identify groups of duplicates among relaxed structures 4) For each group of duplicates: a) Loop over the inidividual files, try to add a property instance until one is successfully added: i) Use :func:`util.aflow_util.AFLOW.get_prototype` to check that the prototype label didn't change ii) Get the library prototype and shortname using :func:`util.aflow_util.AFLOW.get_library_prototype_label_and_shortname` iii) Try to add the properties using :func:`util.property_util.add_property_inst` iv) Validate the properties using :func:`util.physics_validator.validate_binding_energy_crystal` and :func:`util.physics_validator.validate_crystal_structure_npt` 5) If at least one property instance was added, write the ``output/results.edn`` file Date: 3/17/2023 Author: Ilia Nikiforov , Ellad Tadmor """ from multiprocessing.pool import RUN from kim_property import kim_property_dump from kim_python_utils.ase.core import get_isolated_energy_per_atom import numpy as np import pathlib import os, shutil from ase.spacegroup.symmetrize import FixSymmetry from ase.constraints import UnitCellFilter from ase.optimize import BFGS from ase.calculators.kim import KIM from util.aflow_util import AFLOW, read_shortnames from util.property_util import add_property_inst from util.physics_validator import validate_binding_energy_crystal, validate_crystal_structure_npt # Set path to driver test_driver_dir = pathlib.Path(__file__).parent MAXSTEP = 0.05 FTOL = 1e-8 ITER = 10000 COMPARE_DIR = "compare" RELAX_POSCAR_PATH = "compare/{0}" ############################################################################### if __name__ == "__main__": # grab from stdin (or a file) species = input("element(s):\n").split() print(species) prototype_label = input( "AFLOW prototype label:\n" ) print(prototype_label) parameter_names = input("Parameter names:\n").split() print(parameter_names) modeltype = input("model type (only 'standard' supported at this time):\n") print(modeltype) num_param_sets = int(input("number of parameter sets:\n")) print(num_param_sets) param_values_superlist = [] for i in range(num_param_sets): param_values = input("Parameter values for parameter set %d:\n"%i).split() print(param_values) param_values_superlist.append(param_values) modelname = input("model name:\n") print(modelname) # Verify input if modeltype not in ["standard"]: raise RuntimeError( "Unknown model type, received modeltype = {}".format(modeltype) ) # Get dictionary of shortnames for prototype labels shortnames = read_shortnames() spacegroup = int(prototype_label.split("_")[2]) if os.path.isdir(COMPARE_DIR): shutil.rmtree(COMPARE_DIR) os.mkdir(COMPARE_DIR) aflow = AFLOW() # Main minimization loop over parameter sets for i in range(num_param_sets): atoms = aflow.build_atoms_from_prototype(species,prototype_label,param_values_superlist[i]) # Apply symmetry constraints symmetry = FixSymmetry(atoms) atoms.set_constraint(symmetry) atoms_wrapped = UnitCellFilter(atoms) # Attach the KIM calculator for the model being run atoms.calc = KIM(modelname) # Optimize opt = BFGS(atoms_wrapped,maxstep=MAXSTEP) try: converged = opt.run(fmax=FTOL,steps=ITER) iteration_limits_reached = not converged minimization_stalled = False except: minimization_stalled = True iteration_limits_reached = False forces = atoms.get_forces() stress = atoms.get_stress() # Compute the average energy per atom after subtracting out the energies of the # isolated atoms energy_isolated = sum( [get_isolated_energy_per_atom(modelname,sym) for sym in atoms.get_chemical_symbols()] ) energy_per_atom = ( atoms.get_potential_energy() - energy_isolated ) / atoms.get_global_number_of_atoms() print("Minimization "+("converged" if not minimization_stalled else "stalled")+" after "+(("hitting the maximum of "+str(ITER)) if iteration_limits_reached else str(opt.nsteps))+" steps.") print("Maximum force component: "+str(np.max(np.abs(forces)))+" eV/Angstrom") print("Maximum stress component: "+str(np.max(np.abs(stress)))+" eV/Angstrom^3") print("==== Minimized structure obtained from ASE ====") print("symbols = ", atoms.get_chemical_symbols()) print("basis = ", atoms.get_scaled_positions()) print("cellpar = ", atoms.get_cell()) print("forces = ", forces) print("stress = ", stress) print("energy per atom = ", energy_per_atom) print("===============================================") atoms.write(RELAX_POSCAR_PATH.format(i),format='vasp') # compare relaxed structures for duplicates comparison=aflow.compare_materials_dir(COMPARE_DIR) # Try to write at most one property per group of duplicates property_inst = None for materials_group in comparison: # get all indices belonging to this group indices = [int(materials_group['structure_representative']['name'].split('/')[-1])] for structure in materials_group['structures_duplicate']: indices.append(int(structure['name'].split('/')[-1])) print ("Parameter sets %s relaxed to duplicate structures, attempting to write only one of them."%str(indices)) # loop over all materials in this group until we find one that successfully adds a property instance for i in indices: relax_proto_des=aflow.get_prototype(RELAX_POSCAR_PATH.format(i)) relax_proto_label=relax_proto_des['aflow_prototype_label'] if relax_proto_label != prototype_label: print("Prototype label changed during relaxation: test template prototype is %s, while relaxed is %s. Skipping parameter set %d." %(prototype_label,relax_proto_label,i)) if i==indices[-1]: print("No parameter sets in this group successfully added a property instance. Skipping this group.") continue (libproto,shortname)=aflow.get_library_prototype_label_and_shortname(RELAX_POSCAR_PATH.format(i),shortnames) # Try to add property instances. Back up original so we can keep going even if one parameter set fails property_inst_new = property_inst # these are just serialized edn strings, so no need for deepcopy or anything try: property_inst_new = add_property_inst(energy_per_atom,species,relax_proto_des,libproto,shortname,property_inst_new) validate_crystal_structure_npt(property_inst_new) validate_binding_energy_crystal(property_inst_new) property_inst = property_inst_new print("Successfully added property instance for parameter set %d"%i) break except Exception as e: print("Skipping parameter set %d because of error while adding or validating property:\n %s"%(i,e)) if i==indices[-1]: print("No parameter sets in this group successfully added a property instance. Skipping this group.") # Dump the results in a file if not (property_inst is None or \ property_inst in ('None', '', '[]')): with open("output/results.edn", "w") as fp: kim_property_dump(property_inst, fp)