interatomic/twobody/
wca.rs

1// Copyright 2024 Mikael Lund
2//
3// Licensed under the Apache license, version 2.0 (the "license");
4// you may not use this file except in compliance with the license.
5// You may obtain a copy of the license at
6//
7//     http://www.apache.org/licenses/license-2.0
8//
9// Unless required by applicable law or agreed to in writing, software
10// distributed under the license is distributed on an "as is" basis,
11// without warranties or conditions of any kind, either express or implied.
12// See the license for the specific language governing permissions and
13// limitations under the license.
14
15use crate::twobody::IsotropicTwobodyEnergy;
16use crate::{twobody::LennardJones, Cutoff};
17
18#[cfg(feature = "serde")]
19use serde::{Deserialize, Serialize};
20
21/// Weeks-Chandler-Andersen potential
22///
23/// This is a Lennard-Jones type potential, cut and shifted to zero:
24///
25/// $$u(r) = 4 \epsilon \left [ (\sigma_{ij}/r)^{12} - (\sigma_{ij}/r)^6 + \frac{1}{4} \right ]$$
26///
27/// for $r < r_{cut} = 2^{1/6} \sigma_{ij}$; zero otherwise.
28///
29/// Effectively, this provides soft repulsion without any attraction.
30/// More information, see <https://doi.org/ct4kh9>.
31///
32/// # Examples:
33/// ~~~
34/// use interatomic::twobody::*;
35/// let (epsilon, sigma) = (1.5, 2.0);
36/// let wca = WeeksChandlerAndersen::new(epsilon, sigma);
37/// assert_eq!(wca.isotropic_twobody_energy(sigma * sigma), epsilon);
38/// ~~~
39#[derive(Debug, Clone, PartialEq, Default)]
40#[cfg_attr(
41    feature = "serde",
42    derive(Deserialize, Serialize),
43    serde(deny_unknown_fields)
44)]
45pub struct WeeksChandlerAndersen {
46    #[cfg_attr(feature = "serde", serde(flatten))]
47    lennard_jones: LennardJones,
48}
49
50impl WeeksChandlerAndersen {
51    const ONEFOURTH: f64 = 0.25;
52    pub(crate) const TWOTOTWOSIXTH: f64 = 1.2599210498948732; // f64::powf(2.0, 2.0/6.0)
53
54    /// New from epsilon and sigma
55    pub const fn new(epsilon: f64, sigma: f64) -> Self {
56        Self {
57            lennard_jones: LennardJones::new(epsilon, sigma),
58        }
59    }
60}
61
62impl Cutoff for WeeksChandlerAndersen {
63    #[inline(always)]
64    fn cutoff_squared(&self) -> f64 {
65        self.lennard_jones.sigma_squared * Self::TWOTOTWOSIXTH
66    }
67    #[inline(always)]
68    fn cutoff(&self) -> f64 {
69        self.cutoff_squared().sqrt()
70    }
71    fn lower_cutoff(&self) -> f64 {
72        self.lennard_jones.lower_cutoff()
73    }
74}
75
76/// Conversion from Lennard-Jones
77impl From<LennardJones> for WeeksChandlerAndersen {
78    fn from(lennard_jones: LennardJones) -> Self {
79        Self { lennard_jones }
80    }
81}
82
83impl IsotropicTwobodyEnergy for WeeksChandlerAndersen {
84    #[inline(always)]
85    fn isotropic_twobody_energy(&self, distance_squared: f64) -> f64 {
86        if distance_squared > self.cutoff_squared() {
87            return 0.0;
88        }
89        let x6 = (self.lennard_jones.sigma_squared / distance_squared).powi(3); // (σ/r)^6
90        self.lennard_jones.four_times_epsilon * (x6 * x6 - x6 + Self::ONEFOURTH)
91    }
92
93    #[inline(always)]
94    fn isotropic_twobody_force(&self, distance_squared: f64) -> f64 {
95        if distance_squared > self.cutoff_squared() {
96            return 0.0;
97        }
98        self.lennard_jones.isotropic_twobody_force(distance_squared)
99    }
100}
101
102#[cfg(test)]
103mod tests {
104    use super::*;
105    use crate::twobody::LennardJones;
106    use approx::assert_relative_eq;
107
108    #[test]
109    fn test_wca_force() {
110        let (epsilon, sigma) = (1.5, 2.0);
111        let wca = WeeksChandlerAndersen::new(epsilon, sigma);
112        let lj = LennardJones::new(epsilon, sigma);
113        // Inside cutoff: matches LJ force
114        assert_relative_eq!(
115            wca.isotropic_twobody_force(4.0),
116            lj.isotropic_twobody_force(4.0)
117        );
118        // Outside cutoff: zero
119        assert_relative_eq!(wca.isotropic_twobody_force(6.0), 0.0);
120    }
121}
interatomic/twobody/wca.rs

interatomic/twobody/
wca.rs
