#!/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, 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)."""
# Python 2-3 compatible code issues
from __future__ import print_function
try:
input = raw_input
except NameError:
pass
import numpy as np
import random
import textwrap
import kimpy
import kim_python_utils.vc as kim_vc_utils
__version__ = "000"
__author__ = "Daniel S. Karls"
# 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))
else:
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, error = kimpy.species_name.get_species_name(i)
check_error(vc, error, "kimpy.species_name.get_species_name")
species_support, code, error = kim_model.get_species_support_and_code(
species_name
)
check_error(vc, error, "kim_model.get_species_support_and_code")
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:
species_support, code, error = kim_model.get_species_support_and_code(
kimpy.species_name.SpeciesName(s)
)
check_error(
vc, error or not species_support, "kim_model.get_species_support_and_code"
)
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, error = compute_arguments.get_argument_support_status(
kimpy.compute_argument_name.partialEnergy
)
check_error(vc, error, "compute_argument.get_argument_support_status")
model_computes_forces, error = compute_arguments.get_argument_support_status(
kimpy.compute_argument_name.partialForces
)
check_error(vc, error, "compute_argument.get_argument_support_status")
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
"""
error = compute_arguments.set_argument_pointer(
kimpy.compute_argument_name.numberOfParticles, num_particles
)
check_error(vc, error, "kimpy.compute_argument.set_argument_pointer")
error = compute_arguments.set_argument_pointer(
kimpy.compute_argument_name.particleSpeciesCodes, species_codes
)
check_error(vc, error, "kimpy.compute_argument.set_argument_pointer")
error = compute_arguments.set_argument_pointer(
kimpy.compute_argument_name.particleContributing, particle_contributing
)
check_error(vc, error, "kimpy.compute_argument.set_argument_pointer")
error = compute_arguments.set_argument_pointer(
kimpy.compute_argument_name.coordinates, coords
)
check_error(vc, error, "kimpy.compute_argument.set_argument_pointer")
error = compute_arguments.set_argument_pointer(
kimpy.compute_argument_name.partialEnergy, energy
)
check_error(vc, error, "kimpy.compute_argument.set_argument_pointer")
error = compute_arguments.set_argument_pointer(
kimpy.compute_argument_name.partialForces, forces
)
check_error(vc, error, "kimpy.compute_argument.set_argument_pointer")
################################################################################
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 check_error(vc, error, func_name):
"""
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.
error : int
Error code returned by kimpy functions indicating whether the operation succeeded
(0) or failed.
func_name : str
Name of the function that returned error code `error`
Returns
-------
None
"""
if error != 0 and error is not None:
formatted_message = 'ERROR: Calling "{}" failed with error {}.'.format(
func_name, error
)
if vc is None:
print(formatted_message)
else:
vc.rwrite(formatted_message)
raise RuntimeError(formatted_message)
################################################################################
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, modelname, 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.
modelname : 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 is the model refuses the requested units
"""
units_accepted, kim_model, error = kimpy.model.create(
kimpy.numbering.zeroBased,
length_unit,
energy_unit,
kimpy.charge_unit.unused,
kimpy.temperature_unit.unused,
kimpy.time_unit.unused,
modelname,
)
check_error(vc, error, "kimpy.model.create")
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, error = kim_model.compute_arguments_create()
check_error(vc, error, "kim_model.compute_arguments_create")
return kim_model, compute_arguments
################################################################################
def destroy_model(vc, kim_model, compute_arguments):
"""
Deallocate memory for KIM API Model Object and ComputeArguments objects
Parameters
----------
vc : VerificationCheck
An instance of the VerificationCheck class defined in kim_python_utils.vc.
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
Returns
-------
None
"""
# Destroy compute arguments object
error = kim_model.compute_arguments_destroy(compute_arguments)
check_error(vc, error, "kim_model.compute_arguments_destroy")
# Destroy model object
error = kimpy.model.destroy(kim_model)
check_error(vc, error, "kim_model.destroy")
################################################################################
def check_cutoff_conversion(vc, modelname, 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
modelname : 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, compute_arguments = create_model(
vc, modelname, 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)
# Destroy this model object
destroy_model(vc, kim_model, compute_arguments)
# 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, compute_arguments = create_model(
vc, modelname, 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)
# Destroy this model object
destroy_model(vc, kim_model, compute_arguments)
# 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(modelname, 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
----------
modelname : 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
df
vc_comment : str
"""
# 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, err = kimpy.collections.create()
check_error(vc, err, "kimpy.collections.create")
model_type, err = col.get_item_type(modelname)
check_error(vc, err, "kimpy.collections.Collections.get_item_type")
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(modelname),
)
kimpy.collections.destroy(col)
# 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, modelname, 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, modelname, 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.initialize()
# Register get neigh callback
error = compute_arguments.set_callback_pointer(
kimpy.compute_callback_name.GetNeighborList,
kimpy.neighlist.get_neigh_kim(),
neigh,
)
check_error(vc, error, "kimpy.compute_argument.set_callback")
# 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(modelname.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 which 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 addtional 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
error = kimpy.neighlist.build(
neigh,
coords,
model_influence_dist,
np.asarray([1.1 * smallest_model_cutoff], dtype=np.double),
need_neigh,
)
check_error(vc, error, "kimpy.neighlist.build")
# 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
elif 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
else:
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
else:
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 None:
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.",
)
# Destroy neigh list
kimpy.neighlist.clean(neigh)
# Clean up model and compute arguments object
destroy_model(vc, kim_model, compute_arguments)
# 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, modelname, 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.initialize()
# Register get neigh callback
error = compute_arguments.set_callback_pointer(
kimpy.compute_callback_name.GetNeighborList,
kimpy.neighlist.get_neigh_kim(),
neigh,
)
check_error(vc, error, "kimpy.compute_argument.set_callback")
# 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)
error = kimpy.neighlist.build(
neigh,
coords,
model_influence_dist,
np.asarray([1.1 * smallest_model_cutoff], dtype=np.double),
need_neigh,
)
check_error(vc, error, "kimpy.neighlist.build")
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("")
# Destroy neigh list
kimpy.neighlist.clean(neigh)
# Clean up model and compute arguments object
destroy_model(vc, kim_model, compute_arguments)
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__":
vcargs = {
"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:
modelname = input("Model Extended KIM ID = ")
# Execute VC
kim_vc_utils.setup_and_run_vc(do_vc, modelname, **vcargs)