#!/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, Regents of the University of Minnesota.
#  All rights reserved.
#
#  Contributor(s):
#     Ellad B. Tadmor
#
################################################################################

# The docstring below is vc_description
"""Determines whether a model has a continuous energy and first
derivative, i.e. belongs to the C^1 continuity class, for all possible dimers.
For a model supporting N species, there are N + N!/(2(N-2)! distinct dimers for
all possible species combinations.  (For example if N=3, there are 3+3!/2=6
dimers.  If the species are A,B,C, the 6 dimers are AA, BB, CC, AB, AC, BC.)
For each dimer, the equilibrium separation and cutoff are determined.
The continuity across the cutoff is assessed.  Then an analysis is performed
to detect any discontinuities from half the equilibrium distance to the cutoff.
Although the verification check only requires C^1 continuity to pass, continuity
up to 3rd order is checked and reported."""

# Python 2-3 compatible code issues
from __future__ import print_function

try:
    input = raw_input
except NameError:
    pass

import itertools
import math

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
import scipy.optimize as opt
import numpy as np
import numdifftools as nd
from numdifftools.step_generators import MaxStepGenerator

__version__ = "004"
__author__ = "Ellad Tadmor"


################################################################################
#
#   FUNCTIONS
#
################################################################################
def create_species_pairs_list(species):
    """Given a list of species, create all the monoatomic and biatomic dimer
    pairs.

    Parameters
    ----------
    species : array_like
        Contains a list of species strings, e.g. ["Si", "C"]

    Returns
    -------
    list of tuples
        Contains tuples which correspond to all binary combinations (include
        self-combinations) of the species given in the ``species`` array.
    """
    species_pairs = []
    for sp in species:
        species_pairs.append((sp, sp))
    if len(species) > 1:
        species_pairs += list(itertools.combinations(species, 2))
    return species_pairs


################################################################################
def energy(a, dimer, large_cell_len, einf):
    """Return the potential energy of a dimer with a separation ``a``,
    subtracting off the total potential energy of the dimer at very large
    separation. The Atoms object argument, ``dimer``, has its positions
    mutated.

    Parameters
    ----------
    a : float
        Desired dimer separation (Angstroms)
    dimer : Atoms
        ASE Atoms object whose atomic positions will be mutated based on the
        values of ``a`` and ``large_cell_len``
    large_cell_len : float
        Side length of the full (non-periodic) simulation box (Angstroms)
    einf : float
        Potential energy computed for the species pair contained in ``dimer``
        when the atoms are not interacting

    Returns
    -------
    float
        Potential energy of the dimer at separation ``a`` in eV, adjusted for
        isolated atom energies
    """
    dimer.set_positions(get_dimer_positions(a, large_cell_len))
    return dimer.get_potential_energy() - einf


################################################################################
def energy_cheat(a, dimer, large_cell_len, offset, einf):
    """Return the potential energy of a dimer with a separation ``a``,
    subtracting off the total potential energy of the dimer at very large
    separation and adding a small value (``offset``) to ensure that the root
    bisection has a positive value past the cutoff. The Atoms object
    argument, ``dimer``, has its positions mutated.

    Parameters
    ----------
    a : float
        Desired dimer separation (Angstroms)
    dimer : Atoms
        ASE Atoms object whose atomic positions will be mutated based on the
        values of ``a`` and ``large_cell_len``
    large_cell_len : float
        Side length of the full (non-periodic) simulation box (Angstroms)
    offset : float
        Small positive value to add onto the energy to ensure that the energy
        is positive when the dimer separation exceeds the model's cutoff
        distance of the model for the constituent species
    einf : float
        Potential energy computed for the species pair contained in ``dimer``
        when the atoms are not interacting

    Returns
    -------
    float
        Potential energy of the dimer at separation ``a`` in eV, adjusted for
        isolated atom energies and with a small positive value (``offset``)
        added on
    """
    dimer.set_positions(get_dimer_positions(a, large_cell_len))
    return (dimer.get_potential_energy() - einf) + offset


