Error when using a new interaction potential

[jackstahl]

I attempted to add a new interaction potential based on results from a DFT code. The best fit of the data is a modified Lennard-Jones potential:

16.0*( (1.87/xi).powf(6.5) - (1.87/xi).powf(6.0) )

After our conversation yesterday, I attempted to convert this V(r) potential to a dimensionless phi(r) potential by multiplying by r, without worrying about any additional scaling factor.

xi*(16.0*( (1.87/xi).powf(6.5) - (1.87/xi).powf(6.0) ))

This code failed to run, seen in the image below, and I haven't been able to find the backtrace variable to obtain a better backtrace. Let me know if there's any additional information I can give that would be helpful when diagnosing this problem.

[dcurreli]

@jackstahl Plot here the potential you are using. Add also a reference BCA potential, eg. Kr-C, for visual comparison. Make sure to use the same identical formula and units you are entering in the code.

[drobnyjt]

Two possibilities here:

1) The quadrature is not failing, but is producing results that break the center-of-mass to lab-frame transformation, resulting in the theta.is_nan() check panicking.

2) The quadrature is failing for this interaction potential - which would need implementation of an alternative scattering. I hope this isn't the case.

I'll spend some time today investigating this; there may be some easy solutions (detecting when the energy transfer implies "bonding" in the scattering integral and avoiding the calculation in that case) but only playing around with it will tell.

[drobnyjt]

Some updates on this:

The biggest issue right now is that the particular LJ potential (6.5-6) is acting up when I try to calculate the distance of closest approach. I've been able to get it working in Maple, but Maple's root finding algorithms are significantly more advanced than the simple Newton-Raphson I've been using so far for the purely repulsive potentials, and I haven't been able to replicate that success by hand in Python.

One possible strategy is an old method I'm looking into for LJ potentials that involves a complicated transformation of the scattering integral.

[jackstahl]

Hi Jon,

Sorry for taking a few days to get back to you. Thanks so much for looking into this issue for me; I'm glad that we were able to identify what's been causing this issue. Is there anything that I can be doing to assist you in finding a fix to this problem?

[drobnyjt]

@jackstahl Hey Jack - no worries. If you're up to it, the most helpful thing would be to go through this paper by Robinson (let me know if you need the pdf) and replicating his steps for the L-J potential in a higher level language like Python or Matlab. I think there's a contradiction in the text, so if you could get his scattering integral formulation working as a prototype, that would be the fastest way for me to get it implemented into rustbca; I've been held up by other work but this would be my next step as well.

[drobnyjt]

@jackstahl Hey Jack - just a quick update. For L-J potentials, you may want to try implementing the Aberth–Ehrlich method for simultaneous polynomial root finding; the distance of closest approach function (making some assumptions about the positivity of all the constants) can be rearranged into a polynomial which should be solvable by that method. There's an open source C implementation called MPSolve but no official or semi-official Python implementation I've been able to find yet.

[drobnyjt]

Here's a graphical example of the issue. For low energies, R (distance of closest approach) < p (impact parameter) and the distance of closest approach function becomes multivalued. The appropriate polynomial for an n-m L-J potential would be something like:

R^(n - m)*(R^m + 2*epsilon/Er*rm^m) - epsilon/Er*rm^n - p^2*R^(n - 2) = 0

