The matrix appears to be singular

[Felip-552]

When trying to use the package to simulate hydrogen adsorption with saftvrqmie, the codes gives the error: "The matrix appears to be singular"

from feos.si import *  
from feos.dft import *
from feos.saftvrqmie import *
import matplotlib.pyplot as plt

parameters= SaftVRQMieParameters.from_json(
    substances=['hydrogen'], 
    pure_path='.../parameters/saftvrqmie/hammer2023.json')

func = HelmholtzEnergyFunctional.saftvrqmie(parameters, FMTVersion.WhiteBear)

temp = 77*KELVIN
psize= 200*ANGSTROM

solver = DFTSolver().picard_iteration(tol=1e-5, beta=0.05).anderson_mixing()

potential = ExternalPotential.LJ93(3.0, 100.0, 0.08)

bulk = State(func, temp, pressure=30*BAR)
porex = Pore1D(Geometry.Cartesian, psize, potential).initialize(bulk).solve(solver)

plt.plot(porex.z/(ANGSTROM), (porex.density*(METER**3/MOL)/1000).T)

The error is:

Cell In [20], line 9
      6 potential = ExternalPotential.LJ93(3.0, 100.0, 0.08)
      8 bulk = State(func, temp, pressure=30*BAR)
----> 9 porex = Pore1D(Geometry.Cartesian, psize, potential).initialize(bulk).solve(solver)
     11 plt.plot(porex.z/(ANGSTROM), (porex.density*(METER**3/MOL)/1000).T)
     13 plt.xlabel(r"$z$ / A")

RuntimeError: The matrix appears to be singular.

When I try with a bigger pore (300 A) it runs with no problems, and if I try an even smaller pore the error changes to:

Cell In [21], line 9
      6 potential = ExternalPotential.LJ93(3.0, 100.0, 0.08)
      8 bulk = State(func, temp, pressure=30*BAR)
----> 9 porex = Pore1D(Geometry.Cartesian, psize, potential).initialize(bulk).solve(solver)
     11 plt.plot(porex.z/(ANGSTROM), (porex.density*(METER**3/MOL)/1000).T)
     13 plt.xlabel(r"$z$ / A")

RuntimeError: `Picard Iteration` encountered illegal values during the iteration.

I also tried the same configuration using ´pcsaft´ and had no issues for any pore sizes. Could you please look into this errors?

[prehner]

Thank you for your interest in FeOs. There is no guaranteed convergence for the DFT calculations but you can adjust the solver settings. In your case, the calculation converges with, e.g.,

solver = DFTSolver().picard_iteration(max_iter=100, damping_coefficient=0.02).newton()

As you built the Python package from source anyways, you can get the newton solver by pulling the current main branch.
In the update, the beta parameter of the Picard iteration and Anderson mixing were renamed to damping_coefficient.

It appears as though some safeguards against numerical issues that are used in the PC-SAFT model are not yet incorporated in the SAFT-VRQ-Mie model. We will have a look into that. Thank you for bringing it to our attention!

[Felip-552]

Thanks for your kind response, sadly I have tried installing the new update with pip3 install git+https://github.com/feos-org/feos@main on a new environment, and also tried installing the branch [saftvrqmie_refinement] with pip3 install git+https://github.com/feos-org/feos@saftvrqmie_refinement but the installation fails to compile.

I tried the same with an older commit (8f0ebccb5fcb458c9b81b4899f8bd9741405fc2b) and had no issues installing it, but this commit does not have saftvrqmie included.

The current version that I'm using is directly from a fork of morteham, but this one had the convergence issue. Thanks again for your response and work, It's really great.

[prehner]

The installation works on my machine and in the GitHub CI. If you provide details on what goes wrong, we can see if there is a solution. And thank you for the kind words.

[g-bauer]

You might have to upgrade your rust compiler to the latest version. If the installation via pip install git+... is not working, cloning the repository and building via pip install . (from within the repository root) or maturin build --release + pip3 install target/wheels/wheel-name.whl should work.

[prehner]

Closing this issue. The underlying problem is tracked in #87

[g-bauer]

@Felip-552 did an update of the rust compiler resolve the installation issue?

[Felip-552]

@g-bauer Hi, yes it was the rust compiler, now I have no issues installing. Thanks for your help. I was wondering if there is a way to know the molar fraction of the adsorbed phase? so one can determine the selectivity of surfaces when adsorbing mixtures.

[prehner]

You can get the number of particles for each component using profile.moles for a single pore or isotherm.adsorption for an entire adsorption isotherm.