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

Generated on Tue Nov 21 01:17:13 2017 for NAMD by  doxygen 1.4.7