Lattice.h

Go to the documentation of this file.
00001 
00007 #ifndef LATTICE_H
00008 #define LATTICE_H
00009 
00010 #include <stdlib.h>
00011 #include "NamdTypes.h"
00012 #include <math.h>
00013 #include "Tensor.h"
00014 
00015 #ifdef POWERPC_TANINT
00016 extern "builtin" double __tanint(double); //IEEE round
00017 #define latticenearbyint(x)  __tanint(x)
00018 #else
00019 #define latticenearbyint(x)  floor((x)+0.5)
00020 #endif
00021 
00022 // Use latticenearbyint(X) instead of rint(X) because rint() is sensitive
00023 // to the current rounding mode and floor() is not.  It's just safer.
00024 
00025 typedef Vector ScaledPosition;
00026 
00027 class Lattice
00028 {
00029 public:
00030   Lattice(void) : a1(0,0,0), a2(0,0,0), a3(0,0,0),
00031                   b1(0,0,0), b2(0,0,0), b3(0,0,0),
00032                   o(0,0,0), p1(0), p2(0), p3(0) {};
00033 
00034   // maps a transformation triplet onto a single integer
00035   static int index(int i=0, int j=0, int k=0)
00036   {
00037     return 9 * (k+1) + 3 * (j+1) + (i+1);
00038   }
00039 
00040   // sets lattice basis vectors but not origin (fixed center)
00041   void set(Vector A, Vector B, Vector C)
00042   {
00043     set(A,B,C,o);
00044   }
00045 
00046   // sets lattice basis vectors and origin (fixed center)
00047   void set(Vector A, Vector B, Vector C, Position Origin)
00048   {
00049     a1 = A; a2 = B; a3 = C; o = Origin;
00050     p1 = ( a1.length2() ? 1 : 0 );
00051     p2 = ( a2.length2() ? 1 : 0 );
00052     p3 = ( a3.length2() ? 1 : 0 );
00053     if ( ! p1 ) a1 = Vector(1.0,0.0,0.0);
00054     if ( ! p2 ) {
00055       Vector u1 = a1 / a1.length();
00056       Vector e_z(0.0,0.0,1.0);
00057       if ( fabs(e_z * u1) < 0.9 ) { a2 = cross(e_z,a1); }
00058       else { a2 = cross(Vector(1.0,0.0,0.0),a1); }
00059       a2 /= a2.length();
00060     }
00061     if ( ! p3 ) {
00062       a3 = cross(a1,a2);
00063       a3 /= a3.length();
00064     }
00065     if ( volume() < 0.0 ) a3 *= -1.0;
00066     recalculate();
00067   }
00068 
00069   // rescale lattice dimensions by factor, origin doesn't move
00070   void rescale(Tensor factor)
00071   {
00072     a1 = factor * a1;
00073     a2 = factor * a2;
00074     a3 = factor * a3;
00075     recalculate();
00076   }
00077 
00078   // rescale a position, keeping origin constant, assume 3D
00079   void rescale(Position &p, Tensor factor) const
00080   {
00081     p -= o;
00082     p = factor * p;
00083     p += o;
00084   }
00085 
00086   // transform scaled position to unscaled position
00087   Position unscale(ScaledPosition s) const
00088   {
00089     return (o + a1*s.x + a2*s.y + a3*s.z);
00090   }
00091 
00092   // transform unscaled position to scaled position
00093   ScaledPosition scale(Position p) const
00094   {
00095     p -= o;
00096     return Vector(b1*p,b2*p,b3*p);
00097   }
00098 
00099   // transforms a position nearest to a SCALED reference position
00100   Position nearest(Position data, ScaledPosition ref) const
00101   {
00102     ScaledPosition sn = scale(data);
00103     if ( p1 ) {
00104       sn.x -= latticenearbyint(sn.x - ref.x);
00105     }
00106     if ( p2 ) {
00107       sn.y -= latticenearbyint(sn.y - ref.y);
00108     }
00109     if ( p3 ) {
00110       sn.z -= latticenearbyint(sn.z - ref.z);
00111     }
00112     return unscale(sn);
00113   }
00114 
00115   // transforms a position nearest to a SCALED reference position
00116   // adds transform for later reversal
00117   Position nearest(Position data, ScaledPosition ref, Transform *t) const
00118   {
00119     ScaledPosition sn = scale(data);
00120     if ( p1 ) {
00121       BigReal tmp = sn.x - ref.x;
00122       BigReal rit = latticenearbyint(tmp);
00123       sn.x -= rit;
00124       t->i -= (int) rit;
00125     }
00126     if ( p2 ) {
00127       BigReal tmp = sn.y - ref.y;
00128       BigReal rit = latticenearbyint(tmp);
00129       sn.y -= rit;
00130       t->j -= (int) rit;
00131     }
00132     if ( p3 ) {
00133       BigReal tmp = sn.z - ref.z;
00134       BigReal rit = latticenearbyint(tmp);
00135       sn.z -= rit;
00136       t->k -= (int) rit;
00137     }
00138     return unscale(sn);
00139   }
00140 
00141   // applies stored transform to original coordinates
00142   Position apply_transform(Position data, const Transform &t) const
00143   {
00144     return ( data + t.i*a1 + t.j*a2 + t.k*a3 );
00145   }
00146 
00147   // reverses cumulative transformations for output
00148   Position reverse_transform(Position data, const Transform &t) const
00149   {
00150     return ( data - t.i*a1 - t.j*a2 - t.k*a3 );
00151   }
00152 
00153   // calculates shortest vector from p2 to p1 (equivalent to p1 - p2)
00154   Vector delta(const Position &pos1, const Position &pos2) const
00155   {
00156     Vector diff = pos1 - pos2;
00157 #ifdef ARCH_POWERPC   //Prevents stack temporaries
00158     Vector result = diff;
00159     if ( p1 ) {
00160       BigReal fval = latticenearbyint(b1*diff); 
00161       result.x -= a1.x *fval;    
00162       result.y -= a1.y *fval;    
00163       result.z -= a1.z *fval;    
00164     }
00165     if ( p2 ) {
00166       BigReal fval = latticenearbyint(b2*diff);
00167       result.x -= a2.x * fval;
00168       result.y -= a2.y * fval;
00169       result.z -= a2.z * fval;
00170     }
00171     if ( p3 ) {
00172       BigReal fval = latticenearbyint(b3*diff);
00173       result.x -= a3.x * fval;
00174       result.y -= a3.y * fval;
00175       result.z -= a3.z * fval;
00176     }
00177     return result;
00178 #else
00179     BigReal f1 = p1 ? latticenearbyint(b1*diff) : 0.;
00180     BigReal f2 = p2 ? latticenearbyint(b2*diff) : 0.;
00181     BigReal f3 = p3 ? latticenearbyint(b3*diff) : 0.;
00182     diff.x -= f1*a1.x + f2*a2.x + f3*a3.x;
00183     diff.y -= f1*a1.y + f2*a2.y + f3*a3.y;
00184     diff.z -= f1*a1.z + f2*a2.z + f3*a3.z;
00185     return diff;
00186 #endif
00187   }
00188 
00189 #ifdef NAMD_CUDA
00190 #ifdef BONDED_CUDA
00191   // calculates scaled vector v such that vector pos1 - pos2 + v*a is the shortest
00192   inline Vector wrap_delta_scaled(const Position &pos1, const Position &pos2) const
00193   {
00194     Vector diff = pos1 - pos2;
00195     Vector result(0.,0.,0.);
00196     if ( p1 ) result.x = -latticenearbyint(b1*diff);
00197     if ( p2 ) result.y = -latticenearbyint(b2*diff);
00198     if ( p3 ) result.z = -latticenearbyint(b3*diff);
00199     return result;
00200   }
00201 #endif
00202 #endif
00203 
00204   // calculates shortest vector from origin to p1 (equivalent to p1 - o)
00205   Vector delta(const Position &pos1) const
00206   {
00207     Vector diff = pos1 - o;
00208     Vector result = diff;
00209     if ( p1 ) result -= a1*latticenearbyint(b1*diff);
00210     if ( p2 ) result -= a2*latticenearbyint(b2*diff);
00211     if ( p3 ) result -= a3*latticenearbyint(b3*diff);
00212     return result;
00213   }
00214 
00215   // calculates vector to bring p1 closest to origin
00216   Vector wrap_delta(const Position &pos1) const
00217   {
00218     Vector diff = pos1 - o;
00219     Vector result(0.,0.,0.);
00220     if ( p1 ) result -= a1*latticenearbyint(b1*diff);
00221     if ( p2 ) result -= a2*latticenearbyint(b2*diff);
00222     if ( p3 ) result -= a3*latticenearbyint(b3*diff);
00223     return result;
00224   }
00225 
00226   // calculates vector to bring p1 closest to origin in unscaled coordinates
00227   Vector wrap_nearest_delta(Position pos1) const
00228   {
00229     Vector diff = pos1 - o;
00230     Vector result0(0.,0.,0.);
00231     if ( p1 ) result0 -= a1*latticenearbyint(b1*diff);
00232     if ( p2 ) result0 -= a2*latticenearbyint(b2*diff);
00233     if ( p3 ) result0 -= a3*latticenearbyint(b3*diff);
00234     diff += result0;
00235     BigReal dist = diff.length2();
00236     Vector result(0.,0.,0.);
00237     for ( int i1=-p1; i1<=p1; ++i1 ) {
00238       for ( int i2 =-p2; i2<=p2; ++i2 ) {
00239         for ( int i3 =-p3; i3<=p3; ++i3 ) {
00240           Vector newresult = i1*a1+i2*a2+i3*a3;
00241           BigReal newdist = (diff+newresult).length2();
00242           if ( newdist < dist ) {
00243             dist = newdist;
00244             result = newresult;
00245           }
00246         }
00247       }
00248     }
00249     return result0 + result;
00250   }
00251 
00252   Vector offset(int i) const
00253   {
00254     return ( (i%3-1) * a1 + ((i/3)%3-1) * a2 + (i/9-1) * a3 );
00255   }
00256 
00257   static int offset_a(int i) { return (i%3-1); }
00258   static int offset_b(int i) { return ((i/3)%3-1); }
00259   static int offset_c(int i) { return (i/9-1); }
00260 
00261   // lattice vectors
00262   Vector a() const { return a1; }
00263   Vector b() const { return a2; }
00264   Vector c() const { return a3; }
00265 
00266   // only if along x y z axes
00267   int orthogonal() const {
00268     return ( ! ( a1.y || a1.z || a2.x || a2.z || a3.x || a3.y ) );
00269   }
00270 
00271   // origin (fixed center of cell)
00272   Vector origin() const
00273   {
00274     return o;
00275   }
00276 
00277   // reciprocal lattice vectors
00278   Vector a_r() const { return b1; }
00279   Vector b_r() const { return b2; }
00280   Vector c_r() const { return b3; }
00281 
00282   // periodic along this direction
00283   int a_p() const { return p1; }
00284   int b_p() const { return p2; }
00285   int c_p() const { return p3; }
00286 
00287   BigReal volume(void) const
00288   {
00289     return ( p1 && p2 && p3 ? cross(a1,a2) * a3 : 0.0 );
00290   }
00291 
00292 private:
00293   Vector a1,a2,a3; // real lattice vectors
00294   Vector b1,b2,b3; // reciprocal lattice vectors (more or less)
00295   Vector o; // origin (fixed center of cell)
00296   int p1, p2, p3; // periodic along this lattice vector?
00297 
00298   // calculate reciprocal lattice vectors
00299   void recalculate(void) {
00300     {
00301       Vector c = cross(a2,a3);
00302       b1 = c / ( a1 * c );
00303     }
00304     {
00305       Vector c = cross(a3,a1);
00306       b2 = c / ( a2 * c );
00307     }
00308     {
00309       Vector c = cross(a1,a2);
00310       b3 = c / ( a3 * c );
00311     }
00312   }
00313 
00314 };
00315 
00316 #include <pup.h>
00317 PUPbytes (Lattice);  
00318 
00319 #endif
00320 

Generated on Tue Sep 19 01:17:12 2017 for NAMD by  doxygen 1.4.7