where R is the distance of closest approach, epsilon, m, n, rm are the parameters of the LJ potential, Er is the relative energy (E0*M2/(M1 + M2) and p is the impact parameter.

Using scipy.optimize.fsolve and x0=p we get:

Which is largely similar to the confusing plot in Robinson (1997) for the simpler Morse potential:

This can be visualized in this sketch of the possibilities for a purely attractive potential, from S. A. Khrapak's Scattering in the Attractive Yukawa Potential in the Limit of Strong Interaction (2003):

According to Robinson, we need only the largest root of the distance of closest approach function - so using a simultaneous root-finding method like Aberth–Ehrlich and then choosing the largest is almost certainly the way to go (for L-J potentials). For Morse-style potentials and other more complicated ones, something else will probably be required.

[drobnyjt]

@jackstahl actually, I just discovered numpy.roots which seems to work perfectly well! We need an equivalent in Rust still, but having a working example in Python is a great step. Below are the roots from numpy.roots (dots) plotted with the value obtained from fsolve (lines) - the fsolve values are not always correct!

[jackstahl]

@drobnyjt Thanks for the update! I've been busy with the qual recently, but next week I'll take a closer look at this. If I'm understanding your plot correctly, it's really interesting that we can see exactly where the calculations are failing. Are p/a and R/a the impact parameter and apsis, same as the plot above?

[drobnyjt]

Yeah, the numerical root finding is failing right when the number of roots drops from 3 to 2, and it oscillates wildly around a false solution for a while when there's really only one root.

And yes - p/a is the impact parameter divided by the Lindhard screening length (0.8853*A0/(sqrt(Z1) + sqrt(Z2))^(2/3) and R/a is the apsis divided by the same screening length.

[drobnyjt]

@jackstahl Some progress on this - apparently, finding all the roots of an analytic function (which I would like for potentials such as the Morse potential, or Lennard-Jones-Coulomb) is a captial-H Hard problem; as far as I can find, the first truly general-purpose algorithm wasn't written about until the early 2000s.

Fortunately, there are open source implementations of the algorithm in Python and Fortran:

https://pypi.org/project/cxroots/ https://www.sciencedirect.com/science/article/pii/S0010465599004294

Now, for the more narrow L-J case in rustbca, the following changes need to happen: 1) add Gauss-Mehler or Gauss-Legendre quadrature, or re-derive Mendenhall-Weller for the case of arbitrary potentials (instead of screened coulomb) 2) add an option in the distance of closest approach function for polynomial roots 3) implement a polynomial root finding algorithm (or find one already implemented - this would be preferable) 4) Lennard-Jones type potentials experience an orbiting condition (see the figure below) for certain energies and impact parameters, where the scattering angle approaches negative infinity - we should verify that there's enough numerical robustness and NaN-tracking to make sure this doesn't break anything

I've done 1,2, and 3 in Python. In fact, Gauss-Mehler and Gauss-Legendre are going to be covered in my first 598 lecture - they're pretty easy. I'll keep you updated on further progress.

[drobnyjt]

Just for informational purposes, here's the real part of all roots of the distance of closest approach function for Lennard-Jones 12-6, computed with the cxroots package (which unfortunately, takes ~100X longer than fsolve()):

You can see that the places where fsolve() is failing are places where the function has a complex root (presumably, with an imaginary part near zero)!

[jackstahl]

It looks like the branching points at a given epsilon are the places where fsolve() is failing. Does that mean just those points have a complex root, or is that where the entire function begins to have complex roots? And why would it be important that the complex component is near zero?

Would you also be able to elaborate a little bit on why we only need to look for the largest root of the distance of closest approach? My understanding from the Khrapak image above is that any time the incoming particle loops completely around the target particle, you would get multiple roots, and in that case would only care about the first. Is this accurate?

Thanks again for your help!

[drobnyjt]

So, what I think is going on is the following:

F(z) = 1 - V(z)/Er - p^2/z^2 is the unmodified distance of closest approach function.

