feos Extending EoS capabilities

Dear all, I am a PhD in the Engineering Thermodynamics group at the TU Delft. We mostly work on molecular simulations with Monte Carlo and molecular dynamics simulations in Fortran and C++. In my free time, I am learning Rust (I just finished coding my own LJ Monte Carlo program for fun) and came across this package because an office mate was using FeOs. He is using FeOs (from the python side) for phase equilibrium predictions of multi-component mixtures, where he computes among others spinodal lines.

Besides FeOs he is using Refprop and the GERG-2008 EoS (doi: 10.1021/je300655b) and asked me if it would be possible to implement this EoS as well into FeOs. I looked into the Peng-Robinson implementation (python example in the tutorial and the Rust code in the feos-core crate) and in the half hour I spend on it, it seems that FeOs requires Helmholtz free energy functions for correct predictions of its state. This is similar to GERG-2008, although that is a multivariable empirical relation instead of something as elegant and simple as the already implemented PR EoS. I propose the following process, as I have mediocre experience in Rust:

I implement and test the cubic equations of states first in python to check if I have proper understanding of what I am doing.
I implement additional cubic equation of states first to get some Rust practice. I went through the code and expect that I would not have too many problems in implementing the VdW and SRK EoSs for single and multi-component fluids.
When that is implemented and merged successfully in the main FeOs branch, I will start on implementing empirical Helmholtz free energy models such as GERG-2008 and maybe IAPWS95.

For this, I would need the following from the FeOs maintainers:

Some instructions on structure, styling, and testing guidelines in FeOs (as to integrate it well).
Are there other dependencies and requirements besides those mentioned on the prerequisite page in the documentation or included in the Cargo.toml files? (maybe there is some C++ Makefile trauma speaking here :) .)
Besides the: feos/crates/feos-core/src/cubic.rs, which other files should be changed to produce all required python bindings and such?
And connected to that. Why are the cubic EoS rust files in the src root directory and not in the equation_of_state subdirectory like the ideal gas law?
While I have some ideas already on how VdW and SRK can be implemented, I expect more issues with the GERG model. Do you foresee pitfalls for such empirical models, and can you give some advises on how to approach it?

I believe that FeOs would get the following gains from this proposal:

Having extensive, fast cubic EoS support might attract CFD crates to use FeOs for the prediction of fluid properties during simulations. Compressible CFD simulations require good and most of all fast EoS predictions for density, heat capacity and viscosity. I believe that FeOs can be attractive for such programs when there are more cubic equations of state to choose from.
When GERG-2008 is implemented, FeOs can claim that it can compute fluid properties compliant with the ISO 20765-2/3 standard. Scientifically, I do not see too much value in this, however, it might make FeOs suitable for usage in industry.

I cannot work on this project professionally, as my PhD is about computing transport properties in complex aqueous salt mixtures and I do not intend to use FeOs for my own projects. However, I think I can learn aplenty from this side project (especially regarding traits) and I am willing to spend some evenings on it.

Would you guys be interested in such extensions to FeOs? Kind regards, Jelle Lagerweij

May 14 '25 21:05 JelleLagerweij

Hi Jelle,

thank you very much for the issue and your interest in FeOs. I'd happily guide you through the steps needed to implement a new model. You can also ask here at any point or join our discord server (link is in the repo-readme).

The general procedure could look as follows:

Having an implementation that you can compare against is very helpful. So first implementing it in Python (or using existing implementations like thermopack, Clapeyron.jl or teqp) is a good first step.
Then, implement the model in rust and test it in rust. The Python bindings come into play once everything is working correctly. The Python stuff adds a layer of complexity so I'd advise against developing the bindings in parallel to the model.
Think about how to structure your model's parameters - we have some traits that can help with I/O.
Add the Python bindings. For this you need to learn PyO3 (which is awesome) but since we already have several models I think adapting an existing example is not too hard.

Helmholtz energy

In FeOs, residual equations of state (without the ideal gas part) are implemented by defining the Helmholtz energy as function of the thermodynamic conditions. These conditions are not encoded as regular floating point numbers but as generic (hyper-) dual numbers. In essence, these numbers enable us to write a single function (Helmholtz energy) and yield all exact partial derivatives without an analytical implementation. So all you have to care about is how to implement the model - properties are calculated via automatic differentiation.

