Interactions

In general, rose supports local, complex interactions. While we expect that most users will take advantage of the hard-coded Koning-Delaroche potential, we have always written rose with the expectation that many users will want to define their own interactions. The classes below make that process more convenient.

The most basic interaction is affectionately and creatively referred to as: Interaction. It supports fixed-energy interactions whose parameter dependece is affine. The corresponding class, InteractionSpace, generates a list of $\ell$-specific Interactions.

For fixed-energy interactions for which the parameter dependence is non-affine, we have InteractionEIM and InteractionEIMSpace. The classes leverage the Empirical Interpolation Method (EIM) to render that which was non-affine affine.

For non-affine, energy-dependent interactions, we have EnergizedInteractionEIM and EnergizedInteractionEIMSpace. rose works with the energy-scaled Schrödinger equation, so one might think that the scaled interaction is linear in $1/E$. However, because we also work in the dimensionless space, $s\equiv kr$, the dependence is more complex. We again rely on EIM to capture these complexities. (You don't need to know all of that. We just wanted you to impress you.)

Affine, Fixed-Energy Interactions

Wraps the user-defined interaction into a class that stores several relevant parameters of the problem.

Interaction(ell=0, spin_orbit_term=None, coordinate_space_potential=None, n_theta=None, mu=None, energy=None, k=None, Z_1=0, Z_2=0, R_C=0.0, is_complex=False)

Defines a local, (possibly) complex, affine, fixed-energy interaction.

Creates a local, (possibly) complex, affine, fixed-energy interaction.

Parameters:

Name	Type	Description	Default
coordinate_space_potential	Callable[[float, ndarray], float]	V(r, theta)	None
n_theta	int	number of parameters	None
mu	float	reduced mass	None
energy	float	center-of-mass, scattering energy	None
ell	int	angular momentum	0
Z_1	int	charge of particle 1	0
Z_2	int	charge of particle 2	0
R_C	float	Coulomb "cutoff" radius	0.0
is_complex	bool	Is the interaction complex?	False
spin_orbit_term	SpinOrbitTerm	See Spin-Orbit section.	None

Returns:

Name	Type	Description
instance	Interaction	instance of Interaction

Attributes:

Name	Type	Description
v_r	Callable[[float, ndarray], float]	coordinate-space potential; $V(r, \alpha)$
n_theta	int	number of parameters
mu	float	reduced mass
ell	int	angular momentum
k_c	float	Coulomb momentum; $k\eta$
is_complex	bool	Is this a complex potential?
spin_orbit_term	SpinOrbitTerm	See Spin-Orbit section

Source code in src/rose/interaction.py

def __init__(
    self,
    ell: int = 0,
    spin_orbit_term: SpinOrbitTerm = None,
    coordinate_space_potential: Callable[[float, np.array], float] = None,
    n_theta: int = None,
    mu: float = None,
    energy: float = None,
    k: float = None,
    Z_1: int = 0,
    Z_2: int = 0,
    R_C: float = 0.0,
    is_complex: bool = False,
):
    r"""Creates a local, (possibly) complex, affine, fixed-energy interaction.

    Parameters:
        coordinate_space_potential (Callable[[float,ndarray],float]): V(r, theta)
        n_theta (int): number of parameters
        mu (float): reduced mass
        energy (float): center-of-mass, scattering energy
        ell (int): angular momentum
        Z_1 (int): charge of particle 1
        Z_2 (int): charge of particle 2
        R_C (float): Coulomb "cutoff" radius
        is_complex (bool): Is the interaction complex?
        spin_orbit_term (SpinOrbitTerm): See [Spin-Orbit section](#spin-orbit).

    Returns:
        instance (Interaction): instance of `Interaction`

    Attributes:
        v_r (Callable[[float,ndarray],float]): coordinate-space potential; $V(r, \alpha)$
        n_theta (int): number of parameters
        mu (float): reduced mass
        ell (int): angular momentum
        k_c (float): Coulomb momentum; $k\eta$
        is_complex (bool): Is this a complex potential?
        spin_orbit_term (SpinOrbitTerm): See [Spin-Orbit section](#spin-orbit)

    """
    assert coordinate_space_potential is not None
    assert n_theta > 0

    assert ell >= 0
    self.Z_1 = Z_1
    self.Z_2 = Z_2
    self.v_r = coordinate_space_potential
    self.n_theta = n_theta
    self.ell = ell
    self.is_complex = is_complex
    self.spin_orbit_term = spin_orbit_term

    if spin_orbit_term is None:
        self.spin_orbit_term = SpinOrbitTerm()
        self.include_spin_orbit = False
    else:
        self.include_spin_orbit = True

    self.k = k
    self.mu = mu
    self.energy = energy
    self.R_C = R_C
    if R_C <= 0.0:
        R_C = 1e-9
    self.sommerfeld = 0.0

    if mu is not None:
        self.k_c = ALPHA * Z_1 * Z_2 * self.mu / HBARC
    else:
        self.k_c = 0  # TODO mu/energy emulation does not support Coulomb
        assert self.Z_1 * self.Z_1 == 0

    if energy is not None:
        # If the energy is specified (not None as it is when subclass
        # EnergizedInteraction instantiates), set up associated attributes.
        if mu is not None and k is None:
            self.k = np.sqrt(2 * self.mu * self.energy) / HBARC
        self.sommerfeld = self.k_c / self.k

E(alpha)

Energy. Implemented as a function to support energy emulation (where the energy could be a part of the parameter vector, alpha).

Parameters:

Name	Type	Description	Default
alpha	ndarray	parameter vector	required

Returns:

Name	Type	Description
Energy	float	in [MeV]

Source code in src/rose/interaction.py

def E(self, alpha: np.array):
    r"""Energy. Implemented as a function to support energy
    emulation (where the energy could be a part of the parameter vector,
    `alpha`).

    Parameters:
        alpha (ndarray): parameter vector

    Returns:
        Energy (float): in [MeV]
    """
    return self.energy

basis_functions(rho_mesh)

In general, we approximate the potential as

$\hat{U} = \sum_{j} \beta_j(\alpha) u_j$

For affine interactions (like those defined in this class) the basis functions (or "pillars), $u_j$, are just the "naked" parts of the potential. As seen below, it is assumed that the $\beta_j(\alpha)$ coefficients are just the affine parameters, $\alpha$, themselves.

Parameters:

Name	Type	Description	Default
rho_mesh	ndarray	discrete $\rho$ values at which the potential is going to be evaluated	required

Returns:

Name	Type	Description
value	ndarray	values of the scaled potential at provided $\rho$ points

Source code in src/rose/interaction.py

def basis_functions(self, rho_mesh: np.array):
    r"""In general, we approximate the potential as

    $\hat{U} = \sum_{j} \beta_j(\alpha) u_j$

    For affine interactions (like those defined in this class) the basis
    functions (or "pillars), $u_j$, are just the "naked" parts of the
    potential. As seen below, it is assumed that the $\beta_j(\alpha)$
    coefficients are just the affine parameters, $\alpha$, themselves.

    Parameters:
        rho_mesh (ndarray): discrete $\rho$ values at which the potential is
            going to be evaluated

    Returns:
        value (ndarray): values of the scaled potential at provided $\rho$ points
    """
    return np.array([self.tilde(rho_mesh, row) for row in np.eye(self.n_theta)]).T

bundle_gcoeff_args(alpha)

Bundles parameters for the Schrödinger equation

Returns:

Type	Description
	args (tuple) : all the arguments to g_coeff except for $s$

Parameters:

Name	Type	Description	Default
alpha	ndarray)	the parameters for the interaction	required

Source code in src/rose/interaction.py

def bundle_gcoeff_args(self, alpha: np.array):
    r"""Bundles parameters for the Schrödinger equation

    Returns:
        args (tuple) : all the arguments to g_coeff except for $s$

    Parameters:
        alpha (ndarray) : the parameters for the interaction
    """
    k = self.momentum(alpha)
    S_C = self.coulomb_cutoff(alpha) * k
    E = self.E(alpha)
    eta = self.eta(alpha)
    l = self.ell
    v_r = self.v_r
    l_dot_s = self.spin_orbit_term.l_dot_s
    v_so = self.spin_orbit_term.v_so

    return (alpha, k, S_C, E, eta, l, v_r, v_so, l_dot_s)

coefficients(alpha)

As noted in the basis_functions documentation, the coefficients for affine interactions are simply the parameter values. The inverse of the momentum is also returned to support energy emulation.

Parameters:

Name	Type	Description	Default
alpha	ndarray	parameter point	required

Returns:

Name	Type	Description
result	tuple	inverse momentum and coefficients

Source code in src/rose/interaction.py

def coefficients(self, alpha: np.array):  # interaction parameters
    r"""As noted in the `basis_functions` documentation, the coefficients
    for affine interactions are simply the parameter values. The inverse of
    the momentum is also returned to support energy emulation.

    Parameters:
        alpha (ndarray): parameter point

    Returns:
        result (tuple): inverse momentum and coefficients

    """
    return 1 / self.k, alpha

coulomb_cutoff(alpha)

Coulomb cutoff. Implemented as a function to support energy emulation (where the energy/momentum could be a part of the parameter vector, alpha).

Parameters:

Name	Type	Description	Default
alpha	ndarray	parameter vector	required

Returns:

Name	Type	Description
R_C	float	Coulomb cutoff

Source code in src/rose/interaction.py

def coulomb_cutoff(self, alpha: np.array):
    r"""Coulomb cutoff. Implemented as a function to support energy emulation
    (where the energy/momentum could be a part of the parameter vector,
    `alpha`).

    Parameters:
        alpha (ndarray): parameter vector

    Returns:
        R_C (float): Coulomb cutoff
    """
    return self.R_C

eta(alpha)

Sommerfeld parameter. Implemented as a function to support energy emulation (where the energy could be a part of the parameter vector, alpha).

Parameters:

Name	Type	Description	Default
alpha	ndarray	parameter vector	required

Returns:

Name	Type	Description
eta	float	Sommerfeld parameter

Source code in src/rose/interaction.py

def eta(self, alpha: np.array):
    r"""Sommerfeld parameter. Implemented as a function to support energy
    emulation (where the energy could be a part of the parameter vector,
    `alpha`).

    Parameters:
        alpha (ndarray): parameter vector

    Returns:
        eta (float): Sommerfeld parameter
    """
    return self.sommerfeld

momentum(alpha)

Momentum. Implemented as a function to support energy emulation (where the energy/momentum could be a part of the parameter vector, alpha).

Parameters:

Name	Type	Description	Default
alpha	ndarray	parameter vector	required

Returns:

Name	Type	Description
k	float	center-of-mass, scattering momentum

Source code in src/rose/interaction.py

def momentum(self, alpha: np.array):
    r"""Momentum. Implemented as a function to support energy emulation
    (where the energy/momentum could be a part of the parameter vector,
    `alpha`).

    Parameters:
        alpha (ndarray): parameter vector

    Returns:
        k (float): center-of-mass, scattering momentum
    """
    return self.k

reduced_mass(alpha)

Mu. Implemented as a function to support energy emulation (where mu could be a part of the parameter vector, alpha).

Parameters:

Name	Type	Description	Default
alpha	ndarray	parameter vector	required

Returns:

Name	Type	Description
Mu	float	in [MeV/c^2]

Source code in src/rose/interaction.py

def reduced_mass(self, alpha: np.array):
    r"""Mu. Implemented as a function to support energy emulation (where mu could be a
    part of the parameter vector, `alpha`).

    Parameters:
        alpha (ndarray): parameter vector

    Returns:
        Mu (float): in [MeV/c^2]
    """
    return self.mu

tilde(s, alpha)

Scaled potential, $\tilde{U}(s, \alpha, E)$.

• Does not include the Coulomb term.
• $E = E_{c.m.}$; [E] = MeV = [v_r]

Parameters:

Name	Type	Description	Default
s	float	mesh point; $s = pr/\hbar$	required
alpha	ndarray	the varied parameters	required

Returns:

Name	Type	Description
u_tilde	float | complex	value of scaled interaction

Source code in src/rose/interaction.py

def tilde(self, s: float, alpha: np.array):
    r"""Scaled potential, $\tilde{U}(s, \alpha, E)$.

    * Does not include the Coulomb term.
    * $E = E_{c.m.}$; `[E] = MeV = [v_r]`

    Parameters:
        s (float): mesh point; $s = pr/\hbar$
        alpha (ndarray): the varied parameters

    Returns:
        u_tilde (float | complex): value of scaled interaction

    """
    vr = self.v_r(s / self.k, alpha) + self.spin_orbit_term.spin_orbit_potential(
        s / self.k, alpha
    )
    return 1.0 / self.energy * vr

InteractionSpace(l_max=15, interaction_type=Interaction, **kwargs)

Generates a list of $\ell$-specific interactions.

Parameters:

Name	Type	Description	Default
l_max	int	maximum angular momentum	15
interaction_type	Type	type of Interaction to construct	Interaction
kwargs	dict	arguments to constructor of interaction_type	{}

Returns:

Name	Type	Description
instance	InteractionSpace	instance of InteractionSpace

Attributes:

Name	Type	Description
interaction	list	list of Interactions
l_max	int	partial wave cutoff
type	Type	interaction type

Source code in src/rose/interaction.py

def __init__(
    self,
    l_max: int = 15,
    interaction_type=Interaction,
    **kwargs,
):
    r"""Generates a list of $\ell$-specific interactions.

    Parameters:
        l_max (int): maximum angular momentum
        interaction_type (Type): type of `Interaction` to construct
        kwargs (dict): arguments to constructor of `interaction_type`

    Returns:
        instance (InteractionSpace): instance of InteractionSpace

    Attributes:
        interaction (list): list of `Interaction`s
        l_max (int): partial wave cutoff
        type (Type): interaction type
    """
    self.l_max = l_max
    self.type = interaction_type
    self.interactions = []

    if "spin_orbit_term" not in kwargs:
        for l in range(self.l_max + 1):
            self.interactions.append([self.type(ell=l, **kwargs)])
    else:
        spin_orbit_potential = kwargs["spin_orbit_term"]
        kwargs.pop("spin_orbit_term")

        for l in range(self.l_max + 1):
            self.interactions.append(
                [
                    self.type(
                        ell=l,
                        spin_orbit_term=SpinOrbitTerm(spin_orbit_potential, lds),
                        **kwargs,
                    )
                    for lds in couplings(l)
                ]
            )

couplings(l)

For a spin-1/2 nucleon scattering off a spin-0 nucleus, there are maximally 2 different total angular momentum couplings: l+1/2 and l-1/2.

Parameters:

Name	Type	Description	Default
l	int	angular momentum	required

Returns:

Name	Type	Description
couplings	list	epectation value of l dot s

Source code in src/rose/interaction.py

def couplings(l):
    r"""For a spin-1/2 nucleon scattering off a spin-0 nucleus, there are
    maximally 2 different total angular momentum couplings: l+1/2 and l-1/2.

    Parameters:
        l (int): angular momentum

    Returns:
        couplings (list): epectation value of l dot s
    """
    js = [l + 1.0 / 2] if l == 0 else [l + 1.0 / 2, l - 1.0 / 2]
    return [(j * (j + 1) - l * (l + 1) - 0.5 * (0.5 + 1)) for j in js]

Non-Affine, Fixed-Energy Interactions

Interactions that leverage the Empirical Interpolation Method (EIM) to allow the emulation of parameters in which the coordinate-space potential is not affine.

InteractionEIM(training_info=None, n_basis=None, expl_var_ratio_cutoff=None, explicit_training=False, n_train=1000, rho_mesh=DEFAULT_RHO_MESH, match_points=None, method='collocation', **kwargs)

Bases: Interaction

Parameters:

Name	Type	Description	Default
training_info	ndarray	Either (1) parameters bounds or (2) explicit training points If (1): This is a 2-column matrix. The first column are the lower bounds. The second are the upper bounds. Each row maps to a single parameter. If (2): This is an MxN matrix. N is the number of parameters. M is the number of samples.	None
n_basis	int	min number of states in the expansion	None
expl_var_ratio_cutoff	float)	the cutoff in sv2/sum(sv2), sv being the singular values, at which the number of kept bases is chosen	None
explicit_training	bool	Is training_info (1) or (2)? (1) is default	False
n_train	int	How many snapshots to generate? Ignored if explicit_training is True.	1000
rho_mesh	ndarray	coordinate-space points at which the interaction is generated (used for training)	DEFAULT_RHO_MESH
match_points	ndarray	$\rho$ points where agreement with the true potential is enforced	None
method	str)	'collocation' or 'least-squares'. If 'collocation', match_points must be the same length as n_basis; otherwise match_points can be any size.	'collocation'
kwargs	dict	kwargs to Interaction.__init__	{}

