#!/usr/bin/env python3 ################################################################################ # # CDDL HEADER START # # The contents of this file are subject to the terms of the Common Development # and Distribution License Version 1.0 (the "License"). # # You can obtain a copy of the license at # http:# www.opensource.org/licenses/CDDL-1.0. See the License for the # specific language governing permissions and limitations under the License. # # When distributing Covered Code, include this CDDL HEADER in each file and # include the License file in a prominent location with the name LICENSE.CDDL. # If applicable, add the following below this CDDL HEADER, with the fields # enclosed by brackets "[]" replaced with your own identifying information: # # Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved. # # CDDL HEADER END # # Copyright (c) 2019,2021, Regents of the University of Minnesota. # All rights reserved. # # Contributor(s): # Daniel S. Karls # ################################################################################ # The docstring below is vc_description """This verification check ensures that a model can properly convert units of length and energy. Working in one unit system, a random cluster of atoms is generated based on the smallest cutoff of the model such that the energy and/or forces (depending on what the model actually reports) are sizeable. Next, the model is reinitialized with a different set of registered units. The coordinates of the atoms in the cluster are rescaled to correspond to the new length units and the energy and/or forces are requested from the model once again. The energy and/or force returned by the model are then compared to a manual conversion of those originally returned (in the first set of units).""" import random import textwrap import functools import numpy as np import kimpy import kim_python_utils.vc as kim_vc_utils __version__ = "001" __author__ = "Daniel S. Karls" class KimpyError(Exception): """ A call to a kimpy function resulted in a RuntimeError being raised """ pass def check_kimpy_call(vc, func, *args): """ Call a kimpy function using its arguments and, if a RuntimeError is raised, catch it and raise a KimpyError with the exception's message. A VC object is optionally passed in so that, if an exception occurs, it can be echoed to the VC report as well. (Starting with kimpy 2.0.0, a RuntimeError is the only exception type raised when something goes wrong.) Parameters ---------- vc : VerificationCheck or None An instance of the VerificationCheck class defined in kim_python_utils.vc. If None, then any error messages are printed to stdout alone without being written to the report.txt file the verification check creates. f : function The kimpy function that is to be called. Raises ------ KimpyError If the kimpy function raises a RuntimeError when called. """ try: return func(*args) except RuntimeError as e: formatted_message = ( f'Calling kimpy function "{func.__name__}" failed:\n {str(e)}' ) if vc is None: print(formatted_message) else: vc.rwrite(formatted_message) raise KimpyError(formatted_message) # Formatting for various text output DASHWIDTH = 120 indent = " " * 4 indented_wrap = textwrap.TextWrapper( initial_indent=indent, subsequent_indent=indent, width=100 ) unindented_wrap = textwrap.TextWrapper(width=100) # Conversion constants used by the KIM API ANGSTROM_TO_BOHR = np.double(1.889726124558930) HARTREE_TO_EV = np.double(27.21138433) ################################################################################ # # FUNCTIONS # ################################################################################ def wrapped(s, indent=False): """ Wrap a paragraph of text contained in a single string """ if indent: return ("\n").join(indented_wrap.wrap(s)) return ("\n").join(unindented_wrap.wrap(s)) def get_kim_model_supported_species_and_codes(vc, kim_model): """ Get all the supported species of a KIM model and the corresponding codes used by the model. The codes are required Parameters ---------- vc : VerificationCheck An instance of the VerificationCheck class defined in kim_python_utils.vc kim_model : str The KIM ID of the model Returns ------- species_map : dict Keys are chemical symbol strings (e.g. "Ar") and the corresponding values are the corresponding integer codes that the model uses internally to represent these """ def get_kim_model_supported_species(vc, kim_model): """ Get all the supported species by a KIM model Parameters ---------- vc : VerificationCheck An instance of the VerificationCheck class defined in kim_python_utils.vc kim_model : str The KIM ID of the model Returns ------- species: list of str a list of chemical symbols (e.g. ["Mo", "S"]) """ species = [] num_kim_species = kimpy.species_name.get_number_of_species_names() for i in range(num_kim_species): species_name = check_kimpy_call(vc, kimpy.species_name.get_species_name, i) species_support, _ = check_kimpy_call( vc, kim_model.get_species_support_and_code, species_name ) if species_support: species.append(str(species_name)) return species supported_species = get_kim_model_supported_species(vc, kim_model) species_map = dict() for s in supported_species: _, code = check_kimpy_call( vc, kim_model.get_species_support_and_code, kimpy.species_name.SpeciesName(s), ) species_map[s] = code return species_map ################################################################################ def check_compute_arguments(vc, compute_arguments): """ Determine whether a KIM Model computes total energy and/or atomic forces Parameters ---------- vc : VerificationCheck An instance of the VerificationCheck class defined in kim_python_utils.vc compute_arguments : ComputeArguments An instance of the ComputeArguments class defined in kimpy.compute_arguments Returns ------- model_computes_energy : bool Indicates whether the model computes total energy (True) or not (False) model_computes_forces : bool Indicates whether the model computes atomic forces (True) or not (False) """ model_computes_energy = check_kimpy_call( vc, compute_arguments.get_argument_support_status, kimpy.compute_argument_name.partialEnergy, ) model_computes_forces = check_kimpy_call( vc, compute_arguments.get_argument_support_status, kimpy.compute_argument_name.partialForces, ) if model_computes_energy in [ kimpy.support_status.required, kimpy.support_status.optional, ]: model_computes_energy = True else: model_computes_energy = False if model_computes_forces in [ kimpy.support_status.required, kimpy.support_status.optional, ]: model_computes_forces = True else: model_computes_forces = False if not model_computes_energy and not model_computes_forces: report_error( vc, "Model does not compute energy or forces. This verification check " "requires that it compute at least one of these.", ) return model_computes_energy, model_computes_forces ################################################################################ def register_compute_arguments_pointers( vc, compute_arguments, num_particles, species_codes, particle_contributing, coords, energy, forces, ): """ Register pointers to local python arrays in the KIM API using kimpy Parameters ---------- vc : VerificationCheck An instance of the VerificationCheck class defined in kim_python_utils.vc compute_arguments : ComputeArguments An instance of the ComputeArguments class defined in kimpy.compute_arguments num_particles : numpy.ndarray Array containing an np.intc that indicates the total number of particles present species_codes : numpy.ndarray Array containing np.intc that indicate, for each atom present, the corresponding species code used by the model particle_contributing : numpy.ndarray Array containing np.intc values that indicate, for each atom present, whether it is contributing (1) or non-contributing (0) coords : numpy.ndarray Two-dimensional array containing np.double values indicating the coordinates of each atom energy : numpy.ndarray Array containing an np.double that indicates the partial energy of the atomic system forces : numpy.ndarray Two-dimensional array containing np.double values indicating the total internal forces acting on each atom Returns ------- None """ set_argument_pointer = functools.partial( check_kimpy_call, vc, compute_arguments.set_argument_pointer ) arg_pairs = ( (kimpy.compute_argument_name.numberOfParticles, num_particles), (kimpy.compute_argument_name.particleSpeciesCodes, species_codes), (kimpy.compute_argument_name.particleContributing, particle_contributing), (kimpy.compute_argument_name.coordinates, coords), (kimpy.compute_argument_name.partialEnergy, energy), (kimpy.compute_argument_name.partialForces, forces), ) for args in arg_pairs: set_argument_pointer(*args) ################################################################################ def get_cluster_positions(N, rball, dmin): """ Randomly generate `N` positions within a sphere of radius 'rball' (centered at the origin) where no two points are closer than 'dmin' Parameters ---------- N : int Number of atoms to generate positions for rball : float Radius of the sphere containing all of the atoms in the cluster dmin : float Minimum permissible distance between any two atoms in the cluster Returns ------- pos : numpy.ndarray Two-dimensional array containing np.double values indicating the coordinates of each atom """ def pos_has_distance_less_than_dmin(pos, dmin): for i in range(len(pos) - 1): for j in range(i + 1, len(pos)): dist = np.linalg.norm(pos[j] - pos[i]) if dist < dmin: return True return False pos = np.zeros((N, 3)) points_too_close = True while points_too_close: # Create N points uniformly distributed within a sphere of radius 'rball' for i in range(N): phi = random.uniform(0, 2 * np.pi) costheta = random.uniform(-1.0, 1.0) theta = np.arccos(costheta) x = np.sin(theta) * np.cos(phi) y = np.sin(theta) * np.sin(phi) z = np.cos(theta) rad = random.uniform(0, rball) pos[i] = (x * rad, y * rad, z * rad) points_too_close = pos_has_distance_less_than_dmin(pos, dmin) return pos ############################################################################### def format2d(arr, indent=False): """ Parameters ---------- arr : numpy.ndarray Two-dimensional array containing np.double values Returns ------- s : str String representation of the 2d array `arr` that's nicely formatted in double precision """ if indent: s = " [" else: s = "[" shape = arr.shape for row in range(shape[0]): if row == 0: s += "[" else: if indent: s += " [" else: s += " [" for col in range(shape[1]): if col < shape[1] - 1: s += "{: 1.16e} ".format(arr[row, col]) else: s += "{: 1.16e}".format(arr[row, col]) if row < shape[0] - 1: s += "]\n" else: s += "]]\n" return s ################################################################################ def report_error(vc, message=None): """ Print an error message both to stdout and to the report.txt file generated by the verification check """ formatted_message = "ERROR: {}.".format(message) vc.rwrite(formatted_message) raise RuntimeError(formatted_message) ################################################################################ def create_model(vc, model_name, length_unit, energy_unit): """ Create a Model Object with our requested units of length and energy. Charge, temperature, and time units are registered as "unused" (see KIM API documentation). Parameters ---------- vc : VerificationCheck An instance of the VerificationCheck class defined in kim_python_utils.vc. model_name : str The KIM ID of the model for which we want to create a KIM API Model Object length_unit : kimpy.length_unit.LengthUnit Units we want to use for length energy_unit : kimpy.energy_unit.EnergyUnit Units we want to use for energy Returns ------- kim_model : kimpy.model.Model KIM API Model object corresponding to the model name passed in compute_arguments : kimpy.compute_arguments.ComputeArguments An instance of the ComputeArguments class defined in kimpy.compute_arguments Raises ------ RuntimeError if the model refuses the requested units """ units_accepted, kim_model = check_kimpy_call( vc, kimpy.model.create, kimpy.numbering.zeroBased, length_unit, energy_unit, kimpy.charge_unit.unused, kimpy.temperature_unit.unused, kimpy.time_unit.unused, model_name, ) if not units_accepted: report_error( vc, "Requested units not accepted in kimpy.model.create. " "This means that the Model does not support unit conversion at " "all.", ) compute_arguments = check_kimpy_call(vc, kim_model.compute_arguments_create) return kim_model, compute_arguments ################################################################################ def check_cutoff_conversion(vc, model_name, length_conversion_tol): """ See if the model can convert its cutoffs to different unit systems correctly. This is relevant because we need to know the cutoff of the model in Angstroms in order to construct the cluster of atoms that will be used to check unit conversion for energy and/or forces. The check is performed by first asking for the cutoffs in Angstroms and manually converting the smallest of them to Bohr. This value is then compared to the smallest cutoff returned by the model when we request units of Bohr. Parameters ---------- vc : VerificationCheck An instance of the VerificationCheck class defined in kim_python_utils.vc model_name : str The KIM ID of the model for which we want to create a KIM API Model Object length_conversion_tol : float Tolerance used to determine if the Units we want to use for length Returns ------- None Raises ------ RuntimeError is the model refuses the requested units or if the model's conversion of its cutoffs does not fall within the specified tolerance """ # First, ask for the cutoffs in Angstroms kim_model, _ = create_model( vc, model_name, kimpy.length_unit.A, kimpy.energy_unit.eV ) model_cutoffs, _ = kim_model.get_neighbor_list_cutoffs_and_hints() smallest_model_cutoff_angstroms = min(model_cutoffs) # Perform a manual conversion of the cutoff we got in Angstroms to Bohr smallest_cutoff_converted = np.multiply( smallest_model_cutoff_angstroms, ANGSTROM_TO_BOHR ) # Now create another model object where we request the cutoffs in Bohr # First, ask for the cutoff in Angstroms kim_model, _ = create_model( vc, model_name, kimpy.length_unit.Bohr, kimpy.energy_unit.eV ) model_cutoffs, _ = kim_model.get_neighbor_list_cutoffs_and_hints() smallest_model_cutoff_bohr = min(model_cutoffs) # Finally, compare our manual conversion of the smallest cutoff to Bohr with what # the model gave us when we requested it in Bohr cutoff_diff = np.abs(smallest_cutoff_converted - smallest_model_cutoff_bohr) if cutoff_diff < length_conversion_tol: vc.rwrite( wrapped( "SUCCESS: Absolute difference in smallest cutoff is {} Bohr, which is " "below tolerance {}.".format(cutoff_diff, length_conversion_tol), indent=True, ) ) else: msg = ( "FAILURE: Absolute difference in smallest cutoff is {} Bohr, which " "exceeds tolerance {}.".format(cutoff_diff, length_conversion_tol) ) vc.rwrite(wrapped(msg, indent=True)) raise RuntimeError(msg) ################################################################################ def do_vc(model_name, vc): """ Primary function for checking unit conversion. A reference to this function is passed to the setup_and_run_vc function defined in kim_python_utils.vc Parameters ---------- model_name : str The KIM ID of the model on which the verification check procedure is run vc : VerificationCheck An instance of the VerificationCheck class defined in kim_python_utils.vc Returns ------- vc_grade : str Either a passing grade ("P") or a failing grade("F") indicating whether the model performs unit conversion to reasonable accuracy or not. vc_comment : str Additional information related to the outcome of the verification check. """ # Check if this is a Model or a Simulator Model. Since Simulator Models are # just wrappers around simulator-native potentials, they don't have any kind # of unit conversion and so in this case we give a grade of "N/A" col = check_kimpy_call(vc, kimpy.collections.create) model_type = check_kimpy_call(vc, col.get_item_type, model_name) if model_type != kimpy.collection_item_type.portableModel: report_error( vc, "Item {} is not a KIM Portable Model. Only KIM Portable Models are supported." "".format(model_name), ) # Tolerance used to determine if the length conversion performed by the model (used, # for example, in communicating its cutoff to the simulator) is valid. length_conversion_tol = 1e-6 # Used to determine if energy conversion performed by the model is valid. # NOTE: The model may use various combinations of parameters that have different # derived units in order to compute its energy/forces. Furthermore, # conversion constants are usually only determined to some finite precision. # Therefore, one shouldn't make the tolerance too small here energy_conversion_tol = 1e-6 # Tolerance used to decide if total energy is large enough in absolute value # for unit conversion to be non-trivial. energy_nontrivial_tol = 1e3 * energy_conversion_tol # Used to determine if forces communicated to the simulator by the model are being # converted correctly. Note that this conversion isn't done explicitly, but is rather a # consequence of the model's conversion of its parameters. # NOTE: The model may use various combinations of parameters that have different # derived units in order to compute its energy/forces. Furthermore, # conversion constants are usually only determined to some finite precision. # Therefore, one shouldn't make the tolerance too small here force_conversion_tol = 1e-6 # Tolerance used on cumulative norm of forces on all atoms to decide if forces are # large enough for unit conversion to be non-trivial. force_nontrivial_tol = 1e3 * force_conversion_tol # Check if the model converts length correctly. If it doesn't, we don't proceed # with checking conversion of energy/forces since we use the model cutoff to # construct our cluster of atoms vc.rwrite("") vc.rwrite("") vc.rwrite("Checking if model properly converts its cutoffs from Angstrom to Bohr") vc.rwrite("") check_cutoff_conversion(vc, model_name, length_conversion_tol) vc.rwrite("") # Create Model object, registering Angstroms and eV as our requested units for # length and energy vc.rwrite("") vc.rwrite( "Creating Model object with Angstrom and eV as simulator registered units" ) kim_model, compute_arguments = create_model( vc, model_name, kimpy.length_unit.A, kimpy.energy_unit.eV ) # Get the model's native units for length and energy (not used -- just outputted to # user) model_length_unit, model_energy_unit, _, _, _ = kim_model.get_units() # Check what the model can compute model_computes_energy, model_computes_forces = check_compute_arguments( vc, compute_arguments ) # Allocate memory to hold our cluster N = 5 seed = 11 random.seed(seed) coords = np.zeros((N, 3), dtype=np.double) forces = np.zeros((N, 3), dtype=np.double) energy = np.array([0.0], dtype=np.double) num_particles = np.array([N], dtype=np.intc) species_codes = np.zeros(num_particles, dtype=np.intc) particle_contributing = np.ones(num_particles, dtype=np.intc) need_neigh = np.ones(N, dtype=np.intc) # Register pointers to our memory in compute arguments object register_compute_arguments_pointers( vc, compute_arguments, num_particles, species_codes, particle_contributing, coords, energy, forces, ) # Neighbor list neigh = kimpy.neighlist.NeighList() # Register get neigh callback check_kimpy_call( vc, compute_arguments.set_callback_pointer, kimpy.compute_callback_name.GetNeighborList, kimpy.neighlist.get_neigh_kim(), neigh, ) # Influence distance and cutoffs of model model_influence_dist = kim_model.get_influence_distance() model_cutoffs, _ = kim_model.get_neighbor_list_cutoffs_and_hints() smallest_model_cutoff = min(model_cutoffs) # Minimum distance allowed for atoms during cluster generation dmin = 0.1 * smallest_model_cutoff # Get species model supports and its codes for them model_species_and_codes = get_kim_model_supported_species_and_codes(vc, kim_model) species_checked = None # Print out model info and initial simulation parameters vc.rwrite("") vc.rwrite("-" * DASHWIDTH) vc.rwrite("Results for KIM Model : {}".format(model_name.strip())) vc.rwrite( "Supported species : {}".format( " ".join(model_species_and_codes.keys()) ) ) vc.rwrite("Model computes energy : {}".format(model_computes_energy)) vc.rwrite("Model computes forces : {}".format(model_computes_forces)) vc.rwrite("Model native length unit : {}".format(str(model_length_unit))) vc.rwrite("Model native energy unit : {}".format(str(model_energy_unit))) vc.rwrite("Influence distance : {} Angstroms".format(model_influence_dist)) vc.rwrite("Model cutoffs : {} Angstroms".format(model_cutoffs)) vc.rwrite("Smallest model cutoff : {} Angstroms".format(smallest_model_cutoff)) vc.rwrite("") vc.rwrite("Parameters for random cluster generation") vc.rwrite("Random seed: {}".format(seed)) vc.rwrite("Min. distance between atoms (dmin) : {:1.16e} Angstroms".format(dmin)) vc.rwrite("-" * DASHWIDTH) vc.rwrite("") # It's sufficient only to check one species that the model supports. # However, we don't know for sure whether some pure species interactions # are trivially zero. So, we start looping over all of the species until # we find one for which we can generate a 5-atom cluster that has sizeable # energy and/or forces (depending on what the model computes) vc.rwrite("") vc.rwrite("+" * DASHWIDTH) vc.rwrite("") vc.rwrite(" STEP I: Cluster generation") vc.rwrite("") step1_description = ( "Attempting to find a species supported by the Model for which a " "random cluster can be generated that has significant energy " "and/or forces (depending on what the Model is able to compute). We " "construct the cluster by first placing a single atom at the origin " "and then generating additional atom positions that fall within a " "distance 'rball' of the origin; no atoms are allowed to fall within " "a distance of 'dmin' of one another. Note that 'rball' is chosen to " "be 0.9 * the smallest cutoff reported by the Model in Angstroms. In " "the event that the Model's native length unit is not Angstroms and " "the cutoffs it reports are improperly converted to Angstroms, it's " "possible that the cluster generated may actually be far too large " "for atomic interactions to occur. Accordingly, if the energy " "and/or forces are not sizable using the initial value of 'rball', " "we regenerate the cluster positions iteratively using smaller and " "smaller values of 'rball' (while keeping 'dmin' fixed); if the " "energy and/or forces still aren't sizeable when generating a " "cluster using rball = 0.2 * the smallest cutoff reported by the " "Model, that species is skipped and the process is repeated for " "the next species." ) vc.rwrite(wrapped(step1_description, indent=True)) vc.rwrite("") vc.rwrite("+" * DASHWIDTH) vc.rwrite("") for spec, spec_code in model_species_and_codes.items(): for rball_factor in np.arange(0.2, 1, 0.1)[::-1]: rball = rball_factor * smallest_model_cutoff vc.rwrite( wrapped( "Creating random cluster of {} {} atoms with a diameter {:6.6f} " "Angstroms and a minimum separation of {:6.6f} Angstroms " "between all atoms".format(N, spec, rball, dmin) ) ) vc.rwrite("") cluster_pos = get_cluster_positions(N, rball, dmin) vc.rwrite(wrapped("Cluster positions (Angstroms):", indent=True)) vc.rwrite(format2d(cluster_pos, indent=True)) species_codes[:] = spec_code np.copyto(coords, cluster_pos) # preserves address of coords # Construct neigh list check_kimpy_call( vc, neigh.build, coords, model_influence_dist, np.asarray([1.1 * smallest_model_cutoff], dtype=np.double), need_neigh, ) # Call Model's compute function to get energy and/or forces vc.rwrite(wrapped("Calling Model compute...", indent=True)) vc.rwrite("") kim_model.compute(compute_arguments) initial_energy = energy[0] if model_computes_energy: vc.rwrite( wrapped( "Model returned energy of {:16.16e} eV" "".format(initial_energy), indent=True, ) ) vc.rwrite("") initial_forces = np.copy(forces) force_norm = np.linalg.norm(forces) if model_computes_forces: vc.rwrite(wrapped("Model returned forces (eV/Angstrom):", indent=True)) vc.rwrite(format2d(forces, indent=True)) vc.rwrite( wrapped( "Cumulative norm of all forces: {:16.16e} " "eV/Angstrom".format(force_norm), indent=True, ) ) vc.rwrite("") if model_computes_energy and model_computes_forces: if ( np.abs(initial_energy) > energy_nontrivial_tol and force_norm > force_nontrivial_tol ): vc.rwrite( wrapped( "SUCCESS: Absolute value of energy exceeds tolerance {} eV and " "cumulative norm of forces on all atoms exceeds tolerance " "{} eV/Angstrom. Proceeding with unit conversion..." "".format(energy_nontrivial_tol, force_nontrivial_tol), indent=True, ) ) vc.rwrite("") species_checked = spec species_code_checked = spec_code break # from rball loop if np.abs(initial_energy) > energy_nontrivial_tol: vc.rwrite( wrapped( "Absolute value of energy exceeds tolerance {} eV, but " "cumulative norm of forces on all atoms is below " "tolerance {} eV/Angstrom.".format( energy_nontrivial_tol, force_nontrivial_tol ), indent=True, ) ) vc.rwrite("") elif force_norm > force_nontrivial_tol: vc.rwrite( wrapped( "Cumulative norm of forces on all atoms is above tolerance " "{} eV/Angstrom but absolute value of energy is below tolerance " "{} eV.".format( force_nontrivial_tol, energy_nontrivial_tol ), indent=True, ) ) vc.rwrite("") elif model_computes_energy: if np.abs(initial_energy) > energy_nontrivial_tol: vc.rwrite( wrapped( "SUCCESS: Absolute value of energy exceeds tolerance {} eV. Proceeding with unit " "conversion...".format(energy_nontrivial_tol), indent=True, ) ) vc.rwrite("") species_checked = spec species_code_checked = spec_code break # from rball loop vc.rwrite( wrapped( "Absolute value of energy is below tolerance {} eV.".format( energy_nontrivial_tol ), indent=True, ) ) vc.rwrite("") elif model_computes_forces: if force_norm > force_nontrivial_tol: vc.rwrite( wrapped( "SUCCESS: Cumulative norm of forces on all atoms exceeds " "tolerance {} eV/Angstrom. Proceeding with unit conversion." "".format(force_nontrivial_tol), indent=True, ) ) vc.rwrite("") species_checked = spec species_code_checked = spec_code break # from rball loop vc.rwrite( wrapped( "Cumulative norm of forces on all atoms is below tolerance " "{} eV/Angstrom.".format(force_nontrivial_tol), indent=True, ) ) vc.rwrite("") if species_checked is not None: break # from species loop # Print informational message vc.rwrite( wrapped( "Could not generate any clusters for species {} for which " "the energy and/or forces were sizeable enough to proceed with " "unit conversion." ) ) # Couldn't find any species that would work if species_checked is None: vc.rwrite("") report_error( vc, "Could not find any species supported by the Model for " "which a random cluster with significant forces on the atoms could be " "generated to check unit conversion.", ) # Now repeat the calculation for a set of units different from what we originally used vc.rwrite("") vc.rwrite("+" * DASHWIDTH) vc.rwrite("") vc.rwrite( " STEP II: Convert units and ask model to recompute energy and/or forces" ) vc.rwrite("") step2_description = ( "Take the set of coordinates that was found to have sizeable " "energy/forces from Step I and change the units from Angstroms to Bohr " "ourselves. Next, create a new Model object wherein we've registered " "Bohr and Hartree as the unit set. Ask the model for its neighbor list " "cutoffs and influence once more and use them to create a new neighbor " "list using the coordinates we converted to Bohr. Finally, have the " "Model compute the energy and/or forces once more and compare the " "results to our own conversion of the energy and/or forces the model " "returned in Step I." ) vc.rwrite(wrapped(step2_description, indent=True)) vc.rwrite("") vc.rwrite("+" * DASHWIDTH) vc.rwrite("") vc.rwrite( "Creating Model object with Bohr and Hartree as simulator registered units" ) kim_model, compute_arguments = create_model( vc, model_name, kimpy.length_unit.Bohr, kimpy.energy_unit.Hartree ) # Register pointers to our memory in compute arguments object register_compute_arguments_pointers( vc, compute_arguments, num_particles, species_codes, particle_contributing, coords, energy, forces, ) # Neighbor list neigh = kimpy.neighlist.NeighList() # Register get neigh callback check_kimpy_call( vc, compute_arguments.set_callback_pointer, kimpy.compute_callback_name.GetNeighborList, kimpy.neighlist.get_neigh_kim(), neigh, ) # Influence distance and cutoffs of model model_influence_dist = kim_model.get_influence_distance() model_cutoffs, _ = kim_model.get_neighbor_list_cutoffs_and_hints() smallest_model_cutoff = min(model_cutoffs) # Convert coords and neighbor list cutoff from Angstrom to Bohr tmp = np.multiply(coords, ANGSTROM_TO_BOHR) np.copyto(coords, tmp) # preserves address of coords # Recreate neighbor list for our new coords using the Model's newly # reported cutoff (if the model works correctly, the neighbor list # should turn out the same) check_kimpy_call( vc, neigh.build, coords, model_influence_dist, np.asarray([1.1 * smallest_model_cutoff], dtype=np.double), need_neigh, ) vc.rwrite("") vc.rwrite(wrapped("Cluster positions converted to Bohr:")) vc.rwrite(format2d(coords)) # Set species codes species_codes[:] = species_code_checked # Call Model's compute function vc.rwrite(wrapped("Calling Model compute...")) vc.rwrite("") kim_model.compute(compute_arguments) vc.rwrite( wrapped("Model returned energy of {:16.16e} " "Hartree".format(energy[0])) ) vc.rwrite("") vc.rwrite("Model returned forces (Hartree/Bohr):") vc.rwrite(format2d(forces)) force_norm = np.linalg.norm(forces) vc.rwrite( wrapped( "Cumulative norm of all forces: {:16.16e} Hartree/Bohr".format(force_norm) ) ) vc.rwrite("") vc.rwrite( wrapped( "Converting energy and forces returned by model back " "to Angstroms and eV..." ) ) vc.rwrite("") # Convert new energy from Hartree to eV converted_energy = np.multiply(energy[0], HARTREE_TO_EV) vc.rwrite(wrapped("Converted energy is {:16.16e} eV".format(converted_energy))) # Convert new forces from Hartree/Bohr to eV/Angstrom converted_forces = np.multiply(forces, HARTREE_TO_EV * ANGSTROM_TO_BOHR) vc.rwrite("") vc.rwrite("Converted forces (eV/Angstrom):") vc.rwrite(format2d(converted_forces)) vc.rwrite( wrapped( "Cumulative norm of converted forces: {:16.16e} eV/Angstrom".format( np.linalg.norm(converted_forces) ) ) ) # Compute grade vc.rwrite("") vc.rwrite("=" * DASHWIDTH) vc.rwrite( wrapped( "To pass this verification check, the Model must be able to compute " "the energy and/or forces for a random cluster of atoms in two different " "sets of units. If energies are computed, the results must match within " "{} eV. If forces are computed, the cumulative norm of the difference " "between the forces computed in one set of units and the other must be " "less than {} eV/Angstrom.".format( energy_conversion_tol, force_conversion_tol ) ) ) energy_diff = np.abs(converted_energy - initial_energy) force_diff_norm = np.linalg.norm(converted_forces - initial_forces) if model_computes_energy: vc.rwrite("") vc.rwrite( wrapped( "Absolute difference between initial energy and converted energy: {} eV" "".format(energy_diff) ) ) vc.rwrite("") if model_computes_forces: vc.rwrite( wrapped( "Cumulative norm of difference in initial forces and converted " "forces: {} eV/Angstrom".format(force_diff_norm) ) ) vc.rwrite("") if model_computes_energy and model_computes_forces: if ( energy_diff < energy_conversion_tol and force_diff_norm < force_conversion_tol ): vc.rwrite( wrapped( "SUCCESS: Absolute difference in energy is below tolerance {} and " "cumulative norm of difference in forces on all atoms\nis below " "tolerance {}.".format(energy_conversion_tol, force_conversion_tol) ) ) vc.rwrite("") vc_grade = "P" vc_comment = "Unit conversion successful" elif energy_diff < energy_conversion_tol: vc.rwrite( wrapped( "FAILURE: Absolute difference in energy is below tolerance {} but " "cumulative norm of difference in forces on all atoms\nexceeds " "tolerance {}.".format(energy_conversion_tol, force_conversion_tol) ) ) vc.rwrite("") vc_grade = "F" vc_comment = "Unit conversion failed" elif force_diff_norm < force_conversion_tol: vc.rwrite( wrapped( "FAILURE: Cumulative norm of difference in forces on all atoms is " "below tolerance {} but absolute difference in energy\nexceeds " "tolerance {}.".format(force_conversion_tol, energy_conversion_tol) ) ) vc.rwrite("") vc_grade = "F" vc_comment = "Unit conversion failed" else: vc.rwrite( wrapped( "FAILURE: Absolute difference in energy exceeds tolerance {} and " "cumulative norm of difference in forces on all atoms\nexceeds " "tolerance {}.".format(energy_conversion_tol, force_conversion_tol) ) ) vc.rwrite("") vc_grade = "F" vc_comment = "Unit conversion failed" elif model_computes_energy: if energy_diff < energy_conversion_tol: vc.rwrite( wrapped( "SUCCESS: Absolute difference in energy is below tolerance " "{}.".format(energy_conversion_tol) ) ) vc.rwrite("") vc_grade = "P" vc_comment = "Unit conversion successful" else: vc.rwrite( wrapped( "FAILURE: Absolute difference in energy exceeds tolerance {}.".format( energy_conversion_tol ) ) ) vc.rwrite("") vc_grade = "F" vc_comment = "Unit conversion failed" elif model_computes_forces: if force_diff_norm < force_conversion_tol: vc.rwrite( wrapped( "SUCCESS: Cumulative norm of difference in forces on all atoms is " "below tolerance {}.".format(force_conversion_tol) ) ) vc.rwrite("") vc_grade = "P" vc_comment = "Unit conversion successful" else: vc.rwrite( wrapped( "FAILURE: Cumulative norm of difference in forces on all atoms " "exceeds tolerance {}.".format(force_conversion_tol) ) ) vc.rwrite("") vc_grade = "F" vc_comment = "Unit conversion failed" return vc_grade, vc_comment ################################################################################ # # MAIN PROGRAM # ############################################################################### if __name__ == "__main__": vc_args = { "vc_name": "vc-unit-conversion", "vc_author": __author__, "vc_description": kim_vc_utils.vc_stripall(__doc__), "vc_category": "mandatory", "vc_grade_basis": "passfail", "vc_files": [], "vc_debug": False, # Set to True to get exception traceback info } # Get the model extended KIM ID: model_name = input("Model Extended KIM ID = ") # Execute VC kim_vc_utils.setup_and_run_vc(do_vc, model_name, **vc_args)