################################################################################
def get_dimer_positions(a, large_cell_len):
    """Generate positions for a dimer of length 'a' centered in a finite
    simulation box with side length 'large_cell_len'.

    Parameters
    ----------
    a : float
        Desired dimer separation (Angstroms)
    large_cell_len : float
        Side length of the full (non-periodic) simulation box (Angstroms)

    Returns
    -------
    list of three-dimensional lists of floats
        Positions of each of the two atoms in the dimer in Angstroms
    """
    half_cell = 0.5 * large_cell_len
    positions = [
        [half_cell - 0.5 * a, half_cell, half_cell],
        [half_cell + 0.5 * a, half_cell, half_cell],
    ]
    return positions


################################################################################
def get_equilibrium_separation(dimer, large_cell_len, init_guess, einf):
    """Compute the equilibrium separation of a dimer (2-atom chain) using
    Nelder-Mead simplex minimization.

    Parameters
    ----------
    dimer : Atoms
        ASE Atoms object whose atomic positions will be mutated based on the
        values of ``a`` and ``large_cell_len``
    large_cell_len : float
        Side length of the full (non-periodic) simulation box (Angstroms)
    init_guess : float
        Starting guess at the dimer separation to use in the minimization
        algorithm
    einf : float
        Potential energy computed for the species pair contained in ``dimer``
        when the atoms are not interacting

    Returns
    -------
    tuple
        The dimer separation distance (Angstroms) and potential energy (eV)
        obtained at the end of the minimization, as well as a dict containing
        some information about the minimization (number of iterations,
        function calls, and whether a warning was flagged). Note that the
        energy already has the isolated atom energies of the dimer subtracted
        off.
    """
    aopt_arr, eopt, iterations, funccalls, warnflag = opt.fmin(
        energy,
        init_guess,
        args=(dimer, large_cell_len, einf),
        full_output=True,
        xtol=1e-8,
        disp=False,
        maxiter=1000,
    )
    info = {"iterations": iterations, "func_calls": funccalls, "warnflag": warnflag}
    return aopt_arr[0], eopt, info


################################################################################
def get_cutoff_radius(dimer, a, b, offset, large_cell_len, einf):
    """Compute the cutoff radius (separation distance at which energy of
    dimer is zero).

    Parameters
    ----------
    dimer : Atoms
        ASE Atoms object whose atomic positions will be mutated based on the
        values of ``a`` and ``large_cell_len``
    a : float
        Lower bound of the search range (inclusive) for the dimer cutoff
        separation (Angstroms)
    b : float
        Upper bound of the search range (inclusive) for the dimer cutoff
        separation (Angstroms)
    offset : float
        Small positive value to add onto the energy to ensure that the energy
        is positive when the dimer separation exceeds the model's cutoff
        distance of the model for the constituent species
    large_cell_len : float
        Side length of the full (non-periodic) simulation box (Angstroms)
    einf : float
        Potential energy computed for the species pair contained in ``dimer``
        when the atoms are not interacting

    Returns
    -------
    tuple
        The cutoff distance of the dimer in Angstroms and a RootResults
        object that contains information about the bisection procedure
        (including whether it converged)
    """
    x0, results = opt.bisect(
        energy_cheat,
        a,
        b,
        args=(dimer, large_cell_len, offset, einf),
        full_output=True,
        xtol=1e-8,
        maxiter=1000,
    )
    return x0, results


