coulomb 0.5.0

Library for electrolytes and electrostatic 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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163

// Copyright 2023 BjΓΆrn Stenqvist and Mikael Lund
//
// Converted to Rust with modification from the C++ library "CoulombGalore":
// https://zenodo.org/doi/10.5281/zenodo.3522058
//
// 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.

use super::ShortRangeFunction;
use crate::{Cutoff, Matrix3, NalgebraMatrix3, NalgebraVector3, Vector3};

/// Field due to electric multipoles.
pub trait MultipoleField: ShortRangeFunction + Cutoff {
    /// Returns the electrostatic field from a point charge.
    ///
    /// The `charge` is a source point charge and `r` is the distance vector from the charge.
    /// The returned value has units of `input charge` / `input length squared`, i.e.
    /// _surface charge density_ also known as _charge areal density_. To get the field,
    /// divide by 4πΡ₀.
    ///
    /// The field is obtained with the formula:
    ///
    /// 𝐄(𝑧, 𝐫) = 𝑧𝐫 / π‘ŸΒ²γƒ»{ (1 + πœ…π‘Ÿ)・𝑆(π‘ž) - π‘žπ‘†ΚΉ(π‘ž) }・exp{-πœ…π‘Ÿ} where π‘ž = π‘Ÿ / π‘Ÿβœ‚οΈŽ
    ///
    fn ion_field(&self, charge: f64, r: impl Into<Vector3>) -> Vector3 {
        let r: NalgebraVector3 = r.into().into();
        let r2 = r.norm_squared();
        if r2 >= self.cutoff_squared() {
            return NalgebraVector3::zeros().into();
        }
        let r1 = r.norm();
        let q = r1 / self.cutoff();
        let srf0 = self.short_range_f0(q);
        let srf1 = self.short_range_f1(q);
        let result = charge * r / (r2 * r1)
            * match self.kappa() {
                Some(kappa) => ((1.0 + kappa * r1) * srf0 - q * srf1) * (-kappa * r1).exp(),
                None => srf0 - q * srf1,
            };
        result.into()
    }

    /// Returns the electrostatic field scalar from a point charge.
    ///
    /// The `charge` is a source point charge and `r` is the distance vector from the charge.
    ///
    /// 𝐸(𝑧, r) = 𝑧 / π‘ŸΒ²γƒ»{ (1 + πœ…π‘Ÿ)・𝑆(π‘ž) - π‘žπ‘†ΚΉ(π‘ž) }・exp{-πœ…π‘Ÿ} where π‘ž = π‘Ÿ / π‘Ÿβœ‚οΈŽ
    ///
    fn ion_field_scalar(&self, charge: f64, r: f64) -> f64 {
        if r >= self.cutoff() {
            return 0.0;
        }
        let q = r / self.cutoff();
        let srf0 = self.short_range_f0(q);
        let srf1 = self.short_range_f1(q);
        charge / r.powi(2)
            * if let Some(kappa) = self.kappa() {
                ((1.0 + kappa * r) * srf0 - q * srf1) * (-kappa * r).exp()
            } else {
                srf0 - q * srf1
            }
    }

    /// Electrostatic field from point dipole.
    ///
    /// Parameters:
    /// - `dipole`: Point dipole (input length x input charge) [UNIT: (input length) x (input charge)]
    /// - `r`: Distance vector from point dipole (input length) [UNIT: input length]
    ///
    /// Returns:
    /// Field from dipole [UNIT: (input charge) / (input length)^2]
    ///
    /// The field from a point dipole is described by the formula:
    /// E(mu, r) = (3 * (mu.dot(r) * r / rΒ²) - mu) / r3 *
    ///             (s(q) - q * s'(q) + qΒ² / 3 * s''(q)) +
    ///             mu / r3 * (s(q) * πœ…rΒ² - 2 * πœ…r * q * s'(q) + qΒ² / 3 * s''(q))
    fn dipole_field(&self, dipole: impl Into<Vector3>, r: impl Into<Vector3>) -> Vector3 {
        let dipole: NalgebraVector3 = dipole.into().into();
        let r: NalgebraVector3 = r.into().into();
        let r2 = r.norm_squared();
        if r2 >= self.cutoff_squared() {
            return NalgebraVector3::zeros().into();
        }
        let r1 = r.norm();
        let r3_inv = (r1 * r2).recip();
        let q = r1 / self.cutoff();
        let srf0 = self.short_range_f0(q);
        let srf1 = self.short_range_f1(q);
        let srf2 = self.short_range_f2(q);
        let mut field = (3.0 * dipole.dot(&r) * r / r2 - dipole) * r3_inv;

        let result = if let Some(kappa) = self.kappa() {
            let kr = kappa * r1;
            let kr2 = kr * kr;
            field *= srf0 * (1.0 + kr + kr2 / 3.0) - q * srf1 * (1.0 + 2.0 / 3.0 * kr)
                + q * q / 3.0 * srf2;
            let field_i = dipole * r3_inv * (srf0 * kr2 - 2.0 * kr * q * srf1 + srf2 * q * q) / 3.0;
            (field + field_i) * (-kr).exp()
        } else {
            field *= srf0 - q * srf1 + q * q / 3.0 * srf2;
            let field_i = dipole * r3_inv * q * q * srf2 / 3.0;
            field + field_i
        };
        result.into()
    }
    /// Electrostatic field from point quadrupole.
    ///
    /// Parameters:
    /// - `quad`: Point quadrupole (input length^2 x input charge) [UNIT: (input length)^2 x (input charge)]
    /// - `r`: Distance vector from point quadrupole (input length) [UNIT: input length]
    ///
    /// Returns:
    /// Field from quadrupole [UNIT: (input charge) / (input length)^2]
    fn quadrupole_field(&self, quad: impl Into<Matrix3>, r: impl Into<Vector3>) -> Vector3 {
        let quad: NalgebraMatrix3 = quad.into().into();
        let r: NalgebraVector3 = r.into().into();
        let r2 = r.norm_squared();
        if r2 >= self.cutoff_squared() {
            return NalgebraVector3::zeros().into();
        }
        let r1 = r.norm();
        let r_hat = r / r1;
        let q = r1 / self.cutoff();
        let q2 = q * q;
        let r4 = r2 * r2;
        let quadrh = quad * r_hat;
        let quad_trh = quad.transpose() * r_hat;

        let s0 = self.short_range_f0(q);
        let s1 = self.short_range_f1(q);
        let s2 = self.short_range_f2(q);
        let s3 = self.short_range_f3(q);

        let f = (1.0 / r2 * r.transpose() * quad * r)[0]; // 1x1 matrix -> f64 by taking first and only element
        let mut field_d = 3.0 * ((5.0 * f - quad.trace()) * r_hat - quadrh - quad_trh) / r4;

        let result = if let Some(kappa) = self.kappa() {
            let kr = kappa * r1;
            let kr2 = kr * kr;
            field_d *=
                s0 * (1.0 + kr + kr2 / 3.0) - q * s1 * (1.0 + 2.0 / 3.0 * kr) + q2 / 3.0 * s2;
            let field_i = f * r_hat / r4
                * (s0 * (1.0 + kr) * kr2 - q * s1 * (3.0 * kr + 2.0) * kr
                    + s2 * (1.0 + 3.0 * kr) * q2
                    - q2 * q * s3);
            0.5 * (field_d + field_i) * (-kr).exp()
        } else {
            field_d *= s0 - q * s1 + q2 / 3.0 * s2;
            let field_i = f * r_hat / r4 * (s2 * q2 - q2 * q * s3);
            0.5 * (field_d + field_i)
        };
        result.into()
    }
}