Metadata-Version: 2.4
Name: pysfbox
Version: 1.0.0
Summary: Scheutjens-Fleer lattice self-consistent-field (SF-SCF) theory in pure Python/NumPy - a readable, extensively validated reimplementation with Namics/sfbox-compatible input and output
Author-email: Mark Vis <m.vis@tue.nl>
License: GPL-3.0-or-later
Project-URL: Repository, https://github.com/mark-vis/pysfbox
Keywords: Scheutjens-Fleer,self-consistent-field,SCF,polymer physics,lattice SCF,sfbox,Namics
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Chemistry
Classifier: Topic :: Scientific/Engineering :: Physics
Classifier: License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
License-File: NOTICE
Requires-Dist: numpy
Dynamic: license-file

# PySFBox — Scheutjens–Fleer SCF lattice theory in pure Python/NumPy

Scheutjens–Fleer lattice SCF in pure Python/NumPy that **reads Namics input
files and writes Namics output files**, so existing input files, analysis and
plotting scripts, and course material keep working — no C++ toolchain required.

**Install** (Python ≥ 3.9; the only dependency is NumPy):

```bash
pip install git+https://github.com/mark-vis/pysfbox.git   # adds a `pysfbox` command to your PATH
pysfbox my_input.in                                        # writes my_input.kal / my_input.pro
```

`pysfbox` runs from any folder and writes its `.kal`/`.pro` next to each input
file. Or, without installing — just clone the repo (or drop the `pysfbox/`
folder next to your inputs) and run the module directly:

```bash
pip install numpy
python -m pysfbox my_input.in        # run from the folder that contains pysfbox/
```

GPL-3.0-or-later (see `LICENSE` and `NOTICE`; `pysfbox/sfnewton.py` carries its
own Wageningen copyright notice and governing terms).

Regression tests: `python tests/run_tests.py`. The `tests/` folder doubles as
an example collection — every supported feature has a small, commented `.in`
file there (brushes, adsorption, charged walls, weak polyelectrolytes,
branched molecules, 2D/3D grids, searches).

Keyword reference: `REFERENCE.md` documents every input keyword — allowed
values, defaults, and what each computes — as a lookup manual.

## Validation

PySFBox is validated against the compiled Namics: dozens of regression cases
agree with the C++ oracle at or near machine precision, a 28-case exercise set
reproduces the compiled Namics output case for case, and thermodynamic
observables (free energy, grand potential) agree to the shared convergence
floor. A depletion profile additionally cross-validates against an independent
SF-SCF implementation to ~4e-13, and the physics has passed two full
first-principles reviews with adversarial verification. Where PySFBox departs
from Namics it does so deliberately — implementing the corrected physics where
Namics has confirmed bugs (the fixed-potential Poisson factor 2, the
refined-lattice bond length, the weak-charge free-energy/chemical-potential
index terms) — each noted in the relevant section below.

## What it supports (v1)

- **Input**: the Namics `key : name : parameter : value` format, `//`
  comments, `start` blocks (settings accumulate, like Namics), `alias`
  (`#N#` substitution), and `var` parameter scans (`scan`/`step`/`end_value`),
  including alias scans (`scan : N-value`).
- **Lattice**: one, **two, or three gradient directions** (`gradients : 1/2/3`),
  as in Namics. 1-gradient: planar, cylindrical, spherical, `simple_cubic`
  (λ=1/6) or `hexagonal` (λ=1/4), lattice refinement `FJC_choices : 5, 7, …`
  (spherical/cylindrical), `offset_first_layer` shifts the radial origin.
  2-gradient (`simple_cubic`, `FJC 3`, neutral/linear): flat (x,y) and
  cylindrical (r,z); 3-gradient: flat (x,y,z). N-D uses `n_layers_x/y/z` and
  per-axis bounds `lowerbound_x`, `upperbound_y`, … (`mirror` / `surface` /
  `periodic`); one unified finite-volume lattice handles every geometry.