Attributes:

Name	Type	Description
s_mesh	ndarray	$s$ points
singular_values	ndarray	S in U, S, Vt = numpy.linalg.svd(...)
snapshots	ndarray	pillars, columns of U
match_indices	ndarray	indices of points in $\rho$ mesh that are matched to the true potential
match_points	ndarray	points in $\rho$ mesh that are matched to the true potential
Ainv	ndarray	inverse of A matrix (Ax = b)

Source code in src/rose/interaction_eim.py

def __init__(
    self,
    training_info: np.array = None,
    n_basis: int = None,
    expl_var_ratio_cutoff: float = None,
    explicit_training: bool = False,
    n_train: int = 1000,
    rho_mesh: np.array = DEFAULT_RHO_MESH,
    match_points: np.array = None,
    method="collocation",
    **kwargs,
):
    r"""
    Parameters:
        training_info (ndarray): Either (1) parameters bounds or (2)
            explicit training points

            If (1):
                This is a 2-column matrix. The first column are the lower
                bounds. The second are the upper bounds. Each row maps to a
                single parameter.

            If (2):
                This is an MxN matrix. N is the number of parameters. M is
                the number of samples.
        n_basis (int): min number of states in the expansion
        expl_var_ratio_cutoff (float) : the cutoff in sv**2/sum(sv**2), sv
            being the singular values, at which the number of kept bases is chosen
        explicit_training (bool): Is training_info (1) or (2)? (1) is
            default
        n_train (int): How many snapshots to generate? Ignored if
            explicit_training is True.
        rho_mesh (ndarray): coordinate-space points at which the interaction
            is generated (used for training)
        match_points (ndarray): $\rho$ points where agreement with the true
            potential is enforced
        method (str) : 'collocation' or 'least-squares'. If 'collocation',
            match_points must be the same length as n_basis; otherwise match_points
            can be any size.
        kwargs (dict): kwargs to `Interaction.__init__`

    Attributes:
        s_mesh (ndarray): $s$ points
        singular_values (ndarray): `S` in `U, S, Vt = numpy.linalg.svd(...)`
        snapshots (ndarray): pillars, columns of `U`
        match_indices (ndarray): indices of points in $\rho$ mesh that are
            matched to the true potential
        match_points (ndarray): points in $\rho$ mesh that are matched to
            the true potential
        Ainv (ndarray): inverse of A matrix (Ax = b)
    """
    assert training_info is not None

    if n_basis is None:
        if "n_theta" in kwargs:
            n_basis = kwargs["n_theta"]
        else:
            n_basis = 8

    super().__init__(**kwargs)

    self.method = method
    self.n_train = n_train
    self.n_basis = n_basis
    self.training_info = training_info
    self.s_mesh = rho_mesh

    # Generate a basis used to approximate the potential.
    # Did the user specify the training points?
    if explicit_training:
        snapshots = np.array(
            [self.tilde(rho_mesh, theta) for theta in training_info]
        ).T
    else:
        train = latin_hypercube_sample(n_train, training_info)
        snapshots = np.array([self.tilde(rho_mesh, theta) for theta in train]).T

    U, S, _ = np.linalg.svd(snapshots, full_matrices=False)

    # keeping at min n_basis PC's, find cutoff
    # avoids singular matrix in MAXVOL when we find a region of param
    # space w/ very similar potentials
    self.singular_values = S
    if expl_var_ratio_cutoff is not None:
        expl_var = self.singular_values**2 / np.sum(self.singular_values**2)
        n_basis_svs = np.sum(expl_var > expl_var_ratio_cutoff)
        self.n_basis = max(n_basis_svs, self.n_basis)
    else:
        self.n_basis = n_basis

    self.snapshots = U[:, : self.n_basis]
    self.match_points = match_points

    if match_points is not None and method == "collocation":
        self.n_basis = match_points.size
        self.match_points = match_points
        self.match_indices = np.array(
            [np.argmin(np.abs(rho_mesh - ri)) for ri in self.match_points]
        )
        self.r_i = rho_mesh[self.match_indices]
        self.Ainv = np.linalg.inv(self.snapshots[self.match_indices])
    elif match_points is None and method == "collocation":
        # random r points between 0 and 2π fm
        i_max = self.snapshots.shape[0] // 4
        di = i_max // (self.n_basis - 1)
        i_init = np.arange(0, i_max + 1, di)
        self.match_indices = max_vol(self.snapshots, i_init)
        self.match_points = rho_mesh[self.match_indices]
        self.r_i = self.match_points
        self.Ainv = np.linalg.inv(self.snapshots[self.match_indices])
    elif method == "least-squares":
        if match_points is None:
            self.match_points = rho_mesh
        else:
            self.match_points = match_points
        self.match_indices = np.array(
            [find_nearest_idx(self.s_mesh, x) for x in self.match_points]
        )
        self.match_points = self.s_mesh[self.match_indices]
        self.r_i = self.match_points
        self.Ainv = np.linalg.pinv(self.snapshots[self.match_indices])
    else:
        raise ValueError(
            "argument 'method' should be one of `collocation` or `least-squares`"
        )

