From db509d5e37393c1161d92407286c4b74ad900ee7 Mon Sep 17 00:00:00 2001
From: Giovanni Bussi <giovanni.bussi@gmail.com>
Date: Wed, 28 Nov 2012 09:55:38 +0100
Subject: [PATCH] Optimized PBCs
I took advantage of the fact that we have the reduced cell. This
allows to avoid the comparison of different shifts for
many values of the scaled coordinates, as well as using the
sign of the scaled coordinates to check which shifts should
be attempted. The average number of shift vector to be evaluated
is approximately 1 (plus the "null-shift")
There is also a Pbc::test() function which I used during debugging.
I leave it there for future reference.
---
src/tools/Pbc.cpp | 190 +++++++++++++++++++++++++++++++++++++++++++---
src/tools/Pbc.h | 9 ++-
2 files changed, 185 insertions(+), 14 deletions(-)
diff --git a/src/tools/Pbc.cpp b/src/tools/Pbc.cpp
index 9a9f8e9a5..ce60029af 100644
--- a/src/tools/Pbc.cpp
+++ b/src/tools/Pbc.cpp
@@ -23,6 +23,9 @@
#include "Tools.h"
#include "Exception.h"
#include "LatticeReduction.h"
+#include <iostream>
+#include "Random.h"
+#include <cmath>
namespace PLMD{
@@ -33,19 +36,55 @@ Pbc::Pbc():
invBox.zero();
}
+void Pbc::buildShifts(std::vector<Vector> shifts[2][2][2])const{
+ const double small=1e-28;
+ 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();
+ for(int l=-1;l<=1;l++) for(int m=-1;m<=1;m++) for(int n=-1;n<=1;n++){
+ int ishift[3]={l,m,n};
+ Vector dshift(l,m,n);
+ unsigned count=0;
+// count how many components are != 0
+ for(int s=0;s<3;s++) if(ishift[s]!=0) count++;
+// skips trivial (0,0,0) and cases with three shifts
+ if(count==0 || count==3) continue;
+// check if that WS face is perpendicular to the axis
+ 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;
+
+ 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;
+ if(skip) continue;
+ skip=true;
+ for(int s=0;s<3;s++){
+// check that the components of cosdir along the non-shifted directions
+// have the proper sign
+ if(((1-ishift[s]*ishift[s])*block[s])*cosdir[s]<-small) skip=false;
+ }
+ if(skip)continue;
+ shifts[i][j][k].push_back(matmul(transpose(reduced),dshift));
+ }
+ }
+}
+
void Pbc::fullSearch(Vector&d)const{
+ if(type==unset) return;
Vector s=matmul(invReduced.transpose(),d);
for(int i=0;i<3;i++) s[i]=Tools::pbc(s[i]);
d=matmul(reduced.transpose(),s);
- const int smax=1;
+ const int smax=4;
Vector a0(reduced.getRow(0));
Vector a1(reduced.getRow(1));
Vector a2(reduced.getRow(2));
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;
+// 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){
@@ -64,7 +103,7 @@ void Pbc::setBox(const Tensor&b){
type=unset;
double det=box.determinant();
- if(det*det<epsilon)return;
+ if(det*det<epsilon) return;
bool cxy=false;
bool cxz=false;
@@ -76,13 +115,11 @@ void Pbc::setBox(const Tensor&b){
invBox=box.inverse();
if(cxy && cxz && cyz) type=orthorombic;
-// NOT IMPLEMENTED YET; WILL FALL TO GENERIC
-// else if(cxy && cxz) type=yz;
-// else if(cxy && cyz) type=xz;
-// else if(cxz && cyz) type=xy;
else type=generic;
-
+
if(type==orthorombic){
+ reduced=box;
+ invReduced=inverse(reduced);
for(unsigned i=0;i<3;i++){
diag[i]=box[i][i];
hdiag[i]=0.5*box[i][i];
@@ -92,7 +129,9 @@ void Pbc::setBox(const Tensor&b){
reduced=box;
LatticeReduction::reduce(reduced);
invReduced=inverse(reduced);
- }
+ buildShifts(shifts);
+}
+
}
double Pbc::distance( const bool pbc, const Vector& v1, const Vector& v2 ) const {
@@ -100,7 +139,7 @@ double Pbc::distance( const bool pbc, const Vector& v1, const Vector& v2 ) const
else{ return ( delta(v1,v2) ).modulo(); }
}
-Vector Pbc::distance(const Vector&v1,const Vector&v2)const{
+Vector Pbc::distance(const Vector&v1,const Vector&v2,int*nshifts)const{
Vector d=delta(v1,v2);
if(type==unset){
} else if(type==orthorombic) {
@@ -111,7 +150,35 @@ Vector Pbc::distance(const Vector&v1,const Vector&v2)const{
// while(d[i]<=mdiag[i]) d[i]+=diag[i];
// }
} else if(type==generic) {
- fullSearch(d);
+ Vector s=matmul(d,invReduced);
+// check if images have to be computed:
+// if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)){
+// NOTICE: the check in the previous line, albeit correct, is breaking many regtest
+// since it does not apply Tools::pbc in many cases. Moreover, it does not
+// introduce a significant gain. I thus leave it out for the moment.
+ if(true){
+// bring to -0.5,+0.5 region in scaled coordinates:
+ for(int i=0;i<3;i++) s[i]=Tools::pbc(s[i]);
+ d=matmul(s,reduced);
+// check if shifts have to be attempted:
+ if((std::fabs(s[0])+std::fabs(s[1])+std::fabs(s[2])>0.5)){
+// list of shifts is specific for that "octant" (depends on signs of s[i]):
+ const std::vector<Vector> & myshifts(shifts[(s[0]>0?1:0)][(s[1]>0?1:0)][(s[2]>0?1:0)]);
+ Vector best(d);
+ double lbest(modulo2(best));
+// loop over possible shifts:
+ if(nshifts) *nshifts+=myshifts.size();
+ for(int i=0;i<myshifts.size();i++){
+ Vector trial=d+myshifts[i];
+ double ltrial=modulo2(trial);
+ if(ltrial<lbest){
+ lbest=ltrial;
+ best=trial;
+ }
+ }
+ d=best;
+ }
+ }
} else plumed_merror("unknown pbc type");
return d;
}
@@ -136,6 +203,105 @@ const Tensor& Pbc::getInvBox()const{
return invBox;
}
+void Pbc::test(){
+ Random r;
+ r.setSeed(-20);
+ for(int i=0;i<1000;i++){
+// random matrix with some zero element
+ Tensor box;
+ for(int j=0;j<3;j++) for(int k=0;k<3;k++) if(r.U01()>0.2){
+ box[j][k]=2.0*r.U01()-1.0;
+ }
+ int boxtype=i%10;
+ switch(boxtype){
+ case 0:
+// cubic
+ for(int j=0;j<3;j++) for(int k=0;k<3;k++) if(j!=k) box[j][k]=0.0;
+ for(int j=1;j<3;j++) box[j][j]=box[0][0];
+ break;
+ case 1:
+// orthorombic
+ for(int j=0;j<3;j++) for(int k=0;k<3;k++) if(j!=k) box[j][k]=0.0;
+ break;
+ case 2:
+// hexagonal
+ {
+ int perm=r.U01()*100;
+ Vector a;
+ a(0)=r.U01()*2-2; a(1)=0.0;a(2)=0.0;
+ double d=r.U01()*2-2;
+ Vector b(0.0,d,0.0);
+ Vector c(0.0,0.5*d,sqrt(3.0)*d*0.5);
+ box.setRow((perm+0)%3,a);
+ box.setRow((perm+1)%3,b);
+ box.setRow((perm+2)%3,c);
+ }
+ break;
+ case 3:
+// bcc
+ {
+ int perm=r.U01()*100;
+ double d=r.U01()*2-2;
+ Vector a(d,d,d);
+ Vector b(d,-d,d);
+ Vector c(d,d,-d);
+ box.setRow((perm+0)%3,a);
+ box.setRow((perm+1)%3,b);
+ box.setRow((perm+2)%3,c);
+ }
+ break;
+ case 4:
+// fcc
+ {
+ int perm=r.U01()*100;
+ double d=r.U01()*2-2;
+ Vector a(d,d,0);
+ Vector b(d,0,d);
+ Vector c(0,d,d);
+ box.setRow((perm+0)%3,a);
+ box.setRow((perm+1)%3,b);
+ box.setRow((perm+2)%3,c);
+ }
+ break;
+ default:
+// triclinic
+ break;
+ }
+
+ Pbc pbc;
+ pbc.setBox(box);
+ std::cerr<<"( "<<boxtype<<" )\n";
+ std::cerr<<"Box:";
+ for(int j=0;j<3;j++) for(int k=0;k<3;k++) std::cerr<<" "<<box[j][k];
+ std::cerr<<"\n";
+ std::cerr<<"Determinant: "<<determinant(box)<<"\n";
+ std::cerr<<"Shifts:";
+ for(int j=0;j<2;j++) for(int k=0;k<2;k++) for(int l=0;l<2;l++) std::cerr<<" "<<pbc.shifts[j][k][l].size();
+ std::cerr<<"\n";
+ int nshifts=0;
+ int ntot=10000;
+ for(int j=0;j<ntot;j++){
+ Vector v(r.U01()-0.5,r.U01()-0.5,r.U01()-0.5);
+ v*=5;
+ for(int j=0;j<3;j++) if(r.U01()>0.2) v(j)=0.0;
+ Vector full(v);
+ Vector fast=pbc.distance(Vector(0,0,0),v,&nshifts);
+ full=fast;
+
+ pbc.fullSearch(full);
+
+ if(modulo2(fast-full)>1e-10) {
+ std::cerr<<"orig "<<v[0]<<" "<<v[1]<<" "<<v[2]<<"\n";
+ std::cerr<<"fast "<<fast[0]<<" "<<fast[1]<<" "<<fast[2]<<"\n";
+ std::cerr<<"full "<<full[0]<<" "<<full[1]<<" "<<full[2]<<"\n";
+ std::cerr<<"diff "<<modulo2(fast)-modulo2(full)<<std::endl;
+ if(std::fabs(modulo2(fast)-modulo2(full))>1e-15) plumed_error();
+ }
+ }
+ std::cerr<<"Average number of shifts: "<<double(nshifts)/double(ntot)<<"\n";
+ }
+}
+
}
diff --git a/src/tools/Pbc.h b/src/tools/Pbc.h
index bb193146d..238d9ee43 100644
--- a/src/tools/Pbc.h
+++ b/src/tools/Pbc.h
@@ -24,21 +24,26 @@
#include "Vector.h"
#include "Tensor.h"
+#include <vector>
+#include <cstddef>
namespace PLMD{
class Pbc{
- enum {unset,orthorombic,xy,xz,yz,generic} type;
+ enum {unset,orthorombic,generic} type;
Tensor box;
Tensor invBox;
Tensor reduced;
Tensor invReduced;
+ std::vector<Vector> shifts[2][2][2];
Vector diag,hdiag,mdiag;
+ void buildShifts(std::vector<Vector> shifts[2][2][2])const;
void fullSearch(Vector&)const;
public:
+ static void test();
Pbc();
double distance( const bool pbc, const Vector& v1, const Vector& v2 ) const;
- Vector distance(const Vector&,const Vector&)const;
+ Vector distance(const Vector&,const Vector&,int*nshifts=NULL)const;
void setBox(const Tensor&);
const Tensor& getBox()const;
const Tensor& getInvBox()const;
--
GitLab