Skip to content

Commit

Permalink
return Vector3f type instead of plain array
Browse files Browse the repository at this point in the history
  • Loading branch information
ybyygu committed Dec 19, 2019
1 parent efcb24d commit b41f805
Show file tree
Hide file tree
Showing 4 changed files with 49 additions and 59 deletions.
66 changes: 29 additions & 37 deletions src/lib.rs
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
// AUTHOR: Wenping Guo <ybyygu@gmail.com>
// LICENCE: GPL version 3
// CREATED: <2018-04-29 14:27>
// UPDATED: <2019-12-18 Wed 21:14>
// UPDATED: <2019-12-19 Thu 09:28>
//===============================================================================#
// header:1 ends here

Expand Down Expand Up @@ -114,25 +114,20 @@ impl Lattice {
}

/// Set cell origin in Cartesian coordinates
pub fn set_origin(&mut self, loc: [f64; 3]) {
self.origin = Vector3f::from(loc);
pub fn set_origin<T: Into<Vector3f>>(&mut self, loc: T) {
self.origin = loc.into()
}

/// Lattice length parameters: a, b, c
pub fn lengths(&self) -> [f64; 3] {
let lengths = get_cell_lengths(self.matrix);

[lengths[0], lengths[1], lengths[2]]
get_cell_lengths(self.matrix).into()
}

/// Lattice angle parameters in degrees
pub fn angles(&self) -> [f64; 3] {
let angles = get_cell_angles(self.matrix);

[angles[0], angles[1], angles[2]]
get_cell_angles(self.matrix).into()
}

// FIXME: cell widths
/// Scale Lattice by a positive constant
pub fn scale_by(&mut self, v: f64) {
debug_assert!(v > 0.);
Expand All @@ -141,43 +136,38 @@ impl Lattice {
}

/// Get cell origin in Cartesian coordinates
pub fn origin(&self) -> [f64; 3] {
self.origin.into()
pub fn origin(&self) -> Vector3f {
self.origin
}

/// Returns the fractional coordinates given cartesian coordinates.
pub fn to_frac(&self, p: [f64; 3]) -> [f64; 3] {
let v = Vector3f::from(p);
let fs = self.inv_matrix * (v - self.origin);
fs.into()
pub fn to_frac<T: Into<Vector3f>>(&self, p: T) -> Vector3f {
self.inv_matrix * (p.into() - self.origin)
}

/// Returns the cartesian coordinates given fractional coordinates.
pub fn to_cart(&self, p: [f64; 3]) -> [f64; 3] {
let v = Vector3f::from(p);
let fs = self.matrix * v + self.origin;

fs.into()
pub fn to_cart<T: Into<Vector3f>>(&self, p: T) -> Vector3f {
self.matrix * p.into() + self.origin
}

/// Lattice vector a
pub fn vector_a(&self) -> [f64; 3] {
self.matrix.column(0).transpose().into()
pub fn vector_a(&self) -> Vector3f {
self.matrix.column(0).into()
}

/// Lattice vector b
pub fn vector_b(&self) -> [f64; 3] {
self.matrix.column(1).transpose().into()
pub fn vector_b(&self) -> Vector3f {
self.matrix.column(1).into()
}

/// Lattice vector c
pub fn vector_c(&self) -> [f64; 3] {
self.matrix.column(2).transpose().into()
pub fn vector_c(&self) -> Vector3f {
self.matrix.column(2).into()
}

/// Lattice vectors
pub fn vectors(&self) -> [[f64; 3]; 3] {
self.matrix.into()
pub fn matrix(&self) -> Matrix3f {
self.matrix
}

/// Check if lattice is orthorhombic
Expand All @@ -187,10 +177,10 @@ impl Lattice {
m == self.matrix
}

/// Wrap a point to unit cell, obeying the periodic boundary conditions.
pub fn wrap(&self, vec: [f64; 3]) -> [f64; 3] {
let [fx, fy, fz] = self.to_frac(vec);
let fcoords_wrapped = [fx - fx.floor(), fy - fy.floor(), fz - fz.floor()];
/// Wrap a point into unit cell, obeying the periodic boundary conditions.
pub fn wrap<T: Into<Vector3f>>(&self, vec: T) -> Vector3f {
let f = self.to_frac(vec);
let fcoords_wrapped = [f.x - f.x.floor(), f.y - f.y.floor(), f.z - f.z.floor()];
self.to_cart(fcoords_wrapped)
}

Expand All @@ -200,13 +190,15 @@ impl Lattice {
/// Parameters
/// ----------
/// * pi, pj: Cartesian coordinates of point i and point j
pub fn distance(&self, pi: [f64; 3], pj: [f64; 3]) -> f64 {
let pmic = self.apply_mic([pj[0] - pi[0], pj[1] - pi[1], pj[2] - pi[2]]);
pub fn distance<T: Into<Vector3f>>(&self, pi: T, pj: T) -> f64 {
let p = pj.into() - pi.into();
let pmic = self.apply_mic(p);
pmic.norm()
}

/// Return the shortest vector by applying the minimum image convention.
pub fn apply_mic(&self, p: [f64; 3]) -> Vector3f {
/// Return the shortest vector obeying the minimum image convention.
pub fn apply_mic<T: Into<[f64; 3]>>(&self, p: T) -> Vector3f {
let p = p.into();
// Tuckerman algorithm works well for Orthorombic cell
let v_naive = self.apply_mic_tuckerman(p);
if self.is_orthorhombic() {
Expand Down
4 changes: 1 addition & 3 deletions src/mic.rs
Original file line number Diff line number Diff line change
Expand Up @@ -31,9 +31,7 @@ impl Lattice {
}

// transform back to cartesian coordinates
let pij = self.to_cart(fcoords);

Vector3f::from(pij)
self.to_cart(fcoords)
}

/// Return the distance between two points computed using the minimum image
Expand Down
18 changes: 9 additions & 9 deletions src/supercell.rs
Original file line number Diff line number Diff line change
Expand Up @@ -3,25 +3,25 @@
// [[file:~/Workspace/Programming/gchemol-rs/lattice/lattice.note::*base][base:1]]
use crate::Lattice;
use guts::itertools::*;

type Point = [f64; 3];
use vecfx::Vector3f;

impl Lattice {
/// Create a supercell along three cell directions.
pub fn replicate(
pub fn replicate<T: Into<Vector3f> + Copy>(
&self,
points: &[Point],
points: &[T],
ra: impl Iterator<Item = isize> + Clone,
rb: impl Iterator<Item = isize> + Clone,
rc: impl Iterator<Item = isize> + Clone,
) -> Vec<Point> {
) -> Vec<Vector3f> {
iproduct!(ra, rb, rc)
.flat_map(|(i, j, k)| {
let v = [i as f64, j as f64, k as f64];
let [tx, ty, tz] = self.to_cart(v);
points
.iter()
.map(move |[x0, y0, z0]| [x0 + tx, y0 + ty, z0 + tz])
// let [tx, ty, tz]: [f64; 3] = self.to_cart(v).into();
let tv = self.to_cart(v);
points.iter().map(move |&p| {
tv + p.into()
})
})
.collect()
}
Expand Down
20 changes: 10 additions & 10 deletions tests/basic.rs
Original file line number Diff line number Diff line change
Expand Up @@ -14,16 +14,13 @@ fn test_lattice_construct() {
let mut lat = Lattice::default();
let loc = [1.0, 2.0, 3.0];
lat.set_origin(loc);
assert_eq!(loc, lat.origin());

let lat = Lattice::new([[18.256, 0., 0.], [0., 20.534, 0.], [0., 0., 15.084]]);
assert_eq!(true, lat.is_orthorhombic());

let lat = Lattice::new([[15.3643, 0., 0.], [4.5807, 15.5026, 0.], [0., 0., 17.4858]]);

let [a, b, c] = lat.lengths();
assert_eq!(false, lat.is_orthorhombic());

assert_relative_eq!(a, 15.3643, epsilon = 1e-4);
assert_relative_eq!(b, 16.1652, epsilon = 1e-4);
assert_relative_eq!(c, 17.4858, epsilon = 1e-4);
Expand All @@ -41,12 +38,15 @@ fn test_lattice_construct() {
#[test]
fn test_lattice_volume() {
let vts = [[5., 0., 0.], [5., 5., 0.], [1., 0., 5.]];

let mut lat = Lattice::new(vts);
assert_eq!(vts, lat.vectors());
assert_eq!(vts[0], lat.vector_a());
assert_eq!(vts[1], lat.vector_b());
assert_eq!(vts[2], lat.vector_c());
let mat: Matrix3f = vts.into();
assert_eq!(mat, lat.matrix());
let va: [f64; 3] = lat.vector_a().into();
let vb: [f64; 3] = lat.vector_b().into();
let vc: [f64; 3] = lat.vector_c().into();
assert_eq!(vts[0], va);
assert_eq!(vts[1], vb);
assert_eq!(vts[2], vc);

assert_relative_eq!(125.0, lat.volume(), epsilon = 1e-4);
lat.scale_by(4.);
Expand Down Expand Up @@ -83,12 +83,12 @@ fn test_wrap() {

// Orthorhombic unit cell
let cell = Lattice::from_params(3.0, 4.0, 5.0, 90.0, 90.0, 90.0);
let wrapped: Vector3f = cell.wrap([1.0, 1.5, 6.0]).into();
let wrapped = cell.wrap([1.0, 1.5, 6.0]);
assert_relative_eq!(wrapped, Vector3f::from([1.0, 1.5, 1.0]), epsilon = 1e-4);

// Triclinic unit cell
let cell = Lattice::from_params(3.0, 4.0, 5.0, 90.0, 90.0, 90.0);
let wrapped: Vector3f = cell.wrap([1.0, 1.5, 6.0]).into();
let wrapped = cell.wrap([1.0, 1.5, 6.0]);
assert_relative_eq!(wrapped, Vector3f::from([1.0, 1.5, 1.0]), epsilon = 1e-4);
}
// basic.rs:1 ends here

0 comments on commit b41f805

Please sign in to comment.