Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

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

Generated on Sat Aug 30 04:07:40 2008 for NAMD by  doxygen 1.3.9.1