diff --git a/.gitignore b/.gitignore index 0c31119..c8ca2b9 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ +**/*.rs.bk +*.html .ccls-cache/ /target -**/*.rs.bk diff --git a/Cargo.lock b/Cargo.lock index df95dfb..04c51a5 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -359,6 +359,7 @@ dependencies = [ [[package]] name = "factorial" version = "0.2.0" +source = "git+https://github.com/berquist/factorial?branch=signed-double-factorial#225e5823c8ecadeb9ed03900cdb31452516c0bab" dependencies = [ "num-traits 0.2.11 (registry+https://github.com/rust-lang/crates.io-index)", ] @@ -708,7 +709,7 @@ dependencies = [ "cmake 0.1.42 (registry+https://github.com/rust-lang/crates.io-index)", "cpython 0.2.1 (registry+https://github.com/rust-lang/crates.io-index)", "criterion 0.2.11 (registry+https://github.com/rust-lang/crates.io-index)", - "factorial 0.2.0", + "factorial 0.2.0 (git+https://github.com/berquist/factorial?branch=signed-double-factorial)", "ndarray 0.12.1 (registry+https://github.com/rust-lang/crates.io-index)", "ndarray-linalg 0.11.1 (registry+https://github.com/rust-lang/crates.io-index)", "num-integer 0.1.42 (registry+https://github.com/rust-lang/crates.io-index)", @@ -977,6 +978,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" "checksum csv-core 0.1.10 (registry+https://github.com/rust-lang/crates.io-index)" = "2b2466559f260f48ad25fe6317b3c8dac77b5bdb5763ac7d9d6103530663bc90" "checksum either 1.5.3 (registry+https://github.com/rust-lang/crates.io-index)" = "bb1f6b1ce1c140482ea30ddd3335fc0024ac7ee112895426e0a629a6c20adfe3" "checksum env_logger 0.7.1 (registry+https://github.com/rust-lang/crates.io-index)" = "44533bbbb3bb3c1fa17d9f2e4e38bbbaf8396ba82193c4cb1b6445d711445d36" +"checksum factorial 0.2.0 (git+https://github.com/berquist/factorial?branch=signed-double-factorial)" = "" "checksum fuchsia-cprng 0.1.1 (registry+https://github.com/rust-lang/crates.io-index)" = "a06f77d526c1a601b7c4cdd98f54b5eaabffc14d5f2f0296febdc7f357c6d3ba" "checksum glob 0.3.0 (registry+https://github.com/rust-lang/crates.io-index)" = "9b919933a397b79c37e33b77bb2aa3dc8eb6e165ad809e58ff75bc7db2e34574" "checksum hermit-abi 0.1.11 (registry+https://github.com/rust-lang/crates.io-index)" = "8a0d737e0f947a1864e93d33fdef4af8445a00d1ed8dc0c8ddb73139ea6abf15" diff --git a/Cargo.toml b/Cargo.toml index ac1aebe..a81f102 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -7,9 +7,8 @@ build = "build.rs" [dependencies] approx = "0.3" -factorial = { version = "0.2", path = "/home/eric/development/forks/rust/factorial" } +factorial = { version = "*", git = "https://github.com/berquist/factorial", branch = "signed-double-factorial" } num-integer = "0.1" -# num-traits = "*" arrayvec = "0.4" GSL = "1.1" cpython = "0.2" diff --git a/README.md b/README.md new file mode 100644 index 0000000..092c41c --- /dev/null +++ b/README.md @@ -0,0 +1,30 @@ +[![Build Status](https://travis-ci.com/berquist/rchem.svg?branch=master)](https://travis-ci.com/berquist/rchem) +[![License](https://img.shields.io/github/license/berquist/rchem.svg)](LICENSE) +[![Issues](https://img.shields.io/github/issues/berquist/rchem.svg)](https://github.com/berquist/rchem/issues) + + + +# rchem + +_Ab initio_ quantum chemistry in Rust. + +## Dependencies + +- Rust toolchain (stable) +- CMake and a C compiler (for `pyquante2` molecular integrals) +- Python and [`basis-set-exchange`](https://pypi.org/project/basis-set-exchange/) (for reading basis sets) + +## Running + +To run the simple restricted Hartree-Fock program, +``` sh +$ cargo run +``` +To run benchmarks for the different molecular integral implementations, +``` sh +$ cargo bench +``` +To view the library documentation, +``` sh +$ cargo doc +``` diff --git a/TODO.md b/TODO.md deleted file mode 100644 index 82c9c2a..0000000 --- a/TODO.md +++ /dev/null @@ -1,22 +0,0 @@ -# quantum chemistry calculation workflow - -Requirements: -- start with molecular coordinates -- calculate quantum mechanical wavefunction - - -1. Input molecular coordinates from file -2. Gather basis set from file and place on atoms -3. Calculate overlap (S), kinetic energy (T), nuclear-electron attraction (V) integrals -4. Calculate initial guess by diagonalizing H = T + V -5. Do SCF iterations: integral direct TEIs, solve HF equations... -6. ... -7. Profit! - -# Plan post-2019-04-20 - -- [x] Find crate containing gamma function: `rgsl`? -- [x] Find crate containing n-dimensional array crate w/ eigendecomposition (`nalgebra`, `ndarray`) -- [ ] What format should input files be in (Conf, TOML, JSON? not YAML) -- [x] Read in JSON-formatted basis set definitions from EMSL -- [x] ~~How to store basis sets? Can we wrap https://github.com/MolSSI-BSE/basis_set_exchange (Python)?~~ Yes diff --git a/TODO.org b/TODO.org new file mode 100644 index 0000000..2202492 --- /dev/null +++ b/TODO.org @@ -0,0 +1,8 @@ +* TODO Integrals +** TODO write interface to libint https://github.com/berquist/rchem/issues/1 +** TODO write interface to libcint https://github.com/berquist/rchem/issues/2 +** TODO write interface to simint https://github.com/berquist/rchem/issues/5 +- The generated source is located at https://github.com/loriab/simint +- There is also https://anaconda.org/psi4/simint; for CI, it would be advantageous to download the tarball directly from Anaconda Cloud, unpack it, and link directly, rather than going through conda. Though, in the future it might be better to use conda so it's easier to access =libint= and =libcint=. +* TODO Python side +** TODO Better interaction with Python for =basis_set_exchange= dependency https://github.com/berquist/rchem/issues/3 diff --git a/benches/bench_rchem.rs b/benches/bench_rchem.rs index b80d74d..76d6ece 100644 --- a/benches/bench_rchem.rs +++ b/benches/bench_rchem.rs @@ -10,20 +10,17 @@ fn criterion_benchmark(c: &mut Criterion) { let rb = [0.5, 0.8, -0.2]; let powers = [0, 0, 0, 0, 0, 0]; - c.bench_function( - "tho66::pyquante2::pyquante2_overlap_0_0_0_0_0_0", - move |b| { - b.iter(|| { - rchem::integrals::tho66::pyquante2::pyquante2_overlap( - za, - zb, - &ra, - &rb, - &[0, 0, 0, 0, 0, 0], // TODO - ) - }) - }, - ); + c.bench_function("tho66::pyquante2::overlap_0_0_0_0_0_0", move |b| { + b.iter(|| { + rchem::integrals::tho66::pyquante2::overlap( + za, + zb, + &ra, + &rb, + &[0, 0, 0, 0, 0, 0], // TODO + ) + }) + }); c.bench_function("os86::get_overlap_0_0_0_0_0_0", move |b| { b.iter(|| rchem::integrals::os86::get_overlap(za, zb, &ra, &rb, &powers)) }); @@ -40,20 +37,17 @@ fn criterion_benchmark(c: &mut Criterion) { b.iter(|| rchem::integrals::os86::get_overlap(za, zb, &ra, &rb, &powers)) }); let powers = [2, 2, 2, 2, 2, 2]; - c.bench_function( - "tho66::pyquante2::pyquante2_overlap_2_2_2_2_2_2", - move |b| { - b.iter(|| { - rchem::integrals::tho66::pyquante2::pyquante2_overlap( - za, - zb, - &ra, - &rb, - &[2, 2, 2, 2, 2, 2], // TODO - ) - }) - }, - ); + c.bench_function("tho66::pyquante2::overlap_2_2_2_2_2_2", move |b| { + b.iter(|| { + rchem::integrals::tho66::pyquante2::overlap( + za, + zb, + &ra, + &rb, + &[2, 2, 2, 2, 2, 2], // TODO + ) + }) + }); c.bench_function("os86::get_overlap_2_2_2_2_2_2", move |b| { b.iter(|| rchem::integrals::os86::get_overlap(za, zb, &ra, &rb, &powers)) }); @@ -114,10 +108,10 @@ fn criterion_benchmark(c: &mut Criterion) { let powers: [usize; 12] = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]; c.bench_function( - "tho66::pyquante2::pyquante2_coulomb_repulsion_0_0_0_0_0_0_0_0_0_0_0_0", + "tho66::pyquante2::coulomb_repulsion_0_0_0_0_0_0_0_0_0_0_0_0", move |b| { b.iter(|| { - rchem::integrals::tho66::pyquante2::pyquante2_coulomb_repulsion( + rchem::integrals::tho66::pyquante2::coulomb_repulsion( za, zb, zc, @@ -140,10 +134,10 @@ fn criterion_benchmark(c: &mut Criterion) { }); let powers = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]; c.bench_function( - "tho66::pyquante2::pyquante2_coulomb_repulsion_1_0_0_0_0_0_0_0_0_0_0_0", + "tho66::pyquante2::coulomb_repulsion_1_0_0_0_0_0_0_0_0_0_0_0", move |b| { b.iter(|| { - rchem::integrals::tho66::pyquante2::pyquante2_coulomb_repulsion( + rchem::integrals::tho66::pyquante2::coulomb_repulsion( za, zb, zc, @@ -166,10 +160,10 @@ fn criterion_benchmark(c: &mut Criterion) { }); let powers = [2, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0]; c.bench_function( - "tho66::pyquante2::pyquante2_coulomb_repulsion_2_1_0_1_0_0_1_0_0_0_1_0", + "tho66::pyquante2::coulomb_repulsion_2_1_0_1_0_0_1_0_0_0_1_0", move |b| { b.iter(|| { - rchem::integrals::tho66::pyquante2::pyquante2_coulomb_repulsion( + rchem::integrals::tho66::pyquante2::coulomb_repulsion( za, zb, zc, diff --git a/build.rs b/build.rs index 4e2ea66..60c7eae 100644 --- a/build.rs +++ b/build.rs @@ -13,8 +13,26 @@ fn main() { .header("libpyquante2/cints.h") .parse_callbacks(Box::new(bindgen::CargoCallbacks)) .generate() - .expect("Unable to generate bindings"); + .expect("Unable to generate libpyquante2 bindings"); bindings_libpyquante2 .write_to_file(out_dir.join("bindings_libpyquante2.rs")) - .expect("Couldn't write bindings!"); + .expect("Couldn't write libpyquante2 bindings!"); + + let simint_base = env::var("SIMINT_BASE_DIR").unwrap(); + println!("cargo:rustc-link-search={}/lib", simint_base); + println!("cargo:rustc-link-lib=static=simint"); + let bindings_simint = bindgen::Builder::default() + .header(format!("{}/include/simint/simint.h", simint_base)) + .whitelist_var("SIMINT_.*") + .whitelist_function("simint_.*") + .whitelist_type("simint_.*") + // Wanted for the inlines in ostei/ostei_config.h, but doesn't work. + // .generate_inline_functions(true) + // .clang_arg("-fkeep-inline-functions") + .parse_callbacks(Box::new(bindgen::CargoCallbacks)) + .generate() + .expect("Unable to generate simint bindings"); + bindings_simint + .write_to_file(out_dir.join("bindings_simint.rs")) + .expect("Couldn't write simint bindings!"); } diff --git a/src/basis.rs b/src/basis.rs index 4d34c35..e4d883e 100644 --- a/src/basis.rs +++ b/src/basis.rs @@ -1,3 +1,5 @@ +//! Representation of a basis set from the [Basis Set Exchange](https://github.com/MolSSI-BSE/basis_set_exchange) and functions for returning matrices of molecular integrals. + #![allow(non_snake_case)] use std::collections::HashMap as Map; @@ -172,6 +174,7 @@ impl CGTO { } } +/// A representation of a basis set suitable for computing molecular integrals over. #[derive(Debug)] pub struct Basis { name: String, @@ -220,7 +223,9 @@ fn overlap_pgto(a: &PGTO, b: &PGTO) -> f64 { b.powers[1], b.powers[2], ]; - a.norm * b.norm * integrals::get_overlap(a.exponent, b.exponent, &a.origin, &b.origin, &powers) + a.norm + * b.norm + * integrals::os86::get_overlap(a.exponent, b.exponent, &a.origin, &b.origin, &powers) } fn overlap_cgto_left(a: &CGTO, b: &PGTO) -> f64 { @@ -231,6 +236,7 @@ fn overlap_cgto_left(a: &CGTO, b: &PGTO) -> f64 { .sum() } +/// Construct the overlap matrix over a contracted basis. pub fn S(basis_set: &Basis) -> Array { let dim = basis_set.cgtos.len(); let mut mat: Array = Array::zeros((dim, dim)); @@ -259,7 +265,9 @@ fn kinetic_pgto(a: &PGTO, b: &PGTO) -> f64 { b.powers[1], b.powers[2], ]; - a.norm * b.norm * integrals::get_kinetic(a.exponent, b.exponent, &a.origin, &b.origin, &powers) + a.norm + * b.norm + * integrals::os86::get_kinetic(a.exponent, b.exponent, &a.origin, &b.origin, &powers) } fn kinetic_cgto_left(a: &CGTO, b: &PGTO) -> f64 { @@ -270,6 +278,7 @@ fn kinetic_cgto_left(a: &CGTO, b: &PGTO) -> f64 { .sum() } +/// Construct the kinetic energy matrix over a contracted basis. pub fn T(basis_set: &Basis) -> Array { let dim = basis_set.cgtos.len(); let mut mat: Array = Array::zeros((dim, dim)); @@ -300,7 +309,7 @@ fn nuclear_pgto(a: &PGTO, b: &PGTO, atomcoords: &[f64; 3]) -> f64 { ]; a.norm * b.norm - * integrals::get_nuclear( + * integrals::os86::get_nuclear( a.exponent, b.exponent, &a.origin, &b.origin, atomcoords, &powers, ) } @@ -313,6 +322,7 @@ fn nuclear_cgto_left(a: &CGTO, b: &PGTO, atomcoords: &[f64; 3]) -> f64 { .sum() } +/// Construct the nuclear-electron attraction matrix over a contracted basis. pub fn V(basis_set: &Basis, atomcoords: &[[f64; 3]], atomnos: &Vec) -> Array { let dim = basis_set.cgtos.len(); let natoms = atomcoords.len(); @@ -366,7 +376,7 @@ fn coulomb_pgto(a: &PGTO, b: &PGTO, c: &PGTO, d: &PGTO) -> f64 { d.powers[1] as i32, d.powers[2] as i32, ]; - integrals::tho66::pyquante2::pyquante2_coulomb_repulsion( + integrals::tho66::pyquante2::coulomb_repulsion( a.exponent, b.exponent, c.exponent, d.exponent, &a.origin, &b.origin, &c.origin, &d.origin, a.norm, b.norm, c.norm, d.norm, &powers, ) @@ -374,12 +384,13 @@ fn coulomb_pgto(a: &PGTO, b: &PGTO, c: &PGTO, d: &PGTO) -> f64 { // * b.norm // * c.norm // * d.norm - // * integrals::get_coulomb( + // * integrals::os86::get_coulomb( // a.exponent, b.exponent, c.exponent, d.exponent, &a.origin, &b.origin, &c.origin, // &d.origin, &powers, // ) } +/// Build the Coulomb and exchange matrices on the fly (integral direct). pub fn JK_direct(basis_set: &Basis, D: &Array) -> (Array, Array) { let dim = basis_set.cgtos.len(); let mut J: Array = Array::zeros((dim, dim)); @@ -423,6 +434,7 @@ pub fn JK_direct(basis_set: &Basis, D: &Array) -> (Array, Ar (J, K) } +/// Build an explicit rank-4 tensor of all electron repulsion integrals. pub fn build_I(basis_set: &Basis) -> Array { let dim = basis_set.cgtos.len(); let mut I: Array = Array::zeros((dim, dim, dim, dim)); @@ -455,6 +467,7 @@ pub fn build_I(basis_set: &Basis) -> Array { I } +/// Build the Coulomb and exchange matrices from precomputed electron repulsion integrals. pub fn JK_inmem(I: &Array, D: &Array) -> (Array, Array) { let dim = I.shape()[0]; let mut J: Array = Array::zeros((dim, dim)); diff --git a/src/bin/rchem.rs b/src/bin/rchem.rs index fb119a6..e3b3774 100644 --- a/src/bin/rchem.rs +++ b/src/bin/rchem.rs @@ -7,6 +7,7 @@ use ndarray_linalg::*; use rchem::basis; +/// The main driver. Currently just an example of computing the HF/STO-2G energy of a water molecule. fn main() { // http://www.patorjk.com/software/taag/#p=display&f=3D%20Diagonal&t=rchem let logo = r#" diff --git a/src/integrals/mod.rs b/src/integrals/mod.rs index 888e336..f196866 100644 --- a/src/integrals/mod.rs +++ b/src/integrals/mod.rs @@ -1,5 +1,5 @@ +//! Access to molecular integrals, either via direct implementations or bindings. + pub mod os86; +pub mod simint; pub mod tho66; - -// TODO get rid of this re-export -pub use os86::*; diff --git a/src/integrals/os86.rs b/src/integrals/os86.rs index f648d05..81614b6 100644 --- a/src/integrals/os86.rs +++ b/src/integrals/os86.rs @@ -1,3 +1,14 @@ +//! One- and two-electron integrals based on the Obara--Saika (OS) algorithm. +//! +//! Currently available integral types are: +//! - Two-center one-electron overlap (S) +//! - Two-center one-electron kinetic energy (T) +//! - Two-center one-electron nuclear electron attraction (V) +//! - Two-center one-electron multipole moment (M) +//! - Four-center two-electron repulsion (ERI) +//! +//! Reference: [Obara, S.; Saika, A. Efficient recursive computation of molecular integrals over Cartesian Gaussian functions. _J. Chem. Phys._ **1986**, _84_, 3963-3974.](https://doi.org/10.1063/1.450106) + use std::convert::TryInto; use std::f64::consts::PI; @@ -13,14 +24,22 @@ fn remove_item(v: &mut Vec, item: &T) { #[derive(Clone, PartialEq)] enum X2kind { - S, // overlap - T, // kinetic energy - V, // nuclear-electron attraction - M, // multipole moment - L, // angular momentum (not implemented) - E, // electric field (not implemented) - J, // spin-orbit (not implemented) - FC, // Fermi contact (not implemented) + /// overlap + S, + /// kinetic energy + T, + /// nuclear-electron attraction + V, + /// multipole moment + M, + /// angular momentum (not implemented yet) + L, + /// electric field (not implemented yet) + E, + /// spin-orbit (not implemented yet) + J, + /// Fermi contact (not implemented yet) + FC, } #[derive(Clone)] diff --git a/src/integrals/simint.rs b/src/integrals/simint.rs new file mode 100644 index 0000000..ba7c637 --- /dev/null +++ b/src/integrals/simint.rs @@ -0,0 +1,588 @@ +//! Two-electron integrals based on the Obara--Saika algorithm, provided by [Simint](http://www.bennyp.org/research/simint/). + +#![allow(clippy::unreadable_literal)] +#![allow(dead_code)] +#![allow(non_snake_case)] +include!(concat!(env!("OUT_DIR"), "/bindings_simint.rs")); + +impl simint_shell { + /// Create a new **uninitialized** internal `simint_shell` with + /// mutable null pointers. + fn new() -> simint_shell { + simint_shell { + am: 0, + nprim: 0, + x: 0.0, + y: 0.0, + z: 0.0, + alpha: std::ptr::null_mut(), + coef: std::ptr::null_mut(), + memsize: 0, + ptr: std::ptr::null_mut(), + } + } +} + +#[derive(Clone, Debug)] +struct SimintShell { + am: i32, + nprim: i32, + x: f64, + y: f64, + z: f64, + alpha: Vec, + coef: Vec, +} + +impl SimintShell { + // fn new() -> SimintShell { + // SimintShell { + // am: 0, + // nprim: 0, + // x: 0.0, + // y: 0.0, + // z: 0.0, + // alpha: Vec::new(), + // coef: Vec::new(), + // } + // } + + fn from(s: &simint_shell) -> SimintShell { + if s.am < 0 { + panic!(); + } + if s.nprim < 0 { + panic!(); + } + SimintShell { + am: s.am, + nprim: s.nprim, + x: s.x, + y: s.y, + z: s.z, + alpha: unsafe { std::slice::from_raw_parts(s.alpha, s.nprim as usize).to_vec() }, + coef: unsafe { std::slice::from_raw_parts(s.coef, s.nprim as usize).to_vec() }, + } + // TODO Since I'm not doing anything with `s` outside this + // function, I'm probably leaking memory. + } + + fn to(self) -> simint_shell { + unsafe { + simint_init(); + }; + // Strategy: initialize an empty simint_shell, allocate its + // memory, then copy over this SimintShell's contents. + let mut internal_shell = simint_shell::new(); + let p_internal_shell = &mut internal_shell as *mut simint_shell; + unsafe { + simint_allocate_shell(self.nprim, p_internal_shell); + }; + if internal_shell.alpha.is_null() { + panic!("internal_shell.alpha is null"); + } + if internal_shell.coef.is_null() { + panic!("internal_shell.coef is null"); + } + internal_shell.am = self.am; + internal_shell.nprim = self.nprim; + internal_shell.x = self.x; + internal_shell.y = self.y; + internal_shell.z = self.z; + unsafe { + self.alpha + .as_ptr() + .copy_to(internal_shell.alpha, self.alpha.len()); + self.coef + .as_ptr() + .copy_to(internal_shell.coef, self.coef.len()); + } + internal_shell + } +} + +impl PartialEq for SimintShell { + fn eq(&self, other: &Self) -> bool { + self.am == other.am + && self.nprim == other.nprim + && self.x == other.x + && self.y == other.y + && self.z == other.z + && self.alpha == other.alpha + && self.coef == other.coef + } +} + +fn normalize_shells(shells: Vec) -> Vec { + let mut internal_shells: Vec<_> = shells.iter().map(|shell| shell.clone().to()).collect(); + unsafe { + simint_normalize_shells(internal_shells.len() as i32, internal_shells.as_mut_ptr()); + }; + internal_shells + .iter() + .map(|shell| SimintShell::from(&shell.clone())) + .collect() +} + +impl simint_multi_shellpair { + /// Create a new **uninitialized** internal `simint_multi_shellpair` with + /// mutable null pointers. + fn new() -> simint_multi_shellpair { + simint_multi_shellpair { + am1: 0, + am2: 0, + nprim: 0, + nshell12: 0, + nshell12_clip: 0, + nprim12: std::ptr::null_mut(), + AB_x: std::ptr::null_mut(), + AB_y: std::ptr::null_mut(), + AB_z: std::ptr::null_mut(), + x: std::ptr::null_mut(), + y: std::ptr::null_mut(), + z: std::ptr::null_mut(), + PA_x: std::ptr::null_mut(), + PA_y: std::ptr::null_mut(), + PA_z: std::ptr::null_mut(), + PB_x: std::ptr::null_mut(), + PB_y: std::ptr::null_mut(), + PB_z: std::ptr::null_mut(), + alpha: std::ptr::null_mut(), + prefac: std::ptr::null_mut(), + screen: std::ptr::null_mut(), + screen_max: 0.0, + memsize: 0, + ptr: std::ptr::null_mut(), + } + } +} + +#[derive(Debug)] +struct SimintMultiShellpair { + am1: i32, + am2: i32, + nprim: i32, + nshell12: i32, + nshell12_clip: i32, + nprim12: Vec, + dist_ab_x: Vec, + dist_ab_y: Vec, + dist_ab_z: Vec, + x: Vec, + y: Vec, + z: Vec, + pa_x: Vec, + pa_y: Vec, + pa_z: Vec, + pb_x: Vec, + pb_y: Vec, + pb_z: Vec, + alpha: Vec, + // alpha2: Vec, + // beta2: Vec, + prefac: Vec, + screen_method: SimintScreenMethod, + screen: Vec, + screen_max: f64, +} + +#[derive(Clone, Copy, Debug)] +enum SimintScreenMethod { + None, + Schwarz, + FastSchwarz, +} + +impl SimintScreenMethod { + fn from(sm: u32) -> SimintScreenMethod { + match sm { + SIMINT_SCREEN_NONE => SimintScreenMethod::None, + SIMINT_SCREEN_SCHWARZ => SimintScreenMethod::Schwarz, + SIMINT_SCREEN_FASTSCHWARZ => SimintScreenMethod::FastSchwarz, + _ => panic!(), + } + } + + fn to(sm: SimintScreenMethod) -> u32 { + match sm { + SimintScreenMethod::None => SIMINT_SCREEN_NONE, + SimintScreenMethod::Schwarz => SIMINT_SCREEN_SCHWARZ, + SimintScreenMethod::FastSchwarz => SIMINT_SCREEN_FASTSCHWARZ, + } + } +} + +impl SimintMultiShellpair { + fn from(s: &simint_multi_shellpair, screen_method: SimintScreenMethod) -> SimintMultiShellpair { + // TODO check null pointers + let nprim = s.nprim as usize; + let npair = s.nshell12 as usize; + SimintMultiShellpair { + am1: s.am1, + am2: s.am2, + nprim: s.nprim, + nshell12: s.nshell12, + nshell12_clip: s.nshell12_clip, + nprim12: unsafe { std::slice::from_raw_parts(s.nprim12, npair).to_vec() }, + dist_ab_x: unsafe { std::slice::from_raw_parts(s.AB_x, npair).to_vec() }, + dist_ab_y: unsafe { std::slice::from_raw_parts(s.AB_y, npair).to_vec() }, + dist_ab_z: unsafe { std::slice::from_raw_parts(s.AB_z, npair).to_vec() }, + x: unsafe { std::slice::from_raw_parts(s.x, nprim).to_vec() }, + y: unsafe { std::slice::from_raw_parts(s.y, nprim).to_vec() }, + z: unsafe { std::slice::from_raw_parts(s.z, nprim).to_vec() }, + pa_x: unsafe { std::slice::from_raw_parts(s.PA_x, nprim).to_vec() }, + pa_y: unsafe { std::slice::from_raw_parts(s.PA_y, nprim).to_vec() }, + pa_z: unsafe { std::slice::from_raw_parts(s.PA_z, nprim).to_vec() }, + pb_x: unsafe { std::slice::from_raw_parts(s.PB_x, nprim).to_vec() }, + pb_y: unsafe { std::slice::from_raw_parts(s.PB_y, nprim).to_vec() }, + pb_z: unsafe { std::slice::from_raw_parts(s.PB_z, nprim).to_vec() }, + alpha: unsafe { std::slice::from_raw_parts(s.alpha, nprim).to_vec() }, + // TODO We make the choice to hard-code the assumption that SIMINT + // cannot handle derivatives. + // alpha2: Vec::new(), + // beta2: Vec::new(), + prefac: unsafe { std::slice::from_raw_parts(s.prefac, nprim).to_vec() }, + screen_method, + screen: match screen_method { + SimintScreenMethod::None => Vec::new(), + _ => unsafe { std::slice::from_raw_parts(s.screen, nprim).to_vec() }, + }, + screen_max: s.screen_max, + } + } + + fn from_shells( + a: &Vec, + b: &Vec, + sm: SimintScreenMethod, + ) -> SimintMultiShellpair { + let mut internal_multi_shellpair = simint_multi_shellpair::new(); + let p_internal_multi_shellpair = + &mut internal_multi_shellpair as *mut simint_multi_shellpair; + let a_internal: Vec<_> = a.iter().map(|shell| shell.clone().to()).collect(); + let b_internal: Vec<_> = b.iter().map(|shell| shell.clone().to()).collect(); + unsafe { + simint_initialize_multi_shellpair(p_internal_multi_shellpair); + simint_create_multi_shellpair( + a.len() as i32, + a_internal.as_ptr(), + b.len() as i32, + b_internal.as_ptr(), + p_internal_multi_shellpair, + SimintScreenMethod::to(sm) as i32, + ); + }; + SimintMultiShellpair::from(&internal_multi_shellpair, sm) + } + + fn to(&self) -> simint_multi_shellpair { + unsafe { + simint_init(); + }; + // Strategy: initialize an empty simint_multi_shellpair, allocate its + // memory, then copy over this SimintMultiShellpair's contents. + let mut internal_multi_shellpair = simint_multi_shellpair::new(); + let p_internal_multi_shellpair = + &mut internal_multi_shellpair as *mut simint_multi_shellpair; + unsafe { + simint_initialize_multi_shellpair(p_internal_multi_shellpair); + }; + internal_multi_shellpair.am1 = self.am1; + internal_multi_shellpair.am2 = self.am2; + internal_multi_shellpair.nprim = self.nprim; + internal_multi_shellpair.nshell12 = self.nshell12; + internal_multi_shellpair.nshell12_clip = self.nshell12_clip; + unsafe { + self.nprim12 + .as_ptr() + .copy_to(internal_multi_shellpair.nprim12, self.nprim12.len()); + self.dist_ab_x + .as_ptr() + .copy_to(internal_multi_shellpair.AB_x, self.dist_ab_x.len()); + self.dist_ab_y + .as_ptr() + .copy_to(internal_multi_shellpair.AB_y, self.dist_ab_y.len()); + self.dist_ab_z + .as_ptr() + .copy_to(internal_multi_shellpair.AB_z, self.dist_ab_z.len()); + self.x + .as_ptr() + .copy_to(internal_multi_shellpair.x, self.x.len()); + self.y + .as_ptr() + .copy_to(internal_multi_shellpair.y, self.y.len()); + self.z + .as_ptr() + .copy_to(internal_multi_shellpair.z, self.z.len()); + self.pa_x + .as_ptr() + .copy_to(internal_multi_shellpair.PA_x, self.pa_x.len()); + self.pa_y + .as_ptr() + .copy_to(internal_multi_shellpair.PA_y, self.pa_y.len()); + self.pa_z + .as_ptr() + .copy_to(internal_multi_shellpair.PA_z, self.pa_z.len()); + self.pb_x + .as_ptr() + .copy_to(internal_multi_shellpair.PB_x, self.pb_x.len()); + self.pb_y + .as_ptr() + .copy_to(internal_multi_shellpair.PB_y, self.pb_y.len()); + self.pb_z + .as_ptr() + .copy_to(internal_multi_shellpair.PB_z, self.pb_z.len()); + self.alpha + .as_ptr() + .copy_to(internal_multi_shellpair.alpha, self.alpha.len()); + self.prefac + .as_ptr() + .copy_to(internal_multi_shellpair.prefac, self.prefac.len()); + self.screen + .as_ptr() + .copy_to(internal_multi_shellpair.screen, self.screen.len()); + } + internal_multi_shellpair.screen_max = self.screen_max; + internal_multi_shellpair + } +} + +const fn simd_round(x: usize, len: usize) -> usize { + (x + (len - 1)) & (!(len - 1)) +} + +const SIMD_LEN: usize = SIMINT_SIMD_LEN as usize; +const NSHELL_SIMD: usize = SIMINT_NSHELL_SIMD as usize; +const OSTEI_MAXAM: usize = SIMINT_OSTEI_MAXAM as usize; +const NELEMENTS: [usize; OSTEI_MAXAM + 1] = [ + simd_round(NSHELL_SIMD * 1, SIMD_LEN) + simd_round(0, SIMD_LEN) + SIMD_LEN * 1, + simd_round(NSHELL_SIMD * 81, SIMD_LEN) + simd_round(81, SIMD_LEN) + SIMD_LEN * 149, + simd_round(NSHELL_SIMD * 961, SIMD_LEN) + simd_round(4332, SIMD_LEN) + SIMD_LEN * 2405, + simd_round(NSHELL_SIMD * 5476, SIMD_LEN) + simd_round(57512, SIMD_LEN) + SIMD_LEN * 17273, + simd_round(NSHELL_SIMD * 21025, SIMD_LEN) + simd_round(418905, SIMD_LEN) + SIMD_LEN * 79965, + simd_round(NSHELL_SIMD * 63001, SIMD_LEN) + simd_round(2131331, SIMD_LEN) + SIMD_LEN * 280425, + simd_round(NSHELL_SIMD * 159201, SIMD_LEN) + simd_round(8502725, SIMD_LEN) + SIMD_LEN * 811633, + simd_round(NSHELL_SIMD * 355216, SIMD_LEN) + + simd_round(28410752, SIMD_LEN) + + SIMD_LEN * 2040610, +]; + +fn ostei_worksize(derorder: usize, maxam: usize) -> usize { + assert_eq!(derorder, 0); + NELEMENTS[maxam] +} + +fn ostei_workmem(derorder: usize, maxam: usize) -> usize { + // TODO get this machine's size of C double + ostei_worksize(derorder, maxam) * 8 +} + +// ncomputed = simint_compute_eri(&ss_pair, &ss_pair, 0.0, work, targets+ntotal); +fn compute_eri( + p: &SimintMultiShellpair, + q: &SimintMultiShellpair, + screen_tol: f64, + work: *mut f64, + integrals: *mut f64, +) -> isize { + let internal_p = p.to(); + let internal_q = q.to(); + let pr = &internal_p as *const simint_multi_shellpair; + let qr = &internal_q as *const simint_multi_shellpair; + unsafe { simint_compute_eri(pr, qr, screen_tol, work, integrals) as isize } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_simd_round() { + assert_eq!(simd_round(3, 4), 4); + assert_eq!(simd_round(4, 4), 4); + assert_eq!(simd_round(5, 4), 8); + assert_eq!(simd_round(6, 4), 8); + assert_eq!(simd_round(7, 4), 8); + assert_eq!(simd_round(8, 4), 8); + assert_eq!(simd_round(9, 4), 12); + } + + /// A reimplementation of [simint/examples/example1.c](https://github.com/simint-chem/simint-generator/blob/c589bd70e53bbdde1753df093bea6db2eb63c971/skel/examples/example1.c) + #[test] + fn example1() { + let s_shells = vec![ + SimintShell { + am: 0, + nprim: 3, + x: 0.00000, + y: -0.14322, + z: 0.0, + alpha: vec![130.7093200, 23.8088610, 6.4436083], + coef: vec![0.15432897, 0.53532814, 0.44463454], + }, + SimintShell { + am: 0, + nprim: 3, + x: 0.00000, + y: -0.14322, + z: 0.0, + alpha: vec![5.0331513, 1.1695961, 0.3803890], + coef: vec![-0.09996723, 0.39951283, 0.70011547], + }, + SimintShell { + am: 0, + nprim: 3, + x: 1.63804, + y: 1.13654, + z: 0.0, + alpha: vec![3.42525091, 0.62391373, 0.16885540], + coef: vec![0.15432897, 0.53532814, 0.44463454], + }, + SimintShell { + am: 0, + nprim: 3, + x: -1.63804, + y: 1.13654, + z: 0.0, + alpha: vec![3.42525091, 0.62391373, 0.16885540], + coef: vec![0.15432897, 0.53532814, 0.44463454], + }, + ]; + let p_shells = vec![SimintShell { + am: 1, + nprim: 3, + x: 0.0, + y: -0.14322, + z: 0.0, + alpha: vec![5.0331513, 1.1695961, 0.3803890], + coef: vec![0.15591627, 0.60768372, 0.39195739], + }]; + let s_shells_simint: Vec<_> = s_shells.iter().map(|shell| shell.clone().to()).collect(); + let p_shells_simint: Vec<_> = p_shells.iter().map(|shell| shell.clone().to()).collect(); + let s_shells_roundtrip: Vec<_> = s_shells_simint + .iter() + .map(|shell| SimintShell::from(&shell.clone())) + .collect(); + let p_shells_roundtrip: Vec<_> = p_shells_simint + .iter() + .map(|shell| SimintShell::from(&shell.clone())) + .collect(); + assert_eq!(s_shells, s_shells_roundtrip); + assert_eq!(p_shells, p_shells_roundtrip); + let s_shells_normalized = normalize_shells(s_shells_roundtrip); + let p_shells_normalized = normalize_shells(p_shells_roundtrip); + let ss_pair = SimintMultiShellpair::from_shells( + &s_shells_normalized, + &s_shells_normalized, + SimintScreenMethod::None, + ); + let ps_pair = SimintMultiShellpair::from_shells( + &p_shells_normalized, + &s_shells_normalized, + SimintScreenMethod::None, + ); + let pp_pair = SimintMultiShellpair::from_shells( + &p_shells_normalized, + &p_shells_normalized, + SimintScreenMethod::None, + ); + + let work_layout = std::alloc::Layout::from_size_align( + 1000 * ostei_workmem(0, 1), + SIMINT_SIMD_ALIGN as usize, + ) + .unwrap(); + let targets_layout = std::alloc::Layout::from_size_align( + 1000 * 7 * 7 * 7 * 7 * std::mem::size_of::(), + SIMINT_SIMD_ALIGN as usize, + ) + .unwrap(); + let work = unsafe { std::alloc::alloc(work_layout) as *mut f64 }; + let targets = unsafe { std::alloc::alloc(targets_layout) as *mut f64 }; + let mut ncomputed = 0; + let mut ntotal = 0; + + // ncomputed = compute_eri(&ss_pair, &ss_pair, 0.0, work, unsafe { + // targets.offset(ntotal) + // }); + ncomputed = compute_eri(&ss_pair, &ss_pair, 0.0, work, targets); + println!("Computed {} contracted ssss integrals", ncomputed); + ntotal += ncomputed; + + // ncomputed = compute_eri(&ps_pair, &ss_pair, 0.0, work, unsafe { + // targets.offset(ntotal) + // }); + // ncomputed *= 3; + // println!("Computed {} contracted psss integrals", ncomputed); + // ntotal += ncomputed; + + // ncomputed = compute_eri(&ps_pair, &ps_pair, 0.0, work, unsafe { + // targets.offset(ntotal) + // }); + // ncomputed *= 9; + // println!("Computed {} contracted psps integrals", ncomputed); + // ntotal += ncomputed; + + // ncomputed = compute_eri(&pp_pair, &ps_pair, 0.0, work, unsafe { + // targets.offset(ntotal) + // }); + // ncomputed *= 27; + // println!("Computed {} contracted ppps integrals", ncomputed); + // ntotal += ncomputed; + + // ncomputed = compute_eri(&pp_pair, &pp_pair, 0.0, work, unsafe { + // targets.offset(ntotal) + // }); + // ncomputed *= 81; + // println!("Computed {} contracted pppp integrals", ncomputed); + // ntotal += ncomputed; + + println!("++Computed {} contracted integrals", ntotal); + + unsafe { + std::alloc::dealloc(targets as *mut u8, targets_layout); + std::alloc::dealloc(work as *mut u8, work_layout); + simint_finalize(); + }; + } +} + +/* +==67476== Thread 2 integrals::simin: +==67476== Invalid write of size 8 +==67476== at 0x484055B: memmove (vg_replace_strmem.c:1270) +==67476== by 0x124ADC: core::intrinsics::copy (intrinsics.rs:2063) +==67476== by 0x13EBF6: core::ptr::const_ptr::::copy_to (const_ptr.rs:638) +==67476== by 0x15C833: rchem::integrals::simint::SimintMultiShellpair::to (simint.rs:296) +==67476== by 0x15D1EE: rchem::integrals::simint::compute_eri (simint.rs:388) +==67476== by 0x166216: rchem::integrals::simint::tests::example1 (simint.rs:508) +==67476== by 0x14B789: rchem::integrals::simint::tests::example1::{{closure}} (simint.rs:412) +==67476== by 0x14DEBD: core::ops::function::FnOnce::call_once (function.rs:232) +==67476== by 0xE9B25E: as core::ops::function::FnOnce>::call_once (boxed.rs:1017) +==67476== by 0xEF03C6: __rust_maybe_catch_panic (lib.rs:86) +==67476== by 0xEB66BB: try<(),std::panic::AssertUnwindSafe>>> (panicking.rs:281) +==67476== by 0xEB66BB: catch_unwind>>,()> (panic.rs:394) +==67476== by 0xEB66BB: run_test_in_process (lib.rs:542) +==67476== by 0xEB66BB: test::run_test::run_test_inner::{{closure}} (lib.rs:451) +==67476== by 0xE8E485: std::sys_common::backtrace::__rust_begin_short_backtrace (backtrace.rs:130) +==67476== Address 0x0 is not stack'd, malloc'd or (recently) free'd +==67476== +==67476== +==67476== Process terminating with default action of signal 11 (SIGSEGV): dumping core +==67476== Access not within mapped region at address 0x0 +==67476== at 0x484055B: memmove (vg_replace_strmem.c:1270) +==67476== by 0x124ADC: core::intrinsics::copy (intrinsics.rs:2063) +==67476== by 0x13EBF6: core::ptr::const_ptr::::copy_to (const_ptr.rs:638) +==67476== by 0x15C833: rchem::integrals::simint::SimintMultiShellpair::to (simint.rs:296) +==67476== by 0x15D1EE: rchem::integrals::simint::compute_eri (simint.rs:388) +==67476== by 0x166216: rchem::integrals::simint::tests::example1 (simint.rs:508) +==67476== by 0x14B789: rchem::integrals::simint::tests::example1::{{closure}} (simint.rs:412) +==67476== by 0x14DEBD: core::ops::function::FnOnce::call_once (function.rs:232) +==67476== by 0xE9B25E: as core::ops::function::FnOnce>::call_once (boxed.rs:1017) +==67476== by 0xEF03C6: __rust_maybe_catch_panic (lib.rs:86) +==67476== by 0xEB66BB: try<(),std::panic::AssertUnwindSafe>>> (panicking.rs:281) +==67476== by 0xEB66BB: catch_unwind>>,()> (panic.rs:394) +==67476== by 0xEB66BB: run_test_in_process (lib.rs:542) +==67476== by 0xEB66BB: test::run_test::run_test_inner::{{closure}} (lib.rs:451) +==67476== by 0xE8E485: std::sys_common::backtrace::__rust_begin_short_backtrace (backtrace.rs:130) +*/ diff --git a/src/integrals/tho66.rs b/src/integrals/tho66.rs index 6681c43..40d6d9d 100644 --- a/src/integrals/tho66.rs +++ b/src/integrals/tho66.rs @@ -1,5 +1,15 @@ -#![allow(clippy::too_many_arguments)] -#![allow(clippy::many_single_char_names)] +//! One- and two-electron integrals based on the Taketa--Huzinaga--O-ohata (THO) algorithm. +//! +//! Currently available integral types are: +//! - Two-center one-electron overlap (S) +//! - Four-center two-electron repulsion (ERI) +//! +//! Reference: [Taketa, Hiroshi; Huzinaga, Sigeru; O-ohata, Kiyosi. Gaussian-Expansion Methods for Molecular Integrals. _J. Phys. Soc. Jpn._ **1966**, _21_, 2313-2324.](https://doi.org/10.1143/JPSJ.21.2313) +//! +//! Most of this file is commented out because there is an issue with the pure Rust implementation. The implementation used originates from bindings to the pure C implementation taken from [pyquante2](https://github.com/rpmuller/pyquante2/tree/6e34cb4480ae7dbd8c5e44d221d8b27584890c83/cython), stored in [the crate root](https://github.com/berquist/rchem/tree/276a5d44505ac8b817219426e3eced96d2555950/libpyquante2). + +// #![allow(clippy::too_many_arguments)] +// #![allow(clippy::many_single_char_names)] // use factorial::SignedFactorial; // use num_integer; // use std::convert::TryInto; @@ -223,18 +233,23 @@ // } pub mod pyquante2 { - include!(concat!(env!("OUT_DIR"), "/bindings_libpyquante2.rs")); - pub fn pyquante2_overlap(za: f64, zb: f64, ra: &[f64; 3], rb: &[f64; 3], c: &[i32; 6]) -> f64 { + /// The autogenerated bindings, which have public visibility and should be + /// hidden due to their unsafe nature. + mod bindings { + include!(concat!(env!("OUT_DIR"), "/bindings_libpyquante2.rs")); + } + + pub fn overlap(za: f64, zb: f64, ra: &[f64; 3], rb: &[f64; 3], c: &[i32; 6]) -> f64 { unsafe { - overlap( + bindings::overlap( za, c[0], c[1], c[2], ra[0], ra[1], ra[2], zb, c[3], c[4], c[5], rb[0], rb[1], rb[2], ) } } - pub fn pyquante2_coulomb_repulsion( + pub fn coulomb_repulsion( za: f64, zb: f64, zc: f64, @@ -250,7 +265,7 @@ pub mod pyquante2 { c: &[i32; 12], ) -> f64 { unsafe { - coulomb_repulsion( + bindings::coulomb_repulsion( ra[0], ra[1], ra[2], norma, c[0], c[1], c[2], za, rb[0], rb[1], rb[2], normb, c[3], c[4], c[5], zb, rc[0], rc[1], rc[2], normc, c[6], c[7], c[8], zc, rd[0], rd[1], rd[2], normd, c[9], c[10], c[11], zd, @@ -272,16 +287,16 @@ mod tests { let thresh = 1.0e-16; - let integral = pyquante2::pyquante2_overlap(za, zb, &ra, &rb, &[0, 0, 0, 0, 0, 0]); + let integral = pyquante2::overlap(za, zb, &ra, &rb, &[0, 0, 0, 0, 0, 0]); assert!((integral - 0.20373275913014607).abs() < thresh); - let integral = pyquante2::pyquante2_overlap(za, zb, &ra, &rb, &[1, 0, 0, 0, 0, 0]); + let integral = pyquante2::overlap(za, zb, &ra, &rb, &[1, 0, 0, 0, 0, 0]); assert!((integral - 0.062005622343957505).abs() < thresh); - let integral = pyquante2::pyquante2_overlap(za, zb, &ra, &rb, &[1, 1, 0, 1, 1, 0]); + let integral = pyquante2::overlap(za, zb, &ra, &rb, &[1, 1, 0, 1, 1, 0]); assert!((integral - -0.00043801221837779696).abs() < thresh); - let integral = pyquante2::pyquante2_overlap(za, zb, &ra, &rb, &[2, 1, 0, 1, 1, 0]); + let integral = pyquante2::overlap(za, zb, &ra, &rb, &[2, 1, 0, 1, 1, 0]); assert!((integral - -0.0002385994651113168).abs() < thresh); } @@ -300,7 +315,7 @@ mod tests { let thresh = 1.0e-12; let reference = 0.08608517834596989; - let integral = pyquante2::pyquante2_coulomb_repulsion( + let integral = pyquante2::coulomb_repulsion( za, zb, zc,