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 253 by tim, Thu Jan 13 19:40:37 2005 UTC vs.
Revision 1424 by gezelter, Tue Mar 30 15:05:38 2010 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# 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 < namespace oopse {
46 <
47 < RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData), inertiaTensor_(0.0){
48 <
49 < }
50 <
51 < void RigidBody::setPrevA(const RotMat3x3d& a) {
45 > #include "utils/NumericConstant.hpp"
46 > namespace OpenMD {
47 >  
48 >  RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData),
49 >                           inertiaTensor_(0.0){    
50 >  }
51 >  
52 >  void RigidBody::setPrevA(const RotMat3x3d& a) {
53      ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
54 <    //((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
54 <
54 >    
55      for (int i =0 ; i < atoms_.size(); ++i){
56 <        if (atoms_[i]->isDirectional()) {
57 <            atoms_[i]->setPrevA(a * refOrients_[i]);
58 <        }
56 >      if (atoms_[i]->isDirectional()) {
57 >        atoms_[i]->setPrevA(refOrients_[i].transpose() * a);
58 >      }
59      }
60 <
61 < }
62 <
63 <      
64 < void RigidBody::setA(const RotMat3x3d& a) {
60 >    
61 >  }
62 >  
63 >  
64 >  void RigidBody::setA(const RotMat3x3d& a) {
65      ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
66 –    //((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]);
70 <        }
68 >      if (atoms_[i]->isDirectional()) {
69 >        atoms_[i]->setA(refOrients_[i].transpose() * a);
70 >      }
71      }
72 < }    
73 <    
74 < void RigidBody::setA(const RotMat3x3d& a, int snapshotNo) {
72 >  }    
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);
82 <        }
80 >      if (atoms_[i]->isDirectional()) {
81 >        atoms_[i]->setA(refOrients_[i].transpose() * a, snapshotNo);
82 >      }
83      }
84 <
85 < }  
86 <
87 < Mat3x3d RigidBody::getI() {
84 >    
85 >  }  
86 >  
87 >  Mat3x3d RigidBody::getI() {
88      return inertiaTensor_;
89 < }    
90 <
91 < std::vector<double> RigidBody::getGrad() {
92 <     std::vector<double> grad(6, 0.0);
89 >  }    
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 <
135 >      grad[j] = -force[j];
136 >    
137      for (int j = 0; j < 3; j++ ) {
138 <
139 <        grad[3] += torque[j]*ephi[j];
140 <        grad[4] += torque[j]*etheta[j];
141 <        grad[5] += torque[j]*epsi[j];
142 <
138 >      
139 >      grad[3] += torque[j]*ephi[j];
140 >      grad[4] += torque[j]*etheta[j];
141 >      grad[5] += torque[j]*epsi[j];
142 >      
143      }
144      
145      return grad;
146 < }    
147 <
148 < void RigidBody::accept(BaseVisitor* v) {
146 >  }    
147 >  
148 >  void RigidBody::accept(BaseVisitor* v) {
149      v->visit(this);
150 < }    
150 >  }    
151  
152 < /**@todo need modification */
153 < void  RigidBody::calcRefCoords() {
154 <    double mtmp;
152 >  /**@todo need modification */
153 >  void  RigidBody::calcRefCoords() {
154 >    RealType mtmp;
155      Vector3d refCOM(0.0);
156      mass_ = 0.0;
157      for (std::size_t i = 0; i < atoms_.size(); ++i) {
158 <        mtmp = atoms_[i]->getMass();
159 <        mass_ += mtmp;
160 <        refCOM += refCoords_[i]*mtmp;
158 >      mtmp = atoms_[i]->getMass();
159 >      mass_ += mtmp;
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;
166 >      refCoords_[i] -= refCOM;
167      }
168  
169 < // Moment of Inertia calculation
170 <    Mat3x3d Itmp(0.0);
167 <  
169 >    // Moment of Inertia calculation
170 >    Mat3x3d Itmp(0.0);    
171      for (std::size_t i = 0; i < atoms_.size(); i++) {
172 <        mtmp = atoms_[i]->getMass();
173 <        Itmp -= outProduct(refCoords_[i], refCoords_[i]) * mtmp;
174 <        double r2 = refCoords_[i].lengthSquare();
175 <        Itmp(0, 0) += mtmp * r2;
176 <        Itmp(1, 1) += mtmp * r2;
177 <        Itmp(2, 2) += mtmp * r2;
172 >      Mat3x3d IAtom(0.0);  
173 >      mtmp = atoms_[i]->getMass();
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 >        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 185 | Line 197 | void  RigidBody::calcRefCoords() {
197          
198      int nLinearAxis = 0;
199      for (int i = 0; i < 3; i++) {    
200 <        if (fabs(evals[i]) < oopse::epsilon) {
201 <            linear_ = true;
202 <            linearAxis_ = i;
203 <            ++ nLinearAxis;
204 <        }
200 >      if (fabs(evals[i]) < OpenMD::epsilon) {
201 >        linear_ = true;
202 >        linearAxis_ = i;
203 >        ++ nLinearAxis;
204 >      }
205      }
206  
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"
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"
214 <            "\t 3) The programmers did something stupid.\n"
215 <            "\tIt is silly to use a rigid body to describe this situation.  Be smarter.\n"
216 <            );
217 <        painCave.isFatal = 1;
218 <        simError();
208 >      sprintf( painCave.errMsg,
209 >               "RigidBody error.\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"
214 >               "\t 3) The programmers did something stupid.\n"
215 >               "\tIt is silly to use a rigid body to describe this situation.  Be smarter.\n"
216 >               );
217 >      painCave.isFatal = 1;
218 >      simError();
219      }
220    
221 < }
221 >  }
222  
223 < void  RigidBody::calcForcesAndTorques() {
223 >  void  RigidBody::calcForcesAndTorques() {
224      Vector3d afrc;
225      Vector3d atrq;
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  
233 <        afrc = atoms_[i]->getFrc();
234 <        apos = atoms_[i]->getPos();
235 <        rpos = apos - pos;
233 >      afrc = atoms_[i]->getFrc();
234 >      apos = atoms_[i]->getPos();
235 >      rpos = apos - pos;
236          
237 <        frc += afrc;
237 >      frc += afrc;
238  
239 <        trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1];
240 <        trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2];
241 <        trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0];
239 >      trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1];
240 >      trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2];
241 >      trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0];
242  
243 <        // If the atom has a torque associated with it, then we also need to
244 <        // migrate the torques onto the center of mass:
243 >      // If the atom has a torque associated with it, then we also need to
244 >      // migrate the torques onto the center of mass:
245  
246 <        if (atoms_[i]->isDirectional()) {
247 <            atrq = atoms_[i]->getTrq();
248 <            trq += atrq;
249 <        }
246 >      if (atoms_[i]->isDirectional()) {
247 >        atrq = atoms_[i]->getTrq();
248 >        trq += atrq;
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);
300 <    
244 < }
297 >    addFrc(frc);
298 >    addTrq(trq);
299 >    return tau_;
300 >  }
301  
302 < void  RigidBody::updateAtoms() {
302 >  void  RigidBody::updateAtoms() {
303      unsigned int i;
304      Vector3d ref;
305      Vector3d apos;
# Line 253 | Line 309 | void  RigidBody::updateAtoms() {
309      
310      for (i = 0; i < atoms_.size(); i++) {
311      
312 <        ref = body2Lab(refCoords_[i]);
312 >      ref = body2Lab(refCoords_[i]);
313  
314 <        apos = pos + ref;
314 >      apos = pos + ref;
315  
316 <        atoms_[i]->setPos(apos);
316 >      atoms_[i]->setPos(apos);
317  
318 <        if (atoms_[i]->isDirectional()) {
318 >      if (atoms_[i]->isDirectional()) {
319            
320 <          dAtom = (DirectionalAtom *) atoms_[i];
321 <          dAtom->setA(a * refOrients_[i]);
322 <          //dAtom->rotateBy( A );      
267 <        }
320 >        dAtom = (DirectionalAtom *) atoms_[i];
321 >        dAtom->setA(refOrients_[i].transpose() * a);
322 >      }
323  
324      }
325    
326 < }
326 >  }
327  
328  
329 < bool RigidBody::getAtomPos(Vector3d& pos, unsigned int index) {
329 >  void  RigidBody::updateAtoms(int frame) {
330 >    unsigned int i;
331 >    Vector3d ref;
332 >    Vector3d apos;
333 >    DirectionalAtom* dAtom;
334 >    Vector3d pos = getPos(frame);
335 >    RotMat3x3d a = getA(frame);
336 >    
337 >    for (i = 0; i < atoms_.size(); i++) {
338 >    
339 >      ref = body2Lab(refCoords_[i], frame);
340 >
341 >      apos = pos + ref;
342 >
343 >      atoms_[i]->setPos(apos, frame);
344 >
345 >      if (atoms_[i]->isDirectional()) {
346 >          
347 >        dAtom = (DirectionalAtom *) atoms_[i];
348 >        dAtom->setA(refOrients_[i].transpose() * a, frame);
349 >      }
350 >
351 >    }
352 >  
353 >  }
354 >
355 >  void RigidBody::updateAtomVel() {
356 >    Mat3x3d skewMat;;
357 >
358 >    Vector3d ji = getJ();
359 >    Mat3x3d I =  getI();
360 >
361 >    skewMat(0, 0) =0;
362 >    skewMat(0, 1) = ji[2] /I(2, 2);
363 >    skewMat(0, 2) = -ji[1] /I(1, 1);
364 >
365 >    skewMat(1, 0) = -ji[2] /I(2, 2);
366 >    skewMat(1, 1) = 0;
367 >    skewMat(1, 2) = ji[0]/I(0, 0);
368 >
369 >    skewMat(2, 0) =ji[1] /I(1, 1);
370 >    skewMat(2, 1) = -ji[0]/I(0, 0);
371 >    skewMat(2, 2) = 0;
372 >
373 >    Mat3x3d mat = (getA() * skewMat).transpose();
374 >    Vector3d rbVel = getVel();
375 >
376 >
377 >    Vector3d velRot;        
378 >    for (int i =0 ; i < refCoords_.size(); ++i) {
379 >      atoms_[i]->setVel(rbVel + mat * refCoords_[i]);
380 >    }
381 >
382 >  }
383 >
384 >  void RigidBody::updateAtomVel(int frame) {
385 >    Mat3x3d skewMat;;
386 >
387 >    Vector3d ji = getJ(frame);
388 >    Mat3x3d I =  getI();
389 >
390 >    skewMat(0, 0) =0;
391 >    skewMat(0, 1) = ji[2] /I(2, 2);
392 >    skewMat(0, 2) = -ji[1] /I(1, 1);
393 >
394 >    skewMat(1, 0) = -ji[2] /I(2, 2);
395 >    skewMat(1, 1) = 0;
396 >    skewMat(1, 2) = ji[0]/I(0, 0);
397 >
398 >    skewMat(2, 0) =ji[1] /I(1, 1);
399 >    skewMat(2, 1) = -ji[0]/I(0, 0);
400 >    skewMat(2, 2) = 0;
401 >
402 >    Mat3x3d mat = (getA(frame) * skewMat).transpose();
403 >    Vector3d rbVel = getVel(frame);
404 >
405 >
406 >    Vector3d velRot;        
407 >    for (int i =0 ; i < refCoords_.size(); ++i) {
408 >      atoms_[i]->setVel(rbVel + mat * refCoords_[i], frame);
409 >    }
410 >
411 >  }
412 >
413 >        
414 >
415 >  bool RigidBody::getAtomPos(Vector3d& pos, unsigned int index) {
416      if (index < atoms_.size()) {
417  
418 <        Vector3d ref = body2Lab(refCoords_[index]);
419 <        pos = getPos() + ref;
420 <        return true;
418 >      Vector3d ref = body2Lab(refCoords_[index]);
419 >      pos = getPos() + ref;
420 >      return true;
421      } else {
422 <        std::cerr << index << " is an invalid index, current rigid body contains "
423 <                      << atoms_.size() << "atoms" << std::endl;
424 <        return false;
422 >      std::cerr << index << " is an invalid index, current rigid body contains "
423 >                << atoms_.size() << "atoms" << std::endl;
424 >      return false;
425      }    
426 < }
426 >  }
427  
428 < bool RigidBody::getAtomPos(Vector3d& pos, Atom* atom) {
428 >  bool RigidBody::getAtomPos(Vector3d& pos, Atom* atom) {
429      std::vector<Atom*>::iterator i;
430      i = std::find(atoms_.begin(), atoms_.end(), atom);
431      if (i != atoms_.end()) {
432 <        //RigidBody class makes sure refCoords_ and atoms_ match each other
433 <        Vector3d ref = body2Lab(refCoords_[i - atoms_.begin()]);
434 <        pos = getPos() + ref;
435 <        return true;
432 >      //RigidBody class makes sure refCoords_ and atoms_ match each other
433 >      Vector3d ref = body2Lab(refCoords_[i - atoms_.begin()]);
434 >      pos = getPos() + ref;
435 >      return true;
436      } else {
437 <        std::cerr << "Atom " << atom->getGlobalIndex()
438 <                      <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;
439 <        return false;
437 >      std::cerr << "Atom " << atom->getGlobalIndex()
438 >                <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;
439 >      return false;
440      }
441 < }
442 < bool RigidBody::getAtomVel(Vector3d& vel, unsigned int index) {
441 >  }
442 >  bool RigidBody::getAtomVel(Vector3d& vel, unsigned int index) {
443  
444      //velRot = $(A\cdot skew(I^{-1}j))^{T}refCoor$
445  
446      if (index < atoms_.size()) {
447  
448 <        Vector3d velRot;
449 <        Mat3x3d skewMat;;
450 <        Vector3d ref = refCoords_[index];
451 <        Vector3d ji = getJ();
452 <        Mat3x3d I =  getI();
448 >      Vector3d velRot;
449 >      Mat3x3d skewMat;;
450 >      Vector3d ref = refCoords_[index];
451 >      Vector3d ji = getJ();
452 >      Mat3x3d I =  getI();
453  
454 <        skewMat(0, 0) =0;
455 <        skewMat(0, 1) = ji[2] /I(2, 2);
456 <        skewMat(0, 2) = -ji[1] /I(1, 1);
454 >      skewMat(0, 0) =0;
455 >      skewMat(0, 1) = ji[2] /I(2, 2);
456 >      skewMat(0, 2) = -ji[1] /I(1, 1);
457  
458 <        skewMat(1, 0) = -ji[2] /I(2, 2);
459 <        skewMat(1, 1) = 0;
460 <        skewMat(1, 2) = ji[0]/I(0, 0);
458 >      skewMat(1, 0) = -ji[2] /I(2, 2);
459 >      skewMat(1, 1) = 0;
460 >      skewMat(1, 2) = ji[0]/I(0, 0);
461  
462 <        skewMat(2, 0) =ji[1] /I(1, 1);
463 <        skewMat(2, 1) = -ji[0]/I(0, 0);
464 <        skewMat(2, 2) = 0;
462 >      skewMat(2, 0) =ji[1] /I(1, 1);
463 >      skewMat(2, 1) = -ji[0]/I(0, 0);
464 >      skewMat(2, 2) = 0;
465  
466 <        velRot = (getA() * skewMat).transpose() * ref;
466 >      velRot = (getA() * skewMat).transpose() * ref;
467  
468 <        vel =getVel() + velRot;
469 <        return true;
468 >      vel =getVel() + velRot;
469 >      return true;
470          
471      } else {
472 <        std::cerr << index << " is an invalid index, current rigid body contains "
473 <                      << atoms_.size() << "atoms" << std::endl;
474 <        return false;
472 >      std::cerr << index << " is an invalid index, current rigid body contains "
473 >                << atoms_.size() << "atoms" << std::endl;
474 >      return false;
475      }
476 < }
476 >  }
477  
478 < bool RigidBody::getAtomVel(Vector3d& vel, Atom* atom) {
478 >  bool RigidBody::getAtomVel(Vector3d& vel, Atom* atom) {
479  
480      std::vector<Atom*>::iterator i;
481      i = std::find(atoms_.begin(), atoms_.end(), atom);
482      if (i != atoms_.end()) {
483 <        return getAtomVel(vel, i - atoms_.begin());
483 >      return getAtomVel(vel, i - atoms_.begin());
484      } else {
485 <        std::cerr << "Atom " << atom->getGlobalIndex()
486 <                      <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;    
487 <        return false;
485 >      std::cerr << "Atom " << atom->getGlobalIndex()
486 >                <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;    
487 >      return false;
488      }    
489 < }
489 >  }
490  
491 < bool RigidBody::getAtomRefCoor(Vector3d& coor, unsigned int index) {
491 >  bool RigidBody::getAtomRefCoor(Vector3d& coor, unsigned int index) {
492      if (index < atoms_.size()) {
493  
494 <        coor = refCoords_[index];
495 <        return true;
494 >      coor = refCoords_[index];
495 >      return true;
496      } else {
497 <        std::cerr << index << " is an invalid index, current rigid body contains "
498 <                      << atoms_.size() << "atoms" << std::endl;
499 <        return false;
497 >      std::cerr << index << " is an invalid index, current rigid body contains "
498 >                << atoms_.size() << "atoms" << std::endl;
499 >      return false;
500      }
501  
502 < }
502 >  }
503  
504 < bool RigidBody::getAtomRefCoor(Vector3d& coor, Atom* atom) {
504 >  bool RigidBody::getAtomRefCoor(Vector3d& coor, Atom* atom) {
505      std::vector<Atom*>::iterator i;
506      i = std::find(atoms_.begin(), atoms_.end(), atom);
507      if (i != atoms_.end()) {
508 <        //RigidBody class makes sure refCoords_ and atoms_ match each other
509 <        coor = refCoords_[i - atoms_.begin()];
510 <        return true;
508 >      //RigidBody class makes sure refCoords_ and atoms_ match each other
509 >      coor = refCoords_[i - atoms_.begin()];
510 >      return true;
511      } else {
512 <        std::cerr << "Atom " << atom->getGlobalIndex()
513 <                      <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;    
514 <        return false;
512 >      std::cerr << "Atom " << atom->getGlobalIndex()
513 >                <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;    
514 >      return false;
515      }
516  
517 < }
517 >  }
518  
519  
520 < void RigidBody::addAtom(Atom* at, AtomStamp* ats) {
520 >  void RigidBody::addAtom(Atom* at, AtomStamp* ats) {
521  
522 <  Vector3d coords;
523 <  Vector3d euler;
522 >    Vector3d coords;
523 >    Vector3d euler;
524    
525  
526 <  atoms_.push_back(at);
526 >    atoms_.push_back(at);
527  
528 <  if( !ats->havePosition() ){
529 <    sprintf( painCave.errMsg,
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() );
534 <    painCave.isFatal = 1;
535 <    simError();
536 <  }
528 >    if( !ats->havePosition() ){
529 >      sprintf( painCave.errMsg,
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().c_str() );
534 >      painCave.isFatal = 1;
535 >      simError();
536 >    }
537    
538 <  coords[0] = ats->getPosX();
539 <  coords[1] = ats->getPosY();
540 <  coords[2] = ats->getPosZ();
538 >    coords[0] = ats->getPosX();
539 >    coords[1] = ats->getPosY();
540 >    coords[2] = ats->getPosZ();
541  
542 <  refCoords_.push_back(coords);
542 >    refCoords_.push_back(coords);
543  
544 <  RotMat3x3d identMat = RotMat3x3d::identity();
544 >    RotMat3x3d identMat = RotMat3x3d::identity();
545    
546 <  if (at->isDirectional()) {  
546 >    if (at->isDirectional()) {  
547  
548 <    if( !ats->haveOrientation() ){
549 <      sprintf( painCave.errMsg,
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() );
554 <      painCave.isFatal = 1;
555 <      simError();
556 <    }    
548 >      if( !ats->haveOrientation() ){
549 >        sprintf( painCave.errMsg,
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().c_str() );
554 >        painCave.isFatal = 1;
555 >        simError();
556 >      }    
557      
558 <    euler[0] = ats->getEulerPhi();
559 <    euler[1] = ats->getEulerTheta();
560 <    euler[2] = ats->getEulerPsi();
558 >      euler[0] = ats->getEulerPhi() * NumericConstant::PI /180.0;
559 >      euler[1] = ats->getEulerTheta() * NumericConstant::PI /180.0;
560 >      euler[2] = ats->getEulerPsi() * NumericConstant::PI /180.0;
561  
562 <    RotMat3x3d Atmp(euler);
563 <    refOrients_.push_back(Atmp);
562 >      RotMat3x3d Atmp(euler);
563 >      refOrients_.push_back(Atmp);
564      
565 <  }else {
566 <    refOrients_.push_back(identMat);
567 <  }
565 >    }else {
566 >      refOrients_.push_back(identMat);
567 >    }
568    
569    
570 < }
570 >  }
571  
572   }
573  

Diff Legend

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