Create binding to SIMINT

.gitignore

@@ -1,3 +1,4 @@
+**/*.rs.bk
+*.html
 .ccls-cache/
 /target
-**/*.rs.bk

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"

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)
+<!-- [![crates.io](https://img.shields.io/crates/v/rchem.svg)](https://crates.io/crates/rchem) -->
+<!-- [![Docs.rs](https://docs.rs/rchem/badge.svg)](https://docs.rs/rchem) -->
+
+# 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
+```

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

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,

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!");
 }

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<f64, Ix2> {
     let dim = basis_set.cgtos.len();
     let mut mat: Array<f64, _> = 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<f64, Ix2> {
     let dim = basis_set.cgtos.len();
     let mut mat: Array<f64, _> = 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<u64>) -> Array<f64, Ix2> {
     let dim = basis_set.cgtos.len();
     let natoms = atomcoords.len();
@@ -366,20 +376,21 @@ 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,
     )
     // a.norm
     //     * 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<f64, Ix2>) -> (Array<f64, Ix2>, Array<f64, Ix2>) {
     let dim = basis_set.cgtos.len();
     let mut J: Array<f64, _> = Array::zeros((dim, dim));
@@ -423,6 +434,7 @@ pub fn JK_direct(basis_set: &Basis, D: &Array<f64, Ix2>) -> (Array<f64, Ix2>, Ar
     (J, K)
 }
 
+/// Build an explicit rank-4 tensor of all electron repulsion integrals.
 pub fn build_I(basis_set: &Basis) -> Array<f64, Ix4> {
     let dim = basis_set.cgtos.len();
     let mut I: Array<f64, _> = Array::zeros((dim, dim, dim, dim));
@@ -455,6 +467,7 @@ pub fn build_I(basis_set: &Basis) -> Array<f64, Ix4> {
     I
 }
 
+/// Build the Coulomb and exchange matrices from precomputed electron repulsion integrals.
 pub fn JK_inmem(I: &Array<f64, Ix4>, D: &Array<f64, Ix2>) -> (Array<f64, Ix2>, Array<f64, Ix2>) {
     let dim = I.shape()[0];
     let mut J: Array<f64, _> = Array::zeros((dim, dim));

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#"

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::*;