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

Comparing trunk/src/primitives/DirectionalAtom.cpp (file contents):
Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 507 by gezelter, Fri Apr 15 22:04:00 2005 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 40 | Line 40
40   */
41  
42   #include "primitives/DirectionalAtom.hpp"
43 <
43 > #include "utils/simError.h"
44   namespace oopse {
45  
46 < DirectionalAtom::DirectionalAtom(DirectionalAtomType* dAtomType)
47 <                         : Atom(dAtomType){
48 <    objType_= otDAtom;
49 <    if (dAtomType->isMultipole()) {
46 >  DirectionalAtom::DirectionalAtom(DirectionalAtomType* dAtomType)
47 >    : Atom(dAtomType){
48 >      objType_= otDAtom;
49 >      if (dAtomType->isMultipole()) {
50          electroBodyFrame_ = dAtomType->getElectroBodyFrame();
51 +      }
52 +
53 +      //check if one of the diagonal inertia tensor of this directional atom  is zero
54 +      int nLinearAxis = 0;
55 +      Mat3x3d inertiaTensor = getI();
56 +      for (int i = 0; i < 3; i++) {    
57 +        if (fabs(inertiaTensor(i, i)) < oopse::epsilon) {
58 +          linear_ = true;
59 +          linearAxis_ = i;
60 +          ++ nLinearAxis;
61 +        }
62 +      }
63 +
64 +      if (nLinearAxis > 1) {
65 +        sprintf( painCave.errMsg,
66 +                 "Directional Atom error.\n"
67 +                 "\tOOPSE found more than one axis in this directional atom with a vanishing \n"
68 +                 "\tmoment of inertia.");
69 +        painCave.isFatal = 1;
70 +        simError();
71 +      }
72 +      
73      }
52 – }
74  
75 < Mat3x3d DirectionalAtom::getI() {
75 >  Mat3x3d DirectionalAtom::getI() {
76      return static_cast<DirectionalAtomType*>(getAtomType())->getI();
77 < }    
77 >  }    
78  
79 < void DirectionalAtom::setPrevA(const RotMat3x3d& a) {
79 >  void DirectionalAtom::setPrevA(const RotMat3x3d& a) {
80      ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
81      if (atomType_->isMultipole()) {
82 <        ((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;
82 >      ((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;
83      }
84 < }
84 >  }
85  
86        
87 < void DirectionalAtom::setA(const RotMat3x3d& a) {
87 >  void DirectionalAtom::setA(const RotMat3x3d& a) {
88      ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
89  
90      if (atomType_->isMultipole()) {
91 <        ((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;
91 >      ((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;
92      }
93 < }    
93 >  }    
94      
95 < void DirectionalAtom::setA(const RotMat3x3d& a, int snapshotNo) {
95 >  void DirectionalAtom::setA(const RotMat3x3d& a, int snapshotNo) {
96      ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
97  
98      if (atomType_->isMultipole()) {
99 <        ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;    
99 >      ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;    
100      }
101 < }    
101 >  }    
102  
103 < void DirectionalAtom::rotateBy(const RotMat3x3d& m) {
103 >  void DirectionalAtom::rotateBy(const RotMat3x3d& m) {
104      setA(m *getA());
105 < }
105 >  }
106  
107 < std::vector<double> DirectionalAtom::getGrad() {
107 >  std::vector<double> DirectionalAtom::getGrad() {
108      std::vector<double> grad(6, 0.0);
109      Vector3d force;
110      Vector3d torque;
# Line 123 | Line 144 | std::vector<double> DirectionalAtom::getGrad() {
144  
145      //gradient is equal to -force
146      for (int j = 0 ; j<3; j++)
147 <        grad[j] = -force[j];
147 >      grad[j] = -force[j];
148  
149      for (int j = 0; j < 3; j++ ) {
150  
151 <        grad[3] += torque[j]*ephi[j];
152 <        grad[4] += torque[j]*etheta[j];
153 <        grad[5] += torque[j]*epsi[j];
151 >      grad[3] += torque[j]*ephi[j];
152 >      grad[4] += torque[j]*etheta[j];
153 >      grad[5] += torque[j]*epsi[j];
154  
155      }
156      
157      return grad;
158 < }    
158 >  }    
159  
160 < void DirectionalAtom::accept(BaseVisitor* v) {
160 >  void DirectionalAtom::accept(BaseVisitor* v) {
161      v->visit(this);
162 < }    
162 >  }    
163  
164   }
165  

Diff Legend

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