Added functions to calculate the minimal distance between the points

In two Vec3List, with and without periodic boundary conditions

modules/geom/pymod/export_vecmat3_op.cc

@@ -38,7 +38,8 @@ Real (*Mat3Comp)(const Mat3& m, unsigned int i, unsigned int j)       = &Comp;
 Real (*Mat3Minor)(const Mat3& m, unsigned int i, unsigned int j)      = &Minor;
 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;
 
 void export_VecMat3_op()
 {
@@ -64,5 +65,8 @@ void export_VecMat3_op()
   def("OrthogonalVector",OrthogonalVector);
   def("Min",Vec3Min);
   def("Max",Vec3Max);
+  def("Distance2WithPBC",Vec3Distance2WithPBC);
+  def("DistanceWithPBC",Vec3DistanceWithPBC);
   def("MinDistance",MinDistance);
+  def("MinDistanceWithPBC",MinDistanceWithPBC);
 }

modules/geom/src/vecmat3_op.cc

@@ -192,6 +192,7 @@ Real DihedralAngle(const Vec3& p1, const Vec3& p2, const Vec3& p3,
                Dot(r12cross, r23cross));
 }
 
+  
 Real MinDistance(const Vec3List& l1, const Vec3List& l2)
 { 
   // returns the minimal distance between two sets of points (Vec3List)
@@ -207,4 +208,24 @@ Real MinDistance(const Vec3List& l1, const Vec3List& l2)
   return std::sqrt(min);
 }
 
+Real MinDistanceWithPBC(const Vec3List& l1, const Vec3List& l2, Vec3& basis_vec)
+{ 
+  // returns the minimal distance between two sets of points (Vec3List)
+  // given the periodic boundary condition along x,y,z given in the basis_vec
+  if (l1.size()==0 || l2.size()==0){throw std::runtime_error("cannot calculate minimal distance: empty Vec3List");}
+  Real min=Length2(*l1.begin()-*l2.begin());
+  Real d;
+  Vec3 v;
+  for (int i=0; i<3; i++) {
+    basis_vec[i]=std::fabs(basis_vec[i]);
+  }
+  for (Vec3List::const_iterator p1=l1.begin(),e1=l1.end(); p1!=e1; p1++) {
+    for (Vec3List::const_iterator p2=l2.begin(),e2=l2.end(); p2!=e2; p2++) {
+      d=Distance2WithPBC(*p1, *p2, basis_vec);
+      if (d<min) min=d;
+    }
+  }
+  return std::sqrt(min);
+}  
+  
 } // ns

modules/geom/src/vecmat3_op.hh

@@ -194,9 +194,27 @@ inline Real Distance(const Vec3& p1, const Vec3& p2)
     return Length(p1-p2);
 }
 
+
+//! return the squared distance between two points with periodic boundaries in x,y,z given in basis_vec
+inline Real Distance2WithPBC(const Vec3& v1, const Vec3& v2, const Vec3& basis_vec){
+  Vec3 v;
+  v=v1-v2;
+  for (int i=0; i<3; i++) {
+    if (std::fabs(v[i])>basis_vec[i]/2.){ 
+      v[i]=std::fabs(v[i])-basis_vec[i]*int(std::fabs(v[i])/basis_vec[i]+0.5);
+    }
+  }
+  return Length2(v);
+}
+//! return the distance between two points with periodic boundaries in x,y,z given in basis_vec
+inline Real DistanceWithPBC(const Vec3& v1, const Vec3& v2, const Vec3& basis_vec){
+  return sqrt(Distance2WithPBC(v1, v2, basis_vec));
+}
 //! returns the minimal distance between the points in two Vec3List
 Real MinDistance(const Vec3List& l1, const Vec3List& l2);
-
+//! returns the minimal distance between the points in two Vec3List 
+//  with periodic boundaries in x,y,z given in basis_vec
+Real MinDistanceWithPBC(const Vec3List& l1, const Vec3List& l2, Vec3& basis_vec);
 } // ns
 
 #endif