For arbitrary z, F(z) = x + y*i (i.e., F(z) is a complex function). When fsolve is iterating, it's finding a value z0 = x0 + 0.0*i which is very close to an actual complex root ofF(z1 + y1*i) = 0 where y1 is very small (within the tolerance of fsolve, but not within machine tolerance). Since the function F(z) is smooth, if y1 is small enough, it might appear to the solver to be numerically close to a real root (even though it's a complex root) since the error function is the L2 norm (i.e., error = (r0 - r0_new)^2, which would wash out any complex information. I never took an analysis course, so my reasoning may not be 100% correct, but I think it's something close to this.

Robinson (1997) claims the largest root is simply the correct lower bound of the scattering integral, which corroborates a lot of older papers I've read, but I am not 100% sure on the physical reasoning there - I bet a classical mechanics text could elaborate. But integrating from the largest root makes the trajectory symmetric and also avoids integrating around the atom, instead just integrating on either side of it, so that may be it.

[drobnyjt]

Some more food for thought:

John Boyd proposes (and proves!) that Chebyshev-proxy root finding is the optimum for finding real roots of a univariate function on the interval [a,b], in the paper Finding the Zeros of a Univariate Equation: Proxy Rootfinders, Chebyshev Interpolation, and the Companion Matrix and in his (very easy to read) text Solving Transcendental Equations .

I think this may be a strategy worth pursuing. Luckily, we have bounds a and b for the distance of closest approach - a is zero, obviously, but b we can reason out knowing that the distance of closest approach for two interacting atoms will always be less than that for their fully-stripped, coulombic interaction (well, unless you have two interacting negative ions, but that's a special case that probably requires some additional care).

Thankfully, the Rutherford distance-of-closest approach is analytic, giving us b: (Z1*Z2*e^2/4/pi/eps0/Er + sqrt((Z1*Z2*e^2/4/pi/eps0/Er )^2 - 4*p^2))/2. With a and b known, the Chebyshev proxy algorithm should work for any interaction potential of the form P_n(x)*exp(-a*x) + Q_m(x), where P(x) and Q(x) are inverse polynomials of x, i.e.P(x) = sum(a_i/x^i, i=0..n), Q(x) = sum(b_i/x^i, i=0..m) since we can always get rid of the singularity at zero by multiplying the distance of closest approach function through by x^(max(m,n)) without changing the roots we care about.

I have a version working in Maple, I just have to code it up in Python and then Rust and it should be a nice, guaranteed-to-find-all-real-roots replacement of Newton's method (although best-practice does still use Newton's method to refine the results of the Chebyshev-proxy method and to reject any spurious roots introduced by the approximation).

A simpler choice, outlined in Section 6.6.2 of Boyd's book (which you should be able to access online through the library), is an upgrade to Newton's method called the Newton-Armijo method, where the traditional Newton update is multiplied by a small parameter gamma such that x_n+1 = x_n - gamma*f(x)/f'(x) where gamma is searched for by gamma_n,j = 1/2^(j - 1)/2 for j=0,1,2...j_max, and the gamma that results in the smallest residual is used for the n+1 iteration. This process is referred to as a "line search" for gamma. However, there's still no guarantee that it won't skip over the desired largest root, it's just more likely. I'd prefer to have a root finder that works universally, but maybe it can be proved analytically that for specifically the distance of closest approach function that there's an idea initial guess that always converges to the outermost root? If that's the case, then we could get away with just Newton-Armijo or just Newton.

[drobnyjt]

I've implemented a very basic Chebyshev-proxy rootfinding algorithm here:

https://github.com/drobnyjt/chebyshev_proxy_rootfinding

I think it's robust enough to be good for this application, with the following caveats:

1) the distance of closest approach function must be transformed (e.g., substitute y=1/r in the case of a polynomial V(r)) to remove the singularity at r=0. This probably has to be done manually, as I don't really want to implement a Computer Algebra System (e.g., Maple, sympy) into rustbca. 2) the distance of closest approach function must not grow too rapidly as r->+inf, as this introduces numerical issues in the Chebyshev approximation (Boyd calls this the issue of "dynamic range"). This can be handled in two ways: first, with subdivision (i.e., when the degree N of the Chebyshev series gets larger than a user-specified Nmax, split the interval [a,b] in half and find Chebyhev series of each interval [a, b/2] [b/2, b], but I haven't implemented that yet in this Python prototype; second, with a smooth, zero-free scaling function S(r) giving us G(r) = S(r)F(r), where the zeros of G(r) are the same but it grows less quickly. One scaling function that works okay in my testing is S(r) = exp(-C*r), where C>1. Unfortunately, the more the scaling function squishes F(r), the worse convergence of the Chebyshev series becomes. After further experimentation I'll decide if I want to implement automatic subdivision or if we can get away with just scaling F(r). 3) the upper bound b must be as small as possible since F(r) grows as r->+inf. It'll require experimentation to see if the Rutherford upper bound is satisfactory for this purpose.

[drobnyjt]

Well, I went ahead and implemented subdivision. The issues with scaling functions as I played around just got worse and worse - with a strong scaling function, if it approached zero too quickly, the actual function F(x)*S(x) over some interval became numerically equivalent to zero and thus the entire interval was erroneously being identified as roots... But, with a mild scaling function and subdivision, the root-finder is now really quite robust.

Here's an example, for F(x) = r^12 - 1 + 2*r^6 - r^10 (an LJ type potential) and S(x) = exp(-r):

As you can see, the rootfinder automatically subdivided the domain into 3 intervals. It matched two of the intervals with degree 12 Chebyshev series, and one with a degree 24 Chebyshev series. Calculating the roots over each interval is equivalent to finding the eigenvalues of a sparse N-1XN-1 matrix, which is trivial at that size.

