ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/primitives/RigidBody.cpp
Revision: 2058
Committed: Tue Feb 22 18:56:25 2005 UTC (19 years, 4 months ago) by tim
File size: 14826 byte(s)
Log Message:
reactionfield get fixed

File Contents

# Content
1 /*
2 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 *
4 * The University of Notre Dame grants you ("Licensee") a
5 * non-exclusive, royalty free, license to use, modify and
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
19 * notice, this list of conditions and the following disclaimer.
20 *
21 * 3. Redistributions in binary form must reproduce the above copyright
22 * notice, this list of conditions and the following disclaimer in the
23 * documentation and/or other materials provided with the
24 * distribution.
25 *
26 * This software is provided "AS IS," without a warranty of any
27 * kind. All express or implied conditions, representations and
28 * warranties, including any implied warranty of merchantability,
29 * fitness for a particular purpose or non-infringement, are hereby
30 * excluded. The University of Notre Dame and its licensors shall not
31 * be liable for any damages suffered by licensee as a result of
32 * using, modifying or distributing the software or its
33 * derivatives. In no event will the University of Notre Dame or its
34 * licensors be liable for any lost revenue, profit or data, or for
35 * direct, indirect, special, consequential, incidental or punitive
36 * damages, however caused and regardless of the theory of liability,
37 * arising out of the use of or inability to use software, even if the
38 * University of Notre Dame has been advised of the possibility of
39 * such damages.
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
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_;
55
56 for (int i =0 ; i < atoms_.size(); ++i){
57 if (atoms_[i]->isDirectional()) {
58 atoms_[i]->setPrevA(a * refOrients_[i]);
59 }
60 }
61
62 }
63
64
65 void RigidBody::setA(const RotMat3x3d& a) {
66 ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
67 //((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
68
69 for (int i =0 ; i < atoms_.size(); ++i){
70 if (atoms_[i]->isDirectional()) {
71 atoms_[i]->setA(a * refOrients_[i]);
72 }
73 }
74 }
75
76 void RigidBody::setA(const RotMat3x3d& a, int snapshotNo) {
77 ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
78 //((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
79
80 for (int i =0 ; i < atoms_.size(); ++i){
81 if (atoms_[i]->isDirectional()) {
82 atoms_[i]->setA(a * refOrients_[i], snapshotNo);
83 }
84 }
85
86 }
87
88 Mat3x3d RigidBody::getI() {
89 return inertiaTensor_;
90 }
91
92 std::vector<double> RigidBody::getGrad() {
93 std::vector<double> grad(6, 0.0);
94 Vector3d force;
95 Vector3d torque;
96 Vector3d myEuler;
97 double phi, theta, psi;
98 double cphi, sphi, ctheta, stheta;
99 Vector3d ephi;
100 Vector3d etheta;
101 Vector3d epsi;
102
103 force = getFrc();
104 torque =getTrq();
105 myEuler = getA().toEulerAngles();
106
107 phi = myEuler[0];
108 theta = myEuler[1];
109 psi = myEuler[2];
110
111 cphi = cos(phi);
112 sphi = sin(phi);
113 ctheta = cos(theta);
114 stheta = sin(theta);
115
116 // get unit vectors along the phi, theta and psi rotation axes
117
118 ephi[0] = 0.0;
119 ephi[1] = 0.0;
120 ephi[2] = 1.0;
121
122 etheta[0] = cphi;
123 etheta[1] = sphi;
124 etheta[2] = 0.0;
125
126 epsi[0] = stheta * cphi;
127 epsi[1] = stheta * sphi;
128 epsi[2] = ctheta;
129
130 //gradient is equal to -force
131 for (int j = 0 ; j<3; j++)
132 grad[j] = -force[j];
133
134 for (int j = 0; j < 3; j++ ) {
135
136 grad[3] += torque[j]*ephi[j];
137 grad[4] += torque[j]*etheta[j];
138 grad[5] += torque[j]*epsi[j];
139
140 }
141
142 return grad;
143 }
144
145 void RigidBody::accept(BaseVisitor* v) {
146 v->visit(this);
147 }
148
149 /**@todo need modification */
150 void RigidBody::calcRefCoords() {
151 double mtmp;
152 Vector3d refCOM(0.0);
153 mass_ = 0.0;
154 for (std::size_t i = 0; i < atoms_.size(); ++i) {
155 mtmp = atoms_[i]->getMass();
156 mass_ += mtmp;
157 refCOM += refCoords_[i]*mtmp;
158 }
159 refCOM /= mass_;
160
161 // Next, move the origin of the reference coordinate system to the COM:
162 for (std::size_t i = 0; i < atoms_.size(); ++i) {
163 refCoords_[i] -= refCOM;
164 }
165
166 // Moment of Inertia calculation
167 Mat3x3d Itmp(0.0);
168
169 for (std::size_t i = 0; i < atoms_.size(); i++) {
170 mtmp = atoms_[i]->getMass();
171 Itmp -= outProduct(refCoords_[i], refCoords_[i]) * mtmp;
172 double r2 = refCoords_[i].lengthSquare();
173 Itmp(0, 0) += mtmp * r2;
174 Itmp(1, 1) += mtmp * r2;
175 Itmp(2, 2) += mtmp * r2;
176 }
177
178 //project the inertial moment of directional atoms into this rigid body
179 for (std::size_t i = 0; i < atoms_.size(); i++) {
180 if (atoms_[i]->isDirectional()) {
181 RectMatrix<double, 3, 3> Iproject = refOrients_[i].transpose() * atoms_[i]->getI();
182 Itmp(0, 0) += Iproject(0, 0);
183 Itmp(1, 1) += Iproject(1, 1);
184 Itmp(2, 2) += Iproject(2, 2);
185 }
186 }
187
188 //diagonalize
189 Vector3d evals;
190 Mat3x3d::diagonalize(Itmp, evals, sU_);
191
192 // zero out I and then fill the diagonals with the moments of inertia:
193 inertiaTensor_(0, 0) = evals[0];
194 inertiaTensor_(1, 1) = evals[1];
195 inertiaTensor_(2, 2) = evals[2];
196
197 int nLinearAxis = 0;
198 for (int i = 0; i < 3; i++) {
199 if (fabs(evals[i]) < oopse::epsilon) {
200 linear_ = true;
201 linearAxis_ = i;
202 ++ nLinearAxis;
203 }
204 }
205
206 if (nLinearAxis > 1) {
207 sprintf( painCave.errMsg,
208 "RigidBody error.\n"
209 "\tOOPSE found more than one axis in this rigid body with a vanishing \n"
210 "\tmoment of inertia. This can happen in one of three ways:\n"
211 "\t 1) Only one atom was specified, or \n"
212 "\t 2) All atoms were specified at the same location, or\n"
213 "\t 3) The programmers did something stupid.\n"
214 "\tIt is silly to use a rigid body to describe this situation. Be smarter.\n"
215 );
216 painCave.isFatal = 1;
217 simError();
218 }
219
220 }
221
222 void RigidBody::calcForcesAndTorques() {
223 Vector3d afrc;
224 Vector3d atrq;
225 Vector3d apos;
226 Vector3d rpos;
227 Vector3d frc(0.0);
228 Vector3d trq(0.0);
229 Vector3d pos = this->getPos();
230 for (int i = 0; i < atoms_.size(); i++) {
231
232 afrc = atoms_[i]->getFrc();
233 apos = atoms_[i]->getPos();
234 rpos = apos - pos;
235
236 frc += afrc;
237
238 trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1];
239 trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2];
240 trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0];
241
242 // If the atom has a torque associated with it, then we also need to
243 // migrate the torques onto the center of mass:
244
245 if (atoms_[i]->isDirectional()) {
246 atrq = atoms_[i]->getTrq();
247 trq += atrq;
248 }
249
250 }
251
252 setFrc(frc);
253 setTrq(trq);
254
255 }
256
257 void RigidBody::updateAtoms() {
258 unsigned int i;
259 Vector3d ref;
260 Vector3d apos;
261 DirectionalAtom* dAtom;
262 Vector3d pos = getPos();
263 RotMat3x3d a = getA();
264
265 for (i = 0; i < atoms_.size(); i++) {
266
267 ref = body2Lab(refCoords_[i]);
268
269 apos = pos + ref;
270
271 atoms_[i]->setPos(apos);
272
273 if (atoms_[i]->isDirectional()) {
274
275 dAtom = (DirectionalAtom *) atoms_[i];
276 dAtom->setA(a * refOrients_[i]);
277 //dAtom->rotateBy( A );
278 }
279
280 }
281
282 }
283
284
285 void RigidBody::updateAtoms(int frame) {
286 unsigned int i;
287 Vector3d ref;
288 Vector3d apos;
289 DirectionalAtom* dAtom;
290 Vector3d pos = getPos(frame);
291 RotMat3x3d a = getA(frame);
292
293 for (i = 0; i < atoms_.size(); i++) {
294
295 ref = body2Lab(refCoords_[i], frame);
296
297 apos = pos + ref;
298
299 atoms_[i]->setPos(apos, frame);
300
301 if (atoms_[i]->isDirectional()) {
302
303 dAtom = (DirectionalAtom *) atoms_[i];
304 dAtom->setA(a * refOrients_[i], frame);
305 }
306
307 }
308
309 }
310
311 void RigidBody::updateAtomVel() {
312 Mat3x3d skewMat;;
313
314 Vector3d ji = getJ();
315 Mat3x3d I = getI();
316
317 skewMat(0, 0) =0;
318 skewMat(0, 1) = ji[2] /I(2, 2);
319 skewMat(0, 2) = -ji[1] /I(1, 1);
320
321 skewMat(1, 0) = -ji[2] /I(2, 2);
322 skewMat(1, 1) = 0;
323 skewMat(1, 2) = ji[0]/I(0, 0);
324
325 skewMat(2, 0) =ji[1] /I(1, 1);
326 skewMat(2, 1) = -ji[0]/I(0, 0);
327 skewMat(2, 2) = 0;
328
329 Mat3x3d mat = (getA() * skewMat).transpose();
330 Vector3d rbVel = getVel();
331
332
333 Vector3d velRot;
334 for (int i =0 ; i < refCoords_.size(); ++i) {
335 atoms_[i]->setVel(rbVel + mat * refCoords_[i]);
336 }
337
338 }
339
340 void RigidBody::updateAtomVel(int frame) {
341 Mat3x3d skewMat;;
342
343 Vector3d ji = getJ(frame);
344 Mat3x3d I = getI();
345
346 skewMat(0, 0) =0;
347 skewMat(0, 1) = ji[2] /I(2, 2);
348 skewMat(0, 2) = -ji[1] /I(1, 1);
349
350 skewMat(1, 0) = -ji[2] /I(2, 2);
351 skewMat(1, 1) = 0;
352 skewMat(1, 2) = ji[0]/I(0, 0);
353
354 skewMat(2, 0) =ji[1] /I(1, 1);
355 skewMat(2, 1) = -ji[0]/I(0, 0);
356 skewMat(2, 2) = 0;
357
358 Mat3x3d mat = (getA(frame) * skewMat).transpose();
359 Vector3d rbVel = getVel(frame);
360
361
362 Vector3d velRot;
363 for (int i =0 ; i < refCoords_.size(); ++i) {
364 atoms_[i]->setVel(rbVel + mat * refCoords_[i], frame);
365 }
366
367 }
368
369
370
371 bool RigidBody::getAtomPos(Vector3d& pos, unsigned int index) {
372 if (index < atoms_.size()) {
373
374 Vector3d ref = body2Lab(refCoords_[index]);
375 pos = getPos() + ref;
376 return true;
377 } else {
378 std::cerr << index << " is an invalid index, current rigid body contains "
379 << atoms_.size() << "atoms" << std::endl;
380 return false;
381 }
382 }
383
384 bool RigidBody::getAtomPos(Vector3d& pos, Atom* atom) {
385 std::vector<Atom*>::iterator i;
386 i = std::find(atoms_.begin(), atoms_.end(), atom);
387 if (i != atoms_.end()) {
388 //RigidBody class makes sure refCoords_ and atoms_ match each other
389 Vector3d ref = body2Lab(refCoords_[i - atoms_.begin()]);
390 pos = getPos() + ref;
391 return true;
392 } else {
393 std::cerr << "Atom " << atom->getGlobalIndex()
394 <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;
395 return false;
396 }
397 }
398 bool RigidBody::getAtomVel(Vector3d& vel, unsigned int index) {
399
400 //velRot = $(A\cdot skew(I^{-1}j))^{T}refCoor$
401
402 if (index < atoms_.size()) {
403
404 Vector3d velRot;
405 Mat3x3d skewMat;;
406 Vector3d ref = refCoords_[index];
407 Vector3d ji = getJ();
408 Mat3x3d I = getI();
409
410 skewMat(0, 0) =0;
411 skewMat(0, 1) = ji[2] /I(2, 2);
412 skewMat(0, 2) = -ji[1] /I(1, 1);
413
414 skewMat(1, 0) = -ji[2] /I(2, 2);
415 skewMat(1, 1) = 0;
416 skewMat(1, 2) = ji[0]/I(0, 0);
417
418 skewMat(2, 0) =ji[1] /I(1, 1);
419 skewMat(2, 1) = -ji[0]/I(0, 0);
420 skewMat(2, 2) = 0;
421
422 velRot = (getA() * skewMat).transpose() * ref;
423
424 vel =getVel() + velRot;
425 return true;
426
427 } else {
428 std::cerr << index << " is an invalid index, current rigid body contains "
429 << atoms_.size() << "atoms" << std::endl;
430 return false;
431 }
432 }
433
434 bool RigidBody::getAtomVel(Vector3d& vel, Atom* atom) {
435
436 std::vector<Atom*>::iterator i;
437 i = std::find(atoms_.begin(), atoms_.end(), atom);
438 if (i != atoms_.end()) {
439 return getAtomVel(vel, i - atoms_.begin());
440 } else {
441 std::cerr << "Atom " << atom->getGlobalIndex()
442 <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;
443 return false;
444 }
445 }
446
447 bool RigidBody::getAtomRefCoor(Vector3d& coor, unsigned int index) {
448 if (index < atoms_.size()) {
449
450 coor = refCoords_[index];
451 return true;
452 } else {
453 std::cerr << index << " is an invalid index, current rigid body contains "
454 << atoms_.size() << "atoms" << std::endl;
455 return false;
456 }
457
458 }
459
460 bool RigidBody::getAtomRefCoor(Vector3d& coor, Atom* atom) {
461 std::vector<Atom*>::iterator i;
462 i = std::find(atoms_.begin(), atoms_.end(), atom);
463 if (i != atoms_.end()) {
464 //RigidBody class makes sure refCoords_ and atoms_ match each other
465 coor = refCoords_[i - atoms_.begin()];
466 return true;
467 } else {
468 std::cerr << "Atom " << atom->getGlobalIndex()
469 <<" does not belong to Rigid body "<< getGlobalIndex() << std::endl;
470 return false;
471 }
472
473 }
474
475
476 void RigidBody::addAtom(Atom* at, AtomStamp* ats) {
477
478 Vector3d coords;
479 Vector3d euler;
480
481
482 atoms_.push_back(at);
483
484 if( !ats->havePosition() ){
485 sprintf( painCave.errMsg,
486 "RigidBody error.\n"
487 "\tAtom %s does not have a position specified.\n"
488 "\tThis means RigidBody cannot set up reference coordinates.\n",
489 ats->getType() );
490 painCave.isFatal = 1;
491 simError();
492 }
493
494 coords[0] = ats->getPosX();
495 coords[1] = ats->getPosY();
496 coords[2] = ats->getPosZ();
497
498 refCoords_.push_back(coords);
499
500 RotMat3x3d identMat = RotMat3x3d::identity();
501
502 if (at->isDirectional()) {
503
504 if( !ats->haveOrientation() ){
505 sprintf( painCave.errMsg,
506 "RigidBody error.\n"
507 "\tAtom %s does not have an orientation specified.\n"
508 "\tThis means RigidBody cannot set up reference orientations.\n",
509 ats->getType() );
510 painCave.isFatal = 1;
511 simError();
512 }
513
514 euler[0] = ats->getEulerPhi() * NumericConstant::PI /180.0;
515 euler[1] = ats->getEulerTheta() * NumericConstant::PI /180.0;
516 euler[2] = ats->getEulerPsi() * NumericConstant::PI /180.0;
517
518 RotMat3x3d Atmp(euler);
519 refOrients_.push_back(Atmp);
520
521 }else {
522 refOrients_.push_back(identMat);
523 }
524
525
526 }
527
528 }
529