basis_functions(s_mesh)

$u_j$ in $\tilde{U} \approx \hat{U} \equiv \sum_j \beta_j(\alpha) u_j$

Parameters:

Name	Type	Description	Default
s_mesh	ndarray	$s$ mesh points	required

Returns:

Name	Type	Description
u_j	ndarray	"pillars" (MxN matrix; M = number of mesh points; N = number of pillars)

Source code in src/rose/interaction_eim.py

def basis_functions(self, s_mesh: np.array):
    r"""$u_j$ in $\tilde{U} \approx \hat{U} \equiv \sum_j \beta_j(\alpha) u_j$

    Parameters:
        s_mesh (ndarray): $s$ mesh points

    Returns:
        u_j (ndarray): "pillars" (MxN matrix; M = number of mesh points; N = number of pillars)

    """
    return np.copy(self.snapshots)

coefficients(alpha)

Computes the EIM expansion coefficients.

Parameters:

Name	Type	Description	Default
alpha	ndarray	interaction parameters	required

Returns:

Name	Type	Description
coefficients	ndarray	EIM expansion coefficients

Source code in src/rose/interaction_eim.py

def coefficients(self, alpha: np.array):
    r"""Computes the EIM expansion coefficients.

    Parameters:
        alpha (ndarray): interaction parameters

    Returns:
        coefficients (ndarray): EIM expansion coefficients

    """
    u_true = self.tilde(self.r_i, alpha)
    return 1 / self.k, self.Ainv @ u_true

