feat(shielding): Implement Exact STO Integrals and Accuracy Benchmarks

.github/workflows/release.yml

@@ -68,7 +68,7 @@ jobs:
           SOURCE_PATH="target/${{ matrix.target }}/release/$BINARY_NAME"
 
           cp "$SOURCE_PATH" "$ARTIFACT_DIR/$BINARY_NAME"
-          cp LICENSE README.md USAGE.md "$ARTIFACT_DIR/"
+          cp LICENSE README.md USAGE.md BENCHMARKS.md "$ARTIFACT_DIR/"
 
           cd $ARTIFACT_DIR
 

BENCHMARKS.md

@@ -0,0 +1,89 @@
+# Cheq Benchmarks
+
+This benchmark compares the accuracy of **Cheq** (using Slater-Type Orbitals vs Gaussian-Type Orbitals) against **OpenBabel**'s QEq implementation.
+
+The reference values are taken directly from the original QEq paper (Rappe & Goddard, 1991).
+
+## Methodology
+
+- **Cheq STO**: Uses the default Slater-Type Orbital overlap integrals (analytical).
+- **Cheq GTO**: Uses the optional Gaussian-Type Orbital approximation (analytical).
+- **OpenBabel**: Re-implementation of OpenBabel's `qeq.cpp` logic in Rust, using the exact same parameters (`data/qeq.txt`) and GTO integral approximations.
+
+## Results
+
+**Note**: The best result (lowest error) for each row is highlighted in **bold**.
+
+| Molecule             | Paper Result | Cheq STO: Error (Result) | Cheq GTO: Error (Result) | OpenBabel: Error (Result) |
+| -------------------- | ------------ | ------------------------ | ------------------------ | ------------------------- |
+| **Alkali Halides**   |              | **0.0108**               | 0.1902                   | 0.3258                    |
+| LiF (Li)             | _0.7910_     | **0.0117** (0.8027)      | 0.1646 (0.6264)          | 0.1535 (0.9445)           |
+| LiCl (Li)            | _0.9390_     | **0.0190** (0.9580)      | 0.2801 (0.6589)          | 0.3030 (1.2420)           |
+| LiBr (Li)            | _0.9020_     | **0.0194** (0.9214)      | 0.2954 (0.6066)          | 0.3484 (1.2504)           |
+| LiI (Li)             | _0.8410_     | **0.0194** (0.8604)      | 0.3054 (0.5356)          | 0.3937 (1.2347)           |
+| NaF (Na)             | _0.6650_     | **0.0061** (0.6711)      | 0.1166 (0.5484)          | 0.2230 (0.8880)           |
+| NaCl (Na)            | _0.7660_     | **0.0103** (0.7763)      | 0.1896 (0.5764)          | 0.3682 (1.1342)           |
+| NaBr (Na)            | _0.7450_     | **0.0111** (0.7561)      | 0.2007 (0.5443)          | 0.3945 (1.1395)           |
+| NaI (Na)             | _0.7090_     | **0.0112** (0.7202)      | 0.2133 (0.4957)          | 0.4212 (1.1302)           |
+| KF (K)               | _0.6620_     | **0.0050** (0.6670)      | 0.0988 (0.5632)          | 0.2031 (0.8651)           |
+| KCl (K)              | _0.7750_     | **0.0089** (0.7839)      | 0.1682 (0.6068)          | 0.3370 (1.1120)           |
+| KBr (K)              | _0.7680_     | **0.0092** (0.7772)      | 0.1837 (0.5843)          | 0.3647 (1.1327)           |
+| KI (K)               | _0.7540_     | **0.0103** (0.7643)      | 0.2044 (0.5496)          | 0.4043 (1.1583)           |
+| RbF (Rb)             | _0.6530_     | **0.0043** (0.6573)      | 0.0931 (0.5599)          | 0.2039 (0.8569)           |
+| RbCl (Rb)            | _0.7630_     | **0.0078** (0.7708)      | 0.1595 (0.6035)          | 0.3338 (1.0968)           |
+| RbBr (Rb)            | _0.7570_     | **0.0089** (0.7659)      | 0.1742 (0.5828)          | 0.3610 (1.1180)           |
+| RbI (Rb)             | _0.7470_     | **0.0095** (0.7565)      | 0.1960 (0.5510)          | 0.3992 (1.1462)           |
+| **Clusters**         |              | **0.0403**               | 0.2567                   | 0.8057                    |
+| NaCl Monomer (Na)    | _0.7490_     | **0.0273** (0.7763)      | 0.1726 (0.5764)          | 0.3852 (1.1342)           |
+| NaCl Monomer (Cl)    | _-0.7490_    | **0.0273** (-0.7763)     | 0.1726 (-0.5764)         | 0.3852 (-1.1342)          |
+| (NaCl)2 Square (Na)  | _0.8260_     | **0.0187** (0.8073)      | 0.2726 (0.5534)          | 0.8434 (1.6694)           |
+| (NaCl)2 Square (Cl)  | _-0.8260_    | **0.0187** (-0.8073)     | 0.2726 (-0.5534)         | 0.8434 (-1.6694)          |
+| (NaCl)4 Cube (Na)    | _0.8230_     | **0.0750** (0.7480)      | 0.3250 (0.4980)          | 1.1887 (2.0117)           |
+| (NaCl)4 Cube (Cl)    | _-0.8230_    | **0.0750** (-0.7480)     | 0.3250 (-0.4980)         | 1.1887 (-2.0117)          |
+| **Hydrogen**         |              | **0.0329**               | 0.1178                   | 0.6305                    |
+| HF (H)               | _0.4600_     | **0.0193** (0.4407)      | 0.1024 (0.3576)          | 0.7358 (1.1958)           |
+| H2O (H)              | _0.3500_     | **0.0249** (0.3251)      | 0.1404 (0.2096)          | 0.5350 (0.8850)           |
+| NH3 (H)              | _0.2400_     | **0.0211** (0.2189)      | 0.1321 (0.1079)          | 0.5121 (0.7521)           |
+| CH4 (H)              | _0.1500_     | **0.0190** (0.1310)      | 0.1112 (0.0388)          | 0.6887 (0.8387)           |
+| C2H6 (H)             | _0.1600_     | **0.0485** (0.1115)      | 0.1222 (0.0378)          | 0.1553 (0.3153)           |
+| C2H2 (H)             | _0.1300_     | **0.0079** (0.1221)      | 0.0421 (0.0879)          | 0.0323 (0.1623)           |
+| C2H4 (H)             | _0.1500_     | 0.0715 (0.0785)          | 0.1100 (0.0400)          | **0.0610** (0.2110)       |
+| C6H6 (H)             | _0.1000_     | 0.0058 (0.0942)          | 0.0139 (0.0861)          | **0.0025** (0.0975)       |
+| H2CO (O)             | _-0.4300_    | **0.0077** (-0.4377)     | 0.0625 (-0.3675)         | 0.3019 (-0.7319)          |
+| H2CO (C)             | _0.1900_     | **0.0125** (0.2025)      | 0.0998 (0.0902)          | 0.7572 (0.9472)           |
+| H2CO (H)             | _0.1200_     | **0.0024** (0.1176)      | 0.0187 (0.1387)          | 0.2277 (-0.1077)          |
+| H3COH (H on O)       | _0.3600_     | **0.0331** (0.3269)      | 0.1453 (0.2147)          | 0.4772 (0.8372)           |
+| H3COH (O)            | _-0.6600_    | **0.0220** (-0.6380)     | 0.1764 (-0.4836)         | 0.6368 (-1.2968)          |
+| H3COH (C)            | _-0.1500_    | **0.0404** (-0.1096)     | 0.1093 (-0.0407)         | 0.7177 (0.5677)           |
+| H3COH (H gauche)     | _0.1400_     | **0.0126** (0.1274)      | 0.0465 (0.0935)          | 0.2287 (-0.0887)          |
+| H3COH (H trans)      | _0.1800_     | **0.0140** (0.1660)      | 0.0576 (0.1224)          | 0.1108 (0.0692)           |
+| HOC(O)H (O=)         | _-0.4400_    | **0.0479** (-0.3921)     | 0.1633 (-0.2767)         | 0.7590 (-1.1990)          |
+| HOC(O)H (C)          | _0.5600_     | **0.0436** (0.6036)      | 0.3158 (0.2442)          | 2.3439 (2.9039)           |
+| HOC(O)H (H-C)        | _0.1600_     | 0.0774 (0.0826)          | **0.0188** (0.1788)      | 0.9445 (-0.7845)          |
+| HOC(O)H (O-)         | _-0.6500_    | **0.0621** (-0.5879)     | 0.2839 (-0.3661)         | 1.2194 (-1.8694)          |
+| HOC(O)H (H-O)        | _0.3800_     | **0.0861** (0.2939)      | 0.1602 (0.2198)          | 0.5690 (0.9490)           |
+| H3CCN (N)            | _-0.2400_    | **0.0069** (-0.2469)     | 0.0526 (-0.1874)         | 0.6218 (-0.8618)          |
+| H3CCN (C cyano)      | _0.2200_     | **0.0093** (0.2107)      | 0.1742 (0.0458)          | 1.2634 (1.4834)           |
+| H3CCN (C methyl)     | _-0.3700_    | **0.0328** (-0.3372)     | 0.3107 (-0.0593)         | 2.3473 (-2.7173)          |
+| H3CCN (H)            | _0.1300_     | **0.0055** (0.1245)      | 0.0630 (0.0670)          | 0.5685 (0.6985)           |
+| SiH4 (H)             | _0.1300_     | 0.1764 (-0.0464)         | **0.1445** (-0.0145)     | 0.3413 (-0.2113)          |
+| PH3 (H)              | _0.0800_     | **0.0029** (0.0771)      | 0.0422 (0.0378)          | 0.1255 (0.2055)           |
+| HCl (H)              | _0.3200_     | **0.0082** (0.3118)      | 0.0793 (0.2407)          | 0.3701 (0.6901)           |
+| **Polymers**         |              | **0.0703**               | 0.1816                   | 1.4796                    |
+| PE (C)               | _-0.2840_    | **0.0398** (-0.2442)     | 0.1163 (-0.1677)         | 0.4383 (-0.7223)          |
+| PE (H)               | _0.1430_     | **0.0024** (0.1454)      | 0.0615 (0.0815)          | 0.1099 (0.2529)           |
+| PVDF (C-H)           | _-0.0820_    | **0.0195** (-0.0625)     | 0.2179 (0.1359)          | 5.4096 (-5.4916)          |
+| PVDF (H)             | _0.0500_     | **0.0143** (0.0643)      | 0.0414 (0.0914)          | 1.2110 (1.2610)           |
+| PVDF (C-F)           | _0.7760_     | **0.0221** (0.7539)      | 0.4687 (0.3073)          | 4.7492 (5.5252)           |
+| PVDF (F)             | _-0.4150_    | 0.0708 (-0.4858)         | **0.0127** (-0.4023)     | 1.0431 (-1.4581)          |
+| Ala-His-Ala (N-term) | _-0.8100_    | 0.3110 (-0.4990)         | 0.5178 (-0.2922)         | **0.2173** (-1.0273)      |
+| Ala-His-Ala (OH)     | _-0.4500_    | 0.2246 (-0.6746)         | **0.0157** (-0.4343)     | 1.4970 (-1.9470)          |
+| Ala-His-Ala (C=O)    | _-0.4600_    | **0.0602** (-0.5202)     | 0.0830 (-0.3770)         | 0.9053 (-1.3653)          |
+| Ala-His-Ala (N-pep)  | _-0.4600_    | **0.0097** (-0.4503)     | 0.2138 (-0.2462)         | 1.6410 (-2.1010)          |
+| Ala-His-Ala (H-pep)  | _0.2700_     | **0.0164** (0.2536)      | 0.0609 (0.2091)          | 0.6388 (0.9088)           |
+| Ala-His-Ala (C-pep)  | _0.4700_     | **0.0121** (0.4579)      | 0.4069 (0.0631)          | 2.1904 (2.6604)           |
+| Ala-His-Ala (O-pep)  | _-0.5100_    | **0.0348** (-0.5448)     | 0.0869 (-0.4231)         | 0.8353 (-1.3453)          |
+| Ala-His-Ala+ (ND1)   | _-0.5000_    | **0.1599** (-0.3401)     | 0.2686 (-0.2314)         | 1.4420 (-1.9420)          |
+| Ala-His-Ala+ (NE2)   | _-0.5000_    | **0.0178** (-0.4822)     | 0.2304 (-0.2696)         | 1.7639 (-2.2639)          |
+| Ala-His-Ala+ (HD1)   | _0.3300_     | **0.0692** (0.2608)      | 0.1125 (0.2175)          | 0.5334 (0.8634)           |
+| Ala-His-Ala+ (HE2)   | _0.3300_     | **0.1108** (0.2192)      | 0.1728 (0.1572)          | 0.5273 (0.8573)           |

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "cheq"
-version = "0.3.0"
+version = "0.4.0"
 authors = [
     "Tony Kan <tonykan@caltech.edu>",
     "William A. Goddard III <wag@caltech.edu>",
@@ -22,9 +22,10 @@ readme = "README.md"
 [dependencies]
 rayon = "1.11.0"
 libm = "0.2.8"
+sto-ns = "0.1.1"
 faer = "0.23.2"
 thiserror = "2.0.17"
-toml = "0.9.7"
+toml = "0.9.8"
 serde = { version = "1.0.188", features = ["derive"] }
 clap = { version = "4.5.53", features = ["derive"], optional = true }
 prettytable-rs = { version = "0.10.0", optional = true }

README.md

@@ -11,13 +11,28 @@ The core mission of Cheq is to provide a reliable, predictable, and easy-to-inte
 - **Decoupled Architecture**: Agnostic to basic data structures via the flexible `AtomView` trait.
 - **Memory Safe & Fast**: Built in Rust with optimized linear algebra for high performance.
 - **Configurable Parameters**: Includes standard parameters with support for custom TOML-based sets.
-- **Lossy Optimization Knobs**: Optional hard cutoff radius and hydrogen inner iterations (off by default) trade exactness for speed on large systems.
+- **Exact STO Integrals**: Uses exact Slater-Type Orbital integrals for maximum accuracy (default), with optional GTO approximation.
+
+## Benchmarks
+
+Cheq's exact Slater-Type Orbital (STO) implementation delivers significantly higher accuracy compared to traditional Gaussian-Type Orbital (GTO) approximations used in other software like OpenBabel.
+
+| Category           | Cheq (STO) Error | Cheq (GTO) Error | OpenBabel Error |
+| ------------------ | ---------------- | ---------------- | --------------- |
+| Alkali Halides     | **0.0108**       | 0.1902           | 0.3258          |
+| Clusters           | **0.0403**       | 0.2567           | 0.8057          |
+| Hydrogen Molecules | **0.0329**       | 0.1178           | 0.6305          |
+| Polymers           | **0.0703**       | 0.1816           | 1.4796          |
+
+_Values represent the mean absolute error (MAE) in elementary charge units (e) compared to the reference values from Rappé & Goddard (1991)._
+
+See [BENCHMARKS.md](BENCHMARKS.md) for the full detailed breakdown.
 
 ## Getting Started
 
 Cheq ships as both a command-line tool for end users and a library crate for developers.
 
-### For CLI users
+### For CLI Users
 
 Install the binary from crates.io:
 
@@ -35,13 +50,13 @@ cheq water.xyz
 
 The `cheq` command accepts additional options for output formatting, solver tolerances, and custom parameter files. See the [CLI User Manual](USAGE.md) for the full reference and usage patterns.
 
-### For library developers
+### For Library Developers
 
 Add Cheq to your `Cargo.toml`:
 
 ```toml
 [dependencies]
-cheq = "0.3.0"
+cheq = "0.4.0"
 ```
 
 Then, you can use it in your Rust code as follows:
@@ -96,6 +111,7 @@ fn main() {
 ## Tech Stack
 
 - **Core Language**: Rust
+- **STO Integrals**: [`sto-ns`](https://crates.io/crates/sto-ns)
 - **Linear Algebra**: `faer`
 - **Serialization**: `serde` & `toml`
 

USAGE.md

@@ -87,8 +87,9 @@ Cheq exposes a single command with one positional argument and several option gr
 - `--max-iterations <INT>`: Maximum number of solver iterations before aborting. Default: `100`.
 - `--lambda-scale <FLOAT>`: Multiplier applied to the screening length in the Coulomb operator. Default: `0.5`.
 - `--hydrogen-scf <BOOL>`: Enable (true) or disable (false) the hydrogen hardness self-consistency update each iteration. Default: `true`.
-- `--cutoff <FLOAT>`: Hard cutoff radius (Å) for pair interactions; when set, only pairs within the radius are included. Default: unset (uses all pairs).
-- `--hydrogen-inner-iters <INT>`: Extra hydrogen-focused inner iterations before each global solve; `0` disables. Default: `0`.
+- `--basis <sto|gto>`: Basis functions to use for Coulomb integrals. `sto` uses exact Slater-Type Orbitals, `gto` uses approximate Gaussian-Type Orbitals. Default: `sto`.
+- `--damping <auto|fixed|none>`: Damping strategy for the SCF iteration. `auto` adjusts damping based on convergence, `fixed` uses a constant factor, `none` disables damping. Default: `auto`.
+- `--damping-factor <FLOAT>`: Damping factor (0.0-1.0). Used as the fixed value for `fixed` strategy or the initial value for `auto` strategy. Default: `0.4`.
 
 #### General Arguments
 
@@ -181,19 +182,20 @@ cat trajectories/frame_050.xyz | \
 
 ## Argument Reference Table
 
-| CLI Argument (Short) | CLI Argument (Long)      | Value Type | Default      | Description                                       |
-| :------------------- | :----------------------- | :--------- | :----------- | :------------------------------------------------ |
-| _(positional)_       | `INPUT`                  | File Path  | **Required** | XYZ file to process, or `-` for stdin.            |
-| `-o`                 | `--output`               | File Path  | stdout       | Destination for formatted results.                |
-| `-f`                 | `--format`               | Enum       | `pretty`     | Output style: `pretty`, `xyz`, `csv`, or `json`.  |
-| `-p`                 | `--precision`            | Integer    | `6`          | Decimal digits for floating-point values.         |
-| `-P`                 | `--params`               | File Path  | bundled set  | Optional TOML parameter file.                     |
-| `-q`                 | `--total-charge`         | Float      | `0.0`        | Target total charge applied during equilibration. |
-| _(none)_             | `--tolerance`            | Float      | `1e-6`       | RMS change threshold for solver convergence.      |
-| _(none)_             | `--max-iterations`       | Integer    | `100`        | Hard cap on solver iterations.                    |
-| _(none)_             | `--lambda-scale`         | Float      | `0.5`        | Screening length multiplier in Coulomb term.      |
-| _(none)_             | `--hydrogen-scf`         | Bool       | `true`       | Toggle hydrogen hardness self-consistency.        |
-| _(none)_             | `--cutoff`               | Float      | _(unset)_    | Hard cutoff radius (Å) for pair interactions.     |
-| _(none)_             | `--hydrogen-inner-iters` | Integer    | `0`          | Extra hydrogen-focused inner iterations.          |
-| `-h`                 | `--help`                 | Flag       | (N/A)        | Display contextual help and exit.                 |
-| `-V`                 | `--version`              | Flag       | (N/A)        | Show version information and exit.                |
+| CLI Argument (Short) | CLI Argument (Long) | Value Type | Default      | Description                                                   |
+| :------------------- | :------------------ | :--------- | :----------- | :------------------------------------------------------------ |
+| _(positional)_       | `INPUT`             | File Path  | **Required** | XYZ file to process, or `-` for stdin.                        |
+| `-o`                 | `--output`          | File Path  | stdout       | Destination for formatted results.                            |
+| `-f`                 | `--format`          | Enum       | `pretty`     | Output style: `pretty`, `xyz`, `csv`, or `json`.              |
+| `-p`                 | `--precision`       | Integer    | `6`          | Decimal digits for floating-point values.                     |
+| `-P`                 | `--params`          | File Path  | bundled set  | Optional TOML parameter file.                                 |
+| `-q`                 | `--total-charge`    | Float      | `0.0`        | Target total charge applied during equilibration.             |
+| _(none)_             | `--tolerance`       | Float      | `1e-6`       | RMS change threshold for solver convergence.                  |
+| _(none)_             | `--max-iterations`  | Integer    | `100`        | Hard cap on solver iterations.                                |
+| _(none)_             | `--lambda-scale`    | Float      | `0.5`        | Screening length multiplier in Coulomb term.                  |
+| _(none)_             | `--hydrogen-scf`    | Bool       | `true`       | Toggle hydrogen hardness self-consistency.                    |
+| _(none)_             | `--basis`           | Enum       | `sto`        | Coulomb integral basis: `sto` (Slater) or `gto` (Gaussian).   |
+| _(none)_             | `--damping`         | Enum       | `auto`       | SCF damping strategy: `auto`, `fixed`, or `none`.             |
+| _(none)_             | `--damping-factor`  | Float      | `0.4`        | Damping factor (0.0–1.0) for `fixed` or initial `auto` value. |
+| `-h`                 | `--help`            | Flag       | (N/A)        | Display contextual help and exit.                             |
+| `-V`                 | `--version`         | Flag       | (N/A)        | Show version information and exit.                            |

resources/qeq.data.toml

@@ -23,7 +23,7 @@
 "He" = { chi = 9.66, j = 29.84, radius = 0.28, n = 1 }
 
 # --- Row 2 ---
-"Li" = { chi = 3.006, j = 4.772, radius = 1.28, n = 2 }
+"Li" = { chi = 3.006, j = 4.772, radius = 1.557, n = 2 }
 "Be" = { chi = 4.877, j = 8.886, radius = 0.96, n = 2 }
 "B" = { chi = 5.11, j = 9.50, radius = 0.85, n = 2 }
 "C" = { chi = 5.343, j = 10.126, radius = 0.759, n = 2 }