conspire 0.7.7

The Rust interface to conspire.
Documentation
use crate::math::assert::Assert;
use crate::{
    EPSILON,
    math::{Scalar, assert::AssertionError},
    physics::molecular::{
        potential::{Harmonic, Morse},
        single_chain::{ArbitraryPotentialFreelyJointedChain, Ensemble, Thermodynamics},
    },
    units::{BOLTZMANN_CONSTANT, ROOM_TEMPERATURE},
};

const NUM: usize = 333;

#[test]
fn finite_difference() -> Result<(), AssertionError> {
    //
    // Held as the nondimensional stiffness, so the greatest nondimensional
    // force below is the same wherever the chain sits.
    //
    const NONDIMENSIONAL_STIFFNESS: Scalar = 4.403_161_451_317_08e1;
    let a = 1.0;
    let x0 = 1.0;
    let e = NONDIMENSIONAL_STIFFNESS * BOLTZMANN_CONSTANT.value() * ROOM_TEMPERATURE.value()
        / (x0 * x0);
    let eta_max = 0.5 * a * x0 * e / BOLTZMANN_CONSTANT.value() / ROOM_TEMPERATURE.value();
    [Ensemble::Isotensional(ROOM_TEMPERATURE.value())]
        .into_iter()
        .try_for_each(|ensemble| {
            (3..16).into_iter().try_for_each(|number_of_links| {
                let model = ArbitraryPotentialFreelyJointedChain {
                    link_potential: Harmonic {
                        rest_length: x0,
                        stiffness: e,
                    },
                    number_of_links,
                    ensemble,
                };
                (1..NUM)
                    .map(|k| k as Scalar / NUM as Scalar * eta_max)
                    .into_iter()
                    .try_for_each(|mut nondimensional_force| {
                        nondimensional_force += 0.5 * EPSILON;
                        let mut finite_difference_3 = -model
                            .nondimensional_gibbs_free_energy_per_link(nondimensional_force)?;
                        let mut finite_difference_4 =
                            model.nondimensional_extension(nondimensional_force)?;
                        nondimensional_force -= EPSILON;
                        finite_difference_3 -= -model
                            .nondimensional_gibbs_free_energy_per_link(nondimensional_force)?;
                        finite_difference_4 -=
                            model.nondimensional_extension(nondimensional_force)?;
                        nondimensional_force += 0.5 * EPSILON;
                        let nondimensional_extension =
                            model.nondimensional_extension(nondimensional_force)?;
                        let nondimensional_compliance =
                            model.nondimensional_compliance(nondimensional_force)?;
                        Assert::default().eq_within_fd_tol(
                            nondimensional_extension,
                            &(finite_difference_3 / EPSILON),
                        )?;
                        Assert::default().eq_within_fd_tol(
                            nondimensional_compliance,
                            &(finite_difference_4 / EPSILON),
                        )
                    })
            })
        })?;
    [Ensemble::Isotensional(ROOM_TEMPERATURE.value())]
        .into_iter()
        .try_for_each(|ensemble| {
            (3..16).into_iter().try_for_each(|number_of_links| {
                let model = ArbitraryPotentialFreelyJointedChain {
                    link_potential: Morse {
                        rest_length: x0,
                        depth: e,
                        parameter: a,
                    },
                    number_of_links,
                    ensemble,
                };
                (1..NUM)
                    .map(|k| k as Scalar / NUM as Scalar * eta_max)
                    .into_iter()
                    .try_for_each(|mut nondimensional_force| {
                        nondimensional_force += 0.5 * EPSILON;
                        let mut finite_difference_3 = -model
                            .nondimensional_gibbs_free_energy_per_link(nondimensional_force)?;
                        let mut finite_difference_4 =
                            model.nondimensional_extension(nondimensional_force)?;
                        nondimensional_force -= EPSILON;
                        finite_difference_3 -= -model
                            .nondimensional_gibbs_free_energy_per_link(nondimensional_force)?;
                        finite_difference_4 -=
                            model.nondimensional_extension(nondimensional_force)?;
                        nondimensional_force += 0.5 * EPSILON;
                        let nondimensional_extension =
                            model.nondimensional_extension(nondimensional_force)?;
                        let nondimensional_compliance =
                            model.nondimensional_compliance(nondimensional_force)?;
                        Assert::default().eq_within_fd_tol(
                            nondimensional_extension,
                            &(finite_difference_3 / EPSILON),
                        )?;
                        Assert::default().eq_within_fd_tol(
                            nondimensional_compliance,
                            &(finite_difference_4 / EPSILON),
                        )
                    })
            })
        })
}