NAMD
ComputeStir.C
Go to the documentation of this file.
1 
7 /* Barry Isralewitz, July 2001 */
8 
9 #include "InfoStream.h"
10 #include "ComputeStir.h"
11 #include "Node.h"
12 #include "SimParameters.h"
13 #include "HomePatch.h"
14 #include "DataStream.h"
15 #include "Molecule.h"
16 #include "CollectionMgr.h" //for error corr of buffer problem
17 
18 
20  : ComputeHomePatch(c,pid)
21 {
22 
23 
24 
25 
27  axisUnit = simParams->stirAxis.unit();
28  pivot = simParams->stirPivot;
29 
30 
32 
33  //nitVars(p, axisDir);
34  initVars();
35 }
36 /* END OF FUNCTION ComputeStir */
37 
38 
40 
41 {
42 delete reduction;
43 }
44 /* END OF FUNCTION ~ComputeStir */
45 
46 
47 Tensor ComputeStir::matMult (Tensor a, Tensor b) {
48  // should really use operator*(Tensor, Tensor)(, but this should be rarely
49  // used, maybe inline function. Or extendor Tensor class.
50  // Also abusing Tensor math concept to mere 3x3 matrix
51  // just so I can borrow a method from Tensor class.
52  Tensor tmp;
53 
54  //multiply two 3x3 matrices...
55 
56  tmp.xx = a.xx*b.xx + a.xy*b.yx + a.xz*b.zx;
57  tmp.xy = a.xx*b.xy + a.xy*b.yy + a.xz*b.zy;
58  tmp.xz = a.xx*b.xz + a.xy*b.yz + a.xz*b.zz;
59 
60  tmp.yx = a.yx*b.xx + a.yy*b.yx + a.yz*b.zx;
61  tmp.yy = a.yx*b.xy + a.yy*b.yy + a.yz*b.zy;
62  tmp.yz = a.yx*b.xz + a.yy*b.yz + a.yz*b.zz;
63 
64  tmp.zx = a.zx*b.xx + a.zy*b.yx + a.zz*b.zx;
65  tmp.zy = a.zx*b.xy + a.zy*b.yy + a.zz*b.zy;
66  tmp.zz = a.zx*b.xz + a.zy*b.yz + a.zz*b.zz;
67 
68  return tmp;
69 }
70 
71 
72 void ComputeStir::initVars () {
73  //find leftMat and rightMat, using 3x3 arrrays for speed,
74  //so stil must provide translation before/after matrix multuplies
75  BigReal a,b,c,d;
76  Tensor testA, testB;
77  Tensor rotXMat, invRotXMat, rotYMat, invRotYMat;
78  Tensor temp;
79  Vector refPos;
81  Molecule *molecule = Node::Object()->molecule;
82  char statbuf[1024];
83  const int globalNumAtoms = molecule->numAtoms;
84 
85 
86 
87 
88 
89  //iout << "DEBUG: Starting Init ComputeStir Vars...\n" << endi;
90  omega = ( (simParams->stirVel) / 360) * (2 * PI);
91  startingTheta = ( (simParams->stirStartingTheta) / 360) * (2 * PI);
92  //iout <<"DEBUG: omega(deg/ts,stirVel)= "<< simParams->stirVel <<" omega(rad/ts)= "<<omega << " startingTheta(deg,stirStartingTheta)= "<< simParams->stirStartingTheta <<" startingTheta(rad)= "<< startingTheta << "\n"<<endi;
93 
94  //rotation speed in radians/timestep
95 
96  a = axisUnit.x; b = axisUnit.y; c = axisUnit.z;
97  d = sqrt( (b * b) + (c * c) );
98  //iout <<"DEBUG: a="<<a<<" b= "<<b<< "c= "<<c<<" d= "<<d<<"\n"<<endi;
99 
100  testA.xx = 1.1; testA.xy = 3; testA.xz = 3;
101  //printTensor(testA);
102  //iout << "DEBUG: axis unit vector = (" << endi;
103  //iout << axisUnit.x << "," << axisUnit.y << "," << axisUnit.z << ")\n" << endi;
105  //yes, could pre-mupltiply these matrices,
106  //but makes code easier to check
107  //Besides, only done at startup, so speed not an issue
108  rotXMat.xx = 1; rotXMat.xy = 0; rotXMat.xz = 0;
109  rotXMat.yx = 0; rotXMat.yy = c/d; rotXMat.yz = -b/d;
110  rotXMat.zx = 0; rotXMat.zy = b/d; rotXMat.zz = c/d;
111 
112 
113  invRotXMat.xx = 1; invRotXMat.xy = 0; invRotXMat.xz = 0;
114  invRotXMat.yx = 0; invRotXMat.yy = c/d; invRotXMat.yz = b/d;
115  invRotXMat.zx = 0; invRotXMat.zy = -b/d; invRotXMat.zz = c/d;
116 
117 
118  rotYMat.xx = d; rotYMat.xy = 0; rotYMat.xz = -a;
119  rotYMat.yx = 0; rotYMat.yy = 1; rotYMat.yz = 0;
120  rotYMat.zx = a; rotYMat.zy = 0; rotYMat.zz = d;
121 
122 
123  invRotYMat.xx = d; invRotYMat.xy = 0; invRotYMat.xz = a;
124  invRotYMat.yx = 0; invRotYMat.yy = 1; invRotYMat.yz = 0;
125  invRotYMat.zx =-a; invRotYMat.zy = 0; invRotYMat.zz = d;
126 
127  temp = rotYMat;
128  //iout << "DEBUG: Mat check rotYMat\n" << endi;
129  //printTensor (temp);
130 
131  temp = invRotYMat;
132  //iout << "DEBUG: Mat check invRotYMat\n" << endi;
133  //printTensor (temp);
134 
135  temp = matMult (invRotYMat, rotYMat);
136  //iout << "DEBUG: Mat check Y^-1 Y\n" << endi;
137  //printTensor (temp);
138 
139  temp = matMult (invRotXMat, rotXMat);
140  //iout << "DEBUG: Mat check X^-1 X\n" << endi;
141  //printTensor (temp);
142 
143 
144  //tmp.xx = ; tmp.xy = ; tmp.xz =
145  //tmp.yx = ; tmp.yy = ; tmp.yz =
146  //tmp.zx = ; tmp.zy = ; tmp.zz =
147 
148  leftMat = matMult (invRotXMat, invRotYMat);
149  //translate after use
150  //e.g. : matVec(leftMat,p) ) - offset
151  rightMat = matMult ( rotYMat, rotXMat) ;
152  //translate before use
153  //e.g. : matVec(rightMat,p+offset))
154 
155  temp = leftMat;
156  //iout << "DEBUG: Mat check leftMat \n" << endi;
157  //printTensor (temp);
158 
159  temp = rightMat;
160  //iout << "DEBUG: Mat check rightMat \n" << endi;
161  //printTensor (temp);
162 
163  temp = matMult (leftMat,rightMat);
164  //iout << "DEBUG: Mat check matMult (leftMat,rightMat) \n" << endi;
165  //printTensor (temp);
166 
167  temp = arbRotMat (0.34);
168  //iout << "DEBUG: Mat check arbRotMat (0.34) \n" << endi;
169  //printTensor (temp);
170 
171  //now make table of initial thetas
172  // (perhaps could do this back in Molecule.C)
173  for (int i=0; i<globalNumAtoms; i++) {
174  if (molecule->is_atom_stirred(i))
175  {
176  //CkPrintf ("DEBUG: now to get stir params");
177 
178  //use same index as params
179  molecule->get_stir_refPos (refPos, i);
180 
181  molecule->put_stir_startTheta (findTheta(refPos), i);
182  //CkPrintf ("DEBUG: now to get sprintf x y z's");
183  //iout << "DEBUG: For atom i="<<i<<" refPos.x= " << refPos.x << " refPos.y= "<<refPos.y<<" refPos.z "<< refPos.z << " theta= " << molecule->get_stir_startTheta(i)<<"\n"<<endi;
184  // CollectionMgr::Object()->sendDataStream(statbuf);
185 
186  }
187 
188  }
189 }
190 
191 
192 BigReal ComputeStir::findTheta (Vector pos) {
193  BigReal theta;
194  //translate to origin and align to z axis
195  Vector aligned_pos = rightMat * (pos - pivot);
196 
197  if ( ! ((aligned_pos.x == 0 ) && (aligned_pos.y == 0) )) {
198  //we have our rotation counter-clockwise around z-axis (right hand rule
199  //along increasing z
200  theta = atan (aligned_pos.y/aligned_pos.x);
201  if (aligned_pos.x < 0) {
202  theta = theta + PI;
203  }
204  } else {
205  //iout << iWARN << "stirred atom exactly along stir axis, setting theta of this atom to 0.0" << endi;
206  theta = 0.0;
207  }
208 
209 
210  return theta;
211 }
212 
213 Tensor ComputeStir::arbRotMat (BigReal theta) {
214  //needed only to show graphics in output
215  //not called by main doForce
216  //must use +pivot, -pivot after/before
217  Tensor rotZMat;
218  Tensor temp;
219  rotZMat.xx = cos (theta) ; rotZMat.xy = -sin (theta); rotZMat.xz = 0;
220  rotZMat.yx = sin (theta); rotZMat.yy = cos(theta); rotZMat.yz = 0;
221  rotZMat.zx = 0; rotZMat.zy = 0; rotZMat.zz = 1;
222  //iout <<"DEBUG:: in arbRotMat, print Tensor rotZMat\n"<<endi;
223  //printTensor (rotZMat);
224  //iout <<"DEBUG:: in arbRotMat, print Tensor rotZMat*rightMat\n"<<endi;
225  temp =matMult (rotZMat, rightMat);
226  //printTensor (temp);
227 
228  temp = matMult (leftMat,rightMat);
229  //iout << "DEBUG: Mat check matMult (leftMat,rightMat) \n" << endi;
230  //printTensor (temp);
231 
232  temp = matMult (leftMat, (matMult (rotZMat, rightMat) ) ) ;
233 ;
234  //iout <<"DEBUG:: in arbRotMat, print Tensor temp(should be whats returned)\n"<<endi;
235  //printTensor (temp);
236 
237 
238 
239  return matMult (leftMat, (matMult (rotZMat, rightMat) ) );
240 }
241 
242 void ComputeStir::printTensor (Tensor &t) {
243  //just for debugging
244  iout << "DEBUG: Tensor =" << "\n" << endi;
245  iout <<"DEBUG: " << t.xx << " " <<t.xy << " " << t.xz << \
246  "\nDEBUG: " << t.yx << " " << t.yy << " " << t.yz << \
247  "\nDEBUG: " << t.zx << " " << t.zy << " " <<t.zz << "\n" << endi;
248 }
249 
250 BigReal ComputeStir::distanceToRay(Vector dir,Vector origin,Vector x) {
251  return (x - projectionOnRay(dir, origin,x)).length();
252 }
253 BigReal ComputeStir::findHeight(Vector x) {
254  //rotate so axisUnitVec is along z-axis
255  //we have our rotation theta counter-clockwise around z (right hand rule along increasing z)
256  Vector pos =rightMat * (x - pivot);
257  //so, just read height as z
258  return pos.z;
259 };
260 
261 Vector ComputeStir::projectionOnRay (Vector dir, Vector origin, Vector x) {
262  return origin + ( ( ((x - origin ).dot(dir) )/ \
263  (dir.dot(dir)) )
264  * dir);
265 };
266 
267 Vector ComputeStir::placeThetaRadius (BigReal theta, BigReal height, BigReal radius) {
268  Vector a;
269  Tensor rotZMat;
271  rotZMat.xx = cos (theta) ; rotZMat.xy = -sin (theta); rotZMat.xz = 0;
272  rotZMat.yx = sin (theta); rotZMat.yy = cos(theta); rotZMat.yz = 0;
273  rotZMat.zx = 0; rotZMat.zy = 0; rotZMat.zz = 1;
277  //set position to (radius,0,height)
278  a.x = radius; a.y =0; a.z = height;
282  //rotate it by theta, then rotate it to axis system, then transl. to axis sytem
283  a = (matMult(leftMat, rotZMat) * a) +pivot;
284 
286  return a;
287 }
289 
290 
291  Molecule *molecule = Node::Object()->molecule;
292  SimParameters *simParams = Node::Object()->simParameters;
293  Lattice &lattice = homePatch->lattice;
294  char statbuf[1024];
295  Vector curPos; //current atom position
296  Vector proj; //projection on the axis ray
297  Vector forceDir; //direction to apply force in for torque
298  Vector theForce; // the force to apply
299  Vector targPos; //target Position
300  BigReal rayDist, height,targTheta,forceMag;
301  //intermediate calc. results for readability
302 
303  //convert from degrees/timestep to radians/timestep
304  int currentStep = patch->flags.step;
305 
306  //Vector = simParams->eField;
308 
309  int GlobalId;
310  Force *forces = r->f[Results::normal];
311  BigReal energy = 0;
312  Force extForce = 0.;
313  Tensor extVirial;
314 
315  //CkPrintf("DEBUG: In ComputeStir::doForce");
316  // Loop through and check each atom
317  //CkPrintf ("DEBUG: now to loop atoms, numAtoms = %d\n",numAtoms);
318  for (int i=0; i<numAtoms; i++) {
319  if (molecule->is_atom_stirred(p[i].id))
320  {
321  //CkPrintf ("DEBUG: now to get stir params");
322 
323 
324  molecule->get_stir_startTheta (p[i].id);
325 
326  curPos = lattice.reverse_transform(p[i].position,p[i].transform);
327 
328  rayDist = distanceToRay(axisUnit,pivot,curPos);
329  proj = projectionOnRay (axisUnit,pivot,curPos);
330  height = findHeight (proj);
331  targTheta = (omega * currentStep) + \
332  startingTheta + \
333  molecule->get_stir_startTheta (p[i].id);
334 
335  targPos = placeThetaRadius (targTheta,height,rayDist);
336 
337  //project the actual displacement ontoh a unit vector along ( (axis) cross (moment arm) ),
338  //so we are applying only a torque
339 
340  forceDir = (cross(axisUnit, (curPos - proj))).unit();
341  forceMag = forceDir.dot(targPos - curPos);
342 
343  //unitDisp = disp.unit();
344 
345  theForce = (simParams->stirK) * forceMag * forceDir;
346 
347  forces[i] += theForce;
348  extForce += theForce;
349  extVirial += outer(theForce,curPos);
350 
351 
352  //CkPrintf ("DEBUG: now to get sprintf x y z's");
353  sprintf (statbuf, "DEBUG: step= %d atom i= %d globID= %d T= %g %g %g C= %g %g %g F= %g %g %g startTheta= %g targTheta= %g forceMag= %g F.len= %g FC.len= %g\n",currentStep,i,p[i].id, targPos.x, targPos.y, targPos.z, curPos.x, curPos.y, curPos.z,theForce.x, theForce.y, theForce.z, molecule->get_stir_startTheta (p[i].id), targTheta, forceMag, theForce.length(), forces[i].length());
354 
355 
357 
359 
360  }
361 
362  //CkPrintf ("DEBUG: got past check for stirred\n");
363  //GlobalId = p[i].id;
364  //sprintf (statbuf, "DEBUG: i= %d Global= %d\n",i,GlobalId);
365  //dout << "DEBUG: i= " << i << " Global=" << GlobalId <<"\n"<<endd;
366  //CollectionMgr::Object()->sendDataStream(statbuf);
367 
368  }
369 
371  ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,extForce);
372  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,extVirial);
373  reduction->submit();
374 
375 }
376 /* END OF FUNCTION force */
377 
static Node * Object()
Definition: Node.h:86
static CollectionMgr * Object()
Definition: CollectionMgr.h:30
BigReal zy
Definition: Tensor.h:19
void get_stir_refPos(Vector &refPos, int atomnum) const
Definition: Molecule.h:1302
BigReal xz
Definition: Tensor.h:17
int ComputeID
Definition: NamdTypes.h:183
Lattice & lattice
Definition: Patch.h:126
Definition: Vector.h:64
#define ADD_TENSOR_OBJECT(R, RL, D)
Definition: ReductionMgr.h:43
SimParameters * simParameters
Definition: Node.h:178
virtual void doForce(FullAtom *p, Results *r)
Definition: ComputeStir.C:288
static __thread float4 * forces
BigReal & item(int i)
Definition: ReductionMgr.h:312
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
BigReal z
Definition: Vector.h:66
Real get_stir_startTheta(int atomnum) const
Definition: Molecule.h:1314
BigReal yz
Definition: Tensor.h:18
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
zVector stirPivot
SubmitReduction * reduction
Definition: ComputeStir.h:27
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
#define iout
Definition: InfoStream.h:87
BigReal length(void) const
Definition: Vector.h:169
Flags flags
Definition: Patch.h:127
zVector stirAxis
#define PI
Definition: common.h:81
BigReal dot(const Vector &v2)
Definition: Vector.h:222
Bool is_atom_stirred(int atomnum) const
Definition: Molecule.h:1428
BigReal yx
Definition: Tensor.h:18
Force * f[maxNumForces]
Definition: PatchTypes.h:67
ComputeStir(ComputeID c, PatchID pid)
Definition: ComputeStir.C:19
BigReal x
Definition: Vector.h:66
int numAtoms
Definition: Molecule.h:556
int PatchID
Definition: NamdTypes.h:182
void put_stir_startTheta(Real theta, int atomnum) const
Definition: Molecule.h:1308
Tensor outer(const Vector &v1, const Vector &v2)
Definition: Tensor.h:241
BigReal xx
Definition: Tensor.h:17
__device__ __forceinline__ float dot(const float3 v1, const float3 v2)
BigReal stirStartingTheta
BigReal zz
Definition: Tensor.h:19
virtual ~ComputeStir()
Definition: ComputeStir.C:39
#define simParams
Definition: Output.C:127
Position reverse_transform(Position data, const Transform &t) const
Definition: Lattice.h:138
Definition: Tensor.h:15
BigReal xy
Definition: Tensor.h:17
BigReal y
Definition: Vector.h:66
#define ADD_VECTOR_OBJECT(R, RL, D)
Definition: ReductionMgr.h:27
BigReal yy
Definition: Tensor.h:18
void submit(void)
Definition: ReductionMgr.h:323
infostream & endi(infostream &s)
Definition: InfoStream.C:38
gridSize x
BigReal zx
Definition: Tensor.h:19
Molecule * molecule
Definition: Node.h:176
void sendDataStream(const char *)
Vector unit(void) const
Definition: Vector.h:182
HomePatch * homePatch
double BigReal
Definition: common.h:112
int step
Definition: PatchTypes.h:16