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 1442 by gezelter, Mon May 10 17:28:26 2010 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

# Line 35 | Line 35
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).                        
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42   #include <algorithm>
43   #include <math.h>
# Line 52 | Line 53 | namespace OpenMD {
53    void RigidBody::setPrevA(const RotMat3x3d& a) {
54      ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
55      
56 <    for (int i =0 ; i < atoms_.size(); ++i){
56 >    for (unsigned int i = 0 ; i < atoms_.size(); ++i){
57        if (atoms_[i]->isDirectional()) {
58          atoms_[i]->setPrevA(refOrients_[i].transpose() * a);
59        }
# Line 64 | Line 65 | namespace OpenMD {
65    void RigidBody::setA(const RotMat3x3d& a) {
66      ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
67  
68 <    for (int i =0 ; i < atoms_.size(); ++i){
68 >    for (unsigned int i = 0 ; i < atoms_.size(); ++i){
69        if (atoms_[i]->isDirectional()) {
70          atoms_[i]->setA(refOrients_[i].transpose() * a);
71        }
# Line 73 | Line 74 | namespace OpenMD {
74    
75    void RigidBody::setA(const RotMat3x3d& a, int snapshotNo) {
76      ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
77 <    
78 <    //((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;    
78 <    
79 <    for (int i =0 ; i < atoms_.size(); ++i){
77 >        
78 >    for (unsigned int i = 0 ; i < atoms_.size(); ++i){
79        if (atoms_[i]->isDirectional()) {
80          atoms_[i]->setA(refOrients_[i].transpose() * a, snapshotNo);
81        }
# Line 93 | Line 92 | namespace OpenMD {
92      Vector3d force;
93      Vector3d torque;
94      Vector3d myEuler;
95 <    RealType phi, theta, psi;
95 >    RealType phi, theta;
96 >    // RealType psi;
97      RealType cphi, sphi, ctheta, stheta;
98      Vector3d ephi;
99      Vector3d etheta;
# Line 105 | Line 105 | namespace OpenMD {
105      
106      phi = myEuler[0];
107      theta = myEuler[1];
108 <    psi = myEuler[2];
108 >    // psi = myEuler[2];
109      
110      cphi = cos(phi);
111      sphi = sin(phi);
# Line 226 | Line 226 | namespace OpenMD {
226      Vector3d apos;
227      Vector3d rpos;
228      Vector3d frc(0.0);
229 <    Vector3d trq(0.0);    
229 >    Vector3d trq(0.0);
230 >    Vector3d ef(0.0);
231      Vector3d pos = this->getPos();
232 <    for (int i = 0; i < atoms_.size(); i++) {
232 >    AtomType* atype;
233 >    int eCount = 0;
234 >    
235 >    int sl = ((snapshotMan_->getCurrentSnapshot())->*storage_).getStorageLayout();
236 >    
237 >    for (unsigned int i = 0; i < atoms_.size(); i++) {
238  
239 +      atype = atoms_[i]->getAtomType();
240 +
241        afrc = atoms_[i]->getFrc();
242        apos = atoms_[i]->getPos();
243        rpos = apos - pos;
# Line 246 | Line 254 | namespace OpenMD {
254        if (atoms_[i]->isDirectional()) {
255          atrq = atoms_[i]->getTrq();
256          trq += atrq;
257 <      }      
257 >      }
258 >
259 >      if ((sl & DataStorage::dslElectricField) && (atype->isElectrostatic())) {
260 >        ef += atoms_[i]->getElectricField();
261 >        eCount++;
262 >      }
263      }        
264      addFrc(frc);
265      addTrq(trq);    
266 +
267 +    if (sl & DataStorage::dslElectricField)  {
268 +      ef /= eCount;
269 +      setElectricField(ef);
270 +    }
271 +
272    }
273  
274    Mat3x3d RigidBody::calcForcesAndTorquesAndVirial() {
# Line 259 | Line 278 | namespace OpenMD {
278      Vector3d rpos;
279      Vector3d dfrc;
280      Vector3d frc(0.0);
281 <    Vector3d trq(0.0);    
281 >    Vector3d trq(0.0);
282 >    Vector3d ef(0.0);
283 >    AtomType* atype;
284 >    int eCount = 0;
285 >
286      Vector3d pos = this->getPos();
287      Mat3x3d tau_(0.0);
288  
289 <    for (int i = 0; i < atoms_.size(); i++) {
289 >    int sl = ((snapshotMan_->getCurrentSnapshot())->*storage_).getStorageLayout();
290 >
291 >    for (unsigned int i = 0; i < atoms_.size(); i++) {
292        
293 +      atype = atoms_[i]->getAtomType();
294 +
295        afrc = atoms_[i]->getFrc();
296        apos = atoms_[i]->getPos();
297        rpos = apos - pos;
# Line 282 | Line 309 | namespace OpenMD {
309          atrq = atoms_[i]->getTrq();
310          trq += atrq;
311        }
312 +
313 +      if ((sl & DataStorage::dslElectricField) && (atype->isElectrostatic())) {
314 +        ef += atoms_[i]->getElectricField();
315 +        eCount++;
316 +      }
317        
318        tau_(0,0) -= rpos[0]*afrc[0];
319        tau_(0,1) -= rpos[0]*afrc[1];
# Line 296 | Line 328 | namespace OpenMD {
328      }
329      addFrc(frc);
330      addTrq(trq);
331 +
332 +    if (sl & DataStorage::dslElectricField) {
333 +      ef /= eCount;
334 +      setElectricField(ef);
335 +    }
336 +
337      return tau_;
338    }
339  
# Line 317 | Line 355 | namespace OpenMD {
355  
356        if (atoms_[i]->isDirectional()) {
357            
358 <        dAtom = (DirectionalAtom *) atoms_[i];
358 >        dAtom = dynamic_cast<DirectionalAtom *>(atoms_[i]);
359          dAtom->setA(refOrients_[i].transpose() * a);
360        }
361  
# Line 344 | Line 382 | namespace OpenMD {
382  
383        if (atoms_[i]->isDirectional()) {
384            
385 <        dAtom = (DirectionalAtom *) atoms_[i];
385 >        dAtom = dynamic_cast<DirectionalAtom *>(atoms_[i]);
386          dAtom->setA(refOrients_[i].transpose() * a, frame);
387        }
388  
# Line 353 | Line 391 | namespace OpenMD {
391    }
392  
393    void RigidBody::updateAtomVel() {
394 <    Mat3x3d skewMat;;
394 >    Mat3x3d skewMat;
395  
396      Vector3d ji = getJ();
397      Mat3x3d I =  getI();
# Line 375 | Line 413 | namespace OpenMD {
413  
414  
415      Vector3d velRot;        
416 <    for (int i =0 ; i < refCoords_.size(); ++i) {
416 >    for (unsigned int i = 0 ; i < refCoords_.size(); ++i) {
417        atoms_[i]->setVel(rbVel + mat * refCoords_[i]);
418      }
419  
# Line 404 | Line 442 | namespace OpenMD {
442  
443  
444      Vector3d velRot;        
445 <    for (int i =0 ; i < refCoords_.size(); ++i) {
445 >    for (unsigned int i = 0 ; i < refCoords_.size(); ++i) {
446        atoms_[i]->setVel(rbVel + mat * refCoords_[i], frame);
447      }
448  
# Line 451 | Line 489 | namespace OpenMD {
489        Vector3d ji = getJ();
490        Mat3x3d I =  getI();
491  
492 <      skewMat(0, 0) =0;
493 <      skewMat(0, 1) = ji[2] /I(2, 2);
494 <      skewMat(0, 2) = -ji[1] /I(1, 1);
492 >      skewMat(0, 0) =  0;
493 >      skewMat(0, 1) =  ji[2] / I(2, 2);
494 >      skewMat(0, 2) = -ji[1] / I(1, 1);
495  
496 <      skewMat(1, 0) = -ji[2] /I(2, 2);
497 <      skewMat(1, 1) = 0;
498 <      skewMat(1, 2) = ji[0]/I(0, 0);
496 >      skewMat(1, 0) = -ji[2] / I(2, 2);
497 >      skewMat(1, 1) =  0;
498 >      skewMat(1, 2) =  ji[0]/ I(0, 0);
499  
500 <      skewMat(2, 0) =ji[1] /I(1, 1);
501 <      skewMat(2, 1) = -ji[0]/I(0, 0);
502 <      skewMat(2, 2) = 0;
500 >      skewMat(2, 0) =  ji[1] / I(1, 1);
501 >      skewMat(2, 1) = -ji[0] / I(0, 0);
502 >      skewMat(2, 2) =  0;
503  
504        velRot = (getA() * skewMat).transpose() * ref;
505  
506 <      vel =getVel() + velRot;
506 >      vel = getVel() + velRot;
507        return true;
508          
509      } else {
510 <      std::cerr << index << " is an invalid index, current rigid body contains "
510 >      std::cerr << index
511 >                << " is an invalid index, current rigid body contains "
512                  << atoms_.size() << "atoms" << std::endl;
513        return false;
514      }

Diff Legend

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