#!/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)