percent_explained_variance(n=None)

Returns:

Type	Description
	(float) : percent of variance explained in the training set by the first n_basis principal
	components

Source code in src/rose/interaction_eim.py

def percent_explained_variance(self, n=None):
    r"""
    Returns:
        (float) : percent of variance explained in the training set by the first n_basis principal
        components
    """
    if n is None:
        n = self.n_basis
    return (
        100
        * np.sum(self.singular_values[:n] ** 2)
        / np.sum(self.singular_values**2)
    )

tilde_emu(alpha)

Emulated interaction = $\hat{U}(s, \alpha, E)$

Parameters:

Name	Type	Description	Default
alpha	ndarray	interaction parameters	required

Returns:

Name	Type	Description
u_hat	ndarray	emulated interaction

Source code in src/rose/interaction_eim.py

def tilde_emu(self, alpha: np.array):
    r"""Emulated interaction = $\hat{U}(s, \alpha, E)$

    Parameters:
        alpha (ndarray): interaction parameters

    Returns:
        u_hat (ndarray): emulated interaction

    """
    _, x = self.coefficients(alpha)
    emu = np.sum(x * self.snapshots, axis=1)
    return emu

InteractionEIMSpace(l_max=15, interaction_type=InteractionEIM, **kwargs)

Bases: InteractionSpace

Generates a list of $\ell$-specific, EIMed interactions.

Parameters:

Name	Type	Description	Default
interaction_args	list	positional arguments for constructor of interaction_type	required
interaction_kwargs	dict	arguments to constructor of interaction_type	required
l_max	int	maximum angular momentum	15
interaction_type	Type	type of Interaction to construct	InteractionEIM

Returns:

Name	Type	Description
instance	InteractionEIMSpace	instance of InteractionEIMSpace

Attributes:

Name	Type	Description
interaction	list	list of Interactions
l_max	int	partial wave cutoff
type	Type	interaction type

Source code in src/rose/interaction_eim.py

def __init__(
    self,
    l_max: int = 15,
    interaction_type=InteractionEIM,
    **kwargs,
):
    r"""Generates a list of $\ell$-specific, EIMed interactions.

    Parameters:
        interaction_args (list): positional arguments for constructor of `interaction_type`
        interaction_kwargs (dict): arguments to constructor of `interaction_type`
        l_max (int): maximum angular momentum
        interaction_type (Type): type of `Interaction` to construct

    Returns:
        instance (InteractionEIMSpace): instance of InteractionEIMSpace

    Attributes:
        interaction (list): list of `Interaction`s
        l_max (int): partial wave cutoff
        type (Type): interaction type
    """
    super().__init__(interaction_type=interaction_type, l_max=l_max, **kwargs)

Non-Affine, Energy-Emulated Interactions

Defines a class for "affinizing" Interactions using the Empirical Interpolation Method (EIM).

EnergizedInteractionEIM(**kwargs)

Bases: InteractionEIM

Extension of InteractionEIM that supports energy, mu and k as parameters. Expected format for alpha is [energy, mu, k, *rest_of_params]

Parameters:

Name	Type	Description	Default
kwargs	dict	arguments to InteractionEIM. Note; the energy argument should not	{}

Attributes:

Name	Type	Description
singular_values	ndarray	S in U, S, Vt = numpy.linalg.svd(...)
snapshots	ndarray	pillars, columns of U
match_indices	ndarray	indices of points in $\rho$ mesh that are matched to the true potential
match_points	ndarray	points in $\rho$ mesh that are matched to the true potential
Ainv	ndarray	inverse of A matrix (Ax = b)

Source code in src/rose/energized_interaction_eim.py

def __init__(
    self,
    **kwargs,
):
    r"""
    Parameters:
        kwargs (dict): arguments to InteractionEIM. Note; the `energy` argument should not
        be given, it will be ignored. Energy is the first element of alpha. If `mu is None`,
        the reduced mass will expected to be the second element of alpha.

    Attributes:
        singular_values (ndarray): `S` in `U, S, Vt = numpy.linalg.svd(...)`
        snapshots (ndarray): pillars, columns of `U`
        match_indices (ndarray): indices of points in $\rho$ mesh that are
            matched to the true potential
        match_points (ndarray): points in $\rho$ mesh that are matched to
            the true potential
        Ainv (ndarray): inverse of A matrix (Ax = b)
    """
    super().__init__(**kwargs)

E(alpha)

Energy. Implemented as a function to support energy emulation (where the energy could be a part of the parameter vector, alpha).

Parameters:

Name	Type	Description	Default
alpha	ndarray	parameter vector	required

Returns:

Name	Type	Description
Energy	float	in [MeV]

Source code in src/rose/energized_interaction_eim.py

def E(self, alpha: np.array):
    r"""Energy. Implemented as a function to support energy emulation (where the energy
    could be a part of the parameter vector, `alpha`).

    Parameters:
        alpha (ndarray): parameter vector

    Returns:
        Energy (float): in [MeV]
    """
    return alpha[0]

basis_functions(s_mesh)

$u_j$ in $\tilde{U} \approx \hat{U} \equiv \sum_j \beta_j(\alpha) u_j$

Parameters:

Name	Type	Description	Default
s_mesh	ndarray	$s$ mesh points	required

Returns:

Name	Type	Description
u_j	ndarray	"pillars" (MxN matrix; M = number of mesh points; N = number of pillars)

Source code in src/rose/energized_interaction_eim.py

def basis_functions(self, s_mesh: np.array):
    r"""$u_j$ in $\tilde{U} \approx \hat{U} \equiv \sum_j \beta_j(\alpha) u_j$

    Parameters:
        s_mesh (ndarray): $s$ mesh points

    Returns:
        u_j (ndarray): "pillars" (MxN matrix; M = number of mesh points; N = number of pillars)

    """
    return np.copy(self.snapshots)

bundle_gcoeff_args(alpha)

Bundles parameters for the Schrödinger equation

Returns:

Type	Description
	args (tuple) : all the arguments to g_coeff except for $s$

Parameters:

Name	Type	Description	Default
alpha	ndarray)	the parameters for the interaction	required

Source code in src/rose/energized_interaction_eim.py

def bundle_gcoeff_args(self, alpha: np.array):
    r"""Bundles parameters for the Schrödinger equation

    Returns:
        args (tuple) : all the arguments to g_coeff except for $s$

    Parameters:
        alpha (ndarray) : the parameters for the interaction
    """
    k = self.momentum(alpha)
    S_C = self.coulomb_cutoff(alpha) * k
    E = self.E(alpha)
    eta = self.eta(alpha)
    l = self.ell
    v_r = self.v_r
    l_dot_s = self.spin_orbit_term.l_dot_s

    # remove the energy term for alpha, so we return just the parameters that plug into v_r
    return (alpha[3:], k, S_C, E, eta, l, v_r, self.spin_orbit_term.v_so, l_dot_s)

coefficients(alpha)

Computes the EIM expansion coefficients.

Parameters:

Name	Type	Description	Default
alpha	ndarray	interaction parameters	required

Returns:

Name	Type	Description
coefficients	ndarray	EIM expansion coefficients

