#!/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):
# Ellad B. Tadmor
#
################################################################################
# The docstring below is vc_description
"""This verification check ensures that a model supports all of the species that it
claims to. The verification check attempts to compute the energy for all clusters up to
a preset size for all combinations of the species listed by the model. Several random
samples are generated for each cluster. The check is successful if all elements appear
in at least one successful calculation, i.e. the model does not abort during the
calculation. This algorithm accounts for the fact that some models may only support
species in certain combinations. For example a model supporting species A and B may be
for an AB alloy and may only support A-B interactions and not A-A or B-B. The output
provides full information on which interactions are supported and which are not. Note
that this definition of 'support' does not imply that the model has non-trivial energy
and/or force interactions for any of the species."""
# Python 2-3 compatible code issues
from __future__ import print_function
import random
import itertools
import numpy as np
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
from ase import Atoms
__version__ = "002"
__author__ = "Ellad Tadmor"
try:
input = raw_input
except NameError:
pass
################################################################################
#
# FUNCTIONS
#
################################################################################
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
################################################################################
def get_cluster_positions(N, rball, dmin):
"""
Select 'N' positions within a sphere of radius 'rball' where no two points
are closer than 'rmin'.
"""
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 do_vc(model, vc):
"""
Perform Species Supported as Stated VC
"""
# Get supported species
species = get_model_supported_species(model)
species = kim_ase_utils.remove_species_not_supported_by_ASE(list(species))
species.sort()
# Initialize
Nmax = 4 # clusters ranging from 2 to Nmax atoms will be
# tested (as long as 'max_combinations' is not
# exeeded.
max_combinations = 10000 # maximum number of allowed combinations
rball = 1.5 # radius of ball containing clusters
dmin = 1.0 # minimum allowed distance between atoms
Nsamp = 10 # number of samples of each cluster
seed = 13
random.seed(seed)
species_is_supported = {}
for spec in species:
species_is_supported[spec] = False # Assume none of the species
# are supported and now test
yesno = {True: "yes", False: "** NO **"}
# Finite domain in which to embed the finite cluster of atoms
large_cell_len = 7 * rball
# Print VC Header
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("Maximum cluster size considered (#atoms) = %d" % Nmax)
vc.rwrite("Maximum number of allowed combinations = %d" % max_combinations)
vc.rwrite("Radius of ball containing clusters: = %0.3f" % rball)
vc.rwrite("Minimum allowed distance between atoms: = %0.3f" % dmin)
vc.rwrite("Number of random samples for each cluster = %d" % Nsamp)
vc.rwrite("Random number seed: = %d" % seed)
vc.rwrite("-" * dashwidth)
vc.rwrite("")
for spec in species:
# Check if this species has non-trivial force and energy interactions
atoms_interacting_energy, atoms_interacting_force = kim_ase_utils.check_if_atoms_interacting(
model, symbols=[spec, spec]
)
if not atoms_interacting_energy:
vc.rwrite("")
vc.rwrite(
"INFO: The model provided, {}, does not possess a non-trivial energy "
"interaction for species {}, i.e. the energy computed by the model for "
"a given set of atoms of this species is always equal to the energy "
"it computes for an isolated atom of this species multiplied by the "
"number of atoms present. This message is strictly informational and "
"will not affect the results of this Verification Check.".format(
model, spec
)
)
vc.rwrite("")
if not atoms_interacting_force:
vc.rwrite("")
vc.rwrite(
"INFO: The model provided, {}, does not possess a non-trivial force "
"interaction for species {} i.e. the forces computed by the model for "
"a given set of atoms of this species will always be zero. This "
"message is strictly informational, and will not affect the results of "
"this Verification Check.".format(model, spec)
)
vc.rwrite("")
# Loop over all clusters from 2 to Nmax considering all possible
# species combinations without considering order.
# For example for a model supporting 3 species (A,B,C) the loop would
# be over 2-atom and 3-atom clusters with the following combinations:
# A-A, B-B, C-C, A-B, A-C, B-C
# A-A-A, B-B-B, C-C-C, A-A-B, A-A-C, B-B-A, B-B-C, C-C-A, C-C-B, A-B-C
# Test whether the model supports the configuration. The model passes
# the check if all species appear in at least one supported configuration.
# Report this combination
vc.rwrite("%-*s %7s %s" % (2 * Nmax, "Cluster", "Success", "Supported"))
vc.rwrite("-" * (2 * Nmax + 24))
for N in range(2, Nmax + 1):
# Create all combintations of species for a cluster of size N
species_cluster = list(itertools.combinations_with_replacement(species, N))
# Compute number of combinations and skip if too many
num_combinations = len(species_cluster)
if num_combinations > max_combinations:
vc.rwrite(
"...... skipping %d-cluster (%d combinations)" % (N, num_combinations)
)
continue
# Loop over combinations
for i in range(len(species_cluster)):
# Loop over samples
num_success = 0
for n in range(Nsamp):
# Assign positions to atoms in a cluster of N atoms in a ball of
# radius 'rball' with no atoms closer than 'dmin'.
pos = get_cluster_positions(N, rball, dmin)
# Define Atoms object for cluster
atoms = Atoms(
"".join(species_cluster[i]),
positions=pos,
cell=(2 * rball, 2 * rball, 2 * rball),
pbc=(False, False, False),
)
calc = KIM(model)
# Move our finite cluster of atoms to the center of our large cell
atoms.set_cell([large_cell_len, large_cell_len, large_cell_len])
trans = [0.5 * large_cell_len] * 3
atoms.translate(trans)
# Try to compute energy
try:
atoms.set_calculator(calc)
atoms.get_potential_energy()
num_success += 1
for spec in species_cluster[i]:
species_is_supported[spec] = True
except: # noqa: E722
pass
finally:
# Explicitly close calculator to ensure any allocated memory is freed
# (relevant for SMs)
if hasattr(calc, "__del__"):
calc.__del__()
# Report this combination
vc.rwrite(
"%-*s %2d / %2d %s"
% (
2 * Nmax,
"".join(species_cluster[i]),
num_success,
Nsamp,
yesno[num_success > 0],
)
)
vc.rwrite("-" * (2 * Nmax + 24))
vc.rwrite("")
vc.rwrite("=" * dashwidth)
vc.rwrite("To pass this verification check each model element must be supported")
vc.rwrite("in at least one cluster combination.")
vc.rwrite("")
if all(value for value in species_is_supported.values()):
vc_grade = "P"
vc_comment = "All elements claimed by model are supported in at least one of the tested configurations."
else:
vc_grade = "F"
vc_comment = "At least one of the elements claimed by the model is not supported in any tested configuration."
return vc_grade, vc_comment
################################################################################
#
# MAIN PROGRAM
#
###############################################################################
if __name__ == "__main__":
vcargs = {
"vc_name": "vc-species-supported-as-stated",
"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, **vcargs)