Interaction Potentials

Interaction potentials are specified in the input file using the interaction_potential field, which is an NXN array(enum), where N is the number of distinct interactions: atomic interaction potential. Kr-C ("KR_C") is a widely used universal screened Coulomb potential, and is recommended. Other options include:

"MOLIERE": Moliere[6] universal screened Coulomb potential
"KR_C": Kr-C[7] universal screened Coulomb potential
"ZBL": ZBL[8] universal screened Coulomb potential
"LENZ_JENSEN": Lenz-Jensen[9] screeened Coulomb potential
{"LENNARD_JONES_12_6"={sigma = 1E-10, epsilon = 1.6E-19}}, Lennard Jones 12-6 potential with sigma and epsilon adjustable parameters
{"LENNARD_JONES_65_6"={sigma = 1E-10, epsilon = 1.6E-19}} Lennard-Jones 6.5-6 potential with sigma and epsilon adjustable parameters
{"MORSE"={D = 1.6E-19, alpha = 1E10, r0 = 1E-10}} Morse potential with D, alpha and r0 adjustable parameters

Multiple interaction potentials

When using multiple interaction potentials, the root-finder, interaction potentials, and scattering integrals must be defined for every combination of ion and material species for which there is a distinct potential. For example, for He on W, the He-He, He-W, W-He, and W-W interactions must be specified. This takes the form of matrices of interaction potentials, scattering integrals, and root-finders for each combination:

[[He-He, He-W],
[W-He, W-W]]

So, the interaction_potential field may be:

[["KR_C", {"LENNARD_JONES_12_6"={sigma = 1E-10, epsilon = 1.6E-19}}],
[{"LENNARD_JONES_12_6"={sigma = 1E-10, epsilon = 1.6E-19}}, "KR_C",]]

And the corresponding scattering integral and rootfinder fields:

[["MENDENHALL_WELLER", {"GAUSS_MEHLER"={n_points=10}}],
[{"GAUSS_MEHLER"={n_points=10}}, "MENDENHALL_WELLER"]]

[[{"NEWTON"={max_iterations = 100, tolerance = 1E-6}}, {"POLYNOMIAL"={complex_threshold=1E-12}}],
[{"POLYNOMIAL"={complex_threshold=1E-12}}, {"NEWTON"={max_iterations = 100, tolerance = 1E-6}}]]

So He-W uses LJ 12-6 potential, the Polynomial rootfinder, and Gauss-Mehler quadrature, while W-W uses Kr-C, Newton-Raphson, and Mendenhall-Weller quadrature.

For most simulations, however, unless one has specific knowledge of the interaction potential, one should use a "universal" screened Coulomb potential for all species (Kr-C is recommended) with the Newton-Raphson rootfinder and Mendenhall-Weller quadrature:

interaction_potential = [["KR_C"]]

root_finder = [[{"NEWTON"={max_iterations=100, tolerance=1E-6}}]]

scattering_integral = [["MENDENHALL_WELLER"]]

Adding Interaction Potentials

It is possible to add custom interaction potentials to rustBCA. In the future, a special user-defined interaction type will be implemented, that can handle any potential of the form P(r) + Q(r), where P(r) is an inverse power-series in r:

$polynomial in inverse monomial form: P(r) = \sum_{i=0}^N\frac{a_i}{r^i}$

and Q(r) is an exponential inverse-polynomial of the form:

$Exponential polynomial of the form: Q(r)=\sum_{j=0}^M\frac{b_j}{r_j}\sum_{k=0}^Le^{-c_kr}$

For now, implementing a potential in rustBCA consists of the following steps:

in main.rs, add your potential to the InteractionPotential enum using the following syntax: InteractionPotential::INTERACTION_NAME{parameter_1: f64, parameter_2: f64, ...}
add a verbose description of your potential in the impl fmt block so that error messages involving your potential print correctly.
derive the singularity-free distance of closest approach function for your potential by multiplying F(r) = 1 - V(r)/Er - p^2/r^2 by r to the power of the max inverse-degree of r in F(x), r^max(2, max(N, M))
add the interaction potential to interaction_potential() in interactions.rs
add the singularity free distance of closest approach function to interaction.rs
If M=L=0, derive the equivalent monomial-form polynomial coefficients for the polynomial rootfinder
For the CPR rootfinder, implement a scaling function of the form 1/(1 + (r/a)^max(2, max(N, M))), where a is an appropriate scaling distance.
Determine the first derivative of your singularity free distance of closest approach function and add it to interactions.rs
If necessary, add assertions to main.rs to prevent bad combinations of potential/rootfinder/scattering integral.

I recommend a computer algebra system such as Maple, Mathematica, Matlab, or sympy for manipulating the distance of approach function and taking derivatives - tiny errors in the derivative, for example, can break the Newton-Raphson rootfinder and the Newton polishing of the CPR rootfinder.