Source code in src/rose/energized_interaction_eim.py

def coefficients(self, alpha: np.array):  # interaction parameters
    r"""Computes the EIM expansion coefficients.

    Parameters:
        alpha (ndarray): interaction parameters

    Returns:
        coefficients (ndarray): EIM expansion coefficients

    """
    k = self.momentum(alpha)
    u_true = self.tilde(self.r_i, alpha)
    return 1 / k, self.Ainv @ u_true

eta(alpha)

Returns the Sommerfeld parameter.

Parameters:

Name	Type	Description	Default
alpha	ndarray	interaction parameters	required

Returns:

Name	Type	Description
eta	float	Sommerfeld parameter

Source code in src/rose/energized_interaction_eim.py

def eta(self, alpha: np.array):
    r"""Returns the Sommerfeld parameter.

    Parameters:
        alpha (ndarray): interaction parameters

    Returns:
        eta (float): Sommerfeld parameter

    """
    return self.k_c / self.momentum(alpha)

momentum(alpha)

Center-of-mass, scattering momentum. Implemented as a function to support energy emulation (where k could be a part of the parameter vector, alpha).

Parameters:

Name	Type	Description	Default
alpha	ndarray	interaction parameters	required

Returns:

Name	Type	Description
k	float	momentum

Source code in src/rose/energized_interaction_eim.py

def momentum(self, alpha: np.array):
    r"""Center-of-mass, scattering momentum. Implemented as a function to support energy
    emulation (where k could be a part of the parameter vector, `alpha`).

    Parameters:
        alpha (ndarray): interaction parameters

    Returns:
        k (float): momentum
    """
    return alpha[2]

reduced_mass(alpha)

Mu. Implemented as a function to support energy emulation (where mu could be a part of the parameter vector, alpha).

Parameters:

Name	Type	Description	Default
alpha	ndarray	parameter vector	required

Returns:

Name	Type	Description
Mu	float	in [MeV/c^2]

Source code in src/rose/energized_interaction_eim.py

def reduced_mass(self, alpha: np.array):
    r"""Mu. Implemented as a function to support energy emulation (where mu could be a
    part of the parameter vector, `alpha`).

    Parameters:
        alpha (ndarray): parameter vector

    Returns:
        Mu (float): in [MeV/c^2]
    """
    return alpha[1]

tilde(s, alpha)

Computes the energy-scaled interaction.

Parameters:

Name	Type	Description	Default
s	float	mesh point	required
alpha	ndarray	interaction parameters	required

Returns:

Name	Type	Description
u_tilde	float | complex	energy-scaled interaction

Source code in src/rose/energized_interaction_eim.py

def tilde(self, s: float, alpha: np.array):
    r"""Computes the energy-scaled interaction.

    Parameters:
        s (float): mesh point
        alpha (ndarray): interaction parameters

    Returns:
        u_tilde (float | complex): energy-scaled interaction

    """
    energy = self.E(alpha)
    k = self.momentum(alpha)
    alpha_truncated = alpha[3:]
    vr = self.v_r(
        s / k, alpha_truncated
    ) + self.spin_orbit_term.spin_orbit_potential(s / k, alpha_truncated)
    return 1.0 / energy * vr

tilde_emu(s, alpha)

Emulated interaction = $\hat{U}(s, \alpha, E)$

Parameters:

Name	Type	Description	Default
alpha	ndarray	interaction parameters	required

Returns:

Name	Type	Description
u_hat	ndarray	emulated interaction

Source code in src/rose/energized_interaction_eim.py

def tilde_emu(self, s: float, alpha: np.array):
    r"""Emulated interaction = $\hat{U}(s, \alpha, E)$

    Parameters:
        alpha (ndarray): interaction parameters

    Returns:
        u_hat (ndarray): emulated interaction

    """
    _, x = self.coefficients(alpha)
    emu = np.sum(x * self.snapshots, axis=1)
    return emu

EnergizedInteractionEIMSpace(l_max=15, interaction_type=EnergizedInteractionEIM, **kwargs)

Bases: InteractionEIMSpace

Generates a list of $\ell$-specific interactions.

Parameters:

Name	Type	Description	Default
interaction_args	list	positional arguments for constructor of interaction_type	required
interaction_kwargs	dict	arguments to constructor of interaction_type	required
l_max	int	maximum angular momentum	15
interaction_type	Type	type of Interaction to construct	EnergizedInteractionEIM

Returns:

Name	Type	Description
instance	EnergizedInteractionEIMSpaceInteractionSpace	instance of EnergizedInteractionEIMSpace

Attributes:

Name	Type	Description
interaction	list	list of Interactions
l_max	int	partial wave cutoff
type	Type	interaction type

Source code in src/rose/energized_interaction_eim.py

def __init__(
    self,
    l_max: int = 15,
    interaction_type=EnergizedInteractionEIM,
    **kwargs,
):
    r"""Generates a list of $\ell$-specific interactions.

    Parameters:
        interaction_args (list): positional arguments for constructor of `interaction_type`
        interaction_kwargs (dict): arguments to constructor of `interaction_type`
        l_max (int): maximum angular momentum
        interaction_type (Type): type of `Interaction` to construct

    Returns:
        instance (EnergizedInteractionEIMSpaceInteractionSpace):
            instance of EnergizedInteractionEIMSpace

    Attributes:
        interaction (list): list of `Interaction`s
        l_max (int): partial wave cutoff
        type (Type): interaction type
    """
    super().__init__(
        interaction_type=interaction_type,
        l_max=l_max,
        **kwargs,
    )

Spin-Orbit

Defines a class to package the spin-orbit term.

SpinOrbitTerm(spin_orbit_potential=None, l_dot_s=None)

Spin-orbit interaction

Parameters:

Name	Type	Description	Default
spin_orbit_potential	Callable[[float, ndarray, float], float]	coordinate-space, spin-orbit potential	None
l_dot_s	float	$2\ell\cdot s$ matrix elements, $+\ell$ or $-\ell-1$	None

Attributes:

Name	Type	Description
l_dot_s	float	$2\ell\cdot s$ matrix elements, $+\ell$ or $-\ell-1$
spin_orbit_potential	float	(Callable[[float, ndarray, float],float]): coordinate-space, spin-orbit potential

Source code in src/rose/spin_orbit.py

def __init__(
    self,
    spin_orbit_potential: Callable[[float, np.array, float], float] = None,
    l_dot_s: float = None,
):
    r"""Spin-orbit interaction

    Parameters:
        spin_orbit_potential (Callable[[float, ndarray, float],float]):
            coordinate-space, spin-orbit potential
        l_dot_s (float): $2\ell\cdot s$ matrix elements, $+\ell$ or $-\ell-1$

    Attributes:
        l_dot_s (float): $2\ell\cdot s$ matrix elements, $+\ell$ or $-\ell-1$
        spin_orbit_potential: (Callable[[float, ndarray, float],float]):
            coordinate-space, spin-orbit potential

    """
    self.l_dot_s = l_dot_s
    self.v_so = spin_orbit_potential

    if spin_orbit_potential is None:
        self.l_dot_s = 0
        self.v_so = v_so_default

Koning-Delaroche

The Koning-Delaroche potential is a common optical potential for nuclear scattering. It is provided here in simplified form specifically to address this need.

See the Koning-Delaroche paper for details. Equation references are with respect to (w.r.t.) this paper.

EnergizedKoningDelaroche(training_info, l_max=20, n_basis=8, explicit_training=False, n_train=1000, rho_mesh=DEFAULT_RHO_MESH, match_points=None, method='collocation', **kwargs)

Bases: EnergizedInteractionEIMSpace