################################################################################
def do_vc(model, vc):
    """Perform DimerContinuityC1 Verification Check

    Parameters
    ----------
    model : str
        Extended KIM ID of the model on which the verification check is to be
        performed
    vc : VerificationCheck
        Instance of the VerificationCheck class defined in kim_python_utils.vc

    Returns
    -------
    tuple of strings
        The grade assigned to the model by the verification check and a
        comment containing a summary of the result

    Raises
    ------
    RuntimeError
        If the static optimizations to find the equilibrium dimer separation
        fail for every initial separation tried
    RuntimeError
        If the check fails to complete for any species pair supported by the model
    """
    # Initialize
    max_deriv = 3  # largest derivative to be investigated
    # (NOTE: If this is increased, increase nth below)

    continuous = [True] * (1 + max_deriv)  # assume function and all
    # derivs are continuous

    etol_fine = 1e-15  # This tolerance is used to:
    # (1) Determine if a given species pair has energy interactions by checking
    #     if the bond length of a dimer can be reduced so that the total energy
    #     differs from the sum of the corresponding isolated atomic energies for
    #     each species by more than this value.
    # (2) Determine an upper bracket of the search range when attempting to
    #     calculate the cutoff radius.  In order to ensure that a given dimer
    #     separation is above the cutoff, its total energy is compared to the
    #     sum of the corresponding isolated energies; if the absolute value of
    #     the difference is less than this amount, we assert that the atoms are
    #     no longer interacting.

    etol_coarse = 1e-6  # This tolerance is used to:
    # (1) Determine the lower bracket of the search range when attempting to
    #     calculate the cutoff radius by checking if the total energy differs
    #     from the sum of the corresponding isolated atomic energies by more
    #     than this value.
    # (2) Judge whether the energy and its derivatives are truly continuous at
    #     the cutoff.  The energy must be within this value of the sum of the
    #     corresponding isolated atomic energies; the absolute value of the
    #     derivatives must not exceed this value.

    led_tolerance = 1.0  # value considered an indicator of a
    # discontinuity for the local edge detection
    # algorithm

    # Side length of large finite simulation domain
    large_cell_len = 50

    got_atleast_one_set_of_results = False

    nth = {0: "Energy", 1: "1st-Deriv", 2: "2nd-Deriv", 3: "3rd-Deriv"}
    yesno = {True: "yes", False: "** NO **"}

    # Get supported species
    species = get_model_supported_species(model)
    species = kim_ase_utils.remove_species_not_supported_by_ASE(list(species))
    species.sort()

    # Print information specific to this VC
    dashwidth = 101
    vc.rwrite("")
    vc.rwrite("-" * dashwidth)
    vc.rwrite("Results for KIM Model      : %s" % model.strip())
    vc.rwrite("Supported species          : %s" % " ".join(species))

    # Get isolated energy of each species for this model
    isolated_energy_per_atom = {}
    for spec in species:
        isolated_energy_per_atom[spec] = kim_ase_utils.get_isolated_energy_per_atom(
            model, spec
        )

    # Create list of dimers species
    species_pairs = create_species_pairs_list(species)

    # Loop over all dimers and determine the smoothness of each
    for i in range(0, len(species_pairs)):

        atoms_interacting_energy = kim_ase_utils.check_if_atoms_interacting(
            model, symbols=species_pairs[i], check_force=False, etol=etol_fine
        )
        if not atoms_interacting_energy:
            vc.rwrite("")
            vc.rwrite(
                "WARNING: The model provided, {}, does not possess a "
                "non-trivial energy interaction for species {} as required "
                "by this Verification Check. Skipping...".format(model,
                species_pairs[i])
            )
            vc.rwrite("")
            continue

        # Infinite separation energy for this species pair
        einf = (
            isolated_energy_per_atom[species_pairs[i][0]]
            + isolated_energy_per_atom[species_pairs[i][1]]
        )

        # Print header for this dimer
        vc.rwrite("")
        vc.rwrite("-" * dashwidth)
        vc.rwrite("DIMER {0}--{1}".format(*species_pairs[i]))
        vc.rwrite("-" * dashwidth)

        # Build ASE Atoms object for the dimer
        a = 1.0
        dimer = Atoms(
            "".join(species_pairs[i]),
            positions=get_dimer_positions(a, large_cell_len),
            cell=(large_cell_len, large_cell_len, large_cell_len),
            pbc=(False, False, False),
        )

        # Create calculator
        calc = KIM(model)
        dimer.calc = calc

        vc.rwrite("")
        vc.rwrite("E n e r g y  E x t r e m a")

        # Find equilibrium separation for this dimer
        max_equilibrium_separation = 10.0
        got_equil_sep = False
        while not got_equil_sep:
            try:
                aopt, eopt, info = get_equilibrium_separation(
                    dimer, large_cell_len, init_guess=a, einf=einf
                )
                got_equil_sep = True
            except:  # noqa: E722
                a += 0.25
                if a > max_equilibrium_separation:
                    raise RuntimeError(
                        "Cannot find a working configuration within a "
                        "reasonable dimer separation range."
                    )

        warn = "none"
        if info["warnflag"]:
            warn = "** DID NOT CONVERGE **"
        vc.rwrite(
            "{0:25}   {1:>15}   {2:>15}   {3:>9}   {4:>9}    {5}".format(
                "", "distance", "energy", "#iter", "#fn-calls", "warnings"
            )
        )
        vc.rwrite(
            "{0:25}   {1: 11.8e}   {2: 11.8e}   {3:9d}   {4:9d}    {5}".format(
                "Equilibrium separation:",
                aopt,
                eopt,
                info["iterations"],
                info["func_calls"],
                warn,
            )
        )
        if info["warnflag"]:
            vc.rwrite("")
            vc.rwrite(
                "WARNING: NOT CHECKING CONTINUITY FOR "
                + "-".join(species_pairs[i])
                + " DIMER."
            )
            vc.rwrite("")

            # Explicitly close calculator to ensure any allocated memory
            # is freed (relevant for SMs)
            if hasattr(calc, "__del__"):
                calc.__del__()

            continue  # failed to converge, so skip this dimer

        # Check if this dimer collapsed while trying to find the equilibrium distance
        print_aopt_warning = False
        if aopt < np.finfo(float).eps:
            aopt = abs(aopt)  # In case aopt is a very small negative number
            print_aopt_warning = True

        # Find dimer cutoff radius
        #
        # Step 1: Find upper end of bracket delimiting cutoff
        #         Start from equilibrium separation, increase separation until
        #         energy reaches infinite separation energy.
        warn = "none"
        failed = False
        max_upper_cutoff_bracket = 2 * max_equilibrium_separation
        b = aopt
        db = max(aopt, 0.5)
        still_interacting = True
        while still_interacting and not failed:
            b += db
            if b > max_upper_cutoff_bracket:
                failed = True
                warn = "** Unable to obtain upper bracket for cutoff radius **"
            else:
                eb = energy(b, dimer, large_cell_len, einf)
                if abs(eb) < etol_fine:
                    still_interacting = False

        #
        # Step 2: Find lower end of bracket delimiting cutoff
        #         Start from upper end of bracket, and reduce separation until
        #         energy changes by a discernable amount.
        if not failed:
            a = b
            da = 0.01
            not_interacting = True
            while not_interacting:
                a -= da
                if a < 0:
                    failed = True
                    warn = "** Unable to obtain lower bracket for cutoff radius **"
                else:
                    ea = energy(a, dimer, large_cell_len, einf)
                    if abs(ea) > etol_coarse:
                        not_interacting = False

        #
        # Step 3: Use bisection algorithm to get cutoff radius
        #
        #         NOTE: Some Simulator Models have a history dependence due to
        #         them maintaining charges from the previous energy evaluation
        #         to use as an initial guess for the next charge equilibration.
        #         We therefore have to treat them not as single-valued
        #         functions but as distributions, i.e. for a given
        #         configuration you might get any of a range of energy values
        #         depending on the history of your previous energy evaluations.
        #         This is particularly problematic for this step, where we set
        #         up a bisection problem in order to determine the cutoff
        #         radius of the model.  Our solution for this specific case is
        #         to make a very crude estimate of the variance of that
        #         distribution with a 10% factor of safety on it.
        ran_bisection = False
        eb_new = energy(b, dimer, large_cell_len, einf)
        eb_error = abs(eb_new - eb)
        if not failed:
            # compute offset to ensure that energy before and after cutoff have
            # different signs
            if ea < eb:
                offset = -eb + 1.1 * eb_error + np.finfo(float).eps
            else:
                offset = -eb - 1.1 * eb_error - np.finfo(float).eps

            rcut, results = get_cutoff_radius(dimer, a, b, offset, large_cell_len, einf)
            ran_bisection = True
            if not results.converged:
                failed = True
                warn = "** DID NOT CONVERGE **"

        #
        # Step 4: Report to user
        if ran_bisection:
            vc.rwrite(
                "{0:25}   {1: 11.8e}   {2: 11.8e}   {3:9d}   {4:9d}    {5}".format(
                    "Cutoff separation:",
                    rcut,
                    energy(rcut, dimer, large_cell_len, einf),
                    results.iterations,
                    results.function_calls,
                    warn,
                )
            )
        else:
            vc.rwrite(warn)

        if print_aopt_warning:
            vc.rwrite("")
            vc.rwrite(
                "WARNING: Equilibrium distance ({} Angstroms) is smaller than "
                "machine epsilon ({}).".format(aopt, np.finfo(float).eps)
            )

        if failed:
            vc.rwrite("")
            vc.rwrite(
                "WARNING: NOT CHECKING CONTINUITY FOR "
                + "-".join(species_pairs[i])
                + " DIMER."
            )
            vc.rwrite("")

            # Explicitly close calculator to ensure any allocated memory
            # (particularly for Simulator Models) is freed
            if hasattr(calc, "__del__"):
                calc.__del__()
            continue  # failed to converge, so skip this dimer

        vc.rwrite("")
        vc.rwrite("C u t o f f  S m o o t h n e s s")
        vc.rwrite("")

        vc.rwrite("{0:10}   {1:>15}   {2:10}".format("", "value", "continuous"))

        # Check smoothness at cutoff
        sg = MaxStepGenerator(
            base_step=0.5 * (rcut - aopt),
            num_steps=14,
            use_exact_steps=True,
            step_ratio=1.6,
            offset=0,
        )
        for n in range(0, max_deriv + 1):
            Denergy = nd.Derivative(
                energy, step=sg, full_output=True, method="backward", n=n
            )
            val, info = Denergy(rcut, dimer, large_cell_len, einf)
            is_continuous = abs(val) <= etol_coarse
            if not is_continuous:
                continuous[n] = False
            vc.rwrite(
                "{0:10}   {1: 11.8e}   {2:>10}".format(
                    nth[n], np.asscalar(val), yesno[is_continuous]
                )
            )

        vc.rwrite("")
        vc.rwrite("C o n t i n u i t y")
        vc.rwrite("")

        # Set the range and increment for exploring internal discontinuities
        if aopt < rcut / 1000.0:
            # Start with amin=rcut/5 and keep increasing until the energy can be
            # evaluated until amin=rcut/2
            afrac = 0.6
            max_lower_limit = rcut / 2

        else:
            # Start with amin=aopt/2 and keep increasing until the energy can be
            # evaluated until amin=aopt
            afrac = 0.5
            max_lower_limit = aopt

        got_amin = False
        while not got_amin:
            try:
                amin = (1.0 - afrac) * max_lower_limit
                # See that energy can be computed without error.
                energy(amin, dimer, large_cell_len, einf)
                got_amin = True
            except:  # noqa: E722
                afrac *= 0.5

        amax = rcut
        del_a = 0.01
        na = int(math.ceil((amax - amin) / del_a))
        vc.rwrite(
            "Checking continuity for r = [{0:.5f},{1:.5f}] at {2:d} points "
            "(Delta r = {3:.5f})".format( amin, amax, na, del_a)
        )

        vc.rwrite("")
        vc.rwrite(
            "Local edge detection based on a normalized 5th-order local "
            "difference formula T^5 is used to determine the presence of "
            "discontinuities. The tolerance is |T^5|>{0:.5f}.".format(
                led_tolerance)
        )
        vc.rwrite(
            "(For details see Anne Gelb and Eitan Tadmor, J. Sci. Comp., "
            "28:279-306, 2006.)"
        )
        vc.rwrite("")

        for n in range(0, max_deriv + 1):
            # set numerical derivative to constant step to prevent the algorithm
            # going haywire in some cases.
            Denergy = nd.Derivative(energy, full_output=True, n=n, step=0.1 * del_a)
            if n == 0:
                aux_file = "dimer-energy-" + "".join(species_pairs[i]) + ".dat"
            else:
                aux_file = (
                    "dimer-energy-deriv-"
                    + str(n)
                    + "-"
                    + "".join(species_pairs[i])
                    + ".dat"
                )
            vc.vc_files.append(aux_file)
            vc.rwrite("")
            vc.rwrite(
                'Checking {0} ({1} vs. distance in file "{2}")'.format(
                    nth[n], nth[n].lower(), aux_file
                )
            )

            # Generate energy curve
            r = []
            e = []
            for j in range(0, na + 1):
                a = amin + j * del_a
                val, info = Denergy(a, dimer, large_cell_len, einf)
                r.append(a)
                e.append(np.asscalar(val))
            vc.write_aux_x_y(aux_file, r, e)

            # Apply local edge detection algorithm to identify discontinuities
            # Based on: A. Gelb and E. Tadmor, J. Sci. Comp., 28:279-306, 2006.
            fact = 1.0 / 6.0
            is_continuous = True
            for j in range(2, na - 3):
                # use 5-th order local difference formula
                led = fact * (
                    -e[j - 2]
                    + 5 * e[j - 1]
                    - 10 * e[j]
                    + 10 * e[j + 1]
                    - 5 * e[j + 2]
                    + e[j + 3]
                )
                if abs(led) > led_tolerance:
                    continuous[n] = False
                    is_continuous = False
                    vc.rwrite(
                        "==> Suspected discontinuity encountered at r={:11.8e} "
                        "(|T^5| = {:11.8e})".format( r[j], abs(led))
                    )
            if is_continuous:
                vc.rwrite("... No discontinuities found.")

        # Explicitly close calculator to ensure any allocated memory
        # (particularly for leaky Simulator Models) is freed
        if hasattr(calc, "__del__"):
            calc.__del__()

        got_atleast_one_set_of_results = True

    if got_atleast_one_set_of_results:
        # Summary of results
        vc.rwrite("")
        vc.rwrite("=" * dashwidth)
        vc.rwrite("")
        vc.rwrite("SUMMARY of Model Continuity Results Across All Dimers:")
        vc.rwrite("")

        vc.rwrite("{0:10}   {1:10}".format("", "continuous"))
        vc.rwrite("-" * 29)
        for n in range(0, max_deriv + 1):
            vc.rwrite("{0:10}   {1:>10}".format(nth[n], yesno[continuous[n]]))

        # Determine continuity class and write out properties
        k = -1
        n = 0
        while n <= max_deriv and continuous[n]:
            k = n
            n += 1
        vc_comment = "The model is C^{0:d} continuous. ".format(k)
        if k == -1:
            vc_comment += "This means that the model has discontinuous energy."
        if k == 0:
            vc_comment += ("This means that the model has continuous energy, "
                           "but a discontinuous first derivative.")
        if k == 1:
            vc_comment += ("This means that the model has continuous energy "
                           "and continuous first derivative.")
        if k == 2:
            vc_comment += ("This means that the model has continuous energy "
                           "and continuous derivatives up to order 2.")
        if k == 3:
            vc_comment += ("This means that the model has continuous energy "
                           "and continuous derivatives at least up to order "
                           "3. (Derivatives beyond this order were not tested.)")

        vc.rwrite("")
        vc.rwrite("=" * dashwidth)
        vc.rwrite("Continuity must be C^1 or higher to pass this verification "
                  "check."
        )
        vc.rwrite("")
        if k >= 1:
            vc_grade = "P"
        else:
            vc_grade = "F"

        return vc_grade, vc_comment

    else:
        msg = (
            "ERROR: Failed to find any pairs in the species set supported "
            "by the model that have a non-zero energy interaction. This "
            "likely means that the model only computes forces and does "
            "not compute energy; this Verification Check currently "
            "requires that models report energy."
        )
        vc.rwrite("")
        vc.rwrite(msg)
        vc.rwrite("")
        raise RuntimeError(msg)


################################################################################
#
#   MAIN PROGRAM
#
###############################################################################
if __name__ == "__main__":

    vcargs = {
        "vc_name": "vc-dimer-continuity-c1",
        "vc_author": __author__,
        "vc_description": kim_vc_utils.vc_stripall(__doc__),
        "vc_category": "informational",
        "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)