- **mon**: `freedom` free / frozen / pinned, `frozen_range` / `pinned_range`
  (1-gradient `1;3`, `lowerbound`, `upperbound`, `firstlayer`, `lastlayer`;
  2D/3D box grammar `xlo,ylo[,zlo];xhi,yhi[,zhi]`), `chi_X` interaction
  parameters.
- **mol**: linear multi-block compositions `(X)1(A)200(G)1` AND **branched
  architectures** — bracket side chains `(C)18[(O)1](O)1...` (nestable;
  several `[..][..]` on one node), symmetric dendrimers/stars
  `@dend(X,(A)200(G)1,5)`, combs `@comb((B)10;B,(A)25,(B)9,100;(B)1)`, and
  group-repeat sugar `(C)12((O)1(C)2)5(O)1`, on any 1-gradient geometry;
  `freedom` free (`phibulk`), restricted (`n` or `theta`), solvent.
- **search / super-iteration**: `var : mol-X : search : theta` (or `n` /
  `phibulk`) paired with a target — `var : sys-NN : grand_potential : 0`,
  `... : free_energy : V`, or `var : mol-Y : mu / theta / n / phibulk : V` —
  runs an outer scalar root-find (regula falsi) that adjusts the searched
  quantity until the target observable hits its value, each step a full
  warm-started SCF solve. (The `constraint:delta` / `Laplace_pressure`
  membrane-balance search of the Namics examples is not supported.)
- **Electrostatics**: `mon : valence` (fixed charges, strong electrolytes),
  `mon : epsilon` (relative permittivity, default 80), `mon : e.psi0/kT`
  (fixed surface potential on a frozen electrode layer), `mol : freedom :
  neutralizer` (auto-adjusts its bulk fraction for an electroneutral bulk),
  `lat : bondlength` (meters). The dimensionless potential psi joins the
  iteration; profiles via `pro : sys : X : psi` and `q`. Works on **planar,
  cylindrical and spherical** lattices (the curved Poisson uses
  face-area-weighted fluxes; `tests/charged_sphere.in`,
  `tests/charged_cylinder.in` match the compiled Namics at the convergence
  floor, and the spherical potential reproduces analytic Debye–Hückel
  screening). Validated against the compiled Namics at machine precision
  (`tests/polyelectrolyte_wall.in`); fixed-potential systems
  (`tests/edl_fixed_psi.in`) validate against Debye theory instead, because
  PySFBox deliberately corrects a factor-2 bug in the Namics fixed-potential
  Poisson branch (Debye lengths come out √2 too short there). Refined curved
  charged lattices (`FJC_choices > 3`) are supported and refinement-consistent:
  the Poisson equation keeps the physical bond length, where Namics' Debye
  lengths scale spuriously with the refinement — so expect a deviation from an
  unpatched Namics on that combination. Curved fixed-potential electrodes
  raise a clear message.
- **Weak (multistate) charges**: `state : AH : mon : A` + `state : AH :
  valence : 0` attach annealed internal states to a segment;
  `reaction : weak : equation : 1(AH) + 1(H2O) = 1(AM) + 1(H3O)` +
  `reaction : weak : pK : 7` set their equilibria; `state : H3O :
  alphabulk : 1e-7` anchors the pH (site-fraction based, as in Namics:
  "pH" = −log10 alphabulk_H3O). The bulk ionisation is solved to machine
  precision (Henderson–Hasselbalch exactly, in site fractions with the
  water activity included); locally every state equilibrates with the
  fields — charge regulation comes out automatically. Titration scans:
  `var : state-H3O : scan : alphabulk` with `scale : exponential`
  (`steps` = steps per decade, as in Namics). Per-state outputs: kal
  `theta_AM`, `alphabulk_AM`, `phibulk_AM`, `theta_exc_AM` (on the parent
  mon), `mol : X : mu-AM`; pro `phi-AM`, `alpha-AM`, `u-AM`. States
  inherit the mon's chi (state-level chi overrides); the mon-level
  valence is ignored once states exist. Validated against the compiled
  Namics at the convergence floor on the Namics `polE.in` example and a
  73-step titration (`tests/weak_*.in`); PySFBox implements corrected
  physics in a few spots where Namics has known bugs, so F/mu comparisons
  with an unpatched Namics differ whenever state–state chi is nonzero.
