#!/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) 2020-2021, Regents of the University of Minnesota.
# All rights reserved.
#
# Contributor(s):
# Ellad B. Tadmor
# Daniel S. Karls
#
################################################################################
# The docstring below is vc_description
"""This verification check examines whether a model called in parallel by
multiple threads gives the same results as when called sequentially. A number
`num_configs` (preset in the code) of configurations is generated, each
containing a different number of atoms based on a randomly distorted,
periodic, face-centered cubic (fcc) structure containing a random distribution
of all atoms supported by the model. Configurations used for testing are
provided as auxiliary files. The energy and forces for each configuration are
computed in sequence. Then the calculations are repeated using Python
multithreading with `num_configs` threads and the Global Interpreter
Lock (GIL) released to ensure true parallelism. Each thread possesses its
own copy of memory related to the calculation. The threaded calculations
are repeated `num_cycles` times (preset in the code) with the configurations
randomly distributed to the threads to create many varying opportunities
for collisions. To pass the test, the total energy and the forces on all atoms
for a given configuration obtained in all cycles must be be identical to those
obtained in the sequential calculation. Failure of this verification check
implies that the model is inappropriately storing information in persistent
(static) memory during or between calls. Passing this verification check
provides some assurance the the model is thread-safe, but it is NOT a guarantee
due to the inherent randomness of race conditions in unsafe code. To be certain,
the model code would need to be studied and analyzed carefully."""
# Python 2-3 compatible code issues
from __future__ import print_function
from multiprocessing import cpu_count
import random
import threading
import numpy as np
import scipy.optimize
from ase.lattice.cubic import FaceCenteredCubic
from ase.calculators.kim import KIM, get_model_supported_species
import kim_python_utils.ase as kim_ase_utils
import kim_python_utils.vc as kim_vc_utils
import kimpy
try:
input = raw_input
except NameError:
pass
__version__ = "005"
__author__ = "Ellad Tadmor and Daniel S. Karls"
class KimpyError(Exception):
"""
A call to a kimpy function resulted in a RuntimeError being raised
"""
pass
def check_kimpy_call(f, *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.
(Starting with kimpy 2.0.0, a RuntimeError is the only exception type raised
when something goes wrong.)
"""
try:
return f(*args)
except RuntimeError as exception:
raise KimpyError(f'Calling kimpy function "{f.__name__}" failed:\n {str(exception)}')
################################################################################
#
# FUNCTIONS
#
################################################################################
def cubic_cell_energy(alat, atoms, ncells_per_side):
"""
Calculate the energy of the passed 'atoms' structure containing a
cubic structure with 'ncells_per_side'. Scale to lattice constant
'alat' (passed as a nd array of length 1) and return the energy.
"""
acell = alat[0] * ncells_per_side
atoms.set_cell([acell, acell, acell], scale_atoms=True)
return atoms.get_potential_energy()
################################################################################
def compute_config(atoms):
"""
Calculate the energy and forces for the passed 'atoms' structure.
"""
# This forces a calculation of the model
# if anything has changed and the results
# energy, forces, etc. are stored in atoms.
atoms.get_potential_energy()
################################################################################
def do_vc(model, vc):
"""
Perform Thread Safety VC
"""
# Get supported species
species = get_model_supported_species(model)
species = kim_ase_utils.remove_species_not_supported_by_ASE(list(species))
species.sort()
# Basic VC parameters
num_cores = cpu_count()
num_configs = 10
num_cycles = 10
lattice_constant = 2.5
max_lattice_constant = 10.0
mincells_per_side = 2
maxcells_per_side = 5
seed = 13
random.seed(seed)
# Print VC info
agree_text = {True: "OK", False: "Failed"}
dashwidth = 80
vc.rwrite("")
vc.rwrite("-" * dashwidth)
vc.rwrite("Results for KIM Model : %s" % model.strip())
vc.rwrite("Supported species : %s" % " ".join(species))
vc.rwrite("")
vc.rwrite("Number of cores = %d" % num_cores)
vc.rwrite("number of configurations = %d" % num_configs)
vc.rwrite("number of cycles = %d" % num_cycles)
vc.rwrite("random seed = %d" % seed)
vc.rwrite("minimum unit cells per side = %d" % mincells_per_side)
vc.rwrite("maximum unit cells per side = %d" % maxcells_per_side)
vc.rwrite("-" * dashwidth)
vc.rwrite("")
# Check the model type. VC can only be performed for KIM Models
col = check_kimpy_call(kimpy.collections.create)
model_type = check_kimpy_call(col.get_item_type, model)
if model_type != kimpy.collection_item_type.portableModel:
vc.rwrite(model.strip() + " is not a KIM Model.")
vc.rwrite("")
vc.rwrite(
"Thread safety can only be checked for models compatible with the KIM API."
)
vc.rwrite("Check not performed.")
vc_grade = "N/A"
vc_comment = "Thread safety can only be checked for KIM Models."
return vc_grade, vc_comment
# Make sure there are more than one processor available
if num_cores == 1:
vc.rwrite(
"Only one core available on system. Unable to perform thread safety verification check."
)
raise RuntimeError(
"Only one core available on system. Unable to perform thread safety verification check."
)
# Create num_config randomized configurations containing all
# supported species
vc.rwrite("SINGLE THREAD REFERENCE CALCULATIONS")
vc.rwrite("")
vc.rwrite(
"{0:>6} {1:>6} {2:>13} {3:>13}".format(
"Config", "#Atoms", "Energy", "Ave Norm"
)
)
vc.rwrite("-" * 47)
alat = lattice_constant
config_atoms = []
config_energy = []
config_forces = []
for n in range(num_configs):
# Randomly set system size based on allowed range
ncells_per_side = random.randint(mincells_per_side, maxcells_per_side)
# Since this VC only applies to KIM portable models, we know that the
# calculator returned by function "KIM" will be a KIMCalculator instance,
# which has context manager functions (__enter__ and __exit__) defined.
# This will ensure that the destroy function of the model will be called,
# which is helpful since some models may have global variables affected
# by their create and destroy functions which could impact thread safety.
with KIM(model) as calc:
# Compute equilibrium lattice constant
not_done = True
got_equil = False
while not_done:
while True:
atoms = FaceCenteredCubic(
size=(ncells_per_side, ncells_per_side, ncells_per_side),
latticeconstant=alat,
symbol="H",
pbc=False,
)
# Change to the cell size in the following block should
# be run at most once when a selected size is too small
if len(atoms) < len(species):
ncells_per_side += 1
else:
break
# Randomize the species in the cluster using all of the species supported by
# the Model
kim_ase_utils.randomize_species(atoms, species)
# Attach calculator
atoms.set_calculator(calc)
# Find equilibrium lattice constant
try:
res = scipy.optimize.minimize(
cubic_cell_energy,
alat,
args=(atoms, ncells_per_side),
method="Nelder-Mead",
tol=1e-6,
)
alat = res.x[0]
not_done = False
got_equil = True
except: # noqa: E722
# failed for some reason (assume it's because of KIM error)
alat += 0.5
if alat > max_lattice_constant:
not_done = False
if not got_equil:
raise RuntimeError(
"Unable to compute equilibrium lattice constant for one of the configurations."
)
# Rescale and perturb crystal
kim_ase_utils.rescale_to_get_nonzero_forces(atoms, 0.01)
pert_amp = 0.01 * alat
kim_ase_utils.randomize_positions(atoms, pert_amp)
kim_ase_utils.perturb_until_all_forces_sizeable(atoms, pert_amp)
# Write configuration files
aux_file = (
"config-" + str(n).strip() + "-config-" + "".join(species) + ".xyz"
)
vc.vc_files.append(aux_file)
vc.write_aux_ase_atoms(aux_file, atoms, "xyz")
# Store configuration, energy, and forces
config_atoms.append(atoms.copy())
config_energy.append(atoms.get_potential_energy())
config_forces.append(atoms.get_forces())
# Instantiate new calculator for this configuration that releases the GIL so the
# threads will be able to execute in parallel in the next stage of this VC.
#
# NOTE: While defining a different calculator for each configuration means that
# there will be multiple KIM API Model Objects (one for each configuration), only
# a *single* copy of the model's shared library will be loaded into the memory
# space of this python process. This means that if the model is doing something
# thread-unsafe, we'll be able to detect it.
config_atoms[n].set_calculator(KIM(model, options={"release_GIL": True}))
# Report to user
natoms = len(config_atoms[n])
energy = config_energy[n]
avenorm = np.linalg.norm(config_forces[n]) / natoms
vc.rwrite(
"{0:6d} {1:6d} {2:13.6e} {3:13.6e}".format(n, natoms, energy, avenorm)
)
# Close table
vc.rwrite("-" * 47)
vc.rwrite("")
# Cycle through configuration calculations num_cycles times
# using num_config threads to compute all configurations.
vc.rwrite("THREADED CALCULATIONS")
vc.rwrite("")
for cycle in range(num_cycles):
vc.rwrite("### Cycle {} ###".format(cycle + 1))
vc.rwrite("")
# Randomize order of configurations for thread calculations.
# This way the configurations are visited in different order
# during each cycle.
order = list(range(num_configs))
random.shuffle(order)
# Change positions, request energy which will issue a calculate
# operation, then restore positions. This will forcce the threaded
# calculations to recalculate and not use the energy and forces
# stored in the config_atoms[] objects from previous cycle
for n in range(num_configs):
posns = config_atoms[n].get_positions()
config_atoms[n].positions += 1.0
config_atoms[n].get_potential_energy()
config_atoms[n].set_positions(posns)
# Create threads. Each thread will compute one configuration.
threads = []
num_threads = num_configs
for n in range(num_threads):
thr = threading.Thread(
name="thread-" + str(n),
target=compute_config,
args=(config_atoms[order[n]],),
)
threads.append(thr)
thr.start()
# Start all threads and wait until they're all done executing
for thr in threads:
thr.join()
# All threads are complete at this points. Each configuration
# contains the energy and forces computed by a thread. Compare
# these with the original value computed in the earlier reference
# sequential calculation and make sure there are no discrepancies.
vc.rwrite(
"{0:>6} {1:>6} {2:>13} {3:>13} {4:>6} {5:>6}".format(
"Config", "#Atoms", "Energy", "Ave Norm", "Thread", "Status"
)
)
vc.rwrite("-" * 65)
all_agree = True
for n in range(num_configs):
natoms = len(config_atoms[n])
energy = config_atoms[n].get_potential_energy()
forces = config_atoms[n].get_forces()
avenorm = np.linalg.norm(forces) / natoms
thrnum = order.index(n) # number of thread that processed config n
agree = np.allclose(energy, config_energy[n]) and np.allclose(
forces, config_forces[n]
)
agree_s = agree_text[agree]
vc.rwrite(
"{0:6d} {1:6d} {2:13.6e} {3:13.6e} {4:6d} {5:<6}".format(
n, natoms, energy, avenorm, thrnum, agree_s
)
)
all_agree = all_agree and agree
# Close table
vc.rwrite("-" * 65)
vc.rwrite("")
# Stop cycling if an error is detected
if not all_agree:
break
# Report grade
vc.rwrite("=" * dashwidth)
vc.rwrite(
"To pass this verification check all threads must give identical results,"
)
vc.rwrite(
"i.e. the same total energy and the same average norm (force norm divided"
)
vc.rwrite(
'by the number of atoms). This is indicated by the "Status" column in the'
)
vc.rwrite("above table.")
vc.rwrite("")
if all_agree:
vc_grade = "P"
vc_comment = (
"All threads give identical results for tested case. Model "
"appears to be thread-safe."
)
else:
vc_grade = "F"
vc_comment = (
"One or more threads gave different results than a single thread "
"calculation. The model is not thread-safe."
)
return vc_grade, vc_comment
################################################################################
#
# MAIN PROGRAM
#
###############################################################################
if __name__ == "__main__":
vc_args = {
"vc_name": "vc-thread-safe",
"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 = input("Model Extended KIM ID = ")
# Execute VC
kim_vc_utils.setup_and_run_vc(do_vc, model, **vc_args)