Even without additional Newton iterations to polish the root, it still matched 6+decimal places with the fsolve-calculated root. fsolve() is a compiled function, so one would expect it to be much faster than a rootfinder implemented in nearly pure Python. However, this Python CPR rootfinder takes only 0.006 s (Adaptive Chebyshev Interpolation with Subdivision) + (0.001) s (Eigenvalue calculation of roots) versus ~0.0009 for fsolve() - only about a factor of 10 worse!

It's going to be much slower than just Newton, but it has the massive benefit of always identifying real roots and always finding all of them - which, for an arbitrary attractive-repulsive potential, is absolutely necessary to get the right answer for the scattering angle.

A recent paper* on the LJ/s potential (Lennard-Jones-spline) verified their distance of closest approach calculation by running full two-body physics simulations - this is almost guaranteed to be faster!

*https://www.frontiersin.org/articles/10.3389/fphy.2020.00271/full

[drobnyjt]

And just for fun, here's a slightly pathological example - F(x)*sin(x^2):

Evaluation of all the roots took ~0.1 s.

[drobnyjt]

And the CPR rootfinder for a realistic LJ 12-6 potential (epsilon = 0.343 eV, sigma = 2 A):

fsolve() still finds the complex roots with imaginary parts near zero; CPR has no problems with them, and is much, much faster than the complex rootfinding algorithm in cxroots.

[drobnyjt]

And, just to drive the point home on the inherent problems with iterative methods, here's the same plot for the LJ 6.5-6 potential, even if I feed fsolve with the CPR root:

[drobnyjt]

There's one more thing I want to investigate with the Python version before I implement this in Rust. On page 36 of Boyd's Solving Transcendental Equations, he describes a way to map an interval [y0, +inf] to [-1, 1] and produce the so-called Rational Chebyshev basis, which is defined on [y0, +inf].

Practically, this means that we could take F(r), the doca function, transform it with r = L(1 + x)/(1 - x) to G(x), find the roots of G(x) on the interval [-1, 1] using the CPR, and transform the roots back with r0 = L(1 + x0)/(1 - x0) - all without any upper bound b.

Trying this naively, I've run into two difficulties - first, at x = 1, the coordinate r -> +inf and makes the computer unhappy, meaning we have to truncate the interval to x=0.999 or some other number close to 1. Second, around x=1, a lot of spurious zeros are created. These are, in principle, able to be excluded by running Newton's method on them and seeing if they converge, but with a scaling function S(x) = exp(-x) they're so close to zero that Newton's method is perfectly happy to converge to them...

The algorithms in this text were designed with arbitrary precision systems in mind (Maple, MATLAB, etc), so some edge-cases are bound to fail when using 64 bit floating point numbers.

[jackstahl]

Wow, thanks for the detailed updates! I appreciate the glimpse into the process of solving issues like these. It's good to see encouraging results so far, please keep me updated with any more progress!

[drobnyjt]

Yeah, absolutely - these posts are as much for my (future) benefit as yours, since this is a convenient place to write down everything I try in case I forget later.

Another update - I have the Chebyshev proxy rootfinder working in Rust. I caught a small error in my implementation that was forcing the algorithm to subdivide too early, but a throwaway comment in Boyd's book led me to a corrected version. Essentially, sometimes the Chebyshev series sum (a_i * T_i(x), i=0..N) has trailing coefficients that are ~0. The eigenvalue problem uses the last coefficient, a_N, as the denominator of the last row of elements. So, when it's near zero, the eigenvalue problem breaks. The solution is to just keep chugging along with the adaptive interpolation, and when the matrix is constructed, to toss out the trailing terms of a_i below some user-defined threshold. Using this correction, the Python version can find all 300 or so roots of this monster:

(x - 2)(x + 3)(x - 8)(x + 1E-4)(x - 1E-5)(x + 1)sin(x^2)

in ~3 seconds! The Rust version is even faster - typically less than 0.2 seconds, even while running in debug mode with backtracing and no optimizations. Running under release mode, it calculated all the roots in ~0.05 seconds!

Now unfortunately, the eigenvalue calculation relies on BLAS, which is a pain to install on Windows. For now, I have shifted development back to my Linux machine, but I'd like the final code to remain portable so I'm searching for a solution to that.

[drobnyjt]

Good news - I have successfully imported the new root-finding routine into rustbca!

It took some fine-tuning to get working (for whatever reason, the Rust library ndarray-linalg, needed to compute the eigenvalues of the companion matrix to find the roots, refused to work with my openblas installation, so I had to switch to netlib), but now that it's there it's pretty straightforward.