Because we need to formulate the model as dual numbers, you'll find that the relevant functions in the Residual trait, which you'll have to implement, are parameterized over the DualNum trait. DualNum implements all relevant arithmetic operations so you can use dual numbers almost like regular floats (there are some quirks to get used to).

So in a nutshell: Implementing a new equation of state means implementing the Residual trait (and the IdealGas trait if your model needs a specific ideal gas model).

Your questions

Structure/style/testing: Adding a model is done in the feos crate (crates/feos directory). Each model has its own module in feos and is gated behind a cargo-feature. We use rustfmt to format everything and clippy as linter. Tests are usually added inside the module where the functions are defined (see existing models for examples).
Dependencies: No, everything is defined in Cargo.toml :D
Existing cubics module: The implementation of Peng Robinson is purely pedagogical and for having a model to test in feos-core. New cubic models will be in the feos crate.
Model complexity cubics vs GERG: nothing to worry about I think. GERG may require a dedicated ideal gas model iirc? Besides that, only parameter handling will probably a bit more involved but implementing the equations should be fine.

In my free time, I am learning Rust (I just finished coding my own LJ Monte Carlo program for fun)

Awesome, that's how I started with Rust as well :D Our molecular simulation lectures use a Rust code with Python bindings nowadays (an updated version of this code).

May 15 '25 08:05 g-bauer

I started adding cubic equations of state (in a generic way) some time ago. The branch is still here, although we since changed the project's structure a bit (see here). Maybe you can use that as a starting point.

Also, when working with a new model in FeOs, configure your IDE/editor so that only your model is active (disable all other feos-features) and use cargo build (develop profile, not release or even lto). That way, the crates compile way faster.

May 15 '25 08:05 g-bauer

Thank you @JelleLagerweij for the interest in contributing to feos and thank you @g-bauer for the detailed tutorial.

May 16 '25 09:05 prehner

Dear Gernot, Thanks for you quick reply. I see that your impl_cubics branch is indeed quite different (not only in project structure, but it also has a more general cubic eos setup compared to the cubic.rs file in the main branch). To start, I will do the following:

Create a jupyter notebook where I can prepare some test data for PR, VdW, SRK for single a single and two component mixture. I can indeed use Clapeyron.jl for this.
I test the current PR implementation in main with it (that should obviously pass :) ).
I correctly update the impl_cubics branch into the current structure. From your explanation, I understand that this shall be in /crates/feos/src/cubic/ and not /crates/feos-core/src). I will then retest PR and (S)RK using python and Clapeyron created integrated tests in the Rust files (I recognized some examples in your code).
If I then have no conflicts with the main branch and complete all tests we can merge it to the main branch.
I add the VdW EoS completeness and add tests for it. This should be quite trivial if the other steps are already working.

If that work out, then I will see how we can continue with the empirical models. Which, I expect, shall be their own module within the feos crate.

Then regarding sargas, I had some fun with reading through it. We of course have our own molecular simulation course, for which I have been the TA a few times already. We just give the students the book of Frenkel and Smit and tell them you have 3 weeks for coding your own MC simulation for a LJ system. Then give them 3 weeks for their own MD code (with Nosé-Hoover as thermostat) and lastly 3 weeks for transport properties post-processing script to get diffusion and viscosity. They can use the programming language of their choice (all choose python). We noticed that students struggle with the freedom and do not know how to structure their code. Additionally, they have problems with computational speed (#python_and_for_loops), especially for when they need to perform longer simulations to get the transport properties.

The way you structured Sargas might actually avoid these issues, as the structure is already there (included tests for their assignments) and the Rust code will probably perform quite a bit better. We will probably keep our own assignments, but sargas gives us some ideas on how we could make the student experience a bit better.

May 16 '25 09:05 JelleLagerweij

Correct, new models should be in crates/feos/src/model_name.

Feel free to take what I've been building, although some renaming (e.g. of the mixing rule) might be needed before merging. Refactoring things like that is quite easy though. If you want to start with less generic implementations - that's also fine. Even starting with the empirical models wouldn't be necessarily harder to do.

Regarding teaching molecular simulations - that's certainly a fascinating approach. But not keep this issue focused, let's move this to a different conversation. :)

May 16 '25 10:05 g-bauer