Wraps the Koning-Delaroche potential into a rose-friendly class. Saves system-specific information. Allows the user to emulate across energies.

• Does not (yet) support Coulomb.

Parameters:

Name	Type	Description	Default
mu	float	reduced mass of the 2-body system	required
ell	int	angular momentum	required
training_info	ndarray	either (1) parameters bounds or (2) explicit training points If (1): This is a 2-column matrix. The first column are the lower bounds. The second are the upper bounds. Each row maps to a single parameter. If (2): This is an MxN matrix. N is the number of parameters. M is the number of samples.	required
n_basis	int	number of basis states to use for EIM	8
explicit_training	bool	True implies training_info case (2); False implies (1)	False
n_train	int	how many training samples to use	1000
rho_mesh	ndarray	$\rho$ (or $s$) grid values	DEFAULT_RHO_MESH
match_points	ndarray	$\rho$ points at which we demand the EIMed potential match the true potential	None

Returns:

Name	Type	Description
instance	EnergizedKoningDelaroche	instance of the class

Source code in src/rose/koning_delaroche.py

def __init__(
    self,
    training_info: np.array,
    l_max=20,
    n_basis: int = 8,
    explicit_training: bool = False,
    n_train: int = 1000,
    rho_mesh: np.array = DEFAULT_RHO_MESH,
    match_points: np.array = None,
    method="collocation",
    **kwargs,
):
    r"""Wraps the Koning-Delaroche potential into a `rose`-friendly class.
    Saves system-specific information. Allows the user to emulate across
    energies.

    * **Does not (yet) support Coulomb.**

    Parameters:
        mu (float): reduced mass of the 2-body system
        ell (int): angular momentum
        training_info (ndarray): either (1) parameters bounds or (2) explicit training points

            If (1):
            This is a 2-column matrix. The first column are the lower
            bounds. The second are the upper bounds. Each row maps to a
            single parameter.

            If (2):
            This is an MxN matrix. N is the number of parameters. M is the
            number of samples.
        n_basis (int): number of basis states to use for EIM
        explicit_training (bool): True implies training_info case (2); False implies (1)
        n_train (int): how many training samples to use
        rho_mesh (ndarray):  $\rho$ (or $s$) grid values
        match_points (ndarray): $\rho$ points at which we demand the EIMed
            potential match the true potential

    Returns:
        instance (EnergizedKoningDelaroche): instance of the class
    """
    n_params = NUM_PARAMS + 3  # include mu, k, and energy

    super().__init__(
        l_max=l_max,
        coordinate_space_potential=KD_simple,
        n_theta=n_params,
        training_info=training_info,
        Z_1=0,
        Z_2=0,
        is_complex=True,
        spin_orbit_term=KD_simple_so,
        n_basis=n_basis,
        explicit_training=explicit_training,
        n_train=n_train,
        rho_mesh=rho_mesh,
        match_points=match_points,
        method=method,
        **kwargs,
    )

KDGlobal(projectile, param_fpath=None)

Global optical potential in Koning-Delaroche form.

Parameters:

Name	Type	Description	Default
projectile		neutron or proton?	required
param_fpath		path to json file encoding parameter values.	None

Source code in src/rose/koning_delaroche.py

def __init__(self, projectile: Projectile, param_fpath: Path = None):
    r"""
    Parameters:
        projectile : neutron or proton?
        param_fpath : path to json file encoding parameter values.
        Defaults to data/KD_default.json
    """
    if param_fpath is None:
        param_fpath = Path(__file__).parent.resolve() / Path(
            "../data/KD_default.json"
        )

    if projectile == Projectile.neutron:
        tag = "_n"
    elif projectile == Projectile.proton:
        tag = "_p"
    else:
        raise RuntimeError(
            "KDGlobal is defined only for neutron and proton projectiles"
        )

    self.projectile = projectile

    self.param_fpath = param_fpath
    with open(self.param_fpath) as f:
        data = json.load(f)

        if "KDHartreeFock" in data:
            # real central depth
            self.v1_0 = data["KDHartreeFock"]["V1_0"]
            self.v1_asymm = data["KDHartreeFock"]["V1_asymm"]
            self.v1_A = data["KDHartreeFock"]["V1_A"]
            self.v2_0 = data["KDHartreeFock"]["V2_0" + tag]
            self.v2_A = data["KDHartreeFock"]["V2_A" + tag]
            self.v3_0 = data["KDHartreeFock"]["V3_0" + tag]
            self.v3_A = data["KDHartreeFock"]["V3_A" + tag]
            self.v4_0 = data["KDHartreeFock"]["V4_0"]

            # real central form
            self.rv_0 = data["KDHartreeFock"]["r_0"]
            self.rv_A = data["KDHartreeFock"]["r_A"]
            self.av_0 = data["KDHartreeFock"]["a_0"]
            self.av_A = data["KDHartreeFock"]["a_A"]

            # imag volume depth
            self.w1_0 = data["KDImagVolume"]["W1_0" + tag]
            self.w1_A = data["KDImagVolume"]["W1_A" + tag]
            self.w2_0 = data["KDImagVolume"]["W2_0"]
            self.w2_A = data["KDImagVolume"]["W2_A"]

            # imag surface depth
            self.d1_0 = data["KDImagSurface"]["D1_0"]
            self.d1_asymm = data["KDImagSurface"]["D1_asymm"]
            self.d2_0 = data["KDImagSurface"]["D2_0"]
            self.d2_A = data["KDImagSurface"]["D2_A"]
            self.d2_A2 = data["KDImagSurface"]["D2_A2"]
            self.d2_A3 = data["KDImagSurface"]["D2_A3"]
            self.d3_0 = data["KDImagSurface"]["D3_0"]

            # imag surface form
            self.rd_0 = data["KDImagSurface"]["r_0"]
            self.rd_A = data["KDImagSurface"]["r_A"]
            self.ad_0 = data["KDImagSurface"]["a_0" + tag]
            self.ad_A = data["KDImagSurface"]["a_A" + tag]

            # real spin orbit depth
            self.Vso1_0 = data["KDRealSpinOrbit"]["V1_0"]
            self.Vso1_A = data["KDRealSpinOrbit"]["V1_A"]
            self.Vso2_0 = data["KDRealSpinOrbit"]["V2_0"]

            # imag spin orbit form
            self.Wso1_0 = data["KDImagSpinOrbit"]["W1_0"]
            self.Wso2_0 = data["KDImagSpinOrbit"]["W2_0"]

            # spin orbit form
            self.rso_0 = data["KDRealSpinOrbit"]["r_0"]
            self.rso_A = data["KDRealSpinOrbit"]["r_A"]
            self.aso_0 = data["KDRealSpinOrbit"]["a_0"]

            # Coulomb
            if self.projectile == Projectile.proton:
                self.rc_0 = data["KDCoulomb"]["r_C_0"]
                self.rc_A = data["KDCoulomb"]["r_C_A"]
                self.rc_A2 = data["KDCoulomb"]["r_C_A2"]

        elif "KDHartreeFock_V1_0" in data:
            # real central depth
            self.v1_0 = data["KDHartreeFock_V1_0"]
            self.v1_asymm = data["KDHartreeFock_V1_asymm"]
            self.v1_A = data["KDHartreeFock_V1_A"]
            self.v2_0 = data["KDHartreeFock_V2_0" + tag]
            self.v2_A = data["KDHartreeFock_V2_A" + tag]
            self.v3_0 = data["KDHartreeFock_V3_0" + tag]
            self.v3_A = data["KDHartreeFock_V3_A" + tag]
            self.v4_0 = data["KDHartreeFock_V4_0"]

            # real central form
            self.rv_0 = data["KDHartreeFock_r_0"]
            self.rv_A = data["KDHartreeFock_r_A"]
            self.av_0 = data["KDHartreeFock_a_0"]
            self.av_A = data["KDHartreeFock_a_A"]

            # imag volume depth
            self.w1_0 = data["KDImagVolume_W1_0" + tag]
            self.w1_A = data["KDImagVolume_W1_A" + tag]
            self.w2_0 = data["KDImagVolume_W2_0"]
            self.w2_A = data["KDImagVolume_W2_A"]

            # imag surface depth
            self.d1_0 = data["KDImagSurface_D1_0"]
            self.d1_asymm = data["KDImagSurface_D1_asymm"]
            self.d2_0 = data["KDImagSurface_D2_0"]
            self.d2_A = data["KDImagSurface_D2_A"]
            self.d2_A2 = data["KDImagSurface_D2_A2"]
            self.d2_A3 = data["KDImagSurface_D2_A3"]
            self.d3_0 = data["KDImagSurface_D3_0"]

            # imag surface form
            self.rd_0 = data["KDImagSurface_r_0"]
            self.rd_A = data["KDImagSurface_r_A"]
            self.ad_0 = data["KDImagSurface_a_0" + tag]
            self.ad_A = data["KDImagSurface_a_A" + tag]

            # real spin orbit depth
            self.Vso1_0 = data["KDRealSpinOrbit_V1_0"]
            self.Vso1_A = data["KDRealSpinOrbit_V1_A"]
            self.Vso2_0 = data["KDRealSpinOrbit_V2_0"]

            # imag spin orbit form
            self.Wso1_0 = data["KDImagSpinOrbit_W1_0"]
            self.Wso2_0 = data["KDImagSpinOrbit_W2_0"]

            # spin orbit form
            self.rso_0 = data["KDRealSpinOrbit_r_0"]
            self.rso_A = data["KDRealSpinOrbit_r_A"]
            self.aso_0 = data["KDRealSpinOrbit_a_0"]

            # Coulomb
            if self.projectile == Projectile.proton:
                self.rc_0 = data["KDCoulomb_r_C_0"]
                self.rc_A = data["KDCoulomb_r_C_A"]
                self.rc_A2 = data["KDCoulomb_r_C_A2"]
        else:
            raise ValueError("Unrecognized parameter file format for KDUQ!")

        # fermi energy
        if self.projectile == Projectile.neutron:
            self.Ef_0 = -11.2814
            self.Ef_A = 0.02646
        else:
            self.Ef_0 = -8.4075
            self.Ef_A = 0.01378

