From 211af9d10edeeb46220db8fe8e4b605ac2948696 Mon Sep 17 00:00:00 2001
From: Niklaus Johner <nij2003@med.cornell.edu>
Date: Tue, 22 Oct 2013 12:48:17 -0400
Subject: [PATCH] Added Function to Wrap a Vec3 and Vec3List in a
 non-orthogonal unit cell

---
 modules/geom/pymod/export_vecmat3_op.cc | 11 ++++++++---
 modules/geom/src/vecmat3_op.cc          | 20 ++++++++++++++++++++
 modules/geom/src/vecmat3_op.hh          |  7 +++++--
 3 files changed, 33 insertions(+), 5 deletions(-)

diff --git a/modules/geom/pymod/export_vecmat3_op.cc b/modules/geom/pymod/export_vecmat3_op.cc
index 5c06a2499..b2f7c17bf 100644
--- a/modules/geom/pymod/export_vecmat3_op.cc
+++ b/modules/geom/pymod/export_vecmat3_op.cc
@@ -40,7 +40,10 @@ Vec3 (*Vec3Min)(const Vec3&, const Vec3&) = &Min;
 Vec3 (*Vec3Max)(const Vec3&, const Vec3&) = &Max;
 Real (*Vec3Distance2WithPBC)(const Vec3&, const Vec3&, const Vec3&)   = &Distance2WithPBC;
 Real (*Vec3DistanceWithPBC)(const Vec3&, const Vec3&, const Vec3&)    = &DistanceWithPBC;
-
+Vec3 (*wrap_vec3_1)(const Vec3&, const Vec3&, const Vec3&)               = &WrapVec3;
+Vec3 (*wrap_vec3_2)(const Vec3&, const Vec3&, const Vec3&, const Vec3&)  = &WrapVec3;
+Vec3List (*wrap_vec3list1)(const Vec3List&, const Vec3&, const Vec3&)               = &WrapVec3List;
+Vec3List (*wrap_vec3list2)(const Vec3List&, const Vec3&, const Vec3&, const Vec3&)  = &WrapVec3List;
 void export_VecMat3_op()
 {
   using namespace geom;
@@ -70,6 +73,8 @@ void export_VecMat3_op()
   def("MinDistance",MinDistance);
   def("MinDistanceIndices",MinDistanceIndices);
   def("MinDistanceWithPBC",MinDistanceWithPBC);
-  def("WrapVec3",WrapVec3);
-  def("WrapVec3List",WrapVec3List);
+  def("WrapVec3",wrap_vec3_1,(arg("vector"),arg("cell_center"),arg("cell_size")));
+  def("WrapVec3",wrap_vec3_2,(arg("vector"),arg("cell_center"),arg("cell_size"),arg("cell_angles")));
+  def("WrapVec3List",wrap_vec3list1,(arg("vector_list"),arg("cell_center"),arg("cell_size")));
+  def("WrapVec3List",wrap_vec3list2,(arg("vector_list"),arg("cell_center"),arg("cell_size"),arg("cell_angles")));  
 }
diff --git a/modules/geom/src/vecmat3_op.cc b/modules/geom/src/vecmat3_op.cc
index 994b74cd0..3c5cbdea2 100644
--- a/modules/geom/src/vecmat3_op.cc
+++ b/modules/geom/src/vecmat3_op.cc
@@ -288,5 +288,25 @@ Vec3List WrapVec3List(const Vec3List& vl, const Vec3& box_center,const Vec3& uce
   return vl_out;
 }
   
+Vec3 WrapVec3(const Vec3& v1,const Vec3& box_center,const Vec3& ucell_size,const Vec3& ucell_angles){
+  Vec3List vl=CalculateUnitCellVectors(ucell_size,ucell_angles);
+  Vec3 v=v1-box_center,v_wrapped=v1,r;
+  for (int i=0; i<3; i++) {
+    r[2-i]=v[2-i]/vl[2-i][2-i];
+    v=v-r[2-i]*vl[2-i];
+    r[2-i]=(r[2-i] > 0.0) ? std::floor(r[2-i] + 0.5) : std::ceil(r[2-i] - 0.5);
+    v_wrapped=v_wrapped-vl[2-i]*r[2-i];
+  }
+  return v_wrapped;
+}
+
+Vec3List WrapVec3List(const Vec3List& vl, const Vec3& box_center,const Vec3& ucell_size,const Vec3& ucell_angles){
+  Vec3List vl_out;
+  vl_out.reserve(vl_out.size());
+  for (Vec3List::const_iterator v1=vl.begin(),e=vl.end();v1!=e ; v1++) {
+    vl_out.push_back(WrapVec3(*v1,box_center,ucell_size,ucell_angles));
+  }
+  return vl_out;
+}  
   
 } // ns
diff --git a/modules/geom/src/vecmat3_op.hh b/modules/geom/src/vecmat3_op.hh
index 084c69907..fbaaec161 100644
--- a/modules/geom/src/vecmat3_op.hh
+++ b/modules/geom/src/vecmat3_op.hh
@@ -222,9 +222,12 @@ std::vector<unsigned int> DLLEXPORT_OST_GEOM MinDistanceIndices(const Vec3List&
 Vec3List DLLEXPORT_OST_GEOM CalculateUnitCellVectors(const Vec3& ucell_size, const Vec3& ucell_angles);
 //!wraps a vector in a box with periodic boundaries
 Vec3 DLLEXPORT_OST_GEOM WrapVec3(const Vec3& v1,const Vec3& box_center,const Vec3& ucell_size);
-//!wraps all the verctors in a Vec3List in a box with periodic boundaries
+//!wraps all the vectors in a Vec3List in a box with periodic boundaries
 Vec3List DLLEXPORT_OST_GEOM WrapVec3List(const Vec3List& vl,const Vec3& box_center,const Vec3& ucell_size);
-
+//!wraps a vector in a non-rothogonal box with periodic boundaries
+Vec3 DLLEXPORT_OST_GEOM WrapVec3(const Vec3& v1,const Vec3& box_center,const Vec3& ucell_size,const Vec3& ucell_angles);
+//!wraps all the vectors in a Vec3List in a non-rothogonal box with periodic boundaries
+Vec3List DLLEXPORT_OST_GEOM WrapVec3List(const Vec3List& vl,const Vec3& box_center,const Vec3& ucell_size,const Vec3& ucell_angles);
   
 } // ns
 
-- 
GitLab