Dear all, I cloned the main branch of your repository and I am currently refactoring the cubic parts of the impl_cubic into the new structure. The current progress can be found in https://github.com/JelleLagerweij/feos_cubics (or in your fork list of course :) ).

I started by creating the mod.rs and the parameters.rs files. I had some difficulties with a few imports: feos_core::{EoSresult} and feos_core::parameter::{ParameterError}. The last one I found, you put all errors together in an error.rs file which matches the error case using enums. Therefore, I changed the parameters.rs file to use feos_core::{FeosError, FeosResult}.

However, I had more issues with understanding feos_core::{EoSresult}. The impl_cubic implementation for the pub fn new is EosResult<Self>. However, I noticed that the other EoSs return just Self. I guess that the difference comes from the generality that is intended for the cubic EoSs behavior (that all take critical point inputs ect.). What type shall be returned from the new() function in cubic?

Examples:

In crates/feos/src/cubic/mod.rs

impl Cubic {
    /// Generic cubic equation of state with adjustable universal constants.
    pub fn new(parameters: Arc<CubicParameters>, options: CubicOptions) -> EosResult<Self> {
        let p = CriticalParameters::new(&parameters, &options.delta);
        options.alpha.validate(&parameters)?;
        Ok(Self {
            parameters,
            options,
            critical_parameters: p,
        })
    }
  ...
}

In for example crates/feos/src/saftvrmie/eos/mod.rs

impl SaftVRMie {
    pub fn new(parameters: Arc<SaftVRMieParameters>) -> Self {
        Self::with_options(parameters, SaftVRMieOptions::default())
    }

  ...
}

May 16 '25 15:05 JelleLagerweij

Yeah we changed that recently. Now there only exists a single error type, FeosError.

There is no hard rule for the return type of the new function. If it makes sense to return a Result, because something can go wrong, that's fine. In this case returning a Result makes sense because the alpha.validate function is fallible.

May 16 '25 19:05 g-bauer

Dear all, I have made little progress regarding expanding the cubic equations of state support. I had a few busy weeks. However, in my free time, I tried to fiddle around a bit more with feos::cubic, mostly regarding the testing of the PengRobinson1976 which fails.

I compared with Clapeyron.jl for a logical liquid state point (where we should deviate strongly from the ideal gas behavior)

using Clapeyron
using StaticArrays

Tc = 369.83    # Critical Temperature in Kelvin
Pc = 4.21e6  # Critical Pressure in Pascal
ω = 0.153    # Acentric factor
Mw = 44.1   # Molar mass in g/mol

# Create a Peng-Robinson EOS object
eos = PR("propane", userlocations=(Tc = Tc, Pc = Pc, Mw=Mw, acentricfactor = ω))

# Get a reasonable compressed liquid statepoint.
T = 300.0  # K
v = 8.7e-5  # m^3/mol
z=SA[1.]  # mol
p = pressure(eos, v, T, z)  # Pa

It returns the following: P = 1.5476632110008812e6 Pa This is very close (3.4458935260772705e-8 Pa difference) to the PengRobinson1976 as implemented here. And quite different from the comparison state we are using (feos-core::cubic::PengRobinson is about 8047.94 Pa off compared to the the feos:cubic::PengRobinson1976). So I guess that the PengRobinson in feos-core::cubic somewhere has an issue. Any ideas where to start when investigating this? Or is dropping the comparision with feos-core::cubic fine if we replace the test with appropriate comparisons with Clapeyron.jl for now?

