interatomic 0.5.0

Library for calculating inter-particle interactions
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129

// Copyright 2023-2024 Mikael Lund
//
// Licensed under the Apache license, version 2.0 (the "license");
// you may not use this file except in compliance with the license.
// You may obtain a copy of the license at
//
//     http://www.apache.org/licenses/license-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the license is distributed on an "as is" basis,
// without warranties or conditions of any kind, either express or implied.
// See the license for the specific language governing permissions and
// limitations under the license.

//! Implementation of the periodic dihedral potential.

use super::FourbodyAngleEnergy;
#[cfg(feature = "serde")]
use serde::{Deserialize, Serialize};

/// Periodic dihedral potential.
#[derive(Debug, Clone, PartialEq)]
#[cfg_attr(
    feature = "serde",
    derive(Deserialize, Serialize),
    serde(deny_unknown_fields)
)]
pub struct PeriodicDihedral {
    #[cfg_attr(feature = "serde", serde(rename = "phi"))]
    phase_angle: f64,
    #[cfg_attr(feature = "serde", serde(rename = "k"))]
    spring_constant: f64,
    #[cfg_attr(feature = "serde", serde(rename = "n"))]
    periodicity: f64,
}

impl PeriodicDihedral {
    /// Create a new periodic dihedral potential with phase angle, spring constant, and periodicity.
    pub const fn new(phase_angle: f64, spring_constant: f64, periodicity: f64) -> Self {
        Self {
            phase_angle,
            spring_constant,
            periodicity,
        }
    }

    /// Phase angle in degrees.
    pub const fn phase_angle(&self) -> f64 {
        self.phase_angle
    }

    /// Spring constant.
    pub const fn spring_constant(&self) -> f64 {
        self.spring_constant
    }

    /// Periodicity (number of minima per full rotation).
    pub const fn periodicity(&self) -> f64 {
        self.periodicity
    }
}

impl FourbodyAngleEnergy for PeriodicDihedral {
    /// Energy: k * (1 + cos(n*θ - φ)), where θ and φ are in degrees.
    #[inline(always)]
    fn fourbody_angle_energy(&self, angle: f64) -> f64 {
        let angle_rad = angle.to_radians();
        let phase_rad = self.phase_angle.to_radians();
        self.spring_constant * (1.0 + (self.periodicity * angle_rad - phase_rad).cos())
    }

    /// Force: -dU/dφ (in degrees). Returns k*n*sin(n*θ - φ₀) * π/180.
    #[inline(always)]
    fn fourbody_angle_force(&self, angle: f64) -> f64 {
        let angle_rad = angle.to_radians();
        let phase_rad = self.phase_angle.to_radians();
        self.spring_constant
            * self.periodicity
            * (self.periodicity * angle_rad - phase_rad).sin()
            * std::f64::consts::PI
            / 180.0
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use approx::assert_relative_eq;

    #[test]
    fn periodic_dihedral_energy() {
        // k=10, n=2, phi=180°, angle=0° → 10*(1+cos(0-π)) = 0
        let pot = PeriodicDihedral::new(180.0, 10.0, 2.0);
        assert_relative_eq!(pot.fourbody_angle_energy(0.0), 0.0, epsilon = 1e-10);

        // k=10, n=2, phi=180°, angle=90° → 10*(1+cos(π-π)) = 20
        assert_relative_eq!(pot.fourbody_angle_energy(90.0), 20.0, epsilon = 1e-10);

        // k=5, n=3, phi=0°, angle=60° → 5*(1+cos(π)) = 0
        let pot2 = PeriodicDihedral::new(0.0, 5.0, 3.0);
        assert_relative_eq!(pot2.fourbody_angle_energy(60.0), 0.0, epsilon = 1e-10);

        // k=5, n=3, phi=0°, angle=0° → 5*(1+cos(0)) = 10
        assert_relative_eq!(pot2.fourbody_angle_energy(0.0), 10.0, epsilon = 1e-10);

        // k=1, n=1, phi=0°, angle=180° → 1*(1+cos(π)) = 0
        let pot3 = PeriodicDihedral::new(0.0, 1.0, 1.0);
        assert_relative_eq!(pot3.fourbody_angle_energy(180.0), 0.0, epsilon = 1e-10);
    }

    #[test]
    fn test_periodic_dihedral_force() {
        let pot = PeriodicDihedral::new(0.0, 100.0, 3.0);
        // Energy cross-checks
        assert_relative_eq!(pot.fourbody_angle_energy(120.0), 200.0, epsilon = 1e-10);
        assert_relative_eq!(pot.fourbody_angle_energy(60.0), 0.0, epsilon = 1e-10);
        assert_relative_eq!(pot.fourbody_angle_energy(90.0), 100.0, epsilon = 1e-10);

        // At φ=90°: -dU/dφ_deg = kn*sin(270°)*π/180 = -300*π/180
        assert_relative_eq!(
            pot.fourbody_angle_force(90.0),
            -300.0 * std::f64::consts::PI / 180.0,
            epsilon = 1e-10
        );
        // At φ=60° and φ=120°: sin(nφ)=sin(180°)=0 and sin(360°)=0 → force=0
        assert_relative_eq!(pot.fourbody_angle_force(60.0), 0.0, epsilon = 1e-10);
        assert_relative_eq!(pot.fourbody_angle_force(120.0), 0.0, epsilon = 1e-10);
    }
}