diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index f18eb85..538c30b 100644 --- a/.github/workflows/release.yml +++ b/.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 diff --git a/BENCHMARKS.md b/BENCHMARKS.md new file mode 100644 index 0000000..8982a1d --- /dev/null +++ b/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) | diff --git a/Cargo.toml b/Cargo.toml index 6c8d432..e8b4d9f 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "cheq" -version = "0.3.0" +version = "0.4.0" authors = [ "Tony Kan ", "William A. Goddard III ", @@ -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 } diff --git a/README.md b/README.md index d435acd..d4d982d 100644 --- a/README.md +++ b/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` diff --git a/USAGE.md b/USAGE.md index aa89f25..b867684 100644 --- a/USAGE.md +++ b/USAGE.md @@ -87,8 +87,9 @@ Cheq exposes a single command with one positional argument and several option gr - `--max-iterations `: Maximum number of solver iterations before aborting. Default: `100`. - `--lambda-scale `: Multiplier applied to the screening length in the Coulomb operator. Default: `0.5`. - `--hydrogen-scf `: Enable (true) or disable (false) the hydrogen hardness self-consistency update each iteration. Default: `true`. -- `--cutoff `: Hard cutoff radius (Å) for pair interactions; when set, only pairs within the radius are included. Default: unset (uses all pairs). -- `--hydrogen-inner-iters `: Extra hydrogen-focused inner iterations before each global solve; `0` disables. Default: `0`. +- `--basis `: Basis functions to use for Coulomb integrals. `sto` uses exact Slater-Type Orbitals, `gto` uses approximate Gaussian-Type Orbitals. Default: `sto`. +- `--damping `: Damping strategy for the SCF iteration. `auto` adjusts damping based on convergence, `fixed` uses a constant factor, `none` disables damping. Default: `auto`. +- `--damping-factor `: 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. | diff --git a/resources/qeq.data.toml b/resources/qeq.data.toml index 7f99535..2b76e87 100644 --- a/resources/qeq.data.toml +++ b/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 } diff --git a/src/bin/modules/app.rs b/src/bin/modules/app.rs index 8f4414c..09e7901 100644 --- a/src/bin/modules/app.rs +++ b/src/bin/modules/app.rs @@ -1,7 +1,7 @@ -use super::cli::Cli; +use super::cli::{Cli, CliBasisType, CliDampingStrategy}; use super::error::CliError; use super::io; -use cheq::{QEqSolver, SolverOptions, get_default_parameters}; +use cheq::{BasisType, DampingStrategy, QEqSolver, SolverOptions, get_default_parameters}; use indicatif::{ProgressBar, ProgressStyle}; use std::fs; @@ -16,13 +16,26 @@ pub fn run(args: Cli) -> Result<(), CliError> { get_default_parameters().clone() }; + let basis_type = match args.solver.basis { + CliBasisType::Sto => BasisType::Sto, + CliBasisType::Gto => BasisType::Gto, + }; + + let damping = match args.solver.damping { + CliDampingStrategy::Auto => DampingStrategy::Auto { + initial: args.solver.damping_factor, + }, + CliDampingStrategy::Fixed => DampingStrategy::Fixed(args.solver.damping_factor), + CliDampingStrategy::None => DampingStrategy::None, + }; + let solver_options = SolverOptions { tolerance: args.solver.tolerance, max_iterations: args.solver.max_iterations, lambda_scale: args.solver.lambda_scale, hydrogen_scf: args.solver.hydrogen_scf, - cutoff_radius: args.solver.cutoff, - hydrogen_inner_iters: args.solver.hydrogen_inner_iters, + basis_type, + damping, }; let solver = QEqSolver::new(¶ms).with_options(solver_options); @@ -63,7 +76,10 @@ pub fn run(args: Cli) -> Result<(), CliError> { #[cfg(test)] mod tests { - use super::super::cli::{CalculationOptions, OutputFormat, OutputOptions, SolverOptions}; + use super::super::cli::{ + CalculationOptions, CliBasisType, CliDampingStrategy, OutputFormat, OutputOptions, + SolverOptions, + }; use super::super::error::CliError; use super::*; use std::fs; @@ -80,8 +96,9 @@ mod tests { max_iterations: 50, lambda_scale: 0.5, hydrogen_scf: true, - cutoff: None, - hydrogen_inner_iters: 0, + basis: CliBasisType::Sto, + damping: CliDampingStrategy::Auto, + damping_factor: 0.4, } } diff --git a/src/bin/modules/cli.rs b/src/bin/modules/cli.rs index 61a52c7..24a5f40 100644 --- a/src/bin/modules/cli.rs +++ b/src/bin/modules/cli.rs @@ -99,13 +99,17 @@ pub struct SolverOptions { #[arg(long, default_value_t = true, action = clap::ArgAction::Set)] pub hydrogen_scf: bool, - /// Hard cutoff radius (Å) for pair interactions. When unset, all pairs are included. - #[arg(long)] - pub cutoff: Option, + /// Basis functions to use for Coulomb integrals. + #[arg(long, value_enum, default_value_t = CliBasisType::Sto)] + pub basis: CliBasisType, - /// Extra hydrogen-focused inner iterations before each global solve (0 disables). - #[arg(long, default_value_t = 0)] - pub hydrogen_inner_iters: u32, + /// Damping strategy for the SCF iteration. + #[arg(long, value_enum, default_value_t = CliDampingStrategy::Auto)] + pub damping: CliDampingStrategy, + + /// Damping factor (0.0-1.0). Used as fixed value or initial value for auto damping. + #[arg(long, default_value_t = 0.4)] + pub damping_factor: f64, } /// Output format for the calculation results. @@ -120,3 +124,23 @@ pub enum OutputFormat { /// JSON object containing atoms array and metadata. Json, } + +/// Basis function types for Coulomb integral calculations. +#[derive(Clone, ValueEnum)] +pub enum CliBasisType { + /// Slater-Type Orbitals (exact). + Sto, + /// Gaussian-Type Orbitals (approximate). + Gto, +} + +/// Damping strategies for the SCF iteration. +#[derive(Clone, ValueEnum)] +pub enum CliDampingStrategy { + /// Auto-adjust damping based on convergence. + Auto, + /// Fixed damping factor. + Fixed, + /// No damping. + None, +} diff --git a/src/lib.rs b/src/lib.rs index 93a77bb..09dd622 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -60,8 +60,8 @@ // Internal modules are private to enforce encapsulation. mod error; -mod math; mod params; +mod shielding; mod solver; mod types; @@ -71,7 +71,7 @@ pub use self::error::CheqError; pub use self::types::{Atom, AtomView, CalculationResult}; -pub use self::solver::{QEqSolver, SolverOptions}; +pub use self::solver::{BasisType, DampingStrategy, QEqSolver, SolverOptions}; pub use self::params::{ElementData, Parameters}; diff --git a/src/math/mod.rs b/src/math/mod.rs deleted file mode 100644 index d309e45..0000000 --- a/src/math/mod.rs +++ /dev/null @@ -1,20 +0,0 @@ -//! This module provides mathematical utilities and physical constants for the cheq library. -//! -//! It contains fundamental constants for unit conversions and numerical thresholds, as well as -//! functions for computing screened electrostatic interactions between atoms. These components -//! support the charge equilibration algorithm by providing the necessary mathematical infrastructure -//! for potential energy calculations and parameter handling. - -/// Physical and mathematical constants used throughout the library. -/// -/// This module defines essential constants for unit conversions (Bohr to angstrom, Hartree to eV) -/// and numerical tolerances, ensuring consistent handling of different energy and length scales -/// in computational chemistry calculations. -pub mod constants; - -/// Functions for computing screened Coulomb potentials. -/// -/// This module implements the orbital screening model used in charge equilibration, providing -/// analytical integrals for electrostatic interactions between atoms with Gaussian orbital -/// approximations. The screening accounts for charge penetration and orbital overlap effects. -pub mod shielding; diff --git a/src/math/shielding.rs b/src/math/shielding.rs deleted file mode 100644 index 4964530..0000000 --- a/src/math/shielding.rs +++ /dev/null @@ -1,325 +0,0 @@ -//! This module implements screened Coulomb potential calculations for charge equilibration. -//! -//! It provides functions to compute the effective electrostatic interaction between atoms, -//! accounting for orbital overlap and screening effects. The implementation uses Gaussian -//! approximations of Slater-type orbitals to enable efficient analytical integration of -//! Coulomb potentials, which is crucial for the linear system construction in the QEq solver. - -use super::constants::DISTANCE_THRESHOLD_BOHR; -use libm::erf; - -/// Computes the screened Coulomb integral between two atoms using Gaussian orbital approximations. -/// -/// This function calculates the effective electrostatic hardness parameter J_AB between atoms A and B, -/// incorporating screening effects due to orbital overlap. It transforms Slater-type orbitals into -/// equivalent Gaussians and evaluates the resulting analytical integral. -/// -/// The screening parameter λ controls the degree of orbital contraction, with smaller values -/// corresponding to more diffuse, strongly screening orbitals. -/// -/// # Arguments -/// -/// * `distance_bohr` - Interatomic distance in Bohr units. -/// * `n1` - Principal quantum number of atom 1. -/// * `r_cov1_bohr` - Covalent radius of atom 1 in Bohr units. -/// * `n2` - Principal quantum number of atom 2. -/// * `r_cov2_bohr` - Covalent radius of atom 2 in Bohr units. -/// * `lambda` - Screening parameter (dimensionless, typically 0.4-0.6). -/// -/// # Returns -/// -/// Returns the screened Coulomb integral in Hartree units. -#[inline] -pub fn gaussian_coulomb_integral( - distance_bohr: f64, - n1: u8, - r_cov1_bohr: f64, - n2: u8, - r_cov2_bohr: f64, - lambda: f64, -) -> f64 { - let zeta1 = slater_exponent(n1, r_cov1_bohr, lambda); - let zeta2 = slater_exponent(n2, r_cov2_bohr, lambda); - - let alpha1 = equivalent_gaussian_exponent(n1, zeta1); - let alpha2 = equivalent_gaussian_exponent(n2, zeta2); - - screened_potential(distance_bohr, alpha1, alpha2) -} - -/// Calculates the Slater orbital exponent from atomic parameters. -/// -/// The exponent ζ is derived from the principal quantum number and covalent radius, -/// scaled by the screening parameter λ. This follows the relationship ζ ∝ (2n+1)/(2R_cov). -/// -/// # Arguments -/// -/// * `n` - Principal quantum number. -/// * `r_cov_bohr` - Covalent radius in Bohr units. -/// * `lambda` - Screening parameter. -/// -/// # Returns -/// -/// Returns the Slater exponent in inverse Bohr units. -#[inline] -fn slater_exponent(n: u8, r_cov_bohr: f64, lambda: f64) -> f64 { - if r_cov_bohr < 1e-8 { - return 0.0; - } - // Slater exponent scales inversely with orbital size - lambda * (2.0 * n as f64 + 1.0) / (2.0 * r_cov_bohr) -} - -/// Converts a Slater orbital to an equivalent Gaussian exponent. -/// -/// This function maps Slater-type orbitals (1s, 2s, etc.) to Gaussian functions with -/// equivalent radial extent, using fitted coefficients that preserve the orbital shape. -/// -/// # Arguments -/// -/// * `n` - Principal quantum number. -/// * `zeta` - Slater exponent. -/// -/// # Returns -/// -/// Returns the equivalent Gaussian exponent α. -#[inline] -fn equivalent_gaussian_exponent(n: u8, zeta: f64) -> f64 { - if n == 0 { - return 0.0; - } - let c_n = get_c_n_fit_coeff(n); - // Gaussian exponent derived from Slater-Gaussian equivalence - c_n * zeta.powi(2) / (n as f64) -} - -/// Retrieves the fitted coefficient for Slater-to-Gaussian conversion. -/// -/// These coefficients are empirically determined to minimize the difference between -/// Slater and Gaussian orbital shapes for different principal quantum numbers. -/// -/// # Arguments -/// -/// * `n` - Principal quantum number. -/// -/// # Returns -/// -/// Returns the fitting coefficient c_n. -#[inline] -fn get_c_n_fit_coeff(n: u8) -> f64 { - match n { - 1 => 0.270917, - 2 => 0.098800, - 3 => 0.055600, - 4 => 0.039100, - 5 => 0.029600, - // Extrapolation formula for higher quantum numbers - _ => 0.27 * (n as f64).powf(-1.35), - } -} - -/// Evaluates the screened Coulomb potential between two Gaussian orbitals. -/// -/// This function computes the analytical integral ∫∫ φ₁(r) φ₂(r') / |r-r'| dr dr' -/// for Gaussian basis functions, using the error function for the distance-dependent term. -/// -/// # Arguments -/// -/// * `distance_bohr` - Interatomic distance in Bohr units. -/// * `alpha1` - Gaussian exponent for atom 1. -/// * `alpha2` - Gaussian exponent for atom 2. -/// -/// # Returns -/// -/// Returns the potential energy in Hartree units. -#[inline] -fn screened_potential(distance_bohr: f64, alpha1: f64, alpha2: f64) -> f64 { - let beta = (2.0 * alpha1 * alpha2 / (alpha1 + alpha2)).sqrt(); - - if distance_bohr < DISTANCE_THRESHOLD_BOHR { - // Limit at zero distance: erf(∞)/R → 1/R, but βR → 0, so use asymptotic form - 2.0 * beta / std::f64::consts::PI.sqrt() - } else { - erf(beta * distance_bohr) / distance_bohr - } -} - -#[cfg(test)] -mod tests { - use super::*; - use crate::math::constants; - use approx::assert_relative_eq; - - #[test] - fn test_slater_exponent_calculation() { - let r_cov_bohr_h = 0.371 / constants::BOHR_TO_ANGSTROM; - let result = slater_exponent(1, r_cov_bohr_h, 0.5); - let expected = (2.0 * 1.0 + 1.0) / (4.0 * r_cov_bohr_h); - assert_relative_eq!(result, expected, epsilon = 1e-12); - - let r_cov_bohr_c = 0.759 / constants::BOHR_TO_ANGSTROM; - let result_c = slater_exponent(2, r_cov_bohr_c, 0.5); - let expected_c = (2.0 * 2.0 + 1.0) / (4.0 * r_cov_bohr_c); - assert_relative_eq!(result_c, expected_c, epsilon = 1e-12); - - let result_lambda_custom = slater_exponent(1, r_cov_bohr_h, 0.4913); - let expected_lambda_custom = 0.4913 * (2.0 * 1.0 + 1.0) / (2.0 * r_cov_bohr_h); - assert_relative_eq!( - result_lambda_custom, - expected_lambda_custom, - epsilon = 1e-12 - ); - - assert_eq!(slater_exponent(1, 0.0, 0.5), 0.0); - } - - #[test] - fn test_get_c_n_fit_coeff_values() { - assert_eq!(get_c_n_fit_coeff(1), 0.270917); - assert_eq!(get_c_n_fit_coeff(2), 0.098800); - assert_eq!(get_c_n_fit_coeff(3), 0.055600); - assert_eq!(get_c_n_fit_coeff(4), 0.039100); - assert_eq!(get_c_n_fit_coeff(5), 0.029600); - - let n6 = 6u8; - let expected_c6 = 0.27 * (n6 as f64).powf(-1.35); - assert_eq!(get_c_n_fit_coeff(n6), expected_c6); - } - - #[test] - fn test_equivalent_gaussian_exponent_calculation() { - assert_relative_eq!( - equivalent_gaussian_exponent(1, 1.0), - get_c_n_fit_coeff(1) / 1.0 - ); - assert_relative_eq!( - equivalent_gaussian_exponent(2, 1.0), - get_c_n_fit_coeff(2) / 2.0 - ); - - let zeta = 2.0; - let alpha = equivalent_gaussian_exponent(1, zeta); - let expected_alpha = get_c_n_fit_coeff(1) * zeta.powi(2) / 1.0; - assert_relative_eq!(alpha, expected_alpha); - - assert_eq!(equivalent_gaussian_exponent(0, 1.0), 0.0); - } - - #[test] - fn test_integral_at_zero_distance_limit() { - let r_cov_bohr = 0.371 / constants::BOHR_TO_ANGSTROM; - let result_at_zero = gaussian_coulomb_integral(0.0, 1, r_cov_bohr, 1, r_cov_bohr, 0.5); - let result_near_zero = gaussian_coulomb_integral(1e-15, 1, r_cov_bohr, 1, r_cov_bohr, 0.5); - - assert!(result_at_zero > 0.0, "Result at R=0 must be positive"); - assert!( - !result_at_zero.is_infinite(), - "Result at R=0 must be finite" - ); - - assert_relative_eq!(result_at_zero, result_near_zero, epsilon = 1e-9); - } - - #[test] - fn test_integral_at_large_distance_asymptote() { - let r_cov_bohr_h = 0.371 / constants::BOHR_TO_ANGSTROM; - let r_cov_bohr_c = 0.759 / constants::BOHR_TO_ANGSTROM; - let large_distance_bohr = 1000.0; - - let result = - gaussian_coulomb_integral(large_distance_bohr, 1, r_cov_bohr_h, 2, r_cov_bohr_c, 0.5); - - let point_charge_result = 1.0 / large_distance_bohr; - - assert_relative_eq!(result, point_charge_result, epsilon = 1e-9); - } - - #[test] - fn test_integral_with_known_values() { - const LAMBDA: f64 = 0.4913; - - let r_h2_bohr = 0.74 / constants::BOHR_TO_ANGSTROM; - let r_cov_h_bohr = 0.371 / constants::BOHR_TO_ANGSTROM; - let j_ev_h2 = - gaussian_coulomb_integral(r_h2_bohr, 1, r_cov_h_bohr, 1, r_cov_h_bohr, LAMBDA) - * constants::HARTREE_TO_EV; - assert_relative_eq!(j_ev_h2, 14.02504752, epsilon = 1e-8); - - let r_co_bohr = 1.128 / constants::BOHR_TO_ANGSTROM; - let r_cov_c_bohr = 0.759 / constants::BOHR_TO_ANGSTROM; - let r_cov_o_bohr = 0.669 / constants::BOHR_TO_ANGSTROM; - let j_ev_co = - gaussian_coulomb_integral(r_co_bohr, 2, r_cov_c_bohr, 2, r_cov_o_bohr, LAMBDA) - * constants::HARTREE_TO_EV; - assert_relative_eq!(j_ev_co, 5.837573337, epsilon = 1e-8); - } - - #[test] - fn test_integral_symmetry_property() { - let r_cov_bohr_h = 0.371 / constants::BOHR_TO_ANGSTROM; - let r_cov_bohr_si = 1.176 / constants::BOHR_TO_ANGSTROM; - let distance_bohr = 3.0; - let lambda = 0.5; - - let result_h_si = - gaussian_coulomb_integral(distance_bohr, 1, r_cov_bohr_h, 3, r_cov_bohr_si, lambda); - let result_si_h = - gaussian_coulomb_integral(distance_bohr, 3, r_cov_bohr_si, 1, r_cov_bohr_h, lambda); - - assert_eq!(result_h_si, result_si_h); - } - - #[test] - fn test_integral_scaling_with_lambda() { - let r_cov_bohr = 0.7 / constants::BOHR_TO_ANGSTROM; - let distance_bohr = 1.5; - - let result_lambda_low = - gaussian_coulomb_integral(distance_bohr, 1, r_cov_bohr, 1, r_cov_bohr, 0.4); - let result_lambda_high = - gaussian_coulomb_integral(distance_bohr, 1, r_cov_bohr, 1, r_cov_bohr, 0.6); - - assert!( - result_lambda_low < result_lambda_high, - "Smaller lambda should lead to stronger screening (smaller J_AB)" - ); - } - - #[test] - fn test_integral_chemical_trends() { - let r_cov_bohr = 0.7 / constants::BOHR_TO_ANGSTROM; - let j_at_1_angstrom = gaussian_coulomb_integral( - 1.0 / constants::BOHR_TO_ANGSTROM, - 1, - r_cov_bohr, - 1, - r_cov_bohr, - 0.5, - ); - let j_at_3_angstroms = gaussian_coulomb_integral( - 3.0 / constants::BOHR_TO_ANGSTROM, - 1, - r_cov_bohr, - 1, - r_cov_bohr, - 0.5, - ); - assert!( - j_at_1_angstrom > j_at_3_angstroms, - "Interaction should decrease with distance" - ); - - let distance_bohr = 3.0 / constants::BOHR_TO_ANGSTROM; - - let r_h_bohr = 0.371 / constants::BOHR_TO_ANGSTROM; - let j_hh = gaussian_coulomb_integral(distance_bohr, 1, r_h_bohr, 1, r_h_bohr, 0.5); - - let r_na_bohr = 2.085 / constants::BOHR_TO_ANGSTROM; - let j_nana = gaussian_coulomb_integral(distance_bohr, 3, r_na_bohr, 3, r_na_bohr, 0.5); - - assert!( - j_hh > j_nana, - "Interaction between diffuse orbitals (Na-Na) should be weaker than between compact ones (H-H) at the same distance" - ); - } -} diff --git a/src/math/constants.rs b/src/shielding/constants.rs similarity index 97% rename from src/math/constants.rs rename to src/shielding/constants.rs index 07178b9..3ce417b 100644 --- a/src/math/constants.rs +++ b/src/shielding/constants.rs @@ -1,4 +1,4 @@ -//! This module defines fundamental physical and mathematical constants used throughout the cheq library. +//! This module defines fundamental physical and mathematical constants used throughout the library. //! //! These constants are essential for unit conversions and numerical thresholds in charge equilibration //! calculations, ensuring consistent handling of atomic units, energy scales, and geometric tolerances. diff --git a/src/shielding/gto.rs b/src/shielding/gto.rs new file mode 100644 index 0000000..92701bf --- /dev/null +++ b/src/shielding/gto.rs @@ -0,0 +1,267 @@ +//! Gaussian-Type Orbital (GTO) approximation for screened Coulomb interactions. +//! +//! This module implements the calculation of the effective electrostatic interaction +//! between two charge distributions approximated by Gaussian functions. This approach +//! allows for efficient analytical evaluation of the Coulomb integral, serving as a +//! faster but approximate alternative to the exact Slater-Type Orbital (STO) method. + +use super::constants::DISTANCE_THRESHOLD_BOHR; +use libm::erf; +use std::f64::consts::PI; + +/// Computes the screened Coulomb integral between two atoms using +/// the Gaussian-Type Orbital (GTO) approximation. +/// +/// # Arguments +/// +/// * `distance_bohr` - The distance between the two atomic centers in Bohr. +/// * `n1` - Principal quantum number of the first atom. +/// * `r_cov1_bohr` - Covalent radius of the first atom in Bohr. +/// * `n2` - Principal quantum number of the second atom. +/// * `r_cov2_bohr` - Covalent radius of the second atom in Bohr. +/// * `lambda` - The global screening scaling parameter. +/// +/// # Returns +/// +/// The Coulomb integral value in Hartree. +#[inline] +pub fn calculate_integral( + distance_bohr: f64, + n1: u8, + r_cov1_bohr: f64, + n2: u8, + r_cov2_bohr: f64, + lambda: f64, +) -> f64 { + let zeta1 = compute_slater_exponent(n1, r_cov1_bohr, lambda); + let zeta2 = compute_slater_exponent(n2, r_cov2_bohr, lambda); + + let alpha1 = compute_gaussian_exponent(n1, zeta1); + let alpha2 = compute_gaussian_exponent(n2, zeta2); + + compute_screened_potential(distance_bohr, alpha1, alpha2) +} + +/// Computes the Slater exponent from atomic parameters. +#[inline] +fn compute_slater_exponent(n: u8, r_cov_bohr: f64, lambda: f64) -> f64 { + if r_cov_bohr < DISTANCE_THRESHOLD_BOHR { + return 0.0; + } + lambda * (2.0 * n as f64 + 1.0) / (2.0 * r_cov_bohr) +} + +/// Computes the equivalent Gaussian exponent alpha from the Slater exponent zeta. +#[inline] +fn compute_gaussian_exponent(n: u8, zeta: f64) -> f64 { + if n == 0 { + return 0.0; + } + let c_n = get_fitting_coefficient(n); + c_n * zeta * zeta / (n as f64) +} + +/// Returns the empirical scaling coefficient for the GTO approximation. +/// +/// These values are derived to minimize the difference between the Slater and +/// Gaussian density profiles for a given principal quantum number n. +#[inline] +fn get_fitting_coefficient(n: u8) -> f64 { + match n { + 1 => 0.270917, + 2 => 0.098800, + 3 => 0.055600, + 4 => 0.039100, + 5 => 0.029600, + _ => 0.27 * (n as f64).powf(-1.35), + } +} + +/// Evaluates the Coulomb integral between two Gaussian charge distributions. +#[inline] +fn compute_screened_potential(r: f64, alpha1: f64, alpha2: f64) -> f64 { + let beta = (2.0 * alpha1 * alpha2 / (alpha1 + alpha2)).sqrt(); + + if r > DISTANCE_THRESHOLD_BOHR { + erf(beta * r) / r + } else { + 2.0 * beta / PI.sqrt() + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::shielding::constants; + use approx::assert_relative_eq; + + #[test] + fn test_slater_exponent_calculation() { + let r_cov_bohr_h = 0.371 / constants::BOHR_TO_ANGSTROM; + let result = compute_slater_exponent(1, r_cov_bohr_h, 0.5); + let expected = (2.0 * 1.0 + 1.0) / (4.0 * r_cov_bohr_h); + assert_relative_eq!(result, expected, epsilon = 1e-12); + + let r_cov_bohr_c = 0.759 / constants::BOHR_TO_ANGSTROM; + let result_c = compute_slater_exponent(2, r_cov_bohr_c, 0.5); + let expected_c = (2.0 * 2.0 + 1.0) / (4.0 * r_cov_bohr_c); + assert_relative_eq!(result_c, expected_c, epsilon = 1e-12); + + let result_lambda_custom = compute_slater_exponent(1, r_cov_bohr_h, 0.4913); + let expected_lambda_custom = 0.4913 * (2.0 * 1.0 + 1.0) / (2.0 * r_cov_bohr_h); + assert_relative_eq!( + result_lambda_custom, + expected_lambda_custom, + epsilon = 1e-12 + ); + + assert_eq!(compute_slater_exponent(1, 0.0, 0.5), 0.0); + } + + #[test] + fn test_get_fitting_coefficient_values() { + assert_eq!(get_fitting_coefficient(1), 0.270917); + assert_eq!(get_fitting_coefficient(2), 0.098800); + assert_eq!(get_fitting_coefficient(3), 0.055600); + assert_eq!(get_fitting_coefficient(4), 0.039100); + assert_eq!(get_fitting_coefficient(5), 0.029600); + + let n6 = 6u8; + let expected_c6 = 0.27 * (n6 as f64).powf(-1.35); + assert_eq!(get_fitting_coefficient(n6), expected_c6); + } + + #[test] + fn test_gaussian_exponent_calculation() { + assert_relative_eq!( + compute_gaussian_exponent(1, 1.0), + get_fitting_coefficient(1) / 1.0 + ); + assert_relative_eq!( + compute_gaussian_exponent(2, 1.0), + get_fitting_coefficient(2) / 2.0 + ); + + let zeta = 2.0; + let alpha = compute_gaussian_exponent(1, zeta); + let expected_alpha = get_fitting_coefficient(1) * zeta.powi(2) / 1.0; + assert_relative_eq!(alpha, expected_alpha); + + assert_eq!(compute_gaussian_exponent(0, 1.0), 0.0); + } + + #[test] + fn test_integral_at_zero_distance_limit() { + let r_cov_bohr = 0.371 / constants::BOHR_TO_ANGSTROM; + let result_at_zero = calculate_integral(0.0, 1, r_cov_bohr, 1, r_cov_bohr, 0.5); + let result_near_zero = calculate_integral(1e-15, 1, r_cov_bohr, 1, r_cov_bohr, 0.5); + + assert!(result_at_zero > 0.0, "Result at R=0 must be positive"); + assert!( + !result_at_zero.is_infinite(), + "Result at R=0 must be finite" + ); + + assert_relative_eq!(result_at_zero, result_near_zero, epsilon = 1e-9); + } + + #[test] + fn test_integral_at_large_distance_asymptote() { + let r_cov_bohr_h = 0.371 / constants::BOHR_TO_ANGSTROM; + let r_cov_bohr_c = 0.759 / constants::BOHR_TO_ANGSTROM; + let large_distance_bohr = 1000.0; + + let result = calculate_integral(large_distance_bohr, 1, r_cov_bohr_h, 2, r_cov_bohr_c, 0.5); + + let point_charge_result = 1.0 / large_distance_bohr; + + assert_relative_eq!(result, point_charge_result, epsilon = 1e-9); + } + + #[test] + fn test_integral_with_known_values() { + const LAMBDA: f64 = 0.4913; + + let r_h2_bohr = 0.74 / constants::BOHR_TO_ANGSTROM; + let r_cov_h_bohr = 0.371 / constants::BOHR_TO_ANGSTROM; + let j_ev_h2 = calculate_integral(r_h2_bohr, 1, r_cov_h_bohr, 1, r_cov_h_bohr, LAMBDA) + * constants::HARTREE_TO_EV; + assert_relative_eq!(j_ev_h2, 14.02504752, epsilon = 1e-8); + + let r_co_bohr = 1.128 / constants::BOHR_TO_ANGSTROM; + let r_cov_c_bohr = 0.759 / constants::BOHR_TO_ANGSTROM; + let r_cov_o_bohr = 0.669 / constants::BOHR_TO_ANGSTROM; + let j_ev_co = calculate_integral(r_co_bohr, 2, r_cov_c_bohr, 2, r_cov_o_bohr, LAMBDA) + * constants::HARTREE_TO_EV; + assert_relative_eq!(j_ev_co, 5.837573337, epsilon = 1e-8); + } + + #[test] + fn test_integral_symmetry_property() { + let r_cov_bohr_h = 0.371 / constants::BOHR_TO_ANGSTROM; + let r_cov_bohr_si = 1.176 / constants::BOHR_TO_ANGSTROM; + let distance_bohr = 3.0; + let lambda = 0.5; + + let result_h_si = + calculate_integral(distance_bohr, 1, r_cov_bohr_h, 3, r_cov_bohr_si, lambda); + let result_si_h = + calculate_integral(distance_bohr, 3, r_cov_bohr_si, 1, r_cov_bohr_h, lambda); + + assert_eq!(result_h_si, result_si_h); + } + + #[test] + fn test_integral_scaling_with_lambda() { + let r_cov_bohr = 0.7 / constants::BOHR_TO_ANGSTROM; + let distance_bohr = 1.5; + + let result_lambda_low = + calculate_integral(distance_bohr, 1, r_cov_bohr, 1, r_cov_bohr, 0.4); + let result_lambda_high = + calculate_integral(distance_bohr, 1, r_cov_bohr, 1, r_cov_bohr, 0.6); + + assert!( + result_lambda_low < result_lambda_high, + "Smaller lambda should lead to stronger screening (smaller J_AB)" + ); + } + + #[test] + fn test_integral_chemical_trends() { + let r_cov_bohr = 0.7 / constants::BOHR_TO_ANGSTROM; + let j_at_1_angstrom = calculate_integral( + 1.0 / constants::BOHR_TO_ANGSTROM, + 1, + r_cov_bohr, + 1, + r_cov_bohr, + 0.5, + ); + let j_at_3_angstroms = calculate_integral( + 3.0 / constants::BOHR_TO_ANGSTROM, + 1, + r_cov_bohr, + 1, + r_cov_bohr, + 0.5, + ); + assert!( + j_at_1_angstrom > j_at_3_angstroms, + "Interaction should decrease with distance" + ); + + let distance_bohr = 3.0 / constants::BOHR_TO_ANGSTROM; + + let r_h_bohr = 0.371 / constants::BOHR_TO_ANGSTROM; + let j_hh = calculate_integral(distance_bohr, 1, r_h_bohr, 1, r_h_bohr, 0.5); + + let r_na_bohr = 2.085 / constants::BOHR_TO_ANGSTROM; + let j_nana = calculate_integral(distance_bohr, 3, r_na_bohr, 3, r_na_bohr, 0.5); + + assert!( + j_hh > j_nana, + "Interaction between diffuse orbitals (Na-Na) should be weaker than between compact ones (H-H) at the same distance" + ); + } +} diff --git a/src/shielding/mod.rs b/src/shielding/mod.rs new file mode 100644 index 0000000..e8cefbf --- /dev/null +++ b/src/shielding/mod.rs @@ -0,0 +1,8 @@ +//! This module provides functionality for calculating screened Coulomb interactions. +//! +//! It includes constants, and implementations for both Gaussian-Type Orbital (GTO) +//! and Slater-Type Orbital (STO) based integrals. + +pub mod constants; +pub mod gto; +pub mod sto; diff --git a/src/shielding/sto.rs b/src/shielding/sto.rs new file mode 100644 index 0000000..55651d2 --- /dev/null +++ b/src/shielding/sto.rs @@ -0,0 +1,110 @@ +//! Exact two-center Coulomb integrals over ns Slater-type orbitals. +//! +//! This module implements the EXACT computation of the Coulomb integral J_AB(R) +//! between two ns Slater orbital charge distributions, as required by the +//! Charge Equilibration (QEq) method of Rappé & Goddard (1991). +//! +//! This implementation delegates to the high-performance `sto-ns` crate. + +use super::constants::DISTANCE_THRESHOLD_BOHR; + +/// Computes the two-center Coulomb integral between two ns Slater orbital +/// charge distributions using the exact expansion. +/// +/// # Arguments +/// +/// * `distance_bohr` - The distance between the two atomic centers in Bohr. +/// * `n1` - Principal quantum number of the first atom. +/// * `r_cov1_bohr` - Covalent radius of the first atom in Bohr. +/// * `n2` - Principal quantum number of the second atom. +/// * `r_cov2_bohr` - Covalent radius of the second atom in Bohr. +/// * `lambda` - The global screening scaling parameter. +/// +/// # Returns +/// +/// The Coulomb integral value in Hartree. +#[inline] +pub fn calculate_integral( + distance_bohr: f64, + n1: u8, + r_cov1_bohr: f64, + n2: u8, + r_cov2_bohr: f64, + lambda: f64, +) -> f64 { + let zeta1 = compute_slater_exponent(n1, r_cov1_bohr, lambda); + let zeta2 = compute_slater_exponent(n2, r_cov2_bohr, lambda); + + sto_ns::sto_coulomb_integral(distance_bohr, n1, zeta1, n2, zeta2) +} + +/// Computes the Slater exponent from atomic parameters. +#[inline] +fn compute_slater_exponent(n: u8, r_cov_bohr: f64, lambda: f64) -> f64 { + if r_cov_bohr < DISTANCE_THRESHOLD_BOHR { + return 0.0; + } + lambda * (2.0 * n as f64 + 1.0) / (2.0 * r_cov_bohr) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::shielding::constants; + use approx::assert_relative_eq; + + #[test] + fn test_sto_integral_symmetry() { + let r_cov_bohr_h = 0.371 / constants::BOHR_TO_ANGSTROM; + let r_cov_bohr_c = 0.759 / constants::BOHR_TO_ANGSTROM; + let distance_bohr = 2.0; + let lambda = 0.5; + + let result_h_c = + calculate_integral(distance_bohr, 1, r_cov_bohr_h, 2, r_cov_bohr_c, lambda); + let result_c_h = + calculate_integral(distance_bohr, 2, r_cov_bohr_c, 1, r_cov_bohr_h, lambda); + + assert_relative_eq!(result_h_c, result_c_h, epsilon = 1e-12); + } + + #[test] + fn test_sto_integral_large_distance() { + let r_cov_bohr = 1.0; + let distance_bohr = 100.0; + let lambda = 0.5; + + let result = calculate_integral(distance_bohr, 1, r_cov_bohr, 1, r_cov_bohr, lambda); + let expected = 1.0 / distance_bohr; + + assert_relative_eq!(result, expected, epsilon = 1e-6); + } + + #[test] + fn test_sto_integral_zero_distance() { + let r_cov_bohr = 1.0; + let lambda = 0.5; + + let result_zero = calculate_integral(0.0, 1, r_cov_bohr, 1, r_cov_bohr, lambda); + let result_near_zero = calculate_integral(1e-10, 1, r_cov_bohr, 1, r_cov_bohr, lambda); + + assert!(result_zero.is_finite()); + assert!(result_zero > 0.0); + assert_relative_eq!(result_zero, result_near_zero, epsilon = 1e-6); + } + + #[test] + fn test_sto_vs_gto_consistency() { + use crate::shielding::gto; + + let r_cov_bohr = 0.7; + let distance_bohr = 1.5; + let lambda = 0.5; + let n = 2; + + let sto_val = calculate_integral(distance_bohr, n, r_cov_bohr, n, r_cov_bohr, lambda); + let gto_val = gto::calculate_integral(distance_bohr, n, r_cov_bohr, n, r_cov_bohr, lambda); + + assert_relative_eq!(sto_val, gto_val, epsilon = 0.25); + } +} diff --git a/src/solver/implementation.rs b/src/solver/implementation.rs index 4864e2a..a37b1c9 100644 --- a/src/solver/implementation.rs +++ b/src/solver/implementation.rs @@ -7,16 +7,15 @@ //! and `Parameters` for element-specific values, enabling decoupled and flexible molecular //! simulations. -use super::options::SolverOptions; +use super::options::{BasisType, DampingStrategy, SolverOptions}; use crate::{ error::CheqError, - math, params::{ElementData, Parameters}, + shielding::{self, constants}, types::{AtomView, CalculationResult}, }; use faer::{Col, Mat, prelude::*}; use rayon::prelude::*; -use std::collections::HashMap; use std::panic::{self, AssertUnwindSafe}; /// Charge dependence factor for hydrogen hardness, derived from empirical fitting. @@ -24,6 +23,35 @@ use std::panic::{self, AssertUnwindSafe}; /// introducing non-linearity into the QEq equations. const H_CHARGE_DEPENDENCE_FACTOR: f64 = 0.93475415965; +/// A thread-safe wrapper for raw matrix access to enable parallel filling. +/// +/// This struct allows multiple threads to write to disjoint parts of a matrix +/// without locking, which is safe because we ensure unique indices in the parallel iterator. +struct UnsafeMatView { + ptr: *mut f64, + row_stride: isize, + col_stride: isize, +} + +unsafe impl Send for UnsafeMatView {} +unsafe impl Sync for UnsafeMatView {} + +impl UnsafeMatView { + /// Writes a value to the matrix at the specified (row, col) index. + /// + /// # Safety + /// + /// The caller must ensure that: + /// 1. The (row, col) indices are within bounds. + /// 2. No other thread is writing to the same address simultaneously. + unsafe fn write(&self, row: usize, col: usize, val: f64) { + let offset = (row as isize) * self.row_stride + (col as isize) * self.col_stride; + unsafe { + *self.ptr.offset(offset) = val; + } + } +} + /// The main solver for charge equilibration calculations. /// /// This struct holds references to atomic parameters and solver options, providing methods @@ -50,10 +78,11 @@ impl<'p> QEqSolver<'p> { /// # Examples /// /// ``` - /// use cheq::{get_default_parameters, QEqSolver}; + /// use cheq::get_default_parameters; + /// use cheq::QEqSolver; /// /// let params = get_default_parameters(); - /// let solver = QEqSolver::new(¶ms); + /// let solver = QEqSolver::new(params); /// ``` pub fn new(parameters: &'p Parameters) -> Self { Self { @@ -78,11 +107,17 @@ impl<'p> QEqSolver<'p> { /// # Examples /// /// ``` - /// use cheq::{get_default_parameters, QEqSolver, SolverOptions}; + /// use cheq::get_default_parameters; + /// use cheq::{QEqSolver, SolverOptions}; /// /// let params = get_default_parameters(); - /// let options = SolverOptions { tolerance: 1e-8, max_iterations: 50, lambda_scale: 1.0, ..SolverOptions::default() }; - /// let solver = QEqSolver::new(¶ms).with_options(options); + /// let options = SolverOptions { + /// max_iterations: 100, + /// tolerance: 1e-6, + /// ..Default::default() + /// }; + /// + /// let solver = QEqSolver::new(params).with_options(options); /// ``` pub fn with_options(mut self, options: SolverOptions) -> Self { self.options = options; @@ -108,25 +143,33 @@ impl<'p> QEqSolver<'p> { /// /// # Errors /// - /// This function can return the following errors: - /// * `CheqError::NoAtoms` if the atom slice is empty. - /// * `CheqError::ParameterNotFound` if parameters for an atom's element are missing. - /// * `CheqError::LinalgError` if the linear system cannot be solved due to numerical issues. - /// * `CheqError::NotConverged` if the SCF procedure does not converge within the maximum iterations. + /// * `CheqError::NoAtoms` - If the input atom list is empty. + /// * `CheqError::ParameterNotFound` - If an atom's parameters are missing. + /// * `CheqError::NotConverged` - If the SCF procedure fails to converge within the maximum iterations. + /// * `CheqError::LinalgError` - If the linear system solver encounters an error. /// /// # Examples /// /// ``` - /// use cheq::{get_default_parameters, QEqSolver, Atom}; + /// use cheq::get_default_parameters; + /// use cheq::QEqSolver; + /// use cheq::Atom; /// + /// // 1. Setup parameters and solver /// let params = get_default_parameters(); - /// let solver = QEqSolver::new(¶ms); + /// let solver = QEqSolver::new(params); + /// + /// // 2. Define a molecule (e.g., H2) /// let atoms = vec![ - /// Atom { atomic_number: 6, position: [0.0, 0.0, 0.0] }, - /// Atom { atomic_number: 8, position: [1.128, 0.0, 0.0] }, + /// Atom { atomic_number: 1, position: [0.0, 0.0, 0.0] }, + /// Atom { atomic_number: 1, position: [0.74, 0.0, 0.0] }, /// ]; + /// + /// // 3. Run calculation /// let result = solver.solve(&atoms, 0.0).unwrap(); + /// /// assert_eq!(result.charges.len(), 2); + /// println!("Charges: {:?}", result.charges); /// ``` pub fn solve( &self, @@ -139,19 +182,18 @@ impl<'p> QEqSolver<'p> { } let element_data = self.fetch_element_data(atoms)?; + let invariant = self.build_invariant_system(atoms, &element_data, total_charge)?; - let has_hydrogen = !invariant.hydrogen_meta.is_empty(); + let mut charges = Col::zeros(n_atoms); + let has_hydrogen = !invariant.hydrogen_meta.is_empty(); let hydrogen_scf = self.options.hydrogen_scf && has_hydrogen; - let inner_iters = if hydrogen_scf { - self.options.hydrogen_inner_iters - } else { - 0 - }; + + let mut work_matrix = invariant.base_matrix.clone(); if !hydrogen_scf { let (_, equilibrated_potential) = - self.run_single_solve(&invariant, &mut charges, hydrogen_scf)?; + self.run_single_solve(&invariant, &mut work_matrix, &mut charges, false, 1.0)?; return Ok(CalculationResult { charges: charges.as_ref().iter().cloned().collect(), equilibrated_potential, @@ -160,29 +202,50 @@ impl<'p> QEqSolver<'p> { } let mut max_charge_delta = 0.0; + let mut prev_delta = f64::MAX; + + let mut current_damping = match self.options.damping { + DampingStrategy::None => 1.0, + DampingStrategy::Fixed(d) => d, + DampingStrategy::Auto { initial } => initial, + }; for iteration in 1..=self.options.max_iterations { - if inner_iters > 0 { - for _ in 0..inner_iters { - let (delta, _) = - self.run_single_solve(&invariant, &mut charges, hydrogen_scf)?; - if delta < self.options.tolerance { - break; - } - } + if iteration > 1 { + work_matrix.copy_from(&invariant.base_matrix); } - let (delta, equilibrated_potential) = - self.run_single_solve(&invariant, &mut charges, hydrogen_scf)?; + let (delta, equilibrated_potential) = self.run_single_solve( + &invariant, + &mut work_matrix, + &mut charges, + true, + current_damping, + )?; + max_charge_delta = delta; - if !hydrogen_scf || max_charge_delta < self.options.tolerance { + if max_charge_delta < self.options.tolerance { return Ok(CalculationResult { charges: charges.as_ref().iter().cloned().collect(), equilibrated_potential, iterations: iteration, }); } + + if let DampingStrategy::Auto { initial: _ } = self.options.damping { + if max_charge_delta > prev_delta { + current_damping *= 0.5; + } else if max_charge_delta < prev_delta * 0.9 { + current_damping = (current_damping * 1.1).min(1.0); + } + + if current_damping < 0.001 { + current_damping = 0.001; + } + } + + prev_delta = max_charge_delta; } Err(CheqError::NotConverged { @@ -192,21 +255,6 @@ impl<'p> QEqSolver<'p> { } /// Retrieves element data for each atom from the parameters. - /// - /// This helper method maps each atom to its corresponding `ElementData` based on atomic number, - /// ensuring all required parameters are available before proceeding with calculations. - /// - /// # Arguments - /// - /// * `atoms` - A slice of atom data implementing the `AtomView` trait. - /// - /// # Returns - /// - /// A vector of references to `ElementData` for each atom. - /// - /// # Errors - /// - /// Returns `CheqError::ParameterNotFound` if data for any atomic number is missing. fn fetch_element_data( &self, atoms: &[A], @@ -239,7 +287,6 @@ impl<'p> QEqSolver<'p> { let mut base_matrix = Mat::zeros(matrix_size, matrix_size); let mut rhs = Col::zeros(matrix_size); - let mut hydrogen_meta = Vec::new(); for i in 0..n_atoms { @@ -254,27 +301,71 @@ impl<'p> QEqSolver<'p> { let positions: Vec<[f64; 3]> = atoms.iter().map(AtomView::position).collect(); - let off_diagonals = compute_off_diagonal_rows( - &positions, - element_data, - self.options.lambda_scale, - self.options.cutoff_radius, - ); - for (i, entries) in off_diagonals { - for (j, val) in entries { - base_matrix[(i, j)] = val; - base_matrix[(j, i)] = val; + let mat_view = UnsafeMatView { + ptr: base_matrix.as_ptr_mut(), + row_stride: base_matrix.row_stride(), + col_stride: base_matrix.col_stride(), + }; + + let lambda = self.options.lambda_scale; + let basis_type = self.options.basis_type; + + (0..n_atoms).into_par_iter().for_each(|i| { + let data_i = element_data[i]; + let pos_i = positions[i]; + let radius_i_bohr = data_i.radius / constants::BOHR_TO_ANGSTROM; + + for j in (i + 1)..n_atoms { + let pos_j = positions[j]; + let diff_sq = (pos_i[0] - pos_j[0]).powi(2) + + (pos_i[1] - pos_j[1]).powi(2) + + (pos_i[2] - pos_j[2]).powi(2); + + let dist_angstrom = diff_sq.sqrt(); + let dist_bohr = dist_angstrom / constants::BOHR_TO_ANGSTROM; + + let data_j = element_data[j]; + let radius_j_bohr = data_j.radius / constants::BOHR_TO_ANGSTROM; + + let integral_hartree = match basis_type { + BasisType::Gto => shielding::gto::calculate_integral( + dist_bohr, + data_i.principal_quantum_number, + radius_i_bohr, + data_j.principal_quantum_number, + radius_j_bohr, + lambda, + ), + BasisType::Sto => shielding::sto::calculate_integral( + dist_bohr, + data_i.principal_quantum_number, + radius_i_bohr, + data_j.principal_quantum_number, + radius_j_bohr, + lambda, + ), + }; + + let val_ev = integral_hartree * constants::HARTREE_TO_EV; + + // SAFETY: Each unordered pair (i, j) with i < j is handled only by the thread for i. + // That thread writes (i, j) and (j, i), so no two threads write the same entries. + unsafe { + mat_view.write(i, j, val_ev); + mat_view.write(j, i, val_ev); + } } - } + }); base_matrix .col_mut(matrix_size - 1) .subrows_mut(0, n_atoms) - .fill(-1.0); + .fill(1.0); base_matrix .row_mut(matrix_size - 1) .subcols_mut(0, n_atoms) .fill(1.0); + rhs[matrix_size - 1] = total_charge; Ok(InvariantSystem { @@ -288,16 +379,17 @@ impl<'p> QEqSolver<'p> { fn run_single_solve( &self, invariant: &InvariantSystem, + work_matrix: &mut Mat, charges: &mut Col, hydrogen_scf: bool, + damping: f64, ) -> Result<(f64, f64), CheqError> { let n_atoms = charges.nrows(); - let mut work_matrix = invariant.base_matrix.clone(); if hydrogen_scf { for &(idx, hardness) in &invariant.hydrogen_meta { - work_matrix[(idx, idx)] = - hardness * (1.0 + charges[idx] * H_CHARGE_DEPENDENCE_FACTOR); + let q_clamped = charges[idx].clamp(-0.95, 0.95); + work_matrix[(idx, idx)] = hardness * (1.0 + q_clamped * H_CHARGE_DEPENDENCE_FACTOR); } } @@ -309,22 +401,26 @@ impl<'p> QEqSolver<'p> { Ok(sol) => sol, Err(_) => { return Err(CheqError::LinalgError( - "Linear system is likely singular or ill-conditioned due to input geometry." - .to_string(), + "Linear system solver panicked. Matrix might be singular.".to_string(), )); } }; let new_charges = solution.as_ref().subrows(0, n_atoms); + let max_charge_delta = new_charges .as_ref() .iter() .zip(charges.as_ref().iter()) .map(|(new, old): (&f64, &f64)| (*new - *old).abs()) .fold(0.0, f64::max); - charges.as_mut().copy_from(&new_charges); - let equilibrated_potential = -solution[n_atoms]; + for i in 0..n_atoms { + charges[i] = (1.0 - damping) * charges[i] + damping * new_charges[i]; + } + + let equilibrated_potential = solution[n_atoms]; + Ok((max_charge_delta, equilibrated_potential)) } } @@ -336,519 +432,206 @@ struct InvariantSystem { hydrogen_meta: Vec<(usize, f64)>, } -/// Computes screened Coulomb off-diagonal entries for the upper triangle. -fn compute_off_diagonal_rows( - positions: &[[f64; 3]], - element_data: &[&ElementData], - lambda_scale: f64, - cutoff_radius: Option, -) -> Vec<(usize, Vec<(usize, f64)>)> { - match cutoff_radius { - None => compute_dense_off_diagonals(positions, element_data, lambda_scale), - Some(radius) => compute_cutoff_off_diagonals(positions, element_data, lambda_scale, radius), - } -} - -/// Computes all off-diagonal entries without cutoff. -fn compute_dense_off_diagonals( - positions: &[[f64; 3]], - element_data: &[&ElementData], - lambda_scale: f64, -) -> Vec<(usize, Vec<(usize, f64)>)> { - let n_atoms = positions.len(); - - (0..n_atoms) - .into_par_iter() - .map(|i| { - let data_i = element_data[i]; - let pos_i = positions[i]; - let radius_i_bohr = data_i.radius / math::constants::BOHR_TO_ANGSTROM; - let mut row_vals = Vec::with_capacity(n_atoms.saturating_sub(i + 1)); - - for j in (i + 1)..n_atoms { - let dist_sq: f64 = pos_i - .iter() - .zip(positions[j].iter()) - .map(|(pi, pj)| { - let diff = pi - pj; - diff * diff - }) - .sum(); - let val = compute_pair_entry( - dist_sq, - data_i, - element_data[j], - radius_i_bohr, - lambda_scale, - ); - row_vals.push((j, val)); - } - - (i, row_vals) - }) - .collect() -} - -/// Computes off-diagonal entries with a hard cutoff radius. -fn compute_cutoff_off_diagonals( - positions: &[[f64; 3]], - element_data: &[&ElementData], - lambda_scale: f64, - cutoff_radius: f64, -) -> Vec<(usize, Vec<(usize, f64)>)> { - let n_atoms = positions.len(); - let cell_size = cutoff_radius; - let mut cells: HashMap<(i64, i64, i64), Vec> = HashMap::new(); - - for (idx, pos) in positions.iter().enumerate() { - let key = ( - (pos[0] / cell_size).floor() as i64, - (pos[1] / cell_size).floor() as i64, - (pos[2] / cell_size).floor() as i64, - ); - cells.entry(key).or_default().push(idx); - } - - let cells = std::sync::Arc::new(cells); - let cutoff_sq = cutoff_radius * cutoff_radius; - - (0..n_atoms) - .into_par_iter() - .map(|i| { - let data_i = element_data[i]; - let pos_i = positions[i]; - let radius_i_bohr = data_i.radius / math::constants::BOHR_TO_ANGSTROM; - let key = ( - (pos_i[0] / cell_size).floor() as i64, - (pos_i[1] / cell_size).floor() as i64, - (pos_i[2] / cell_size).floor() as i64, - ); - - let mut entries = Vec::new(); - - for dx in -1..=1 { - for dy in -1..=1 { - for dz in -1..=1 { - let neighbor_key = (key.0 + dx, key.1 + dy, key.2 + dz); - if let Some(indices) = cells.get(&neighbor_key) { - for &j in indices { - if j <= i { - continue; - } - let dist_sq: f64 = pos_i - .iter() - .zip(positions[j].iter()) - .map(|(pi, pj)| { - let diff = pi - pj; - diff * diff - }) - .sum(); - if dist_sq > cutoff_sq { - continue; - } - let val = compute_pair_entry( - dist_sq, - data_i, - element_data[j], - radius_i_bohr, - lambda_scale, - ); - entries.push((j, val)); - } - } - } - } - } - - entries.sort_by_key(|(j, _)| *j); - (i, entries) - }) - .collect() -} - -/// Computes the screened Coulomb interaction entry between two atoms. -fn compute_pair_entry( - dist_sq: f64, - data_i: &ElementData, - data_j: &ElementData, - radius_i_bohr: f64, - lambda_scale: f64, -) -> f64 { - let distance_angstrom = dist_sq.sqrt(); - let distance_bohr = distance_angstrom / math::constants::BOHR_TO_ANGSTROM; - let radius_j_bohr = data_j.radius / math::constants::BOHR_TO_ANGSTROM; - - let j_ij_hartree = math::shielding::gaussian_coulomb_integral( - distance_bohr, - data_i.principal_quantum_number, - radius_i_bohr, - data_j.principal_quantum_number, - radius_j_bohr, - lambda_scale, - ); - - j_ij_hartree * math::constants::HARTREE_TO_EV -} - #[cfg(test)] mod tests { use super::*; - use crate::get_default_parameters; - use crate::types::Atom; + use crate::{get_default_parameters, types::Atom}; use approx::assert_relative_eq; - use std::panic; - - fn get_test_solver() -> QEqSolver<'static> { - let params = get_default_parameters(); - QEqSolver::new(params) - } #[test] - fn test_linear_system_carbon_monoxide() { - let solver = get_test_solver(); - let atoms = vec![ - Atom { - atomic_number: 6, - position: [0.0, 0.0, 0.0], - }, - Atom { - atomic_number: 8, - position: [1.128, 0.0, 0.0], - }, - ]; - - let result = solver.solve(&atoms, 0.0).unwrap(); - assert_eq!(result.iterations, 1); - let (q_c, q_o) = (result.charges[0], result.charges[1]); - assert!(q_o < 0.0); - assert!(q_c > 0.0); - assert_relative_eq!(q_c + q_o, 0.0, epsilon = 1e-9); - assert!( - q_o < -0.2 && q_o > -0.6, - "Oxygen charge magnitude seems off. Got: {}", - q_o - ); - } + fn test_h2_molecule_symmetry() { + let params = get_default_parameters(); + let solver = QEqSolver::new(params); - #[test] - fn test_nonlinear_system_water_molecule() { - let solver = get_test_solver(); - let bond_angle_rad = 104.45f64.to_radians(); - let bond_length = 0.9575; let atoms = vec![ - Atom { - atomic_number: 8, - position: [0.0, 0.0, 0.0], - }, Atom { atomic_number: 1, - position: [bond_length, 0.0, 0.0], + position: [0.0, 0.0, 0.0], }, Atom { atomic_number: 1, - position: [ - bond_length * bond_angle_rad.cos(), - bond_length * bond_angle_rad.sin(), - 0.0, - ], + position: [0.74, 0.0, 0.0], }, ]; let result = solver.solve(&atoms, 0.0).unwrap(); - assert!(result.iterations > 1); - let (q_o, q_h1, q_h2) = (result.charges[0], result.charges[1], result.charges[2]); - assert!(q_o < 0.0); - assert!(q_h1 > 0.0); - assert_relative_eq!(q_h1, q_h2, epsilon = 1e-7); - assert_relative_eq!(q_o + q_h1 + q_h2, 0.0, epsilon = 1e-9); - - assert!( - q_h1 > 0.15 && q_h1 < 0.4, - "Hydrogen charge is outside a reasonable physical range for water" - ); + assert_relative_eq!(result.charges[0], 0.0, epsilon = 1e-6); + assert_relative_eq!(result.charges[1], 0.0, epsilon = 1e-6); + assert_relative_eq!(result.charges[0], result.charges[1], epsilon = 1e-10); } #[test] - fn test_symmetry_dihydrogen_molecule() { - let solver = get_test_solver(); + fn test_hf_molecule_polarity() { + let params = get_default_parameters(); + let solver = QEqSolver::new(params); + let atoms = vec![ Atom { atomic_number: 1, position: [0.0, 0.0, 0.0], }, Atom { - atomic_number: 1, - position: [0.74, 0.0, 0.0], + atomic_number: 9, + position: [0.917, 0.0, 0.0], }, ]; let result = solver.solve(&atoms, 0.0).unwrap(); - assert_relative_eq!(result.charges[0], 0.0, epsilon = 1e-9); - assert_relative_eq!(result.charges[1], 0.0, epsilon = 1e-9); - } - - #[test] - fn test_ionic_system_lithium_hydride() { - let solver = get_test_solver(); - let atoms = vec![ - Atom { - atomic_number: 3, - position: [0.0, 0.0, 0.0], - }, - Atom { - atomic_number: 1, - position: [1.595, 0.0, 0.0], - }, - ]; - let result = solver.solve(&atoms, 0.0).unwrap(); - assert!(result.iterations > 1); - let (q_li, q_h) = (result.charges[0], result.charges[1]); + assert!(result.charges[0] > 0.0); + assert!(result.charges[1] < 0.0); - assert!(q_h < 0.0, "Hydrogen should be negative in LiH"); - assert!(q_li > 0.0, "Lithium should be positive in LiH"); - assert_relative_eq!(q_li + q_h, 0.0, epsilon = 1e-9); + assert_relative_eq!(result.charges[0] + result.charges[1], 0.0, epsilon = 1e-6); } #[test] - fn test_total_charge_constraint_hydroxide_ion() { - let solver = get_test_solver(); + fn test_water_molecule_sto_vs_gto() { + let params = get_default_parameters(); + let atoms = vec![ Atom { atomic_number: 8, - position: [0.0, 0.0, 0.0], + position: [0.000000, 0.000000, 0.117300], }, Atom { atomic_number: 1, - position: [0.96, 0.0, 0.0], - }, - ]; - - let result = solver.solve(&atoms, -1.0).unwrap(); - let (q_o, q_h) = (result.charges[0], result.charges[1]); - - assert_relative_eq!(q_o + q_h, -1.0, epsilon = 1e-9); - } - - #[test] - fn test_chemical_trend_nacl_vs_nabr() { - let solver = get_test_solver(); - - let nacl_bond_length = 2.36; - let nacl_atoms = vec![ - Atom { - atomic_number: 11, - position: [0.0, 0.0, 0.0], - }, - Atom { - atomic_number: 17, - position: [nacl_bond_length, 0.0, 0.0], - }, - ]; - - let nabr_bond_length = 2.50; - let nabr_atoms = vec![ - Atom { - atomic_number: 11, - position: [0.0, 0.0, 0.0], + position: [0.000000, 0.757200, -0.469200], }, Atom { - atomic_number: 35, - position: [nabr_bond_length, 0.0, 0.0], + atomic_number: 1, + position: [0.000000, -0.757200, -0.469200], }, ]; - let nacl_result = solver.solve(&nacl_atoms, 0.0).unwrap(); - let nabr_result = solver.solve(&nabr_atoms, 0.0).unwrap(); - - let charge_na_in_nacl = nacl_result.charges[0]; - let charge_na_in_nabr = nabr_result.charges[0]; - - assert!( - charge_na_in_nacl > charge_na_in_nabr, - "Chemical trend is incorrect: Charge on Na in NaCl should be greater than in NaBr." - ); + let solver_sto = QEqSolver::new(params); + let res_sto = solver_sto.solve(&atoms, 0.0).unwrap(); - assert!(charge_na_in_nacl > 0.0); - assert!(charge_na_in_nabr > 0.0); - } + assert!(res_sto.charges[0] < -0.1); + assert!(res_sto.charges[1] > 0.05); + assert_relative_eq!(res_sto.charges[1], res_sto.charges[2], epsilon = 1e-6); - #[test] - fn test_error_handling_no_atoms() { - let solver = get_test_solver(); - let atoms: Vec = vec![]; - let result = solver.solve(&atoms, 0.0); - assert!(matches!(result, Err(CheqError::NoAtoms))); - } + let options_gto = SolverOptions { + basis_type: BasisType::Gto, + ..SolverOptions::default() + }; + let solver_gto = QEqSolver::new(params).with_options(options_gto); + let res_gto = solver_gto.solve(&atoms, 0.0).unwrap(); - #[test] - fn test_error_handling_parameter_not_found() { - let solver = get_test_solver(); - let atoms = vec![Atom { - atomic_number: 118, - position: [0.0, 0.0, 0.0], - }]; - let result = solver.solve(&atoms, 0.0); - assert!(matches!(result, Err(CheqError::ParameterNotFound(118)))); + assert_relative_eq!(res_sto.charges[0], res_gto.charges[0], epsilon = 0.3); } #[test] - fn test_not_converged_error() { + fn test_convergence_failure() { let params = get_default_parameters(); let options = SolverOptions { max_iterations: 1, - tolerance: 1e-20, - lambda_scale: 1.0, ..SolverOptions::default() }; let solver = QEqSolver::new(params).with_options(options); + let atoms = vec![ Atom { - atomic_number: 8, + atomic_number: 3, position: [0.0, 0.0, 0.0], }, Atom { atomic_number: 1, - position: [0.9575, 0.0, 0.0], - }, - Atom { - atomic_number: 1, - position: [0.9575 * (-0.5), 0.9575 * (3.0_f64).sqrt() / 2.0, 0.0], + position: [1.595, 0.0, 0.0], }, ]; + let result = solver.solve(&atoms, 0.0); assert!(matches!(result, Err(CheqError::NotConverged { .. }))); } #[test] - fn test_linalg_error_singular_matrix() { - let solver = get_test_solver(); - let atoms = vec![ - Atom { - atomic_number: 6, - position: [0.0, 0.0, 0.0], - }, - Atom { - atomic_number: 6, - position: [0.0, 0.0, 0.0], - }, - ]; + fn test_error_no_atoms() { + let params = get_default_parameters(); + let solver = QEqSolver::new(params); + + let atoms: Vec = vec![]; let result = solver.solve(&atoms, 0.0); - match result { - Err(CheqError::LinalgError(_)) => (), - Ok(_) => (), - _ => panic!("Unexpected error"), - } + assert!(matches!(result, Err(CheqError::NoAtoms))); } #[test] - fn test_with_options() { + fn test_error_parameter_not_found() { let params = get_default_parameters(); - let options = SolverOptions { - max_iterations: 100, - tolerance: 1e-10, - lambda_scale: 1.0, - ..SolverOptions::default() - }; - let solver = QEqSolver::new(params).with_options(options); + let solver = QEqSolver::new(params); + let atoms = vec![Atom { - atomic_number: 6, + atomic_number: 118, position: [0.0, 0.0, 0.0], }]; - let result = solver.solve(&atoms, 0.0).unwrap(); - assert_eq!(result.charges.len(), 1); + let result = solver.solve(&atoms, 0.0); + + assert!(matches!(result, Err(CheqError::ParameterNotFound(118)))); } #[test] - fn test_cutoff_large_matches_dense() { - let atoms = vec![ - Atom { - atomic_number: 1, - position: [0.0, 0.0, 0.0], - }, - Atom { - atomic_number: 1, - position: [0.74, 0.0, 0.0], - }, - ]; - - let base_solver = get_test_solver(); - let base = base_solver.solve(&atoms, 0.0).unwrap(); - + fn test_damping_strategies() { let params = get_default_parameters(); - let cutoff_solver = QEqSolver::new(params).with_options(SolverOptions { - cutoff_radius: Some(10.0), - ..SolverOptions::default() - }); - let cutoff = cutoff_solver.solve(&atoms, 0.0).unwrap(); - - assert_relative_eq!(cutoff.charges[0], base.charges[0], epsilon = 1e-9); - assert_relative_eq!(cutoff.charges[1], base.charges[1], epsilon = 1e-9); - assert_relative_eq!(cutoff.charges.iter().sum::(), 0.0, epsilon = 1e-9); - } - #[test] - fn test_hydrogen_inner_iters_respects_charge() { let atoms = vec![ Atom { - atomic_number: 8, + atomic_number: 1, position: [0.0, 0.0, 0.0], }, Atom { - atomic_number: 1, - position: [0.96, 0.0, 0.0], + atomic_number: 9, + position: [0.917, 0.0, 0.0], }, ]; - let base = get_test_solver().solve(&atoms, -1.0).unwrap(); - - let params = get_default_parameters(); - let solver = QEqSolver::new(params).with_options(SolverOptions { - hydrogen_inner_iters: 2, + let options_fixed = SolverOptions { + damping: DampingStrategy::Fixed(0.5), ..SolverOptions::default() - }); - let result = solver.solve(&atoms, -1.0).unwrap(); + }; + let solver_fixed = QEqSolver::new(params).with_options(options_fixed); + let result_fixed = solver_fixed.solve(&atoms, 0.0).unwrap(); + assert!(result_fixed.charges[0] > 0.0); - assert_relative_eq!(result.charges.iter().sum::(), -1.0, epsilon = 1e-9); - assert_relative_eq!(result.charges[0], base.charges[0], epsilon = 1e-5); - assert_relative_eq!(result.charges[1], base.charges[1], epsilon = 1e-5); + let options_none = SolverOptions { + damping: DampingStrategy::None, + ..SolverOptions::default() + }; + let solver_none = QEqSolver::new(params).with_options(options_none); + let result_none = solver_none.solve(&atoms, 0.0).unwrap(); + assert!(result_none.charges[0] > 0.0); + + assert_relative_eq!( + result_fixed.charges[0], + result_none.charges[0], + epsilon = 0.01 + ); } #[test] - fn test_hydrogen_scf_off_converges_in_one_iteration() { + fn test_basis_type_gto() { let params = get_default_parameters(); - let solver = QEqSolver::new(params).with_options(SolverOptions { - hydrogen_scf: false, - ..SolverOptions::default() - }); - let bond_length = 0.9575; - let bond_angle_rad = 104.45f64.to_radians(); let atoms = vec![ - Atom { - atomic_number: 8, - position: [0.0, 0.0, 0.0], - }, Atom { atomic_number: 1, - position: [bond_length, 0.0, 0.0], + position: [0.0, 0.0, 0.0], }, Atom { - atomic_number: 1, - position: [ - bond_length * bond_angle_rad.cos(), - bond_length * bond_angle_rad.sin(), - 0.0, - ], + atomic_number: 9, + position: [0.917, 0.0, 0.0], }, ]; + let options = SolverOptions { + basis_type: BasisType::Gto, + ..SolverOptions::default() + }; + let solver = QEqSolver::new(params).with_options(options); let result = solver.solve(&atoms, 0.0).unwrap(); - assert_eq!(result.iterations, 1); - assert_relative_eq!(result.charges.iter().sum::(), 0.0, epsilon = 1e-9); - assert_relative_eq!(result.charges[1], result.charges[2], epsilon = 1e-9); + assert!(result.charges[0] > 0.0); + assert!(result.charges[1] < 0.0); + assert_relative_eq!(result.charges[0] + result.charges[1], 0.0, epsilon = 1e-6); } } diff --git a/src/solver/mod.rs b/src/solver/mod.rs index 1d42e22..ec333ab 100644 --- a/src/solver/mod.rs +++ b/src/solver/mod.rs @@ -6,14 +6,5 @@ mod implementation; mod options; -/// The `QEqSolver` struct for performing charge equilibration calculations. -/// -/// This is the main solver that implements the SCF iterative procedure to compute partial -/// atomic charges based on the QEq method. pub use implementation::QEqSolver; - -/// The `SolverOptions` struct for configuring solver parameters. -/// -/// This struct allows customization of convergence tolerance, maximum iterations, and -/// shielding parameters for the charge equilibration process. -pub use options::SolverOptions; +pub use options::{BasisType, DampingStrategy, SolverOptions}; diff --git a/src/solver/options.rs b/src/solver/options.rs index b22f929..23263d8 100644 --- a/src/solver/options.rs +++ b/src/solver/options.rs @@ -4,6 +4,46 @@ //! criteria, iteration limits, and screening parameters for the QEq algorithm. These options //! control the trade-off between computational cost and solution accuracy. +/// Specifies the type of basis functions used for Coulomb integrals. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] +pub enum BasisType { + /// Gaussian-Type Orbitals (GTO). + /// + /// Uses Gaussian approximations for Slater orbitals. This allows for fast analytical + /// integration but introduces a small approximation error. + Gto, + + /// Slater-Type Orbitals (STO). + /// + /// Uses exact Slater orbitals. This is more accurate but computationally more expensive + /// as it requires evaluating complex analytical expansions. + #[default] + Sto, +} + +/// Strategy for damping charge updates during SCF iterations. +#[derive(Debug, Clone, Copy, PartialEq)] +pub enum DampingStrategy { + /// No damping is applied (equivalent to damping = 1.0). + /// + /// Fastest convergence for simple, well-behaved systems, but prone to oscillation + /// or divergence in complex cases (e.g., LiH, large systems). + None, + + /// A fixed damping factor (0.0 < d <= 1.0). + /// + /// A constant mixing rate. Lower values (e.g., 0.1) are more stable but slower. + /// Higher values (e.g., 0.8) are faster but risk instability. + Fixed(f64), + + /// Automatically adjusts damping based on convergence behavior. + /// + /// Starts with an initial value. If the error increases (divergence), it reduces damping (brakes). + /// If the error decreases significantly (convergence), it increases damping (accelerates). + /// This is the recommended strategy for general use. + Auto { initial: f64 }, +} + /// Configuration parameters for the charge equilibration solver. /// /// This struct encapsulates the numerical settings that control the iterative solution process @@ -36,29 +76,26 @@ pub struct SolverOptions { /// When disabled, hydrogen uses a fixed hardness (lossy simplification). Enabled by default. pub hydrogen_scf: bool, - /// Optional hard cutoff radius (in Å) for pair interactions. + /// The type of basis functions to use for Coulomb integrals. /// - /// If `None`, all pairs are included (lossless). If `Some(r)`, pairs beyond `r` are skipped - /// (lossy, hard cutoff). - pub cutoff_radius: Option, + /// Defaults to `BasisType::Sto` (Slater-Type Orbitals) for maximum accuracy. + pub basis_type: BasisType, - /// Number of optional inner iterations focused on hydrogen before each global solve. + /// The damping strategy for the SCF iteration. /// - /// Defaults to `0` (disabled, lossless). When greater than zero, performs extra hydrogen-focused - /// solves using the same tolerance, then finishes with a full constrained solve to restore total - /// charge conservation (lossy path). - pub hydrogen_inner_iters: u32, + /// Controls how charge updates are mixed between iterations to ensure convergence. + pub damping: DampingStrategy, } impl Default for SolverOptions { fn default() -> Self { Self { tolerance: 1.0e-6, - max_iterations: 100, + max_iterations: 2000, lambda_scale: 0.5, hydrogen_scf: true, - cutoff_radius: None, - hydrogen_inner_iters: 0, + basis_type: BasisType::default(), + damping: DampingStrategy::Auto { initial: 0.4 }, } } } diff --git a/tests/alkali.rs b/tests/alkali.rs new file mode 100644 index 0000000..5228790 --- /dev/null +++ b/tests/alkali.rs @@ -0,0 +1,185 @@ +mod common; + +use cheq::Atom; +use common::{TestCase, run_group_test}; + +fn make_diatomic(z1: u8, z2: u8, dist: f64) -> Vec { + vec![ + Atom { + atomic_number: z1, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: z2, + position: [dist, 0.0, 0.0], + }, + ] +} + +#[test] +fn test_alkali_halides_group() { + let cases = vec![ + TestCase { + name: "LiF", + atoms: make_diatomic(3, 9, 1.564), + expected: vec![(0, 0.791)], + }, + TestCase { + name: "LiCl", + atoms: make_diatomic(3, 17, 2.021), + expected: vec![(0, 0.939)], + }, + TestCase { + name: "LiBr", + atoms: make_diatomic(3, 35, 2.170), + expected: vec![(0, 0.902)], + }, + TestCase { + name: "LiI", + atoms: make_diatomic(3, 53, 2.392), + expected: vec![(0, 0.841)], + }, + TestCase { + name: "NaF", + atoms: make_diatomic(11, 9, 1.926), + expected: vec![(0, 0.665)], + }, + TestCase { + name: "NaCl", + atoms: make_diatomic(11, 17, 2.361), + expected: vec![(0, 0.766)], + }, + TestCase { + name: "NaBr", + atoms: make_diatomic(11, 35, 2.502), + expected: vec![(0, 0.745)], + }, + TestCase { + name: "NaI", + atoms: make_diatomic(11, 53, 2.711), + expected: vec![(0, 0.709)], + }, + TestCase { + name: "KF", + atoms: make_diatomic(19, 9, 2.171), + expected: vec![(0, 0.662)], + }, + TestCase { + name: "KCl", + atoms: make_diatomic(19, 17, 2.667), + expected: vec![(0, 0.775)], + }, + TestCase { + name: "KBr", + atoms: make_diatomic(19, 35, 2.821), + expected: vec![(0, 0.768)], + }, + TestCase { + name: "KI", + atoms: make_diatomic(19, 53, 3.048), + expected: vec![(0, 0.754)], + }, + TestCase { + name: "RbF", + atoms: make_diatomic(37, 9, 2.270), + expected: vec![(0, 0.653)], + }, + TestCase { + name: "RbCl", + atoms: make_diatomic(37, 17, 2.787), + expected: vec![(0, 0.763)], + }, + TestCase { + name: "RbBr", + atoms: make_diatomic(37, 35, 2.945), + expected: vec![(0, 0.757)], + }, + TestCase { + name: "RbI", + atoms: make_diatomic(37, 53, 3.177), + expected: vec![(0, 0.747)], + }, + TestCase { + name: "CsF", + atoms: make_diatomic(55, 9, 2.345), + expected: vec![(0, 0.655)], + }, + TestCase { + name: "CsCl", + atoms: make_diatomic(55, 17, 2.906), + expected: vec![(0, 0.769)], + }, + TestCase { + name: "CsBr", + atoms: make_diatomic(55, 35, 3.072), + expected: vec![(0, 0.767)], + }, + TestCase { + name: "CsI", + atoms: make_diatomic(55, 53, 3.315), + expected: vec![(0, 0.763)], + }, + TestCase { + name: "CO2", + atoms: vec![ + Atom { + atomic_number: 6, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: 8, + position: [1.160, 0.0, 0.0], + }, + Atom { + atomic_number: 8, + position: [-1.160, 0.0, 0.0], + }, + ], + expected: vec![(1, -0.45)], + }, + TestCase { + name: "Ketene", + atoms: { + let r_co = 1.160; + let r_cc = 1.314; + let r_ch = 1.079; + let angle_hch = 122.0f64.to_radians(); + let half_angle = angle_hch / 2.0; + + vec![ + Atom { + atomic_number: 6, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: 8, + position: [r_co, 0.0, 0.0], + }, + Atom { + atomic_number: 6, + position: [-r_cc, 0.0, 0.0], + }, + Atom { + atomic_number: 1, + position: [ + -r_cc - r_ch * half_angle.cos(), + r_ch * half_angle.sin(), + 0.0, + ], + }, + Atom { + atomic_number: 1, + position: [ + -r_cc - r_ch * half_angle.cos(), + -r_ch * half_angle.sin(), + 0.0, + ], + }, + ] + }, + expected: vec![(1, -0.45), (0, 0.42), (2, -0.23)], + }, + ]; + + run_group_test("Alkali Halides & Others", cases, 0.01, 0.20); +} diff --git a/tests/clusters.rs b/tests/clusters.rs new file mode 100644 index 0000000..0e8b32a --- /dev/null +++ b/tests/clusters.rs @@ -0,0 +1,100 @@ +mod common; + +use cheq::Atom; +use common::{TestCase, run_group_test}; + +fn make_nacl_monomer(dist: f64) -> Vec { + vec![ + Atom { + atomic_number: 11, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: 17, + position: [dist, 0.0, 0.0], + }, + ] +} + +fn make_nacl_dimer_square(dist: f64) -> Vec { + vec![ + Atom { + atomic_number: 11, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: 17, + position: [dist, 0.0, 0.0], + }, + Atom { + atomic_number: 17, + position: [0.0, dist, 0.0], + }, + Atom { + atomic_number: 11, + position: [dist, dist, 0.0], + }, + ] +} + +fn make_nacl_tetramer_cube(dist: f64) -> Vec { + vec![ + Atom { + atomic_number: 11, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: 17, + position: [dist, 0.0, 0.0], + }, + Atom { + atomic_number: 17, + position: [0.0, dist, 0.0], + }, + Atom { + atomic_number: 11, + position: [dist, dist, 0.0], + }, + Atom { + atomic_number: 17, + position: [0.0, 0.0, dist], + }, + Atom { + atomic_number: 11, + position: [dist, 0.0, dist], + }, + Atom { + atomic_number: 11, + position: [0.0, dist, dist], + }, + Atom { + atomic_number: 17, + position: [dist, dist, dist], + }, + ] +} + +#[test] +fn test_clusters_group() { + let d = 2.361; + + let cases = vec![ + TestCase { + name: "NaCl Monomer", + atoms: make_nacl_monomer(d), + expected: vec![(0, 0.749), (1, -0.749)], + }, + TestCase { + name: "(NaCl)2 Square", + atoms: make_nacl_dimer_square(d), + expected: vec![(0, 0.826), (1, -0.826)], + }, + TestCase { + name: "(NaCl)4 Cube", + atoms: make_nacl_tetramer_cube(d), + expected: vec![(0, 0.823), (4, -0.823)], + }, + ]; + + run_group_test("Ionic Clusters", cases, 0.05, 0.10); +} diff --git a/tests/common/mod.rs b/tests/common/mod.rs new file mode 100644 index 0000000..abacacf --- /dev/null +++ b/tests/common/mod.rs @@ -0,0 +1,93 @@ +use cheq::{Atom, QEqSolver, get_default_parameters}; + +pub struct TestCase<'a> { + pub name: &'a str, + pub atoms: Vec, + pub expected: Vec<(usize, f64)>, +} + +pub fn run_group_test( + group_name: &str, + cases: Vec, + group_avg_limit: f64, + group_max_limit: f64, +) { + let params = get_default_parameters(); + let solver = QEqSolver::new(params); + + let mut group_total_error = 0.0; + let mut group_max_error = 0.0; + let mut total_data_points = 0; + + println!("\nRunning Group Test: {}", group_name); + println!("{:-<80}", ""); + println!( + "{:<25} | {:<10} | {:<10} | {:<10}", + "Molecule", "Atom Idx", "Expected", "Calculated" + ); + + for case in cases { + let result = match solver.solve(&case.atoms, 0.0) { + Ok(res) => res, + Err(e) => { + println!( + "{:<25} | {:<10} | {:<10} | ERROR: {:?}", + case.name, "-", "-", e + ); + group_total_error += 100.0; + group_max_error = 100.0; + total_data_points += 1; + continue; + } + }; + + for (index, expected_q) in &case.expected { + let calculated_q = result.charges[*index]; + let error = (calculated_q - expected_q).abs(); + + println!( + "{:<25} | {:<10} | {:<10.4} | {:<10.4} (Err: {:.4})", + case.name, index, expected_q, calculated_q, error + ); + + group_total_error += error; + if error > group_max_error { + group_max_error = error; + } + total_data_points += 1; + } + } + + let group_avg_error = if total_data_points > 0 { + group_total_error / total_data_points as f64 + } else { + 0.0 + }; + + println!("{:-<80}", ""); + println!("Group Statistics for '{}':", group_name); + println!(" Total Data Points: {}", total_data_points); + println!( + " Group Avg Error: {:.4} (Limit: {:.4})", + group_avg_error, group_avg_limit + ); + println!( + " Group Max Error: {:.4} (Limit: {:.4})", + group_max_error, group_max_limit + ); + println!("{:-<80}\n", ""); + + assert!( + group_avg_error <= group_avg_limit, + "Group average error {:.4} exceeds limit {:.4}", + group_avg_error, + group_avg_limit + ); + + assert!( + group_max_error <= group_max_limit, + "Group maximum error {:.4} exceeds limit {:.4}", + group_max_error, + group_max_limit + ); +} diff --git a/tests/hydrogen.rs b/tests/hydrogen.rs new file mode 100644 index 0000000..047d9c3 --- /dev/null +++ b/tests/hydrogen.rs @@ -0,0 +1,435 @@ +mod common; + +use cheq::Atom; +use common::{TestCase, run_group_test}; +use std::f64::consts::PI; + +fn make_diatomic(z1: u8, z2: u8, dist: f64) -> Vec { + vec![ + Atom { + atomic_number: z1, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: z2, + position: [dist, 0.0, 0.0], + }, + ] +} + +fn make_water_like(z_center: u8, z_outer: u8, r: f64, angle_deg: f64) -> Vec { + let angle_rad = angle_deg.to_radians(); + let half_angle = angle_rad / 2.0; + vec![ + Atom { + atomic_number: z_center, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: z_outer, + position: [r * half_angle.cos(), r * half_angle.sin(), 0.0], + }, + Atom { + atomic_number: z_outer, + position: [r * half_angle.cos(), -r * half_angle.sin(), 0.0], + }, + ] +} + +fn make_pyramidal_xy3(z_center: u8, z_outer: u8, r: f64, angle_hxh_deg: f64) -> Vec { + let theta_rad = angle_hxh_deg.to_radians(); + let d_hh = 2.0 * r * (theta_rad / 2.0).sin(); + let big_r = d_hh / 3.0f64.sqrt(); + let z_h = (r.powi(2) - big_r.powi(2)).sqrt(); + + let z_coord = if angle_hxh_deg >= 119.99 { 0.0 } else { -z_h }; + + vec![ + Atom { + atomic_number: z_center, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: z_outer, + position: [big_r, 0.0, z_coord], + }, + Atom { + atomic_number: z_outer, + position: [-big_r * 0.5, big_r * 3.0f64.sqrt() * 0.5, z_coord], + }, + Atom { + atomic_number: z_outer, + position: [-big_r * 0.5, -big_r * 3.0f64.sqrt() * 0.5, z_coord], + }, + ] +} + +fn make_tetrahedral_xy4(z_center: u8, z_outer: u8, r: f64) -> Vec { + let s = r / 3.0f64.sqrt(); + vec![ + Atom { + atomic_number: z_center, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: z_outer, + position: [s, s, s], + }, + Atom { + atomic_number: z_outer, + position: [s, -s, -s], + }, + Atom { + atomic_number: z_outer, + position: [-s, s, -s], + }, + Atom { + atomic_number: z_outer, + position: [-s, -s, s], + }, + ] +} + +#[test] +fn test_hydrogen_molecules_group() { + let cases = vec![ + TestCase { + name: "HF", + atoms: make_diatomic(1, 9, 0.917), + expected: vec![(0, 0.46)], + }, + TestCase { + name: "H2O", + atoms: make_water_like(8, 1, 0.958, 104.5), + expected: vec![(1, 0.35)], + }, + TestCase { + name: "NH3", + atoms: make_pyramidal_xy3(7, 1, 1.012, 106.7), + expected: vec![(1, 0.24)], + }, + TestCase { + name: "CH4", + atoms: make_tetrahedral_xy4(6, 1, 1.087), + expected: vec![(1, 0.15)], + }, + TestCase { + name: "C2H6", + atoms: { + let r_cc = 1.535; + let r_ch = 1.094; + let angle_hcc = 111.2f64.to_radians(); + + let z_h1 = -r_ch * angle_hcc.cos(); + let r_h1 = r_ch * angle_hcc.sin(); + + let z_h2 = r_cc + r_ch * angle_hcc.cos(); + + let mut atoms = vec![ + Atom { + atomic_number: 6, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: 6, + position: [0.0, 0.0, r_cc], + }, + ]; + + for i in 0..3 { + let phi = (i as f64 * 120.0).to_radians(); + atoms.push(Atom { + atomic_number: 1, + position: [r_h1 * phi.cos(), r_h1 * phi.sin(), z_h1], + }); + } + + for i in 0..3 { + let phi = (i as f64 * 120.0 + 60.0).to_radians(); + atoms.push(Atom { + atomic_number: 1, + position: [r_h1 * phi.cos(), r_h1 * phi.sin(), z_h2], + }); + } + atoms + }, + expected: vec![(2, 0.16)], + }, + TestCase { + name: "C2H2", + atoms: vec![ + Atom { + atomic_number: 6, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: 6, + position: [1.203, 0.0, 0.0], + }, + Atom { + atomic_number: 1, + position: [-1.061, 0.0, 0.0], + }, + Atom { + atomic_number: 1, + position: [1.203 + 1.061, 0.0, 0.0], + }, + ], + expected: vec![(2, 0.13)], + }, + TestCase { + name: "C2H4", + atoms: { + let r_cc = 1.339; + let r_ch = 1.086; + let _angle_hch = 117.6f64.to_radians(); + let angle_cch = 121.2f64.to_radians(); + + vec![ + Atom { + atomic_number: 6, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: 6, + position: [r_cc, 0.0, 0.0], + }, + Atom { + atomic_number: 1, + position: [-r_ch * angle_cch.cos(), r_ch * angle_cch.sin(), 0.0], + }, + Atom { + atomic_number: 1, + position: [-r_ch * angle_cch.cos(), -r_ch * angle_cch.sin(), 0.0], + }, + Atom { + atomic_number: 1, + position: [r_cc + r_ch * angle_cch.cos(), r_ch * angle_cch.sin(), 0.0], + }, + Atom { + atomic_number: 1, + position: [r_cc + r_ch * angle_cch.cos(), -r_ch * angle_cch.sin(), 0.0], + }, + ] + }, + expected: vec![(2, 0.15)], + }, + TestCase { + name: "C6H6", + atoms: { + let r_cc = 1.397; + let r_ch = 1.084; + let mut atoms = Vec::new(); + for i in 0..6 { + let phi = (i as f64 * 60.0).to_radians(); + atoms.push(Atom { + atomic_number: 6, + position: [r_cc * phi.cos(), r_cc * phi.sin(), 0.0], + }); + atoms.push(Atom { + atomic_number: 1, + position: [(r_cc + r_ch) * phi.cos(), (r_cc + r_ch) * phi.sin(), 0.0], + }); + } + atoms + }, + expected: vec![(1, 0.10)], + }, + TestCase { + name: "H2CO", + atoms: { + let r_co = 1.205; + let r_ch = 1.111; + let _angle_hch = 116.5f64.to_radians(); + let angle_och = 121.75f64.to_radians(); + + vec![ + Atom { + atomic_number: 6, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: 8, + position: [r_co, 0.0, 0.0], + }, + Atom { + atomic_number: 1, + position: [r_ch * angle_och.cos(), r_ch * angle_och.sin(), 0.0], + }, + Atom { + atomic_number: 1, + position: [r_ch * angle_och.cos(), -r_ch * angle_och.sin(), 0.0], + }, + ] + }, + expected: vec![(1, -0.43), (0, 0.19), (2, 0.12)], + }, + TestCase { + name: "H3COH", + atoms: { + let r_co = 1.427; + let r_oh = 0.956; + let r_ch_trans = 1.096; + let r_ch_gauche = 1.094; + let angle_coh = 108.9f64.to_radians(); + + let h_o_x = r_co + r_oh * (PI - angle_coh).cos(); + let h_o_y = r_oh * (PI - angle_coh).sin(); + + let angle_hco = 109.5f64.to_radians(); + + let h_trans_x = r_ch_trans * angle_hco.cos(); + let h_trans_y = -r_ch_trans * angle_hco.sin(); + + let cos120 = -0.5f64; + let sin120 = 3.0f64.sqrt() / 2.0; + + let h_g1_x = r_ch_gauche * angle_hco.cos(); + let h_g1_y = (-r_ch_gauche * angle_hco.sin()) * cos120; + let h_g1_z = (-r_ch_gauche * angle_hco.sin()) * sin120; + + let h_g2_x = r_ch_gauche * angle_hco.cos(); + let h_g2_y = (-r_ch_gauche * angle_hco.sin()) * cos120; + let h_g2_z = -h_g1_z; + + vec![ + Atom { + atomic_number: 6, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: 8, + position: [r_co, 0.0, 0.0], + }, + Atom { + atomic_number: 1, + position: [h_o_x, h_o_y, 0.0], + }, + Atom { + atomic_number: 1, + position: [h_trans_x, h_trans_y, 0.0], + }, + Atom { + atomic_number: 1, + position: [h_g1_x, h_g1_y, h_g1_z], + }, + Atom { + atomic_number: 1, + position: [h_g2_x, h_g2_y, h_g2_z], + }, + ] + }, + expected: vec![(2, 0.36), (1, -0.66), (0, -0.15), (4, 0.14), (3, 0.18)], + }, + TestCase { + name: "HOC(O)H", + atoms: { + let r_co_double = 1.202; + let r_co_single = 1.343; + let r_oh = 0.972; + let r_ch = 1.097; + let angle_oco = 125.0f64.to_radians(); + let angle_hco_double = 124.1f64.to_radians(); + + let o_double_pos = [r_co_double, 0.0, 0.0]; + + let o_single_pos = [ + r_co_single * angle_oco.cos(), + r_co_single * angle_oco.sin(), + 0.0, + ]; + + let h_c_pos = [ + r_ch * angle_hco_double.cos(), + -r_ch * angle_hco_double.sin(), + 0.0, + ]; + + let angle_oh_global = (125.0f64 + 180.0 - 106.3).to_radians(); + let h_o_pos = [ + o_single_pos[0] + r_oh * angle_oh_global.cos(), + o_single_pos[1] + r_oh * angle_oh_global.sin(), + 0.0, + ]; + + vec![ + Atom { + atomic_number: 6, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: 8, + position: o_double_pos, + }, + Atom { + atomic_number: 8, + position: o_single_pos, + }, + Atom { + atomic_number: 1, + position: h_c_pos, + }, + Atom { + atomic_number: 1, + position: h_o_pos, + }, + ] + }, + expected: vec![(1, -0.44), (0, 0.56), (3, 0.16), (2, -0.65), (4, 0.38)], + }, + TestCase { + name: "H3CCN", + atoms: { + let r_cc = 1.458; + let r_cn = 1.157; + let r_ch = 1.104; + let angle_hcc = 109.5f64.to_radians(); + + let h_x = r_ch * angle_hcc.cos(); + let h_r = r_ch * angle_hcc.sin(); + + let mut atoms = vec![ + Atom { + atomic_number: 6, + position: [0.0, 0.0, 0.0], + }, + Atom { + atomic_number: 6, + position: [r_cc, 0.0, 0.0], + }, + Atom { + atomic_number: 7, + position: [r_cc + r_cn, 0.0, 0.0], + }, + ]; + + for i in 0..3 { + let phi = (i as f64 * 120.0).to_radians(); + atoms.push(Atom { + atomic_number: 1, + position: [h_x, h_r * phi.cos(), h_r * phi.sin()], + }); + } + atoms + }, + expected: vec![(2, -0.24), (1, 0.22), (0, -0.37), (3, 0.13)], + }, + TestCase { + name: "SiH4", + atoms: make_tetrahedral_xy4(14, 1, 1.480), + expected: vec![(1, 0.13)], + }, + TestCase { + name: "PH3", + atoms: make_pyramidal_xy3(15, 1, 1.421, 93.3), + expected: vec![(1, 0.08)], + }, + TestCase { + name: "HCl", + atoms: make_diatomic(1, 17, 1.275), + expected: vec![(0, 0.32)], + }, + ]; + + run_group_test("Small Molecules with Hydrogen", cases, 0.035, 0.20); +} diff --git a/tests/polymers.rs b/tests/polymers.rs new file mode 100644 index 0000000..ab7073f --- /dev/null +++ b/tests/polymers.rs @@ -0,0 +1,494 @@ +mod common; + +use cheq::Atom; +use common::{TestCase, run_group_test}; + +struct Turtle { + position: [f64; 3], + heading: [f64; 3], + up: [f64; 3], + left: [f64; 3], +} + +impl Turtle { + fn new() -> Self { + Self { + position: [0.0, 0.0, 0.0], + heading: [1.0, 0.0, 0.0], + up: [0.0, 0.0, 1.0], + left: [0.0, 1.0, 0.0], + } + } + + fn forward(&mut self, dist: f64) { + self.position[0] += self.heading[0] * dist; + self.position[1] += self.heading[1] * dist; + self.position[2] += self.heading[2] * dist; + } + + fn roll(&mut self, angle_deg: f64) { + let rad = angle_deg.to_radians(); + let (s, c) = rad.sin_cos(); + + let new_up = [ + self.up[0] * c + self.left[0] * s, + self.up[1] * c + self.left[1] * s, + self.up[2] * c + self.left[2] * s, + ]; + + let new_left = [ + -self.up[0] * s + self.left[0] * c, + -self.up[1] * s + self.left[1] * c, + -self.up[2] * s + self.left[2] * c, + ]; + + self.up = new_up; + self.left = new_left; + self.normalize(); + } + + fn pitch(&mut self, angle_deg: f64) { + let rad = angle_deg.to_radians(); + let (s, c) = rad.sin_cos(); + + let new_heading = [ + self.heading[0] * c + self.up[0] * s, + self.heading[1] * c + self.up[1] * s, + self.heading[2] * c + self.up[2] * s, + ]; + + let new_up = [ + -self.heading[0] * s + self.up[0] * c, + -self.heading[1] * s + self.up[1] * c, + -self.heading[2] * s + self.up[2] * c, + ]; + + self.heading = new_heading; + self.up = new_up; + self.normalize(); + } + + fn normalize(&mut self) { + let mag_h = + (self.heading[0].powi(2) + self.heading[1].powi(2) + self.heading[2].powi(2)).sqrt(); + self.heading.iter_mut().for_each(|x| *x /= mag_h); + + let mag_u = (self.up[0].powi(2) + self.up[1].powi(2) + self.up[2].powi(2)).sqrt(); + self.up.iter_mut().for_each(|x| *x /= mag_u); + + let mag_l = (self.left[0].powi(2) + self.left[1].powi(2) + self.left[2].powi(2)).sqrt(); + self.left.iter_mut().for_each(|x| *x /= mag_l); + } + + fn extend(&mut self, r: f64, theta: f64, phi: f64) { + self.roll(phi); + self.pitch(180.0 - theta); + self.forward(r); + } + + fn save(&self) -> Self { + Self { + position: self.position, + heading: self.heading, + up: self.up, + left: self.left, + } + } + + fn restore(&mut self, state: &Self) { + self.position = state.position; + self.heading = state.heading; + self.up = state.up; + self.left = state.left; + } +} + +fn make_polymer_chain(n_units: usize, is_pvdf: bool) -> Vec { + let mut atoms = Vec::new(); + let mut turtle = Turtle::new(); + + atoms.push(Atom { + atomic_number: 6, + position: turtle.position, + }); + + let mut c_states = Vec::new(); + c_states.push(turtle.save()); + + let r_cc = 1.54; + let angle_ccc = 109.5; + let torsion_trans = 180.0; + + for _ in 1..n_units { + turtle.extend(r_cc, angle_ccc, torsion_trans); + atoms.push(Atom { + atomic_number: 6, + position: turtle.position, + }); + c_states.push(turtle.save()); + } + + for (i, state) in c_states.iter().enumerate() { + let mut t = Turtle::new(); + t.restore(state); + + let is_cf2 = is_pvdf && (i % 2 != 0); + let z_sub = if is_cf2 { 9 } else { 1 }; + let r_sub = if is_cf2 { 1.35 } else { 1.09 }; + + let mut t1 = t.save(); + t1.extend(r_sub, 109.5, 60.0); + atoms.push(Atom { + atomic_number: z_sub, + position: t1.position, + }); + + let mut t2 = t.save(); + t2.extend(r_sub, 109.5, -60.0); + atoms.push(Atom { + atomic_number: z_sub, + position: t2.position, + }); + + if i == 0 { + let mut t_cap = t.save(); + t_cap.extend(1.09, 0.0, 0.0); + atoms.push(Atom { + atomic_number: 1, + position: t_cap.position, + }); + } + if i == n_units - 1 { + let mut t_cap = t.save(); + t_cap.extend(1.09, 109.5, 180.0); + atoms.push(Atom { + atomic_number: 1, + position: t_cap.position, + }); + } + } + + atoms +} + +fn make_ala_his_ala(protonated: bool) -> Vec { + let mut atoms = Vec::new(); + let mut t = Turtle::new(); + + let r_n_ca = 1.46; + let r_ca_c = 1.51; + let r_c_n = 1.33; + let r_ca_cb = 1.54; + let r_c_o = 1.23; + let r_n_h = 1.01; + + let a_n_ca_c = 111.0; + let a_ca_c_n = 116.0; + let a_c_n_ca = 122.0; + + let phi = -135.0; + let psi = 135.0; + let omega = 180.0; + + atoms.push(Atom { + atomic_number: 7, + position: t.position, + }); + + let mut t_h = t.save(); + t_h.roll(180.0); + t_h.pitch(180.0 - 109.5); + t_h.forward(r_n_h); + atoms.push(Atom { + atomic_number: 1, + position: t_h.position, + }); + + let mut t_h2 = t.save(); + t_h2.roll(180.0 + 120.0); + t_h2.pitch(180.0 - 109.5); + t_h2.forward(r_n_h); + atoms.push(Atom { + atomic_number: 1, + position: t_h2.position, + }); + + t.forward(r_n_ca); + atoms.push(Atom { + atomic_number: 6, + position: t.position, + }); + let ca1_state = t.save(); + + let mut t_ha = t.save(); + t_ha.extend(1.09, 109.5, psi + 120.0); + atoms.push(Atom { + atomic_number: 1, + position: t_ha.position, + }); + + let mut t_cb = t.save(); + t_cb.extend(r_ca_cb, 109.5, psi - 120.0); + atoms.push(Atom { + atomic_number: 6, + position: t_cb.position, + }); + let cb_state = t_cb.save(); + for k in 0..3 { + let mut t_m = cb_state.save(); + t_m.extend(1.09, 109.5, k as f64 * 120.0 + 60.0); + atoms.push(Atom { + atomic_number: 1, + position: t_m.position, + }); + } + + t.restore(&ca1_state); + t.extend(r_ca_c, a_n_ca_c, psi); + atoms.push(Atom { + atomic_number: 6, + position: t.position, + }); + + let mut t_o = t.save(); + t_o.extend(r_c_o, 121.0, 0.0); + atoms.push(Atom { + atomic_number: 8, + position: t_o.position, + }); + + t.extend(r_c_n, a_ca_c_n, omega); + atoms.push(Atom { + atomic_number: 7, + position: t.position, + }); + let mut t_hn = t.save(); + t_hn.extend(r_n_h, 120.0, 0.0); + atoms.push(Atom { + atomic_number: 1, + position: t_hn.position, + }); + t.extend(r_n_ca, a_c_n_ca, phi); + atoms.push(Atom { + atomic_number: 6, + position: t.position, + }); + let ca2_state = t.save(); + + let mut t_ha2 = t.save(); + t_ha2.extend(1.09, 109.5, psi + 120.0); + atoms.push(Atom { + atomic_number: 1, + position: t_ha2.position, + }); + + let mut t_cb2 = t.save(); + t_cb2.extend(r_ca_cb, 109.5, psi - 120.0); + atoms.push(Atom { + atomic_number: 6, + position: t_cb2.position, + }); + + let mut t_hb1 = t_cb2.save(); + t_hb1.extend(1.09, 109.5, 120.0); + atoms.push(Atom { + atomic_number: 1, + position: t_hb1.position, + }); + let mut t_hb2 = t_cb2.save(); + t_hb2.extend(1.09, 109.5, 240.0); + atoms.push(Atom { + atomic_number: 1, + position: t_hb2.position, + }); + + t_cb2.extend(1.50, 109.5, -60.0); + atoms.push(Atom { + atomic_number: 6, + position: t_cb2.position, + }); + let cg_state = t_cb2.save(); + + let mut t_ring = cg_state.save(); + t_ring.extend(1.38, 126.0, 90.0); + atoms.push(Atom { + atomic_number: 7, + position: t_ring.position, + }); + let nd1_state = t_ring.save(); + + t_ring.extend(1.32, 108.0, 180.0); + atoms.push(Atom { + atomic_number: 6, + position: t_ring.position, + }); + let ce1_state = t_ring.save(); + + t_ring.extend(1.32, 108.0, 0.0); + atoms.push(Atom { + atomic_number: 7, + position: t_ring.position, + }); + let ne2_state = t_ring.save(); + + t_ring.extend(1.38, 108.0, 0.0); + atoms.push(Atom { + atomic_number: 6, + position: t_ring.position, + }); + let cd2_state = t_ring.save(); + + let mut t_hd2 = cd2_state.save(); + t_hd2.extend(1.08, 126.0, 180.0); + atoms.push(Atom { + atomic_number: 1, + position: t_hd2.position, + }); + + let mut t_he1 = ce1_state.save(); + t_he1.extend(1.08, 126.0, 180.0); + atoms.push(Atom { + atomic_number: 1, + position: t_he1.position, + }); + + let mut t_hd1 = nd1_state.save(); + t_hd1.extend(1.01, 126.0, 0.0); + atoms.push(Atom { + atomic_number: 1, + position: t_hd1.position, + }); + + if protonated { + let mut t_he2 = ne2_state.save(); + t_he2.extend(1.01, 126.0, 180.0); + atoms.push(Atom { + atomic_number: 1, + position: t_he2.position, + }); + } + + t.restore(&ca2_state); + t.extend(r_ca_c, a_n_ca_c, psi); + atoms.push(Atom { + atomic_number: 6, + position: t.position, + }); + let c2_state = t.save(); + + let mut t_o2 = c2_state.save(); + t_o2.extend(r_c_o, 121.0, 0.0); + atoms.push(Atom { + atomic_number: 8, + position: t_o2.position, + }); + + t.extend(r_c_n, a_ca_c_n, omega); + atoms.push(Atom { + atomic_number: 7, + position: t.position, + }); + let mut t_hn3 = t.save(); + t_hn3.extend(r_n_h, 120.0, 0.0); + atoms.push(Atom { + atomic_number: 1, + position: t_hn3.position, + }); + + t.extend(r_n_ca, a_c_n_ca, phi); + atoms.push(Atom { + atomic_number: 6, + position: t.position, + }); + let ca3_state = t.save(); + + let mut t_ha3 = t.save(); + t_ha3.extend(1.09, 109.5, psi + 120.0); + atoms.push(Atom { + atomic_number: 1, + position: t_ha3.position, + }); + + let mut t_cb3 = t.save(); + t_cb3.extend(r_ca_cb, 109.5, psi - 120.0); + atoms.push(Atom { + atomic_number: 6, + position: t_cb3.position, + }); + let cb3_state = t_cb3.save(); + for k in 0..3 { + let mut t_m = cb3_state.save(); + t_m.extend(1.09, 109.5, k as f64 * 120.0 + 60.0); + atoms.push(Atom { + atomic_number: 1, + position: t_m.position, + }); + } + + t.restore(&ca3_state); + t.extend(r_ca_c, a_n_ca_c, psi); + atoms.push(Atom { + atomic_number: 6, + position: t.position, + }); + let c3_state = t.save(); + + let mut t_o3 = c3_state.save(); + t_o3.extend(r_c_o, 121.0, 0.0); + atoms.push(Atom { + atomic_number: 8, + position: t_o3.position, + }); + + let mut t_oh = c3_state.save(); + t_oh.extend(1.34, 115.0, 180.0); + atoms.push(Atom { + atomic_number: 8, + position: t_oh.position, + }); + + let mut t_ho = t_oh.save(); + t_ho.extend(0.97, 109.0, 180.0); + atoms.push(Atom { + atomic_number: 1, + position: t_ho.position, + }); + + atoms +} + +#[test] +fn test_polymers_group() { + let cases = vec![ + TestCase { + name: "PE (5 units)", + atoms: make_polymer_chain(5, false), + expected: vec![(2, -0.284), (10, 0.143)], + }, + TestCase { + name: "PVDF (5 units)", + atoms: make_polymer_chain(5, true), + expected: vec![(2, -0.082), (10, 0.050), (3, 0.776), (12, -0.415)], + }, + TestCase { + name: "Ala-His-Ala (Neutral)", + atoms: make_ala_his_ala(false), + expected: vec![ + (0, -0.81), + (38, -0.45), + (37, -0.46), + (11, -0.46), + (12, 0.27), + (9, 0.47), + (10, -0.51), + ], + }, + TestCase { + name: "Ala-His-Ala (Protonated)", + atoms: make_ala_his_ala(true), + expected: vec![(19, -0.50), (21, -0.50), (25, 0.33), (26, 0.33)], + }, + ]; + + run_group_test("Polymers & Biomolecules", cases, 0.10, 0.35); +}