Besides that, I intended to use State<Cubic> instead of StateHD<f64> for testing our implementation (as that gives inputs to which we can control the units of the input, and it interacts more similar with the EoS compared with real usecases so we could use it as an example as well. I could use it easily to compute the pressure comparison.

let eos = Arc::new(
    Cubic::peng_robinson(
        parameters,
        Some(Alpha::PengRobinson1976(PengRobinson1976)),
        Some(MixingRule::Quadratic(Quadratic)),
    )
    .unwrap(),
);

// Build the test state
let temp = 300.0 * KELVIN;
let vol = 8.7e-5 * METER.powi::<P3>(); 
let mol = arr1(&[1.0]) * MOL;

let state = StateBuilder::new(&eos)
    .temperature(temp)
    .volume(vol)
    .moles(&mol)
    .build()
    .unwrap();

state.pressure(Total)  // Or IdealGas or Residual

However, if I wanted to get the Helmholtz energy (a = helmholtz_energy) I get a trait error:

error[E0599]: the method `helmholtz_energy` exists for struct `State<Cubic>`, but its trait bounds were not satisfied
   --> crates\feos\src\cubic\mod.rs:388:31
    |
116 | pub struct Cubic {
    | ---------------- doesn't satisfy `Cubic: feos_core::IdealGas`
...
388 |             state_implemented.helmholtz_energy(Total)
    |                               ^^^^^^^^^^^^^^^^
    |
    = note: the following trait bounds were not satisfied:
            `Cubic: feos_core::IdealGas`
note: the trait `feos_core::IdealGas` must be implemented
   --> C:\Users\vlagerweij\Documents\TU jaar 6\Project_FEOS\feos_cubics\crates\feos-core\src\equation_of_state\ideal_gas.rs:7:1
    |
7   | pub trait IdealGas: Components + Sync + Send {
    | ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

This is fully logical, because in feos-core/src/tate/properties.rs we require impl<E: Residual + IdealGas> State<E> { for the entire code block. This is not required for the pressure function in the residual_properties.rs file. Could you explain to me how to make Cubic satisfy feos_core::IdealGas?

Jun 02 '25 08:06 JelleLagerweij

@JelleLagerweij My first attempt at coding in Rust was to expand the cubic EOS to include the temperature-dependent volume translation: c = c0 + c1*(T - 288.15). I have put it on the back-burner because I have not yet had enough time or expertise to tackle testing in Rust. If it's not too much trouble, I would appreciate it if you could add this to your cubic EOS implementation.

This was my attempt at generalizing the cubic to calculate a, b, kappa depending on the specified version. There are undoubtedly better ways to do it that allow more flexibility with the alpha-function. It seems that @g-bauer was on his way to implementing that in his version.

impl Parameter for PengRobinsonParameters {
    type Pure = PengRobinsonRecord;
    type Binary = f64;

    /// Creates parameters from pure component records.
    fn from_records(
        pure_records: Vec<PureRecord<Self::Pure>>,
        binary_records: Option<Array2<Self::Binary>>,
    ) -> FeosResult<Self> {
        let n = pure_records.len();

        let mut tc = Array1::zeros(n);
        let mut a = Array1::zeros(n);
        let mut b = Array1::zeros(n);
        let mut molarweight = Array1::zeros(n);
        let mut kappa = Array1::zeros(n);
        let mut c0_arr = Array1::zeros(n);
        let mut c1_arr = Array1::zeros(n);

        for (i, record) in pure_records.iter().enumerate() {
            molarweight[i] = record.molarweight;
            let r = &record.model_record;
            tc[i] = r.tc;
            // Compute aᵢ, bᵢ in SI. Use r.pc already in Pa.
            let eos = r.eos_type.to_uppercase();
            match eos.as_str() {
                "PR76" => {
                    a[i] = 0.45724 * r.tc.powi(2) * KB_A3 / r.pc;
                    b[i] = 0.07780 * r.tc * KB_A3 / r.pc;
                    kappa[i] = 0.37464 + (1.54226 - 0.26992 * r.acentric_factor) * r.acentric_factor;
                }
                "PR78" => {
                    // aᵢ and bᵢ use PR‐78 coefficients
                    a[i] = 0.45308 * r.tc.powi(2) * KB_A3 / r.pc;
                    b[i] = 0.07780 * r.tc * KB_A3 / r.pc;

                    // kappa: if ω < 0.49, fall back to PR-76 form; else use PR-78 form
                    if r.acentric_factor < 0.49 {
                        kappa[i] = 0.37464 + (1.54226 - 0.26992 * r.acentric_factor) * r.acentric_factor;
                    } else {
                        kappa[i] = 0.379642
                            + (1.48503 - 0.164423 * r.acentric_factor) * r.acentric_factor
                            + 0.016666 * r.acentric_factor.powi(2);
                    }
                }
                "PR19" => {
                    a[i] = 0.45723553 * r.tc.powi(2) * KB_A3 / r.pc;
                    b[i] = 0.07779607 * r.tc * KB_A3 / r.pc;
                    kappa[i] = 0.3919
                        + (1.4996 - 0.2721 * r.acentric_factor) * r.acentric_factor
                        + 0.1063 * r.acentric_factor.powi(2);
                }
                "RKS72" => {
                    a[i] = 0.42747 * r.tc.powi(2) * KB_A3 / r.pc;
                    b[i] = 0.08664 * r.tc * KB_A3 / r.pc;
                    kappa[i] = 0.480 + (1.574 - 0.176 * r.acentric_factor) * r.acentric_factor;
                }
                "RKS19" => {
                    a[i] = 0.427481 * r.tc.powi(2) * KB_A3 / r.pc;
                    b[i] = 0.086641 * r.tc * KB_A3 / r.pc;
                    kappa[i] = 0.481
                        + (1.5963 - 0.2963 * r.acentric_factor) * r.acentric_factor
                        + 0.1223 * r.acentric_factor.powi(2);
                }
                _ => {
                    return Err(FeosError::IncompatibleParameters(
                        format!("Unsupported EOS type: {}", eos),
                    ));
                }
            }

            // Convert volume-shift from mL/mol → m³/mol
            c0_arr[i] = r.volume_shift_c0 * 1e-6;
            c1_arr[i] = r.volume_shift_c1 * 1e-6;
        }

Note that the RKS19 and PR19 are based on the work from Pina-Martinez et al: 10.1016/j.fluid.2018.12.007.

Jun 02 '25 17:06 cjsisco

This is fully logical, because in feos-core/src/tate/properties.rs we require impl<E: Residual + IdealGas> State<E> { for the entire code block. This is not required for the pressure function in the residual_properties.rs file. Could you explain to me how to make Cubic satisfy feos_core::IdealGas?

Hi, in this case, you just want to use state.residual_helmholtz_energy or state.residual_molar_helmholtz_energy, which are independent of the ideal gas contribution and probably the values you want to compare to the results from other implementations.

Regarding the deviations in feos-core::cubic, I don't have the overview over the cubics right now to answer that, maybe @g-bauer has a spontaneous idea, otherwise we can investigate it more thoroughly.

In that context also thanks to @cjsisco for the suggestion! Sounds like a good addition.

Jun 02 '25 21:06 prehner

As Philipp wrote, since we only want to compare the residual models, we don't need the ideal gas model. For completeness, though, you can use the EquationOfState struct from core - but then you have to specify the ideal gas model:

use feos_core::EquationOfState;

let eos = EquationOfState {
    ideal_gas: Arc::new(Dippr...), // or other models
    residual: Arc::new(Cubic...)
};

// eos has methods for <IdealGas + Residual>

Again, not necessary here, but might come into play when you think about GERG.

Regarding the differences between implementations: The difference stems from the two omega values. In our core implementation, we use

//        v omega_a
a[i] = 0.45724 * r.tc.powi(2) * KB_A3 / r.pc;
b[i] = 0.07780 * r.tc * KB_A3 / r.pc;
//        ^ omega_b

whereas in the generic cubic implementation, we follow Privat et al., which allows to generate PR and SRK (and possibly others) from a single expression:

// Calculate universal critical constants from universal cubic parameters.
//
// See https://doi.org/10.1016/j.fluid.2012.05.008
fn critical_constants(&self) -> (f64, f64) {
    let (r1, r2) = (-self.d1, -self.d2);
    let eta_c = 1.0
        / (((1.0 - r1) * (1.0 - r2).powi(2)).cbrt()
            + ((1.0 - r2) * (1.0 - r1).powi(2)).cbrt()
            + 1.0);
    let omega_a = (1.0 - eta_c * r1) * (1.0 - eta_c * r2) / (1.0 - eta_c)
        * (2.0 - eta_c * (r1 + r2))
        / (3.0 - eta_c * (1.0 + r1 + r2)).powi(2);
    let omega_b = eta_c / (3.0 - eta_c * (1.0 + r1 + r2));
    (omega_a, omega_b)
}

which yields

omega_a = 0.45723552892138214
omega_b = 0.07779607390388846

If you simply replace the calculated omega's with the values of the core implementation, you'll find that the values match very well. E.g. for the Helmholtz energy:

implemented: -0.1675017780826407
compare: -0.16750177808264066

So as summary - it is probably more sensible to compare against Clapeyron.jl than against our core implementation. ;)

Jun 02 '25 22:06 g-bauer

@cjsisco Thanks for the message. All alpha functions are implemented in this file (might be a good idea to take a very close look at the used constants). A future optimization could be to extract the m (your kappa) and only compute it once.

And yes - also including volume translation is a good idea. A good place would be to put it into the CubicOptions that also contain the alpha function and the mixing rule. Once everything is properly implemented and tested, it'd make sense to take a close look at the structures and see if we can simplify it a bit to make it more readable.

I appreciate all the input.

Jun 02 '25 22:06 g-bauer

@g-bauer According to Jaubert et al. (10.1016/j.fluid.2016.03.012), the Helmholtz energy is unaffected by volume translation (even the temperature-dependent form), so it seems natural to include a volume_shift_c0 and volume_shift_c1 in the default parameter set and default their values =0. Unless there are other volume translation forms than the linear Péneloux form, then I don't foresee the user needing to specify anything in CubicOptions.

If they don't provide volume_shift_* values, then use the defaults. If they do, then apply the simple Péneloux form.

Jun 02 '25 22:06 cjsisco

Ok perfect, I already noticed that I need to get some more hands-on experience before implementing my own EOS fully, so such smaller changes are welcome to me :).

To implement all the above correctly, I propose to restructure the cubic file tree a bit. This would than allow anyone to pick 1 out of every submodule to combine them till their own EoS. I could implement a default mixing rule and a default no-translation method in their specific subtree, and we can use the CubicOptions to make people switch between them. If as @cjsisco states that volume translations does not affect the residual Helmholtz energy (I went to a conference where Jaubert was talking about this as well), we could also couple it less strictly to the cubic parts of FeOs. But I believe that from a user perspective, it is more logical to combine all in one. CubicOptions even if this isn't strictly required.

 src/cubic/
        ├── alpha/
        │   ├── mod.rs
        │   ├── mathisas_copeman.rs
        │   └── ... (other alpha functions split per family)
        ├── mixing/
        │   ├── mod.rs
        │   ├── quadratic.rs
        │   └── ... (other mixing rules)
        ├── translate/
        │   ├── mod.rs
        │   ├── constant_translation.rs
        │   ├── temperature_depending.rs
        │   └── ... (other translate functions per family)
        ├── mod.rs
        └── parameters.rs

Is this file structure ok?

PS: I will still think about implementing the testing in such a way that we do not depend on a magical external value which is computed through Clapeyron. For the general non-mixed PR and (S)RK eos I can probably compute the pressure directly from their original expressions and compare that to the pressure computed from the implemented residual Helmholtz energy. I will come back to this :).

Jun 03 '25 07:06 JelleLagerweij

Sounds good. I have no strong opinion about having each alpha-function, mixing rule and translation method in its own file.

I will still think about implementing the testing in such a way that we do not depend on a magical external value which is computed through Clapeyron.

If you have a good idea how to circumvent this, cool. But it is also common practice to use values from other libraries for comparison. We could, for example, use two libraries such as Clapeyron and thermopack as a first check and then use their values (if they agree) to test our output. That would be totally fine.

Jun 03 '25 13:06 g-bauer

@g-bauer What would it look like to repurpose elements of SAFT for the cubic in feos?

For example, the CPA model of Kontogeorgis is the same association term originally proposed by Walter Chapman for the SAFT model and simply added to the basic repulsive and attractive Helmholtz terms of the cubic. A few years ago, my group hybridized the chain and association terms from the original SAFT into the cubic to make what we call CPCA (cubic+chain+association, 10.1021/acs.iecr.3c02774). This involved a little bit of a change to the repulsive and attractive Helmholtz terms because of the introduction of the chain length.

How much modification do you think would be required to make generic chain and association files that could be used by other models?

Jun 03 '25 17:06 cjsisco

@cjsisco We already have association as separate module here and reuse it across the different implementations. Now, we'd have to take a look at the specific model to see if and how to extend what we already have.

The code is probably not so easy to read if you're not used to generics in Rust. Importantly, one would need to implement the AssociationStrength trait if the current implementation does not suffice.

The chain term on the other hand, we added in the SAFT model's modules because it is a compact implementation. Generalizing contributions often comes at the cost of more complexity because we'd probably have to introduce a new trait. It is always possible to copy the code into the cubic model as a first step, though.

Jun 03 '25 20:06 g-bauer