ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/primitives/RigidBody.cpp
(Generate patch)

Comparing trunk/src/primitives/RigidBody.cpp (file contents):
Revision 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 1424 by gezelter, Tue Mar 30 15:05:38 2010 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 + * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41   #include <algorithm>
42   #include <math.h>
43   #include "primitives/RigidBody.hpp"
44   #include "utils/simError.h"
45   #include "utils/NumericConstant.hpp"
46 < namespace oopse {
47 <
48 <  RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData), inertiaTensor_(0.0){
49 <
46 > namespace OpenMD {
47 >  
48 >  RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData),
49 >                           inertiaTensor_(0.0){    
50    }
51 <
51 >  
52    void RigidBody::setPrevA(const RotMat3x3d& a) {
53      ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
54 <    //((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
55 <
54 >    
55      for (int i =0 ; i < atoms_.size(); ++i){
56        if (atoms_[i]->isDirectional()) {
57 <        atoms_[i]->setPrevA(a * refOrients_[i]);
57 >        atoms_[i]->setPrevA(refOrients_[i].transpose() * a);
58        }
59      }
60 <
60 >    
61    }
62 <
63 <      
62 >  
63 >  
64    void RigidBody::setA(const RotMat3x3d& a) {
65      ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
67 –    //((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
66  
67      for (int i =0 ; i < atoms_.size(); ++i){
68        if (atoms_[i]->isDirectional()) {
69 <        atoms_[i]->setA(a * refOrients_[i]);
69 >        atoms_[i]->setA(refOrients_[i].transpose() * a);
70        }
71      }
72    }    
73 <    
73 >  
74    void RigidBody::setA(const RotMat3x3d& a, int snapshotNo) {
75      ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
76 +    
77      //((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;    
78 <
78 >    
79      for (int i =0 ; i < atoms_.size(); ++i){
80        if (atoms_[i]->isDirectional()) {
81 <        atoms_[i]->setA(a * refOrients_[i], snapshotNo);
81 >        atoms_[i]->setA(refOrients_[i].transpose() * a, snapshotNo);
82        }
83      }
84 <
84 >    
85    }  
86 <
86 >  
87    Mat3x3d RigidBody::getI() {
88      return inertiaTensor_;
89    }    
90 <
91 <  std::vector<double> RigidBody::getGrad() {
92 <    std::vector<double> grad(6, 0.0);
90 >  
91 >  std::vector<RealType> RigidBody::getGrad() {
92 >    std::vector<RealType> grad(6, 0.0);
93      Vector3d force;
94      Vector3d torque;
95      Vector3d myEuler;
96 <    double phi, theta, psi;
97 <    double cphi, sphi, ctheta, stheta;
96 >    RealType phi, theta, psi;
97 >    RealType cphi, sphi, ctheta, stheta;
98      Vector3d ephi;
99      Vector3d etheta;
100      Vector3d epsi;
101 <
101 >    
102      force = getFrc();
103      torque =getTrq();
104      myEuler = getA().toEulerAngles();
105 <
105 >    
106      phi = myEuler[0];
107      theta = myEuler[1];
108      psi = myEuler[2];
109 <
109 >    
110      cphi = cos(phi);
111      sphi = sin(phi);
112      ctheta = cos(theta);
113      stheta = sin(theta);
114 <
114 >    
115      // get unit vectors along the phi, theta and psi rotation axes
116 <
116 >    
117      ephi[0] = 0.0;
118      ephi[1] = 0.0;
119      ephi[2] = 1.0;
120 <
120 >    
121 >    //etheta[0] = -sphi;
122 >    //etheta[1] =  cphi;
123 >    //etheta[2] =  0.0;
124 >    
125      etheta[0] = cphi;
126      etheta[1] = sphi;
127 <    etheta[2] = 0.0;
128 <
127 >    etheta[2] =  0.0;
128 >    
129      epsi[0] = stheta * cphi;
130      epsi[1] = stheta * sphi;
131      epsi[2] = ctheta;
132 <
132 >    
133      //gradient is equal to -force
134      for (int j = 0 ; j<3; j++)
135        grad[j] = -force[j];
136 <
136 >    
137      for (int j = 0; j < 3; j++ ) {
138 <
138 >      
139        grad[3] += torque[j]*ephi[j];
140        grad[4] += torque[j]*etheta[j];
141        grad[5] += torque[j]*epsi[j];
142 <
142 >      
143      }
144      
145      return grad;
146    }    
147 <
147 >  
148    void RigidBody::accept(BaseVisitor* v) {
149      v->visit(this);
150    }    
151  
152    /**@todo need modification */
153    void  RigidBody::calcRefCoords() {
154 <    double mtmp;
154 >    RealType mtmp;
155      Vector3d refCOM(0.0);
156      mass_ = 0.0;
157      for (std::size_t i = 0; i < atoms_.size(); ++i) {
# Line 157 | Line 160 | namespace oopse {
160        refCOM += refCoords_[i]*mtmp;
161      }
162      refCOM /= mass_;
163 <
163 >    
164      // Next, move the origin of the reference coordinate system to the COM:
165      for (std::size_t i = 0; i < atoms_.size(); ++i) {
166        refCoords_[i] -= refCOM;
167      }
168  
169      // Moment of Inertia calculation
170 <    Mat3x3d Itmp(0.0);
168 <  
170 >    Mat3x3d Itmp(0.0);    
171      for (std::size_t i = 0; i < atoms_.size(); i++) {
172 +      Mat3x3d IAtom(0.0);  
173        mtmp = atoms_[i]->getMass();
174 <      Itmp -= outProduct(refCoords_[i], refCoords_[i]) * mtmp;
175 <      double r2 = refCoords_[i].lengthSquare();
176 <      Itmp(0, 0) += mtmp * r2;
177 <      Itmp(1, 1) += mtmp * r2;
178 <      Itmp(2, 2) += mtmp * r2;
179 <    }
180 <
181 <    //project the inertial moment of directional atoms into this rigid body
179 <    for (std::size_t i = 0; i < atoms_.size(); i++) {
174 >      IAtom -= outProduct(refCoords_[i], refCoords_[i]) * mtmp;
175 >      RealType r2 = refCoords_[i].lengthSquare();
176 >      IAtom(0, 0) += mtmp * r2;
177 >      IAtom(1, 1) += mtmp * r2;
178 >      IAtom(2, 2) += mtmp * r2;
179 >      Itmp += IAtom;
180 >      
181 >      //project the inertial moment of directional atoms into this rigid body
182        if (atoms_[i]->isDirectional()) {
183 <        RectMatrix<double, 3, 3> Iproject = refOrients_[i].transpose() * atoms_[i]->getI();
184 <        Itmp(0, 0) += Iproject(0, 0);
183 <        Itmp(1, 1) += Iproject(1, 1);
184 <        Itmp(2, 2) += Iproject(2, 2);
185 <      }
183 >        Itmp += refOrients_[i].transpose() * atoms_[i]->getI() * refOrients_[i];
184 >      }
185      }
186  
187 +    //    std::cout << Itmp << std::endl;
188 +
189      //diagonalize
190      Vector3d evals;
191      Mat3x3d::diagonalize(Itmp, evals, sU_);
# Line 196 | Line 197 | namespace oopse {
197          
198      int nLinearAxis = 0;
199      for (int i = 0; i < 3; i++) {    
200 <      if (fabs(evals[i]) < oopse::epsilon) {
200 >      if (fabs(evals[i]) < OpenMD::epsilon) {
201          linear_ = true;
202          linearAxis_ = i;
203          ++ nLinearAxis;
# Line 206 | Line 207 | namespace oopse {
207      if (nLinearAxis > 1) {
208        sprintf( painCave.errMsg,
209                 "RigidBody error.\n"
210 <               "\tOOPSE found more than one axis in this rigid body with a vanishing \n"
210 >               "\tOpenMD found more than one axis in this rigid body with a vanishing \n"
211                 "\tmoment of inertia.  This can happen in one of three ways:\n"
212                 "\t 1) Only one atom was specified, or \n"
213                 "\t 2) All atoms were specified at the same location, or\n"
# Line 225 | Line 226 | namespace oopse {
226      Vector3d apos;
227      Vector3d rpos;
228      Vector3d frc(0.0);
229 <    Vector3d trq(0.0);
229 >    Vector3d trq(0.0);    
230      Vector3d pos = this->getPos();
231      for (int i = 0; i < atoms_.size(); i++) {
232  
# Line 245 | Line 246 | namespace oopse {
246        if (atoms_[i]->isDirectional()) {
247          atrq = atoms_[i]->getTrq();
248          trq += atrq;
249 <      }
249 >      }      
250 >    }        
251 >    addFrc(frc);
252 >    addTrq(trq);    
253 >  }
254 >
255 >  Mat3x3d RigidBody::calcForcesAndTorquesAndVirial() {
256 >    Vector3d afrc;
257 >    Vector3d atrq;
258 >    Vector3d apos;
259 >    Vector3d rpos;
260 >    Vector3d dfrc;
261 >    Vector3d frc(0.0);
262 >    Vector3d trq(0.0);    
263 >    Vector3d pos = this->getPos();
264 >    Mat3x3d tau_(0.0);
265 >
266 >    for (int i = 0; i < atoms_.size(); i++) {
267 >      
268 >      afrc = atoms_[i]->getFrc();
269 >      apos = atoms_[i]->getPos();
270 >      rpos = apos - pos;
271          
272 +      frc += afrc;
273 +
274 +      trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1];
275 +      trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2];
276 +      trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0];
277 +
278 +      // If the atom has a torque associated with it, then we also need to
279 +      // migrate the torques onto the center of mass:
280 +
281 +      if (atoms_[i]->isDirectional()) {
282 +        atrq = atoms_[i]->getTrq();
283 +        trq += atrq;
284 +      }
285 +      
286 +      tau_(0,0) -= rpos[0]*afrc[0];
287 +      tau_(0,1) -= rpos[0]*afrc[1];
288 +      tau_(0,2) -= rpos[0]*afrc[2];
289 +      tau_(1,0) -= rpos[1]*afrc[0];
290 +      tau_(1,1) -= rpos[1]*afrc[1];
291 +      tau_(1,2) -= rpos[1]*afrc[2];
292 +      tau_(2,0) -= rpos[2]*afrc[0];
293 +      tau_(2,1) -= rpos[2]*afrc[1];
294 +      tau_(2,2) -= rpos[2]*afrc[2];
295 +
296      }
297 <    
298 <    setFrc(frc);
299 <    setTrq(trq);
254 <    
297 >    addFrc(frc);
298 >    addTrq(trq);
299 >    return tau_;
300    }
301  
302    void  RigidBody::updateAtoms() {
# Line 273 | Line 318 | namespace oopse {
318        if (atoms_[i]->isDirectional()) {
319            
320          dAtom = (DirectionalAtom *) atoms_[i];
321 <        dAtom->setA(a * refOrients_[i]);
277 <        //dAtom->rotateBy( A );      
321 >        dAtom->setA(refOrients_[i].transpose() * a);
322        }
323  
324      }
# Line 301 | Line 345 | namespace oopse {
345        if (atoms_[i]->isDirectional()) {
346            
347          dAtom = (DirectionalAtom *) atoms_[i];
348 <        dAtom->setA(a * refOrients_[i], frame);
348 >        dAtom->setA(refOrients_[i].transpose() * a, frame);
349        }
350  
351      }
# Line 486 | Line 530 | namespace oopse {
530                 "RigidBody error.\n"
531                 "\tAtom %s does not have a position specified.\n"
532                 "\tThis means RigidBody cannot set up reference coordinates.\n",
533 <               ats->getType() );
533 >               ats->getType().c_str() );
534        painCave.isFatal = 1;
535        simError();
536      }
# Line 506 | Line 550 | namespace oopse {
550                   "RigidBody error.\n"
551                   "\tAtom %s does not have an orientation specified.\n"
552                   "\tThis means RigidBody cannot set up reference orientations.\n",
553 <                 ats->getType() );
553 >                 ats->getType().c_str() );
554          painCave.isFatal = 1;
555          simError();
556        }    

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines