Repository navigation
Restructuring plan
This is the plan for some restructuring of the EDRIXS code.
For the backend interface intended for the PETSc implementation, see the Implementation guide.
For a detailed description of the reference SciPy implementation, see the SciPy backend implementation note.
Contains what the user will usually call to define a particular problem or model. It has negligible computational overhead and provides infrastructure to describe the problem, with many calls to the rest of the EDRIXS Python interface. The functions involved will be
-
model_<>where<>is currently1v1c,2v1c,siam, etc. --- setup problem
We had initially used the term setup here, but that is arguably not a good choice given the widespread use of setup and setup.py for package installation.
We need to build in more flexibility here regarding how the basis is defined. We should implement a class that enforces either
- a
subspaces= a list of tuples [(orb, noccu), (orb, noccu), (orb, noccu)] or -
full= np.array([int32])
This can be enforced using a dataclass.
Contains what the user will usually call to solve a particular problem. It should contain code that is independent of the method used to compute the spectra. It is made up of functions
- ``build_op` --- build a many-body Fock-space operator from two-fermion and/or four-fermion terms
-
get_ops--- generate set of initial and intermediate state Hamiltonians and photon operators -
ed--- compute low-energy eigenenergies and eigenvectors. -
xas--- generate final spectra -
rixs--- generate final spectra
get_ops, ed, xas, and rixs orchestrate calls into specific solvers. build_op will be called by get_ops and belongs in solvers.py since it is highly useful for the user. We plan to implement build_op, get_ops, ed, xas, rixs as thin wrapper functions to call into the backend specific functions.
Where <>= scipy, petsc etc. A series of different modules implementing the underlying methods for different backends. The current dense compatibility path delegates to the SciPy implementation, while the PETSc backend is currently an interface stub. Backend implementations provide the following functions.
- ``build_op_<>` --- Build a many-body Fock-space operator from basis objects and orbital-space two-fermion and/or four-fermion terms.
-
ed_<>--- obtain low-energy eigenenergies and eigenvectors. xas_<>rixs_<>
All of these will call into some underlying backend, which will depend on a few parameters. For uniformity, we should pass these as a dictionary backend_kws. solvers.py remains method-independent and dispatches to backend entry points. Method-specific numerical routines live in the backend implementation and shared numerical helper modules such as krylov.py.
The approach previously in manybody_operator_csr.py is backend-specific, so this algorithm appears in scipy_backend/scipy_backend.py as build_op_scipy(). This is because the operator building process requires careful optimization that will probably be backend-dependent.
The class FockBasis belongs in fock_basis.py. manybody_operator_csr.py is removed entirely, while the existing Fock-basis utilities remain in fock_basis.py.
This should contain stuff that is backend independent and typically not called by the user, including backend loading/inference, model helpers, legacy workflow helpers, and functions like _expand_broadening.
Define the problem by calling model_<> where <> is currently 1v1c, 2v1c, siam
emat_i, umat_i, basis_i, emat_n, umat_n, basis_n, trans_mat = model_1v1c(...)
where (...) specifies the problem at hand independent of which backend (scipy, petsc, etc.) is going to be called.
emat (umat) is the two(four)-fermion Hamiltonian in the single-particle basis. Suffixes _i, _n denote the initial and intermediate states. basis describes (but does not build) the basis. trans_mat contains the dipole or quadrupole transition components from the core to valence states, again in the single-particle basis.
hmat_i, hmat_n, trans_ops = get_ops(emat_i, umat_i, basis_i, emat_n, umat_n, basis_n, trans_mat,
*, backend=`scipy`)
hmat_i (hmat_n) is a many-body Hamiltonian for the initial (intermediate) states. trans_ops contains the many-body absorption operators corresponding to the components in trans_mat.
Also within solvers, we should have a wrapper function
op = build_op(emat, umat, lb, rb=None, backend=backend)
where lb and rb are two basis objects. If rb is not included, it will set rb=lb.
get_ops constructs the initial and intermediate Hamiltonians and each component of the photon transition operator through build_op. Under the hood, build_op would call, for example, build_op_scipy from scipy_backend/scipy_backend.py.
From here, the matrices will be owned by the selected backend. The same backend may be specified explicitly for the other solvers. If it is omitted, the current interface infers it from the backend-owned operators.
eval_all, evec_all = ed(hmat_i, num_evals=1, backend=backend, , backend_kws = {"tol": 1e-10, "maxiter": 200})
Obtain low-energy eigenvalues and eigenvectors. num_evals is the number of states. Rather than managing exact details of the diagonalization algorithm, I suggest exposing all the keyword arguments via backend_kws which just gets passed to the underlying backend. Within the specific backends, we can either define our own defaults (which would be overwritten by backend_kws if these are passed in) or leave it to the package defaults.
Ixas = xas(eval_i, evec_i, hmat_n, trans_op, ominc, *, gamma_c=0.1, thin=1.0, phi=0.0, pol_type=None, temperature=1.0, scatter_axis=None, backend=None, backend_kws=None)
This should take a 1d array of ominc and a sequence of pol_type values with a simple for loop and return a 2d array with shape (len(ominc), len(pol_type)). In performance in sensitive situations, higher-level code can parallelize multiple calls into this function. This level of the code we deliberately leave out to minimize complexity and maximize flexibility.
nkrylis somewhat method-dependent and should be passed throughbackend_kws`.
Irixs = rixs(eval_i, evec_i, hmat_i, hmat_n, trans_op, ominc, eloss, *, gamma_c=0.1, gamma_f=0.01, thin=1.0, thout=1.0, phi=0.0, pol_type=None, temperature=1.0, scatter_axis=None, skip_gs=False, return_poles=False, backend=None, backend_kws=None):
This should take a 1d array of incident energies and polarizations equivalent to what was done for xas.
Keep skip_gs since it has already been released.
There are a few different things we want to achieve under the broad umbrella of "testing". Here, try to break down this with clarity of what we want to achieve and where it should live in the repo.
It's been reasonably suggested that the tests folder lives in the root of the repository rather than under edrixs.
We should lean towards having simple fast unit tests for as functions as appropriate in the form test_<module_name> within the tests folder. These, of course, should execute in the CI.
These run one small problem against each other within using different backends/methods in order to check the correctness of calculations. We plan for these to live in tests/solver_consistency_tests/<filename>.py. These should be small enough to execute in CI. We can and should verify correctness without benchmarking against the Fortran. And actually, since the Fortran has a lot of lines and proven very difficult to inspect due to the complexity of the code and most people's limited knowledge of Fortran, it's perhaps a particularly good candidate for bugs.
In the current main branch, the only actual CI tests happening involve executing examples, which are primarily to inform users of how to use the software. Right now these live examples examples
It seems sensible to keep running some of the computationally inexpensive examples as tests in this way?
We probably want to run some tests of speed comparing to the old Fortran implementation to the new code. While there are many sensible options for how one might do this.
Current plan is that these exist in benchmark_checks in the root folder
A suitable methods for installing via conda has been identified:
conda install --channel conda-forge --strict-channel-priority "petsc=*=cuda12_complex_*" "slepc=*=cuda12_complex_*" petsc4py slepc4py mpi4py
-
Note this has cuda12. cuda13 is coming.
-
We need to decide if we make petsc slepc strict requirements to install or not. Not doing so would make quite a lot of sense, as it would open some native functionality on windows.
Work on the petsc_edrixs github branch. Raising PRs to petsc_edrixs as we have been doing.
Considering making a pre release
v0.1.2 _> v1.0.a1, v1.0.b1 etc.
Do not invest heavily in restructuring functions that we anticipate being eventually deleted. Just mark them with a depreciation notice.
At least in the short term, parallelization with from concurrent.futures import ProcessPoolExecutor is removed so that the decision of how to do the parallelization can be deferred until later in our process.
Functions such as rixs_scipy should remain serial and use simple loops (for now). Parallelism, if needed, can be implemented at a higher level by making independent calls to rixs_scipy, for example with a single incident energy per call. We therefore do not expose threading or multiprocessing controls in rixs_scipy itself.
Some notes for tests could be written here...
Any particular problem you see for typing?
A minor note is we should stop EDRIXS spewing printouts to the command line whether you request it or not.
We need to fix potentially confusing name shadowing were the naive command edrixs.build_opers Still exists from the old methods. edrixs.build_ops But is too similar to the new command.
edrixs.cb_op
also presumes numpy dense methods. It would be better if it covered other matrices.
The current F2PY extension does not transfer solver data. It only passes three MPI-derived integers (comm, rank, size) into pyapi.f90; the Fortran routines then read the same files the standalone programs use. Python writes config.in, hopping_.in, coulomb_.in, fock_*.in, and transition matrices; after execution it reads files such as eigvals.dat, denmat.dat, and pole files. Moreover, the repository already builds independent MPI executables:
- ed.x
- xas.x
- rixs.x
- opavg.x
The F2PY routines in src/pyapi.f90 (line 1) largely duplicate the corresponding *_main.f90 programs. The only reason normal pip install edrixs compiles Fortran is that setup.py (line 66) unconditionally enables EDRIXS_PY_INTERFACE, which causes CMake to generate and build the NumPy/F2PY extension. I would recommend:
- Make the base Python package pure Python. Do not require CMake, a Fortran compiler, MPI development headers, LAPACK/ARPACK, or NumPy development headers just to install and use the SciPy/Python functionality.
- Treat the Fortran solver as an optional external backend. Let the existing Python file-writing workflow invoke ed.x, xas.x, or rixs.x with subprocess, normally through a configurable launcher such as mpiexec, srun, or jsrun.
- Distribute the Fortran side separately:
- a separately installable system/Conda package containing the executables, and/or
- documented CMake build instructions for users who want to compile it. The Python layer should accept an executable path or discover it on PATH.
- Deprecate and eventually remove the F2PY module and edrixs/scripts.py (line 1) wrappers that import it.
This has a few clear benefits:
- Pure-Python users get a straightforward install.
- It removes the fragile Python/NumPy/F2PY/Fortran/MPI ABI coupling.
- It avoids passing an MPI communicator handle across language and MPI-library boundaries.
- The Python package becomes easier to wheel and test.
- It removes duplicated Fortran entry-point logic. The main caveat is that this does not eliminate compilation for people who use the Fortran solver—it relocates it to the optional solver distribution, which is the appropriate boundary. Shipping portable prebuilt MPI binaries on PyPI is difficult because MPI runtime and scheduler environments vary, so a separately compiled executable or Conda package is likely the practical delivery mechanism. One important API choice: do not preserve the current comm-based calling convention when moving to subprocesses. A Python process launched under MPI cannot reliably “nest” another mpiexec. Instead, expose a runner configuration (executable, launcher command, process count, working directory), and launch the Fortran job as its own MPI application. I think this can be done with the mpi5py command spawn? The already file-based design makes that natural. I’d frame this as a packaging/backend separation, not a rewrite of the solver protocol: keep the established file format initially, add an external-executable runner, make F2PY optional/off by default, then remove it after a compatibility period.
We need things to be generally well documented, including updating all examples. Adding more description of the numerical methods and possible sources of finite errors will be important.
- Consider performance from big problems that can be memory bound. Should matrices be recomputed during the procedure. Consider matrix free implementation.
- Debshika has implemented shared Krylov space method