get_params(A, Z, mu, E_lab, k)

Calculates Koning-Delaroche global neutron-nucleus OMP parameters for given A, Z, and COM-frame energy, returns params in form useable by EnergizedKoningDelaroche

Source code in src/rose/koning_delaroche.py

def get_params(self, A, Z, mu, E_lab, k):
    """
    Calculates Koning-Delaroche global neutron-nucleus OMP parameters for given A, Z,
    and COM-frame energy, returns params in form useable by EnergizedKoningDelaroche
    """

    N = A - Z
    delta = (N - Z) / A
    factor = 1.0
    if self.projectile == Projectile.proton:
        delta *= -1.0
        factor = -1.0

    # fermi energy
    Ef = self.Ef_0 + self.Ef_A * A

    # real central depth
    v1 = self.v1_0 - self.v1_asymm * delta - self.v1_A * A
    v2 = self.v2_0 - self.v2_A * A * factor
    v3 = self.v3_0 - self.v3_A * A * factor
    v4 = self.v4_0
    vv = Vv(E_lab, v1, v2, v3, v4, Ef)

    # real central form
    rv = self.rv_0 - self.rv_A * A ** (-1.0 / 3.0)
    av = self.av_0 - self.av_A * A

    # imag volume depth
    w1 = self.w1_0 + self.w1_A * A
    w2 = self.w2_0 + self.w2_A * A
    wv = Wv(E_lab, w1, w2, Ef)

    # imag volume form
    rwv = rv
    awv = av

    # imag surface depth
    d1 = self.d1_0 - self.d1_asymm * delta
    d2 = self.d2_0 + self.d2_A / (1 + np.exp((A - self.d2_A3) / self.d2_A2))
    d3 = self.d3_0
    wd = Wd(E_lab, d1, d2, d3, Ef)

    # imag surface form
    rd = self.rd_0 - self.rd_A * A ** (1.0 / 3.0)
    ad = self.ad_0 - self.ad_A * A * factor

    # real spin orbit depth
    vso1 = self.Vso1_0 + self.Vso1_A * A
    vso2 = self.Vso2_0
    vso = Vso(E_lab, vso1, vso2, Ef)

    # real spin orbit form
    rso = self.rso_0 - self.rso_A * A ** (-1.0 / 3.0)
    aso = self.aso_0

    # imag spin orbit form
    wso1 = self.Wso1_0
    wso2 = self.Wso2_0
    wso = Wso(E_lab, wso1, wso2, Ef)

    # imag spin orbit form
    rwso = rso
    awso = aso

    # Coulomb radius
    R_C = 0
    if self.projectile == Projectile.proton:
        # Coulomb radius
        rc0 = (
            self.rc_0
            + self.rc_A * A ** (-2.0 / 3.0)
            + self.rc_A2 * A ** (-5.0 / 3.0)
        )
        R_C = rc0 * A ** (1.0 / 3.0)

        # Coulomb correction
        Vcbar = 1.73 / rc0 * Z * A ** (-1.0 / 3.0)
        Vc = delta_VC(E_lab, Vcbar, v1, v2, v3, v4, Ef)
        vv += Vc

    # 15 params total
    params = np.array(
        [
            vv,
            rv * A ** (1.0 / 3.0),
            av,
            wv,
            rwv * A ** (1.0 / 3.0),
            awv,
            wd,
            rd * A ** (1.0 / 3.0),
            ad,
            vso,
            rso * A ** (1.0 / 3.0),
            aso,
            wso,
            rwso * A ** (1.0 / 3.0),
            awso,
        ]
    )

    return R_C, params

KoningDelaroche(energy, training_info, mu, l_max=20, n_basis=8, explicit_training=False, n_train=1000, rho_mesh=DEFAULT_RHO_MESH, match_points=None, method='collocation', **kwargs)

Bases: InteractionEIMSpace

Koning-Delaroche potential (without energy-dependent strength coefficients) for arbitrary systems defined by mu, energy, ell, Z_1, and Z_2.

Wraps the Koning-Delaroche potential into a rose-friendly class. Saves system-specific information.

Parameters:

Name	Type	Description	Default
mu	float	reduced mass of the 2-body system	required
ell	int	angular momentum	required
energy	float	center-of-mass, scattering energy	required
training_info	ndarray	either (1) parameters bounds or (2) explicit training points If (1): This is a 2-column matrix. The first column are the lower bounds. The second are the upper bounds. Each row maps to a single parameter. If (2): This is an MxN matrix. N is the number of parameters. M is the number of samples.	required
n_basis	int	number of basis states to use for EIM	8
explicit_training	bool	True implies training_info case (2); False implies (1)	False
n_train	int	how many training samples to use	1000
rho_mesh	ndarray	$\rho$ (or $s$) grid values	DEFAULT_RHO_MESH
match_points	ndarray	$\rho$ points at which we demand the EIMed potential match the true potential	None

Returns:

Name	Type	Description
instance	KoningDelaroche	instance of the class

Source code in src/rose/koning_delaroche.py

