ComputeStir.C

Go to the documentation of this file.
00001 
00007 /* Barry Isralewitz, July 2001 */
00008 
00009 #include "InfoStream.h" 
00010 #include "ComputeStir.h"
00011 #include "Node.h"
00012 #include "SimParameters.h"
00013 #include "HomePatch.h"
00014 #include "DataStream.h"
00015 #include "Molecule.h"
00016 #include "CollectionMgr.h" //for error corr of buffer problem
00017 
00018 
00019 ComputeStir::ComputeStir(ComputeID c, PatchID pid)
00020   : ComputeHomePatch(c,pid)
00021 {
00022 
00023   
00024 
00025   
00026   SimParameters* simParams = Node::Object()->simParameters;
00027   axisUnit = simParams->stirAxis.unit();
00028   pivot = simParams->stirPivot;
00029   
00030   
00031   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00032   
00033   //nitVars(p, axisDir);
00034   initVars();
00035 }
00036 /*                      END OF FUNCTION ComputeStir             */
00037 
00038 
00039 ComputeStir::~ComputeStir()
00040 
00041 {
00042 delete reduction;
00043 }
00044 /*                      END OF FUNCTION ~ComputeStir            */
00045 
00046 
00047 Tensor ComputeStir::matMult (Tensor a, Tensor b) {
00048   // should really use operator*(Tensor, Tensor)(, but this should be rarely
00049   // used, maybe inline function.  Or extendor Tensor class.
00050   // Also abusing Tensor math concept to mere 3x3 matrix
00051   // just so I can borrow a method from Tensor class. 
00052   Tensor tmp;
00053 
00054   //multiply two 3x3 matrices...
00055 
00056   tmp.xx = a.xx*b.xx + a.xy*b.yx + a.xz*b.zx;
00057   tmp.xy = a.xx*b.xy + a.xy*b.yy + a.xz*b.zy;
00058   tmp.xz = a.xx*b.xz + a.xy*b.yz + a.xz*b.zz;
00059 
00060   tmp.yx = a.yx*b.xx + a.yy*b.yx + a.yz*b.zx;
00061   tmp.yy = a.yx*b.xy + a.yy*b.yy + a.yz*b.zy;
00062   tmp.yz = a.yx*b.xz + a.yy*b.yz + a.yz*b.zz;
00063 
00064   tmp.zx = a.zx*b.xx + a.zy*b.yx + a.zz*b.zx;
00065   tmp.zy = a.zx*b.xy + a.zy*b.yy + a.zz*b.zy;
00066   tmp.zz = a.zx*b.xz + a.zy*b.yz + a.zz*b.zz;
00067   
00068   return tmp;
00069 }
00070 
00071 
00072 void ComputeStir::initVars () {
00073   //find leftMat and rightMat, using 3x3 arrrays for speed,
00074   //so stil must provide translation before/after matrix multuplies
00075   BigReal a,b,c,d;
00076   Tensor testA, testB;
00077   Tensor rotXMat, invRotXMat, rotYMat, invRotYMat;
00078   Tensor temp;
00079   Vector refPos;
00080   SimParameters *simParams = Node::Object()->simParameters;
00081   Molecule *molecule = Node::Object()->molecule;
00082   char statbuf[1024];
00083   const int globalNumAtoms = molecule->numAtoms;
00084 
00085 
00086 
00087  
00088 
00089   //iout << "DEBUG: Starting Init ComputeStir Vars...\n" << endi;
00090   omega = ( (simParams->stirVel) /  360) * (2 * PI);
00091   startingTheta  = ( (simParams->stirStartingTheta) /  360) * (2 * PI);
00092   //iout <<"DEBUG: omega(deg/ts,stirVel)= "<< simParams->stirVel <<"  omega(rad/ts)= "<<omega << "  startingTheta(deg,stirStartingTheta)= "<< simParams->stirStartingTheta <<"  startingTheta(rad)= "<< startingTheta << "\n"<<endi;
00093 
00094   //rotation speed in radians/timestep
00095 
00096   a = axisUnit.x; b = axisUnit.y; c = axisUnit.z;
00097   d = sqrt( (b * b) + (c * c) );
00098   //iout <<"DEBUG: a="<<a<<" b= "<<b<< "c= "<<c<<" d= "<<d<<"\n"<<endi;
00099 
00100   testA.xx = 1.1; testA.xy = 3; testA.xz = 3;
00101   //printTensor(testA);  
00102   //iout << "DEBUG:   axis unit vector = ("  << endi;
00103   //iout  << axisUnit.x << "," <<    axisUnit.y << "," << axisUnit.z << ")\n" << endi;
00105   //yes, could pre-mupltiply these matrices,
00106   //but makes code easier to check
00107   //Besides, only done at startup, so speed not an issue
00108   rotXMat.xx = 1;  rotXMat.xy =   0;  rotXMat.xz =    0;
00109   rotXMat.yx = 0;  rotXMat.yy = c/d;  rotXMat.yz = -b/d;
00110   rotXMat.zx = 0;  rotXMat.zy = b/d;  rotXMat.zz =  c/d;
00111 
00112   
00113   invRotXMat.xx =  1; invRotXMat.xy =     0;  invRotXMat.xz =   0;
00114   invRotXMat.yx =  0; invRotXMat.yy =   c/d;  invRotXMat.yz = b/d;
00115   invRotXMat.zx =  0; invRotXMat.zy =  -b/d;  invRotXMat.zz = c/d;
00116 
00117 
00118   rotYMat.xx = d; rotYMat.xy = 0; rotYMat.xz = -a;
00119   rotYMat.yx = 0; rotYMat.yy = 1;  rotYMat.yz = 0;
00120   rotYMat.zx = a; rotYMat.zy = 0;  rotYMat.zz = d;
00121 
00122 
00123   invRotYMat.xx = d;   invRotYMat.xy =  0;      invRotYMat.xz = a;
00124   invRotYMat.yx = 0;   invRotYMat.yy =  1;      invRotYMat.yz = 0;
00125   invRotYMat.zx =-a;   invRotYMat.zy =  0;      invRotYMat.zz = d;
00126 
00127   temp = rotYMat;
00128   //iout << "DEBUG: Mat check rotYMat\n" << endi;
00129   //printTensor (temp);
00130  
00131   temp = invRotYMat;
00132   //iout << "DEBUG: Mat check invRotYMat\n" << endi;
00133   //printTensor (temp);
00134 
00135   temp = matMult (invRotYMat, rotYMat);
00136   //iout << "DEBUG: Mat check Y^-1 Y\n" << endi;
00137   //printTensor (temp);
00138 
00139     temp = matMult (invRotXMat, rotXMat);
00140   //iout << "DEBUG: Mat check X^-1 X\n" << endi;
00141   //printTensor (temp);
00142   
00143 
00144   //tmp.xx =  ; tmp.xy =   ;  tmp.xz = 
00145   //tmp.yx = ;  tmp.yy =  ;  tmp.yz =
00146   //tmp.zx =  ; tmp.zy =  ;  tmp.zz = 
00147  
00148   leftMat = matMult (invRotXMat, invRotYMat);
00149   //translate after use
00150   //e.g. :  matVec(leftMat,p) ) - offset
00151   rightMat = matMult ( rotYMat, rotXMat) ;
00152   //translate before use
00153   //e.g. :  matVec(rightMat,p+offset))
00154 
00155   temp =  leftMat;
00156   //iout << "DEBUG: Mat check leftMat \n" << endi;
00157   //printTensor (temp); 
00158 
00159      temp = rightMat;
00160   //iout << "DEBUG: Mat check rightMat \n" << endi;
00161   //printTensor (temp); 
00162 
00163   temp = matMult (leftMat,rightMat);
00164   //iout << "DEBUG: Mat check matMult (leftMat,rightMat) \n" << endi;
00165   //printTensor (temp); 
00166   
00167   temp = arbRotMat (0.34);
00168   //iout << "DEBUG: Mat check arbRotMat (0.34) \n" << endi;
00169   //printTensor (temp);   
00170 
00171   //now make table of initial thetas
00172   // (perhaps could do this back in Molecule.C)
00173   for (int i=0; i<globalNumAtoms; i++) {
00174       if (molecule->is_atom_stirred(i))
00175       {
00176         //CkPrintf ("DEBUG: now to get stir params");   
00177         
00178         //use same index as params
00179         molecule->get_stir_refPos (refPos, i);
00180 
00181         molecule->put_stir_startTheta (findTheta(refPos), i); 
00182         //CkPrintf ("DEBUG: now to get sprintf x y z's");
00183         //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; 
00184           // CollectionMgr::Object()->sendDataStream(statbuf);  
00185       
00186       }
00187 
00188   }
00189 }
00190   
00191 
00192 BigReal ComputeStir::findTheta (Vector pos) {
00193   BigReal theta;
00194   //translate to origin and align to z axis
00195   Vector aligned_pos = rightMat * (pos - pivot);
00196   
00197   if ( ! ((aligned_pos.x == 0 ) && (aligned_pos.y  == 0) )) {
00198     //we have our rotation counter-clockwise around z-axis (right hand rule
00199     //along increasing z
00200     theta = atan (aligned_pos.y/aligned_pos.x);
00201     if (aligned_pos.x < 0) {
00202       theta = theta + PI;
00203       }
00204   } else { 
00205     //iout << iWARN << "stirred atom exactly along stir axis, setting theta of this atom to 0.0" << endi;
00206     theta = 0.0;
00207   }
00208   
00209   
00210   return theta;
00211 }
00212 
00213 Tensor ComputeStir::arbRotMat (BigReal theta) {
00214   //needed only to show graphics in output
00215   //not called by main doForce
00216   //must use +pivot, -pivot after/before 
00217   Tensor rotZMat;
00218   Tensor temp;
00219   rotZMat.xx = cos (theta) ; rotZMat.xy = -sin (theta);  rotZMat.xz = 0;
00220   rotZMat.yx = sin (theta);  rotZMat.yy =  cos(theta);   rotZMat.yz = 0;
00221   rotZMat.zx =  0;           rotZMat.zy =      0;        rotZMat.zz = 1;
00222   //iout <<"DEBUG:: in arbRotMat, print Tensor rotZMat\n"<<endi;
00223   //printTensor (rotZMat);
00224   //iout <<"DEBUG:: in arbRotMat, print Tensor rotZMat*rightMat\n"<<endi;
00225   temp =matMult (rotZMat, rightMat);
00226   //printTensor (temp);
00227 
00228   temp = matMult (leftMat,rightMat);
00229   //iout << "DEBUG: Mat check matMult (leftMat,rightMat) \n" << endi;
00230   //printTensor (temp); 
00231 
00232   temp = matMult (leftMat, (matMult (rotZMat, rightMat) ) ) ;
00233 ;
00234  //iout <<"DEBUG:: in arbRotMat, print Tensor temp(should be whats returned)\n"<<endi;
00235  //printTensor (temp); 
00236 
00237 
00238 
00239   return matMult (leftMat, (matMult (rotZMat, rightMat) ) );
00240 }
00241 
00242 void ComputeStir::printTensor (Tensor &t) {
00243   //just for debugging
00244   iout << "DEBUG: Tensor =" << "\n" << endi;
00245   iout <<"DEBUG:   " << t.xx << " " <<t.xy  << " " << t.xz << \
00246     "\nDEBUG:   " << t.yx << " " << t.yy << " " << t.yz << \
00247     "\nDEBUG:   " << t.zx << " " << t.zy << " " <<t.zz << "\n"  << endi;
00248 }
00249 
00250 BigReal ComputeStir::distanceToRay(Vector dir,Vector origin,Vector x) {
00251           return (x - projectionOnRay(dir, origin,x)).length();
00252 }
00253 BigReal ComputeStir::findHeight(Vector x) {
00254   //rotate so axisUnitVec is along z-axis
00255   //we have our rotation theta counter-clockwise around z (right hand rule along increasing z)
00256   Vector pos =rightMat * (x - pivot);  
00257   //so, just read height as z
00258   return pos.z;
00259 };
00260 
00261 Vector ComputeStir::projectionOnRay (Vector dir, Vector origin, Vector x) {
00262   return origin + ( ( ((x - origin ).dot(dir) )/    \
00263                      (dir.dot(dir))                )
00264                                * dir);
00265 };
00266 
00267 Vector ComputeStir::placeThetaRadius (BigReal theta, BigReal height, BigReal radius) {
00268   Vector a; 
00269   Tensor rotZMat;
00271   rotZMat.xx = cos (theta) ; rotZMat.xy = -sin (theta);  rotZMat.xz = 0;
00272   rotZMat.yx = sin (theta);  rotZMat.yy =  cos(theta);   rotZMat.yz = 0;
00273   rotZMat.zx =  0;           rotZMat.zy =      0;        rotZMat.zz = 1;
00277   //set position to (radius,0,height)  
00278   a.x = radius; a.y =0; a.z = height;
00282   //rotate it by theta, then rotate it to axis system, then transl. to axis sytem
00283   a = (matMult(leftMat, rotZMat) * a) +pivot;
00284   
00286   return a;
00287 }
00288 void ComputeStir::doForce (FullAtom* p, Results* r) {
00289   
00290   
00291   Molecule *molecule = Node::Object()->molecule;
00292   SimParameters *simParams = Node::Object()->simParameters;
00293   Lattice &lattice = homePatch->lattice;
00294   char statbuf[1024];
00295   Vector curPos; //current atom position
00296   Vector proj;  //projection on the axis ray
00297   Vector forceDir;  //direction to apply force in for torque
00298   Vector theForce; // the force to apply
00299   Vector targPos; //target Position
00300   BigReal rayDist, height,targTheta,forceMag;
00301   //intermediate calc. results for readability
00302 
00303   //convert from degrees/timestep to radians/timestep
00304   int currentStep = patch->flags.step;
00305 
00306   //Vector  = simParams->eField;
00308 
00309   int GlobalId;
00310   Force *forces = r->f[Results::normal];
00311   BigReal energy = 0;
00312   Force extForce = 0.;
00313   Tensor extVirial;
00314   
00315   //CkPrintf("DEBUG: In ComputeStir::doForce");
00316   //  Loop through and check each atom
00317   //CkPrintf ("DEBUG: now to loop atoms, numAtoms = %d\n",numAtoms);    
00318     for (int i=0; i<numAtoms; i++) {
00319       if (molecule->is_atom_stirred(p[i].id))
00320       {
00321         //CkPrintf ("DEBUG: now to get stir params");   
00322 
00323 
00324         molecule->get_stir_startTheta (p[i].id);
00325         
00326         curPos = lattice.reverse_transform(p[i].position,p[i].transform);
00327 
00328         rayDist = distanceToRay(axisUnit,pivot,curPos);
00329         proj = projectionOnRay (axisUnit,pivot,curPos);
00330         height = findHeight (proj);
00331         targTheta = (omega * currentStep) + \
00332                      startingTheta + \
00333                       molecule->get_stir_startTheta (p[i].id);
00334 
00335         targPos = placeThetaRadius (targTheta,height,rayDist);
00336 
00337         //project the actual displacement ontoh a unit vector along ( (axis) cross (moment arm) ),
00338         //so we are applying only  a torque
00339       
00340         forceDir = (cross(axisUnit, (curPos - proj))).unit();
00341         forceMag = forceDir.dot(targPos - curPos); 
00342         
00343                        //unitDisp = disp.unit();
00344         
00345         theForce = (simParams->stirK) * forceMag * forceDir;
00346         
00347         forces[i]  += theForce;
00348         extForce += theForce;
00349         extVirial += outer(theForce,curPos);
00350 
00351                              
00352           //CkPrintf ("DEBUG: now to get sprintf x y z's");
00353           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());
00354 
00355         
00356           CollectionMgr::Object()->sendDataStream(statbuf);     
00357           
00359 
00360       }
00361 
00362       //CkPrintf ("DEBUG: got past check for stirred\n");       
00363       //GlobalId = p[i].id;
00364       //sprintf (statbuf, "DEBUG: i= %d Global= %d\n",i,GlobalId);
00365       //dout << "DEBUG: i= " << i << " Global=" << GlobalId <<"\n"<<endd;
00366       //CollectionMgr::Object()->sendDataStream(statbuf);
00367     
00368     }
00369 
00370   reduction->item(REDUCTION_MISC_ENERGY) += energy;
00371   ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,extForce);
00372   ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,extVirial);
00373   reduction->submit();
00374 
00375 }
00376 /*                      END OF FUNCTION force                           */
00377  

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