I'd like to expose those root-finding parameters (I'll note them at the bottom of this comment) to the user in the input file, as well as figure out a scheme for intelligently determining some of them. The downside of this rootfinder is that it is significantly slower than Newton alone - but, on the plus side, it's all-but guaranteed to give the right answer (unlike Newton-only, which only gives the right answer for simple roots and purely repulsive or purely attractive potentials) - when the wrong root is chosen, the quadrature error goes through the roof and the scattering angles are all wrong. On my test run using the CPR rootfinder, the time per ion is up to 0.001-0.006 seconds from ~0.00001 seconds, so it's a significant cost (100X), but for an attractive-repulsive potential the correctness guarantee is likely worth it.

Next, I have to add some new Gaussian quadrature integration into rustBCA. The current version only handles screened potentials, but thankfully Gaussian quadrature is ~6 lines of code, so that should be easy going.

I expect a prototype Lennard-Jones-type BCA run by early next week.

[a, b]: interval to search for zeros on. Currently b is set to an arbitrarily high value of 100 screening lengths N0: initial degree of chebyshev polynomial - 5th degree seems a decent starting point for this problem, as it converges rapidly and for most parameters doesn't have to subdivide the interval. epsilon: absolute error of Chebyshev interpolation. 1E-3 is relatively high, but it speeds things up. N_max: maximum degree of Chebyshev polynomial - solving the roots involves solving an NXN matrix, so if matrix solving is faster than subdividing, this number should be high complex_threshold: roots with an imaginary part below this value are kept truncation_threshold: Chebyshev terms with coefficients below this value will be thrown away interval_limit: the minimum allowed subinterval size - if you hit this, you're on your way to a stack overflow so it panics

[drobnyjt]

Gauss-Mehler and Gauss-Legendre have been implemented and tested - they're outputting scattering angles to within 0.1% of the Mendenhall-Weller quadrature (which is a modified Gauss-Lobatto). Right now everything is still stuck in screened-coulomb mode, but as I move the scattering integrals away from integer flags for selecting interaction potentials to function pointers, that'll get updated too.

[drobnyjt]

@jackstahl Alright, the code is messy and needs to be cleaned up but! I've performed the first run of the bca using a Lennard-Jones 12-6 potential. A clean working version should be ready to go in the next day or two.

[drobnyjt]

Another update, this time a theoretical one: I just found an analytical form for determining when any reasonable LJ potential has multiple roots in the distance of closest approach function - I'll outline the solution (verified by Maple) here so I can reference it later and just in case someone else catches a flaw in my reasoning:

F(r) = 1 - V(r)/Er - b^2/r^2is the distance of closest approach function. Er is the CoM energy, b is the impact parameter, r is the interatomic distance, and V is the interaction potential. For an LJ potential of type m-n:

V(r) = 4*epsilon*((sigma/r)^m - (sigma/r)^n) where m>n, m and n are integers and usually n divides m - that is, there's an integer p>1 such that m = n*p

We can't solve for the values of the roots directly, that's immediately non-analytic (and the whole point of the CPR root finder), but we can solve for b(r), the impact parameter, instead, which has one positive root:

b(r) = sqrt(Er*(V(r) + Er)*r)/Er

This can be plotted (plot below) which is great! It's the same plot as the numeric plots of r(b) above, but rotated 90 degrees. Unfortunately, inversion of this is just the same root-finding problem we started with, so this doesn't immediately help. However, we can see that r(b) has multiple roots when b(r) has extrema - so taking db(r)/dr=0, we can find the value of r that produces extrema, which is the root of this polynomial:

2*epsilon*(sigma/z)^m*m-2*epsilon*(sigma/z)^n*n-4*epsilon*(sigma/z)^m+4*epsilon*(sigma/z)^n+Er

