NAMD
Lattice.h
Go to the documentation of this file.
1 
7 #ifndef LATTICE_H
8 #define LATTICE_H
9 
10 #include <stdlib.h>
11 #include "NamdTypes.h"
12 #include <math.h>
13 #include "Tensor.h"
14 
16 
17 class Lattice
18 {
19 public:
20  Lattice(void) : a1(0,0,0), a2(0,0,0), a3(0,0,0),
21  b1(0,0,0), b2(0,0,0), b3(0,0,0),
22  o(0,0,0), p1(0), p2(0), p3(0) {};
23 
24  // maps a transformation triplet onto a single integer
25  static int index(int i=0, int j=0, int k=0)
26  {
27  return 9 * (k+1) + 3 * (j+1) + (i+1);
28  }
29 
30  // sets lattice basis vectors but not origin (fixed center)
31  void set(Vector A, Vector B, Vector C)
32  {
33  set(A,B,C,o);
34  }
35 
36  // sets lattice basis vectors and origin (fixed center)
38  {
39  a1 = A; a2 = B; a3 = C; o = Origin;
40  p1 = ( a1.length2() ? 1 : 0 );
41  p2 = ( a2.length2() ? 1 : 0 );
42  p3 = ( a3.length2() ? 1 : 0 );
43  if ( ! p1 ) a1 = Vector(1.0,0.0,0.0);
44  if ( ! p2 ) {
45  Vector u1 = a1 / a1.length();
46  Vector e_z(0.0,0.0,1.0);
47  if ( fabs(e_z * u1) < 0.9 ) { a2 = cross(e_z,a1); }
48  else { a2 = cross(Vector(1.0,0.0,0.0),a1); }
49  a2 /= a2.length();
50  }
51  if ( ! p3 ) {
52  a3 = cross(a1,a2);
53  a3 /= a3.length();
54  }
55  if ( volume() < 0.0 ) a3 *= -1.0;
56  recalculate();
57  }
58 
59  // rescale lattice dimensions by factor, origin doesn't move
60  void rescale(Tensor factor)
61  {
62  a1 = factor * a1;
63  a2 = factor * a2;
64  a3 = factor * a3;
65  recalculate();
66  }
67 
68  // rescale a position, keeping origin constant, assume 3D
69  void rescale(Position &p, Tensor factor) const
70  {
71  p -= o;
72  p = factor * p;
73  p += o;
74  }
75 
76  // transform scaled position to unscaled position
78  {
79  return (o + a1*s.x + a2*s.y + a3*s.z);
80  }
81 
82  // transform unscaled position to scaled position
84  {
85  p -= o;
86  return Vector(b1*p,b2*p,b3*p);
87  }
88 
89  // transforms a position nearest to a SCALED reference position
91  {
92  ScaledPosition sn = scale(data);
93  if ( p1 ) {
94  sn.x -= namdnearbyint(sn.x - ref.x);
95  }
96  if ( p2 ) {
97  sn.y -= namdnearbyint(sn.y - ref.y);
98  }
99  if ( p3 ) {
100  sn.z -= namdnearbyint(sn.z - ref.z);
101  }
102  return unscale(sn);
103  }
104 
105  // transforms a position nearest to a SCALED reference position
106  // adds transform for later reversal
108  {
109  ScaledPosition sn = scale(data);
110  if ( p1 ) {
111  BigReal tmp = sn.x - ref.x;
112  BigReal rit = namdnearbyint(tmp);
113  sn.x -= rit;
114  t->i -= (int) rit;
115  }
116  if ( p2 ) {
117  BigReal tmp = sn.y - ref.y;
118  BigReal rit = namdnearbyint(tmp);
119  sn.y -= rit;
120  t->j -= (int) rit;
121  }
122  if ( p3 ) {
123  BigReal tmp = sn.z - ref.z;
124  BigReal rit = namdnearbyint(tmp);
125  sn.z -= rit;
126  t->k -= (int) rit;
127  }
128  return unscale(sn);
129  }
130 
131  // applies stored transform to original coordinates
133  {
134  return ( data + t.i*a1 + t.j*a2 + t.k*a3 );
135  }
136 
137  // reverses cumulative transformations for output
139  {
140  return ( data - t.i*a1 - t.j*a2 - t.k*a3 );
141  }
142 
143  // calculates shortest vector from p2 to p1 (equivalent to p1 - p2)
144  Vector delta(const Position &pos1, const Position &pos2) const
145  {
146  Vector diff = pos1 - pos2;
147 #ifdef ARCH_POWERPC //Prevents stack temporaries
148  Vector result = diff;
149  if ( p1 ) {
150  BigReal fval = namdnearbyint(b1*diff);
151  result.x -= a1.x *fval;
152  result.y -= a1.y *fval;
153  result.z -= a1.z *fval;
154  }
155  if ( p2 ) {
156  BigReal fval = namdnearbyint(b2*diff);
157  result.x -= a2.x * fval;
158  result.y -= a2.y * fval;
159  result.z -= a2.z * fval;
160  }
161  if ( p3 ) {
162  BigReal fval = namdnearbyint(b3*diff);
163  result.x -= a3.x * fval;
164  result.y -= a3.y * fval;
165  result.z -= a3.z * fval;
166  }
167  return result;
168 #else
169  BigReal f1 = p1 ? namdnearbyint(b1*diff) : 0.;
170  BigReal f2 = p2 ? namdnearbyint(b2*diff) : 0.;
171  BigReal f3 = p3 ? namdnearbyint(b3*diff) : 0.;
172  diff.x -= f1*a1.x + f2*a2.x + f3*a3.x;
173  diff.y -= f1*a1.y + f2*a2.y + f3*a3.y;
174  diff.z -= f1*a1.z + f2*a2.z + f3*a3.z;
175  return diff;
176 #endif
177  }
178 
179 #ifdef NAMD_CUDA
180 #ifdef BONDED_CUDA
181  // calculates scaled vector v such that vector pos1 - pos2 + v*a is the shortest
182  inline Vector wrap_delta_scaled(const Position &pos1, const Position &pos2) const
183  {
184  Vector diff = pos1 - pos2;
185  Vector result(0.,0.,0.);
186  if ( p1 ) result.x = -namdnearbyint(b1*diff);
187  if ( p2 ) result.y = -namdnearbyint(b2*diff);
188  if ( p3 ) result.z = -namdnearbyint(b3*diff);
189  return result;
190  }
191 #endif
192 #endif
193 
194  // calculates shortest vector from origin to p1 (equivalent to p1 - o)
195  Vector delta(const Position &pos1) const
196  {
197  Vector diff = pos1 - o;
198  Vector result = diff;
199  if ( p1 ) result -= a1*namdnearbyint(b1*diff);
200  if ( p2 ) result -= a2*namdnearbyint(b2*diff);
201  if ( p3 ) result -= a3*namdnearbyint(b3*diff);
202  return result;
203  }
204 
205  // calculates vector to bring p1 closest to origin
206  Vector wrap_delta(const Position &pos1) const
207  {
208  Vector diff = pos1 - o;
209  Vector result(0.,0.,0.);
210  if ( p1 ) result -= a1*namdnearbyint(b1*diff);
211  if ( p2 ) result -= a2*namdnearbyint(b2*diff);
212  if ( p3 ) result -= a3*namdnearbyint(b3*diff);
213  return result;
214  }
215 
216  // calculates vector to bring p1 closest to origin in unscaled coordinates
218  {
219  Vector diff = pos1 - o;
220  Vector result0(0.,0.,0.);
221  if ( p1 ) result0 -= a1*namdnearbyint(b1*diff);
222  if ( p2 ) result0 -= a2*namdnearbyint(b2*diff);
223  if ( p3 ) result0 -= a3*namdnearbyint(b3*diff);
224  diff += result0;
225  BigReal dist = diff.length2();
226  Vector result(0.,0.,0.);
227  for ( int i1=-p1; i1<=p1; ++i1 ) {
228  for ( int i2 =-p2; i2<=p2; ++i2 ) {
229  for ( int i3 =-p3; i3<=p3; ++i3 ) {
230  Vector newresult = i1*a1+i2*a2+i3*a3;
231  BigReal newdist = (diff+newresult).length2();
232  if ( newdist < dist ) {
233  dist = newdist;
234  result = newresult;
235  }
236  }
237  }
238  }
239  return result0 + result;
240  }
241 
242  Vector offset(int i) const
243  {
244  return ( (i%3-1) * a1 + ((i/3)%3-1) * a2 + (i/9-1) * a3 );
245  }
246 
247  static int offset_a(int i) { return (i%3-1); }
248  static int offset_b(int i) { return ((i/3)%3-1); }
249  static int offset_c(int i) { return (i/9-1); }
250 
251  // lattice vectors
252  Vector a() const { return a1; }
253  Vector b() const { return a2; }
254  Vector c() const { return a3; }
255 
256  // only if along x y z axes
257  int orthogonal() const {
258  return ( ! ( a1.y || a1.z || a2.x || a2.z || a3.x || a3.y ) );
259  }
260 
261  // origin (fixed center of cell)
262  Vector origin() const
263  {
264  return o;
265  }
266 
267  // reciprocal lattice vectors
268  Vector a_r() const { return b1; }
269  Vector b_r() const { return b2; }
270  Vector c_r() const { return b3; }
271 
272  // periodic along this direction
273  int a_p() const { return p1; }
274  int b_p() const { return p2; }
275  int c_p() const { return p3; }
276 
277  BigReal volume(void) const
278  {
279  return ( p1 && p2 && p3 ? cross(a1,a2) * a3 : 0.0 );
280  }
281 
282 private:
283  Vector a1,a2,a3; // real lattice vectors
284  Vector b1,b2,b3; // reciprocal lattice vectors (more or less)
285  Vector o; // origin (fixed center of cell)
286  int p1, p2, p3; // periodic along this lattice vector?
287 
288  // calculate reciprocal lattice vectors
289  void recalculate(void) {
290  {
291  Vector c = cross(a2,a3);
292  b1 = c / ( a1 * c );
293  }
294  {
295  Vector c = cross(a3,a1);
296  b2 = c / ( a2 * c );
297  }
298  {
299  Vector c = cross(a1,a2);
300  b3 = c / ( a3 * c );
301  }
302  }
303 
304 };
305 
306 #include <pup.h>
307 PUPbytes (Lattice);
308 
309 #endif
310 
static int index(int i=0, int j=0, int k=0)
Definition: Lattice.h:25
static int offset_b(int i)
Definition: Lattice.h:248
#define namdnearbyint(x)
Definition: common.h:74
signed char j
Definition: NamdTypes.h:38
static int offset_c(int i)
Definition: Lattice.h:249
Vector a_r() const
Definition: Lattice.h:268
signed char i
Definition: NamdTypes.h:38
Vector delta(const Position &pos1) const
Definition: Lattice.h:195
const BigReal A
Definition: Vector.h:64
Vector c_r() const
Definition: Lattice.h:270
Vector wrap_delta(const Position &pos1) const
Definition: Lattice.h:206
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
BigReal z
Definition: Vector.h:66
int orthogonal() const
Definition: Lattice.h:257
BigReal length(void) const
Definition: Vector.h:169
Vector origin() const
Definition: Lattice.h:262
void set(Vector A, Vector B, Vector C, Position Origin)
Definition: Lattice.h:37
Lattice(void)
Definition: Lattice.h:20
Position nearest(Position data, ScaledPosition ref) const
Definition: Lattice.h:90
Vector b_r() const
Definition: Lattice.h:269
signed char k
Definition: NamdTypes.h:38
Position nearest(Position data, ScaledPosition ref, Transform *t) const
Definition: Lattice.h:107
void rescale(Tensor factor)
Definition: Lattice.h:60
static int offset_a(int i)
Definition: Lattice.h:247
void set(Vector A, Vector B, Vector C)
Definition: Lattice.h:31
AtomID Origin
Definition: ComputeQM.h:110
PUPbytes(Lattice)
Vector delta(const Position &pos1, const Position &pos2) const
Definition: Lattice.h:144
BigReal x
Definition: Vector.h:66
BigReal volume(void) const
Definition: Lattice.h:277
Vector ScaledPosition
Definition: Lattice.h:15
BigReal length2(void) const
Definition: Vector.h:173
Vector offset(int i) const
Definition: Lattice.h:242
Position apply_transform(Position data, const Transform &t) const
Definition: Lattice.h:132
Position unscale(ScaledPosition s) const
Definition: Lattice.h:77
Position reverse_transform(Position data, const Transform &t) const
Definition: Lattice.h:138
Definition: Tensor.h:15
BigReal y
Definition: Vector.h:66
Vector b() const
Definition: Lattice.h:253
const BigReal B
ScaledPosition scale(Position p) const
Definition: Lattice.h:83
void rescale(Position &p, Tensor factor) const
Definition: Lattice.h:69
Vector wrap_nearest_delta(Position pos1) const
Definition: Lattice.h:217
int b_p() const
Definition: Lattice.h:274
int a_p() const
Definition: Lattice.h:273
Vector a() const
Definition: Lattice.h:252
Vector c() const
Definition: Lattice.h:254
double BigReal
Definition: common.h:112
int c_p() const
Definition: Lattice.h:275