- **sys**: `initial_guess` — `polymer_adsorption` / `membrane` / `micelle`
  analytic starting potentials (the Namics formulas, e.g. u = −λ·χ at
  solid-adjacent sites), honored on the first calculation only, exactly
  like Namics; `previous_result`/`none` accepted; guess files and
  `membrane_torus` raise a clear message.
- **Solver**: the Namics pseudohessian quasi-Newton (a licensed Python
  translation of `sfnewton.cpp`) as primary, with Anderson-accelerated
  Picard and two rescue stages as fallbacks (reads `newton : ... :
  tolerance / iterationlimit / deltamax / m / method`; the Namics `method`
  choices `pseudohessian` / `DIIS` / `Picard` / `LBFGS` / `hessian` map onto
  this cascade — `hessian` is the Namics full-Hessian mode, recomputing a
  finite-difference Hessian every iteration). Propagators are renormalised
  with log-domain bookkeeping, so long or strongly adsorbing chains cannot
  overflow. Large 2D/3D grids (beyond the dense pseudohessian) are handled
  by the Anderson stage.
- **Output**: `.kal` (header + one row per calculation/var step, `%.16e`)
  and `.pro` (x at half-integer layers + `%.20g` columns), with Namics file
  numbering. Unknown output properties print `NiN` plus a warning — exactly
  like Namics. Implemented kal properties: `sys` grand_potential,
  free_energy (and `free_energy (po)`), iterations, residual; `lat`
  n_layers / volume; `mol` theta, theta_exc, phibulk, Mu/mu (chemical
  potential), mu-STATE, n, N, chainlength, GN, phiM, phiMax, phiMin,
  `<alias>-value`; `mon` theta, theta_exc, phibulk, chi_X,
  1st_M_phi_z, 2nd_M_phi_z, RMS, fluctuations, and per-state
  alphabulk_S / valence_S / phibulk_S / theta_S / theta_exc_S. Profiles
  (`pro`): `mol : X : phi`, `mon : X : phi`, `mon : X : u`,
  `mon : X : phi-S / alpha-S / u-S`, `sys : X : alpha / psi / q`.

## Beyond this release

Full-Namics features that this release does not implement raise a clear
`NotImplementedError` naming the feature, rather than silently producing
wrong numbers — so an input either runs correctly or tells you exactly what
it needs. Currently outside scope: ring and asymmetric-dendrimer
architectures, the `constraint:delta` / `Laplace_pressure` membrane-balance
search, initial-guess files, and mesodyn/cleng/teng calculations.

The **one-gradient** path (planar / cylindrical / spherical, including
`FJC_choices` refinement) covers the full feature set above. The **two- and
three-gradient** lattices — 2D flat, 2D cylindrical (r,z), 3D flat, as in
Namics — currently cover neutral linear chains on `simple_cubic`; charged and
branched systems there need `gradients : 1`. Curved fixed-potential
electrodes raise. For weak charges, unbalanced reactions and states on frozen
monomers are rejected (the latter also in Namics).

## Extending (where things live)

The module layout mirrors the Namics classes, so the C++ original is the
roadmap for extensions:

| file | contents | Namics counterpart |
|---|---|---|
| `inputreader.py` | input parsing, alias, var | `input.cpp` |
| `lattice.py` | 1-gradient geometry, ghost layers, site averages, moments | `lattice.cpp`, `LGrad1.cpp` |
| `latticend.py` | 2- and 3-gradient finite-volume lattice | `LGrad2.cpp`, `LGrad3.cpp` |
| `model.py` | segments, molecules (linear + branched trees), renormalised propagators | `segment.cpp`, `molecule.cpp`, `mol_linear.cpp`, `mol_branched.cpp`, `mol_comb.cpp`, `mol_dendrimer.cpp` |
| `reactions.py` | internal states, reactions, bulk ionisation solve | `state.cpp`, `reaction.cpp` |
| `system.py` | residual, masks, observables, solver cascade, initial guesses | `system.cpp`, `solve_scf.cpp` |
| `sfnewton.py` | pseudohessian quasi-Newton (licensed translation) | `sfnewton.cpp` |
| `output.py`, `runner.py` | .kal/.pro writers, start/var loops, search | `output.cpp`, `namics.cpp`, `variate.cpp` |

Examples: new output properties = extend the tables in `System.get_value` /
`System.get_profile`; new chain architectures = generalise the tree walk in
`Molecule` (see `mol_branched.cpp`).

## Convergence envelope

The primary solver IS the Namics quasi-Newton (pseudohessian): `sfnewton.py`
is a faithful, licensed Python translation of `sfnewton.cpp`, run from the
warm-start guess with Namics-like iteration counts. Anderson-accelerated
Picard (with a damped-Picard warmup on cold starts) is the fallback, plus
two rescue stages (extended-budget pseudohessian, one-shot full-Hessian
anchor) for the stiffest cases. This converges the entire shipped
validation set, including the extreme corner — very strong adsorption of very long
chains (the shipped `homopolymer_adsorption.in`, N = 10000, chi_Si = -6)
converges via the extended rescue and matches the compiled Namics exactly.
Genuine non-convergence still fails *gracefully* (a clear "no convergence"
message with actionable advice, never a wrong number). `var` scans
warm-start each step from the previous solution — ramping a parameter
(e.g. a co-solvent) is the intended way into hard regimes.

## Speed

The default path is deliberately pure NumPy — one dependency, every equation
inspectable — running the same pseudohessian quasi-Newton as Namics, with
Namics-like iteration counts. On a 28-case benchmark set the compiled Namics
C++ takes ~18 s where PySFBox takes ~2 minutes, and every one of those 28 cases
matches the compiled Namics output: typical single problems solve in seconds,
the stiffest shipped case (`homopolymer_adsorption.in`, N = 10000, chi_Si = −6)
in a few minutes. For throughput-bound parameter sweeps, the same input files
run unchanged on the compiled Namics — compatibility works in both directions.

## Caveats for quantitative use

- The primary solver is a translation of Namics' quasi-Newton; iteration
  counts are Namics-like but not identical (fallback stages differ), and
  converged solutions agree to the shared convergence floor.
- The grand potential and free energy implement the standard uncharged SF
  expressions; cross-validated against compiled Namics on a two-brush case
  (free energy 6e-9, grand potential 5e-8 relative — the shared convergence
  floor; see `tests/two_brushes_quick.in`). For research use, spot-check
  your own geometry the same way. Profiles, theta, theta_exc, and moments
  follow the Namics definitions exactly (`LGrad1::Moment`, `WeightedSum`).
- Default tolerance 1e-7, iterationlimit 1000, deltamax 0.1 (as in Namics).

## Acknowledgements & lineage

PySFBox is a reimplementation of the Scheutjens–Fleer self-consistent-field
lattice theory as realised in **sfbox/Namics** by F.A.M. Leermakers and the
Physical Chemistry group at Wageningen University — the input/output formats,
property names, and numerical conventions all follow Namics so that existing
material keeps working. `pysfbox/sfnewton.py` is a Python translation of
Namics' `sfnewton.cpp` and carries the original Wageningen University
copyright and reproduction notice in its file header, which governs that
file. Many thanks to Frans Leermakers for sfbox and Namics — the theory and
the code that PySFBox builds on — and for discussions along the way.

PySFBox itself is licensed under the GNU General Public License v3.0 or later
(GPL-3.0-or-later) — see `LICENSE` and `NOTICE`.
