trunk/src/unitcell.cc

   1 // Voro++, a 3D cell-based Voronoi library
   2 //
   3 // Author   : Chris H. Rycroft (LBL / UC Berkeley)
   4 // Email    : chr@alum.mit.edu
   5 // Date     : August 30th 2011
   6
   7 /** \file unitcell.cc
   8  * \brief Function implementations for the unitcell class. */
   9
  10 #include <cmath>
  11 #include <queue>
  12
  13 #include "unitcell.hh"
  14 #include "cell.hh"
  15
  16 namespace voro {
  17
  18 /** Initializes the unit cell class for a particular non-orthogonal periodic
  19  * geometry, corresponding to a parallelepiped with sides given by three
  20  * vectors. The class constructs the unit Voronoi cell corresponding to this
  21  * geometry.
  22  * \param[in] (bx_) The x coordinate of the first unit vector.
  23  * \param[in] (bxy_,by_) The x and y coordinates of the second unit vector.
  24  * \param[in] (bxz_,byz_,bz_) The x, y, and z coordinates of the third unit
  25  *                            vector. */
  26 unitcell::unitcell(double bx_,double bxy_,double by_,double bxz_,double byz_,double bz_)
  27         : bx(bx_), bxy(bxy_), by(by_), bxz(bxz_), byz(byz_), bz(bz_),
  28         unit_voro(max_unit_voro_shells*max_unit_voro_shells*4*(bx*bx+by*by+bz*bz)) {
  29         int i,j,l=1;
  30
  31         // Initialize the Voronoi cell to be a very large rectangular box
  32         const double ucx=max_unit_voro_shells*bx,ucy=max_unit_voro_shells*by,ucz=max_unit_voro_shells*bz;
  33         unit_voro.init(-ucx,ucx,-ucy,ucy,-ucz,ucz);
  34
  35         // Repeatedly cut the cell by shells of periodic image particles
  36         while(l<2*max_unit_voro_shells) {
  37
  38                 // Check to see if any of the planes from the current shell
  39                 // will cut the cell
  40                 if(unit_voro_intersect(l)) {
  41
  42                         // If they do, apply the plane cuts from the current
  43                         // shell
  44                         unit_voro_apply(l,0,0);
  45                         for(i=1;i<l;i++) {
  46                                 unit_voro_apply(l,i,0);
  47                                 unit_voro_apply(-l,i,0);
  48                         }
  49                         for(i=-l;i<=l;i++) unit_voro_apply(i,l,0);
  50                         for(i=1;i<l;i++) for(j=-l+1;j<=l;j++) {
  51                                 unit_voro_apply(l,j,i);
  52                                 unit_voro_apply(-j,l,i);
  53                                 unit_voro_apply(-l,-j,i);
  54                                 unit_voro_apply(j,-l,i);
  55                         }
  56                         for(i=-l;i<=l;i++) for(j=-l;j<=l;j++) unit_voro_apply(i,j,l);
  57                 } else {
  58
  59                         // Calculate a bound on the maximum y and z coordinates
  60                         // that could possibly cut the cell. This is based upon
  61                         // a geometric result that particles with z>l can't cut
  62                         // a cell lying within the paraboloid
  63                         // z<=(l*l-x*x-y*y)/(2*l). It is always a tighter bound
  64                         // than the one based on computing the maximum radius
  65                         // of a Voronoi cell vertex.
  66                         max_uv_y=max_uv_z=0;
  67                         double y,z,q,*pts=unit_voro.pts,*pp=pts;
  68                         while(pp<pts+4*unit_voro.p) {
  69                                 q=*(pp++);y=*(pp++);z=*pp;pp+=2;q=sqrt(q*q+y*y+z*z);
  70                                 if(y+q>max_uv_y) max_uv_y=y+q;
  71                                 if(z+q>max_uv_z) max_uv_z=z+q;
  72                         }
  73                         max_uv_z*=0.5;
  74                         max_uv_y*=0.5;
  75                         return;
  76                 }
  77                 l++;
  78         }
  79
  80         // If the routine makes it here, then the unit cell still hasn't been
  81         // completely bounded by the plane cuts. Give the memory error code,
  82         // because this is mainly a case of hitting a safe limit, than any
  83         // inherent problem.
  84         voro_fatal_error("Periodic cell computation failed",VOROPP_MEMORY_ERROR);
  85 }
  86
  87 /** Applies a pair of opposing plane cuts from a periodic image point
  88  * to the unit Voronoi cell.
  89  * \param[in] (i,j,k) the index of the periodic image to consider. */
  90 inline void unitcell::unit_voro_apply(int i,int j,int k) {
  91         double x=i*bx+j*bxy+k*bxz,y=j*by+k*byz,z=k*bz;
  92         unit_voro.plane(x,y,z);
  93         unit_voro.plane(-x,-y,-z);
  94 }
  95
  96 /** Calculates whether the unit Voronoi cell intersects a given periodic image
  97  * of the domain.
  98  * \param[in] (dx,dy,dz) the displacement of the periodic image.
  99  * \param[out] vol the proportion of the unit cell volume within this image,
 100  *                 only computed in the case that the two intersect.
 101  * \return True if they intersect, false otherwise. */
 102 bool unitcell::intersects_image(double dx,double dy,double dz,double &vol) {
 103         const double bxinv=1/bx,byinv=1/by,bzinv=1/bz,ivol=bxinv*byinv*bzinv;
 104         voronoicell c;
 105         c=unit_voro;
 106         dx*=2;dy*=2;dz*=2;
 107         if(!c.plane(0,0,bzinv,dz+1)) return false;
 108         if(!c.plane(0,0,-bzinv,-dz+1)) return false;
 109         if(!c.plane(0,byinv,-byz*byinv*bzinv,dy+1)) return false;
 110         if(!c.plane(0,-byinv,byz*byinv*bzinv,-dy+1)) return false;
 111         if(!c.plane(bxinv,-bxy*bxinv*byinv,(bxy*byz-by*bxz)*ivol,dx+1)) return false;
 112         if(!c.plane(-bxinv,bxy*bxinv*byinv,(-bxy*byz+by*bxz)*ivol,-dx+1)) return false;
 113         vol=c.volume()*ivol;
 114         return true;
 115 }
 116
 117 /** Computes a list of periodic domain images that intersect the unit Voronoi cell.
 118  * \param[out] vi a vector containing triplets (i,j,k) corresponding to domain
 119  *                images that intersect the unit Voronoi cell, when it is
 120  *                centered in the middle of the primary domain.
 121  * \param[out] vd a vector containing the fraction of the Voronoi cell volume
 122  *                within each corresponding image listed in vi. */
 123 void unitcell::images(std::vector<int> &vi,std::vector<double> &vd) {
 124         const int ms2=max_unit_voro_shells*2+1,mss=ms2*ms2*ms2;
 125         bool *a=new bool[mss],*ac=a+max_unit_voro_shells*(1+ms2*(1+ms2)),*ap=a;
 126         int i,j,k;
 127         double vol;
 128
 129         // Initialize mask
 130         while(ap<ac) *(ap++)=true;
 131         *(ap++)=false;
 132         while(ap<a+mss) *(ap++)=true;
 133
 134         // Set up the queue and add (0,0,0) image to it
 135         std::queue<int> q;
 136         q.push(0);q.push(0);q.push(0);
 137
 138         while(!q.empty()) {
 139
 140                 // Read the next entry on the queue
 141                 i=q.front();q.pop();
 142                 j=q.front();q.pop();
 143                 k=q.front();q.pop();
 144
 145                 // Check intersection of this image
 146                 if(intersects_image(i,j,k,vol)) {
 147
 148                         // Add this entry to the output vectors
 149                         vi.push_back(i);
 150                         vi.push_back(j);
 151                         vi.push_back(k);
 152                         vd.push_back(vol);
 153
 154                         // Add neighbors to the queue if they have not been
 155                         // tested
 156                         ap=ac+i+ms2*(j+ms2*k);
 157                         if(k>-max_unit_voro_shells&&*(ap-ms2*ms2)) {q.push(i);q.push(j);q.push(k-1);*(ap-ms2*ms2)=false;}
 158                         if(j>-max_unit_voro_shells&&*(ap-ms2)) {q.push(i);q.push(j-1);q.push(k);*(ap-ms2)=false;}
 159                         if(i>-max_unit_voro_shells&&*(ap-1)) {q.push(i-1);q.push(j);q.push(k);*(ap-1)=false;}
 160                         if(i<max_unit_voro_shells&&*(ap+1)) {q.push(i+1);q.push(j);q.push(k);*(ap+1)=false;}
 161                         if(j<max_unit_voro_shells&&*(ap+ms2)) {q.push(i);q.push(j+1);q.push(k);*(ap+ms2)=false;}
 162                         if(k<max_unit_voro_shells&&*(ap+ms2*ms2)) {q.push(i);q.push(j);q.push(k+1);*(ap+ms2*ms2)=false;}
 163                 }
 164         }
 165
 166         // Remove mask memory
 167         delete [] a;
 168 }
 169
 170 /** Tests to see if a shell of periodic images could possibly cut the periodic
 171  * unit cell.
 172  * \param[in] l the index of the shell to consider.
 173  * \return True if a point in the shell cuts the cell, false otherwise. */
 174 bool unitcell::unit_voro_intersect(int l) {
 175         int i,j;
 176         if(unit_voro_test(l,0,0)) return true;
 177         for(i=1;i<l;i++) {
 178                 if(unit_voro_test(l,i,0)) return true;
 179                 if(unit_voro_test(-l,i,0)) return true;
 180         }
 181         for(i=-l;i<=l;i++) if(unit_voro_test(i,l,0)) return true;
 182         for(i=1;i<l;i++) for(j=-l+1;j<=l;j++) {
 183                 if(unit_voro_test(l,j,i)) return true;
 184                 if(unit_voro_test(-j,l,i)) return true;
 185                 if(unit_voro_test(-l,-j,i)) return true;
 186                 if(unit_voro_test(j,-l,i)) return true;
 187         }
 188         for(i=-l;i<=l;i++) for(j=-l;j<=l;j++) if(unit_voro_test(i,j,l)) return true;
 189         return false;
 190 }
 191
 192 /** Tests to see if a plane cut from a particular periodic image will cut th
 193  * unit Voronoi cell.
 194  * \param[in] (i,j,k) the index of the periodic image to consider.
 195  * \return True if the image cuts the cell, false otherwise. */
 196 inline bool unitcell::unit_voro_test(int i,int j,int k) {
 197         double x=i*bx+j*bxy+k*bxz,y=j*by+k*byz,z=k*bz;
 198         double rsq=x*x+y*y+z*z;
 199         return unit_voro.plane_intersects(x,y,z,rsq);
 200 }
 201
 202 /** Draws the periodic domain in gnuplot format.
 203  * \param[in] fp the file handle to write to. */
 204 void unitcell::draw_domain_gnuplot(FILE *fp) {
 205         fprintf(fp,"0 0 0\n%g 0 0\n%g %g 0\n%g %g 0\n",bx,bx+bxy,by,bxy,by);
 206         fprintf(fp,"%g %g %g\n%g %g %g\n%g %g %g\n%g %g %g\n",bxy+bxz,by+byz,bz,bx+bxy+bxz,by+byz,bz,bx+bxz,byz,bz,bxz,byz,bz);
 207         fprintf(fp,"0 0 0\n%g %g 0\n\n%g %g %g\n%g %g %g\n\n",bxy,by,bxz,byz,bz,bxy+bxz,by+byz,bz);
 208         fprintf(fp,"%g 0 0\n%g %g %g\n\n%g %g 0\n%g %g %g\n\n",bx,bx+bxz,byz,bz,bx+bxy,by,bx+bxy+bxz,by+byz,bz);
 209 }
 210
 211 /** Draws the periodic domain in POV-Ray format.
 212  * \param[in] fp the file handle to write to. */
 213 void unitcell::draw_domain_pov(FILE *fp) {
 214         fprintf(fp,"cylinder{0,0,0>,<%g,0,0>,rr}\n"
 215                    "cylinder{<%g,%g,0>,<%g,%g,0>,rr}\n",bx,bxy,by,bx+bxy,by);
 216         fprintf(fp,"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n"
 217                    "cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n",bxz,byz,bz,bx+bxz,byz,bz,bxy+bxz,by+byz,bz,bx+bxy+bxz,by+byz,bz);
 218         fprintf(fp,"cylinder{<0,0,0>,<%g,%g,0>,rr}\n"
 219                    "cylinder{<%g,0,0>,<%g,%g,0>,rr}\n",bxy,by,bx,bx+bxy,by);
 220         fprintf(fp,"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n"
 221                    "cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n",bxz,byz,bz,bxy+bxz,by+byz,bz,bx+bxz,byz,bz,bx+bxy+bxz,by+byz,bz);
 222         fprintf(fp,"cylinder{<0,0,0>,<%g,%g,%g>,rr}\n"
 223                    "cylinder{<%g,0,0>,<%g,%g,%g>,rr}\n",bxz,byz,bz,bx,bx+bxz,byz,bz);
 224         fprintf(fp,"cylinder{<%g,%g,0>,<%g,%g,%g>,rr}\n"
 225                    "cylinder{<%g,%g,0>,<%g,%g,%g>,rr}\n",bxy,by,bxy+bxz,by+byz,bz,bx+bxy,by,bx+bxy+bxz,by+byz,bz);
 226         fprintf(fp,"sphere{<0,0,0>,rr}\nsphere{<%g,0,0>,rr}\n"
 227                    "sphere{<%g,%g,0>,rr}\nsphere{<%g,%g,0>,rr}\n",bx,bxy,by,bx+bxy,by);
 228         fprintf(fp,"sphere{<%g,%g,%g>,rr}\nsphere{<%g,%g,%g>,rr}\n"
 229                    "sphere{<%g,%g,%g>,rr}\nsphere{<%g,%g,%g>,rr}\n",bxz,byz,bz,bx+bxz,byz,bz,bxy+bxz,by+byz,bz,bx+bxy+bxz,by+byz,bz);
 230 }
 231
 232 }