From d399cbc6ce83ad7576251a130bba3ce8ddb7c17e Mon Sep 17 00:00:00 2001
From: Giovanni Bussi <giovanni.bussi@gmail.com>
Date: Sun, 28 Apr 2013 11:41:49 +0200
Subject: [PATCH] Added some developer documentation for Pbc.
---
src/tools/Pbc.cpp | 25 +++++++++++++++++++----
src/tools/Pbc.h | 51 +++++++++++++++++++++++++++++++++++++++++++++--
2 files changed, 70 insertions(+), 6 deletions(-)
diff --git a/src/tools/Pbc.cpp b/src/tools/Pbc.cpp
index e2fc52bf6..15e1ba401 100644
--- a/src/tools/Pbc.cpp
+++ b/src/tools/Pbc.cpp
@@ -38,23 +38,40 @@ Pbc::Pbc():
void Pbc::buildShifts(std::vector<Vector> shifts[2][2][2])const{
const double small=1e-28;
+
+// clear all shifts
for(int i=0;i<2;i++) for(int j=0;j<2;j++) for(int k=0;k<2;k++) shifts[i][j][k].clear();
+
+// enumerate all possible shifts
+// since box is reduced, only 27 shifts have to be attempted
for(int l=-1;l<=1;l++) for(int m=-1;m<=1;m++) for(int n=-1;n<=1;n++){
+
+// int/double shift vectors
int ishift[3]={l,m,n};
Vector dshift(l,m,n);
- unsigned count=0;
+
// count how many components are != 0
+ unsigned count=0;
for(int s=0;s<3;s++) if(ishift[s]!=0) count++;
+
// skips trivial (0,0,0) and cases with three shifts
+// only 18 shifts survive past this point
if(count==0 || count==3) continue;
-// check if that WS face is perpendicular to the axis
+
+// check if that Wigner-Seitz face is perpendicular to the axis.
+// this allows to eliminate shifts in symmetric cells.
+// e.g., if one lactice vector is orthogonal to the plane spanned
+// by the other two vectors, that shift should never be tried
Vector cosdir=matmul(reduced,transpose(reduced),dshift);
double dp=dotProduct(dshift,cosdir);
double ref=modulo2(dshift)*modulo2(cosdir);
if(std::fabs(ref-dp*dp)<small) continue;
+// here we start pruning depending on the sign of the scaled coordinate
for(int i=0;i<2;i++) for(int j=0;j<2;j++) for(int k=0;k<2;k++){
+
int block[3]={2*i-1,2*j-1,2*k-1};
+
// skip cases where shift would bring too far from origin
bool skip=false;
for(int s=0;s<3;s++) if(ishift[s]*block[s]>0) skip=true;
@@ -66,6 +83,8 @@ void Pbc::buildShifts(std::vector<Vector> shifts[2][2][2])const{
if(((1-ishift[s]*ishift[s])*block[s])*cosdir[s]<-small) skip=false;
}
if(skip)continue;
+
+// if we arrive to this point, shift is eligible and is added to the list
shifts[i][j][k].push_back(matmul(transpose(reduced),dshift));
}
}
@@ -83,8 +102,6 @@ void Pbc::fullSearch(Vector&d)const{
Vector best(d);
double lbest=d.modulo2();
for(int i=-smax;i<=smax;i++) for(int j=-smax;j<=smax;j++) for(int k=-smax;k<=smax;k++){
-// int x=i*i+j*j+k*k;
-// if(x==0 || x==3) continue;
Vector trial=d+i*a0+j*a1+k*a2;
double ltrial=trial.modulo2();
if(ltrial<lbest){
diff --git a/src/tools/Pbc.h b/src/tools/Pbc.h
index 1aa9dde9e..804feae5b 100644
--- a/src/tools/Pbc.h
+++ b/src/tools/Pbc.h
@@ -29,29 +29,76 @@
namespace PLMD{
+/*
+Tool to deal with periodic boundary conditions.
+
+This class is useful to apply periodic boundary conditions on interatomic
+distances. It stores privately information about reduced lattice vectors
+*/
class Pbc{
+/// Type of box
enum {unset,orthorombic,generic} type;
+/// Box
Tensor box;
+/// Inverse box
Tensor invBox;
+/// Reduced box.
+/// This is a set of lattice vectors generating the same lattice
+/// but "minimally skewed". Useful to optimize image search.
Tensor reduced;
+/// Inverse of the reduced box
Tensor invReduced;
+/// List of shifts that should be attempted.
+/// Depending on the sign of the scaled coordinates representing
+/// a distance vector, a different set of shifts must be tried.
std::vector<Vector> shifts[2][2][2];
+/// Alternative representation for orthorombic cells.
+/// Not really used, but could be used to optimize search in
+/// orthorombic cells.
Vector diag,hdiag,mdiag;
+/// Build list of shifts.
+/// This is expensive, and must be called only when box is
+/// reset. It allows building a minimal set of shifts
+/// depending on the sign of the scaled coordinates representing
+/// a distance vector.
void buildShifts(std::vector<Vector> shifts[2][2][2])const;
+/// Full search (for testing)
void fullSearch(Vector&)const;
+/// internal version of distance, also returns the number
+/// of attempted shifts (used in Pbc::test()).
+ Vector distance(const Vector&,const Vector&,int*nshifts)const;
public:
+/// Perform some check. Useful for debugging.
static void test();
+/// Constructor
Pbc();
+/// Compute modulo of (v2-v1), using or not pbc depending on bool pbc.
double distance( const bool pbc, const Vector& v1, const Vector& v2 ) const;
- Vector distance(const Vector&,const Vector&,int*nshifts=NULL)const;
- void setBox(const Tensor&);
+/// Computes v2-v1, using minimal image convention
+ Vector distance(const Vector& v1,const Vector& v2)const;
+/// Set the lattice vectors.
+/// b[i][j] is the j-th component of the i-th vector
+ void setBox(const Tensor&b);
+/// Returns the box
const Tensor& getBox()const;
+/// Returns the inverse matrix of box.
+/// Thus: pbc.getInvBox() == inverse(pbc.getBox()).
const Tensor& getInvBox()const;
+/// Transform a vector in real space to a vector in scaled coordinates.
+/// Thus:pbc.realToScaled(v) == matmul(transpose(inverse(pbc.getBox(),v)));
Vector realToScaled(const Vector&)const;
+/// Transform a vector in scaled coordinates to a vector in real space.
+/// Thus:pbc.scaledToRead(v) == matmul(transpose(pbc.getBox()),v);
Vector scaledToReal(const Vector&)const;
+/// Returns true if the box vectors are orthogonal
bool isOrthorombic()const;
};
+inline
+Vector Pbc::distance(const Vector& v1,const Vector& v2)const{
+ return distance(v1,v2,NULL);
+}
+
}
#endif
--
GitLab