When this polynomial has all real roots, the extrema exist, and r(b) is multivalued for some impact parameters - when it has all imaginary roots, the extrema go away, and r(b) only has one root (and we can use Newton's method).

This polynomial can be cast into this form:

A = (4*epsilon/Er - 2*epsilon/Er*n)*sigma^n
B = (2*epsilon/Er*m - 4*epsilon/Er)*sigma^m
z^m + A*z^(m-n) + B = 0

Since m = p*n, we can u-substitute with u=z^n:

u^p + A*u^(p-1) + B = 0 And now we have 3 cases where this is analytic: p=4, 3 and 2 (p=0 is trivial). Using the discriminants of degree 2, 3, and 4 polynomials, we find, for each case, the condition when multiple roots of r(b) exist:

p=2: A^2 - 4*B
p=3: -4*A^3*B - 27*B^2
p=4: 256*B^3 - 27*A^4*B^2

For the LJ 12-6 potential, this works out to a very, very simple condition for the distance of closest approach function having only one root: Er > 4/5 epsilon! Which is incredibly simple to check for in the BCA and simply run Newton for those cases instead (since now it's guaranteed to converge to the correct answer).

Unfortunately, this only works when m=p*n, p<5. But it covers the classic 12-6 and 9-3 and some other unusual LJ potentials (16-4) that nobody's covered before. As far as I could search, nobody's done this exact analysis before - all the LJ papers that cover scattering theory either use an iterative method with a heuristic, find the root using the eigenvalues of a companion matrix, or use an RK4 two-body solver to find accurate distances of approach. This should speed up the CPR case by a lot, since E < 4/5 epsilon only occurs at the very, very end of the particle trajectories.

Blue: E =0.01 eV Red: E= 4/5 epsilon Blue: E = 10 eV

As you can see, at 4/5 epsilon, the curve entirely flattens out, and the multiple solutions to r(b) vanish!

[drobnyjt]

@jackstahl Hey Jack, I've got a prototype version working on the branch cpr_rootfinder on the repository - could you send me the input file for your specific problem and the values in SI units for the LJ-6.5-6 potential you're using (epsilon, rm/sigma)? I'd like to do some tests on an actual problem. I've done some limited testing, and it appears that the CPR rootfinder parameters need to be chosen carefully for each ion-atom pair - the relative energy is playing a role in the magnitude of the DOCA function at arbitrary r, so the absolute tolerance epsilon needs to be adjusted up or down depending on the reduced mass.

Currently, they're controlled by the user in the options in input.toml: cpr_n0 cpr_nmax cpr_epsilon ... So some playing around will be necessary for each problem.

[drobnyjt]

One more update as I close this issue - the code is working, as far as I can tell, perfectly with the new interaction potentials and rootfinders. the CPR rootfinder isn't happy with the 6.5-6 form, however - I think the Newton iterations used to refine the guesses from the Chebyshev interpolant are converging very poorly, since it's not a nice integer-degree polynomial. However, by using the substitution r = y^2, I can convert it into a 13-degree polynomial and use a polynomial rootfinder, so long as I convert the roots back using the same substitution!

So on the rustBCA wiki, there's an example using a 12-6 potential: https://github.com/lcpp-org/RustBCA/wiki/Example:-Multiple-Interaction-Potentials

Just be sure and use the rootfinder {"POLYNOMIAL"={complex_threshold=1E-12}} for the 6.5-6, neither Newton nor CPR can successfully churn through it. There may be an algebra mistake in my distance of closest approach and derivative of the distance of closest approach function - which would slow down Newton a lot - but in any case, it's not necessary, the polynomial root-finder is faster than CPR so it should be used whenever the potential can be cast into polynomial form, such as is the case here.

[drobnyjt]

Oh, note also that the input file has changed dramatically - I've updated all the pages on the wiki with the new format. It should be a lot more self-explanatory, at the cost of being a little harder to write by hand.

[jackstahl]

Hi Jon, sorry for the radio silence lately-- another big project that I'm on has been finishing up and I haven't had the time to catch up with what you've been doing here. I'll test this out tonight and let you know how it goes!

[drobnyjt]

No problem - let me know if you run into any issues!

On Sep 20, 2020 3:04 PM, Jack Stahl notifications@github.com wrote:

Hi Jon, sorry for the radio silence lately-- another big project that I'm on has been finishing up and I haven't had the time to catch up with what you've been doing here. I'll test this out tonight and let you know how it goes!

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHubhttps://github.com/lcpp-org/RustBCA/issues/31#issuecomment-695829768, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AJBUE2HC5VIKP7F2NDIY3IDSGZN3ZANCNFSM4PEA75BA.

[jackstahl]

Just to check my understanding, it seems that I need to follow the Ubuntu installation instructions with WSL in order to compile the build with rcpr on a Windows machine. Would this also involve reinstalling python and its libraries as well?

[drobnyjt]

@jackstahl Yep, that should be right. Python isn't strictly required (and in fact the rustbca python scripts are probably outdated now anyway) but if you want to use it to run the code directly, I believe you will have to install python on WSL.

[jackstahl]

Hi Jon, I ended up installing Ubuntu on a VM and was following along with your instructions there. Everything was working as expected until I tried compiling rustBCA, where it failed when building netlib-src. I tried this on my Windows installation as well, just to check, and ran into the same error. Could there have been an update to this program causing it to fail to install?

[drobnyjt]

One thing to try would be to git clone rcpr and replace the line in rustBCA's Cargo.toml that downloads rcpr from github with your local version. Then, in rcpr's Cargo.toml, under the ndarray-linalg dependency, replace "netlib" with one of the other lapack/blas backends - first, try "openblas" and see if you get a different error message.

[drobnyjt]

Oh, and this is a silly question, but is cmake installed? It may be silently necessary (and have to be added to the build instructions).

[jackstahl]

I don't believe it was, actually. I'm downloading it to test right now.

[jackstahl]

I installed cmake through the ubuntu storefront. This seems to have made some difference because the error message has changed, though I'm not sure how to troubleshoot the Fortran compiler in the VM.

For now, I'll look into your previous suggestion on replacing the rcpr line.

[drobnyjt]

Ok, I've seen this error message before - it's not happening at the rustBCA or rcpr level, it's actually a problem with ndarray-linalg, one of rcpr's dependencies. That page has some troubleshooting hints that may help. One thing to try would be to make sure gfortran is installed:

sudo apt-get install gfortran

And then changing the environment variable that cmake is looking for, "FC" to gfortran:

export FC=gfortran

[drobnyjt]

I'm giving it a go on my Windows machine running WSL right now.

[jackstahl]

gfortran was not installed; I had used a fresh VM for this. In order to change the cmake environment variable, do I need to first navigate into the cmake files?

[drobnyjt]

Nah, you shouldn't have to do that; they're being built automatically by the rust cmake crate so I don't think they're even manually editable. Try just setting it as a user environment variable by the command line.

If that doesn't work, try installing the following: sudo apt-get install libatlas-base libatlas-dev libblas-dev liblapack-dev

After installing all of those I was able to get it to build on WSL with the cpr_rootfinder - so clearly one of them is necessary! Also, make sure to run cargo update, mine had an old version of ndarray-linalg cached that was giving me some new problems.

[jackstahl]

After running export FC=gfortran in the command line, the installation seemed to stall out. I didn't have to run for all that long, but it attempted to build netlib-src for a few minutes with no errors or progress.

I also tried installing sudo apt-get install libatlas-base libatlas-dev libblas-dev liblapack-dev, and got the following message:

[drobnyjt]

I must have mistyped the command - go ahead and use the suggested replacement!

Also, netlib-src does just take 10-15+ minutes to build for the first time on my machine (which has a previous generation i5)! So it may have been on its way towards success.

[jackstahl]

With the replacement I was still receiving an error:

E: Unable to locate package libatlas-base

For now I checked that cargo was updated, and am going to let the installation run for about 20 minutes. Hopefully the lib packages won't be needed! I'll let you know how things turn out afterwards.

[drobnyjt]

Here's how I got it to work with a fresh Ubuntu install on WSL:

install rustup run the suggested shell command to add cargo to path git clone rustbca Now here, it failed with linker cc not found so I ran the following: sudo apt-get update sudo apt-get install build-essential gcc Now it gave me a generic permission denied os error, but it mentioned cmake files in the error message, so I tried: sudo apt-get install cmake And now I get the "unknown Fortran compiler error." So next: sudo apt-get install gfortran And now netlib-src builds successfully.

And now, actually, it just runs. I know on my Linux Mint machine I had to install at least the liblapack-dev package for it to run, but perhaps that comes with Ubuntu for Windows Subsystem for Linux? I'm going to add this info to the installation page - let me know how your effort goes!

[jackstahl]

Made it much further along with the installation-- hdf5-sys v0.6.1 failed now. It looks like it can't find the HDF5 installation at all, so maybe it hasn't been installed yet?

[drobnyjt]

Oh, go ahead and ignore HDF5 and don't build with it - that's only for running very large (~GB) simulations coming from a particle source, such as a plasma sheath solver, where it would be impractical to have the code write all billion particles to the single line of input.toml that stores particles. It's an experimental feature that's not quite ready for external use yet.

[drobnyjt]

If in the future you'd like to use the HDF5 input, I think all that's required is: sudo apt-get install libhdf5-dev