def __init__(
    self,
    energy: float,
    training_info: np.array,
    mu: float,
    l_max=20,
    n_basis: int = 8,
    explicit_training: bool = False,
    n_train: int = 1000,
    rho_mesh: np.array = DEFAULT_RHO_MESH,
    match_points: np.array = None,
    method="collocation",
    **kwargs,
):
    r"""Wraps the Koning-Delaroche potential into a `rose`-friendly class.
    Saves system-specific information.

    Parameters:
        mu (float): reduced mass of the 2-body system
        ell (int): angular momentum
        energy (float): center-of-mass, scattering energy
        training_info (ndarray): either (1) parameters bounds or (2) explicit training points

            If (1):
            This is a 2-column matrix. The first column are the lower
            bounds. The second are the upper bounds. Each row maps to a
            single parameter.

            If (2):
            This is an MxN matrix. N is the number of parameters. M is the
            number of samples.
        n_basis (int): number of basis states to use for EIM
        explicit_training (bool): True implies training_info case (2); False implies (1)
        n_train (int): how many training samples to use
        rho_mesh (ndarray):  $\rho$ (or $s$) grid values
        match_points (ndarray): $\rho$ points at which we demand the EIMed
            potential match the true potential

    Returns:
        instance (KoningDelaroche): instance of the class

    """
    super().__init__(
        coordinate_space_potential=KD_simple,
        n_theta=NUM_PARAMS,
        mu=mu,
        energy=energy,
        training_info=training_info,
        Z_1=0,
        Z_2=0,
        is_complex=True,
        spin_orbit_term=KD_simple_so,
        n_basis=n_basis,
        explicit_training=explicit_training,
        n_train=n_train,
        rho_mesh=rho_mesh,
        match_points=match_points,
        method=method,
        l_max=l_max,
        **kwargs,
    )

KD(r, E, v1, v2, v3, v4, w1, w2, d1, d2, d3, Ef, Rv, av, Rd, ad)

Koning-Delaroche without the spin-orbit terms - Eq. (1)

Source code in src/rose/koning_delaroche.py

@njit
def KD(r, E, v1, v2, v3, v4, w1, w2, d1, d2, d3, Ef, Rv, av, Rd, ad):
    """Koning-Delaroche without the spin-orbit terms - Eq. (1)"""
    return (
        -Vv(E, v1, v2, v3, v4, Ef) * woods_saxon_safe(r, Rv, av)
        - 1j * Wv(E, w1, w2, Ef) * woods_saxon_safe(r, Rv, av)
        - 1j * (-4 * ad) * Wd(E, d1, d2, d3, Ef) * woods_saxon_prime_safe(r, Rd, ad)
    )

KD_simple(r, alpha)

simplified Koning-Delaroche without the spin-orbit terms

Take Eq. (1) and remove the energy dependence of the coefficients.

Source code in src/rose/koning_delaroche.py

@njit
def KD_simple(r, alpha):
    r"""simplified Koning-Delaroche without the spin-orbit terms

    Take Eq. (1) and remove the energy dependence of the coefficients.
    """
    vv, rv, av, wv, rwv, awv, wd, rd, ad = decompose_alpha(alpha)[0]
    return (
        -vv * woods_saxon_safe(r, rv, av)
        - 1j * wv * woods_saxon_safe(r, rwv, awv)
        - 1j * (-4 * ad) * wd * woods_saxon_prime_safe(r, rd, ad)
    )

KD_simple_so(r, alpha, lds)

simplified Koning-Delaroche spin-orbit terms

Take Eq. (1) and remove the energy dependence of the coefficients.

lds: l • s = 1/2 * (j(j+1) - l(l+1) - s(s+1))

Source code in src/rose/koning_delaroche.py

@njit
def KD_simple_so(r, alpha, lds):
    r"""simplified Koning-Delaroche spin-orbit terms

    Take Eq. (1) and remove the energy dependence of the coefficients.

    lds: l • s = 1/2 * (j(j+1) - l(l+1) - s(s+1))
    """
    vso, rso, aso, wso, rwso, awso = decompose_alpha(alpha)[1]
    return lds * (
        vso / MASS_PION**2 * thomas_safe(r, rso, aso)
        + 1j * wso / MASS_PION**2 * thomas_safe(r, rwso, awso)
    )

Vso(E, vso1, vso2, Ef)

energy-dependent, spin-orbit strength --- real term, Eq. (7)

Source code in src/rose/koning_delaroche.py

@njit
def Vso(E, vso1, vso2, Ef):
    """energy-dependent, spin-orbit strength --- real term, Eq. (7)"""
    return vso1 * np.exp(-vso2 * (E - Ef))

Vv(E, v1, v2, v3, v4, Ef)

energy-dependent, volume-central strength - real term, Eq. (7)

Source code in src/rose/koning_delaroche.py

@njit
def Vv(E, v1, v2, v3, v4, Ef):
    r"""energy-dependent, volume-central strength - real term, Eq. (7)"""
    return v1 * (1 - v2 * (E - Ef) + v3 * (E - Ef) ** 2 - v4 * (E - Ef) ** 3)

Wd(E, d1, d2, d3, Ef)

energy-dependent, surface-central strength - imaginary term (no real term), Eq. (7)

Source code in src/rose/koning_delaroche.py

@njit
def Wd(E, d1, d2, d3, Ef):
    """energy-dependent, surface-central strength - imaginary term (no real
    term), Eq. (7)
    """
    return d1 * (E - Ef) ** 2 / ((E - Ef) ** 2 + d3**2) * np.exp(-d2 * (E - Ef))

Wso(E, wso1, wso2, Ef)

energy-dependent, spin-orbit strength --- imaginary term, Eq. (7)

Source code in src/rose/koning_delaroche.py

@njit
def Wso(E, wso1, wso2, Ef):
    """energy-dependent, spin-orbit strength --- imaginary term, Eq. (7)"""
    return wso1 * (E - Ef) ** 2 / ((E - Ef) ** 2 + wso2**2)

Wv(E, w1, w2, Ef)

energy-dependent, volume-central strength - imaginary term, Eq. (7)

Source code in src/rose/koning_delaroche.py

@njit
def Wv(E, w1, w2, Ef):
    """energy-dependent, volume-central strength - imaginary term, Eq. (7)"""
    return w1 * (E - Ef) ** 2 / ((E - Ef) ** 2 + w2**2)

decompose_alpha(alpha)

Splits the parameter-space vector into non-spin-orbit and spin-orbit parameters.

Parameters:

Name	Type	Description	Default
alpha	ndarray	interaction parameters	required

Returns:

Name	Type	Description
parameters	tuple	2-tuple of non-spin-orbit (parameters[0]) and spin-orbit parameters (parameters[1])

Source code in src/rose/koning_delaroche.py

@njit
def decompose_alpha(alpha):
    r"""Splits the parameter-space vector into non-spin-orbit and spin-orbit
    parameters.

    Parameters:
        alpha (ndarray): interaction parameters

    Returns:
        parameters (tuple): 2-tuple of non-spin-orbit (`parameters[0]`) and
            spin-orbit parameters (`parameters[1]`)

    """
    vv, rv, av, wv, rwv, awv, wd, rd, ad, vso, rso, aso, wso, rwso, awso = alpha
    return (vv, rv, av, wv, rwv, awv, wd, rd, ad), (vso, rso, aso, wso, rwso, awso)

delta_VC(E, Vcbar, v1, v2, v3, v4, Ef)

energy dependent Coulomb correction term, Eq. 23

Source code in src/rose/koning_delaroche.py

@njit
def delta_VC(E, Vcbar, v1, v2, v3, v4, Ef):
    """energy dependent Coulomb correction term, Eq. 23"""
    return v1 * Vcbar * (v2 - 2 * v3 * (E - Ef) + 3 * v4 * (E - Ef) ** 2)