ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/visitors/ReplacementVisitor.cpp
(Generate patch)

Comparing:
trunk/src/visitors/ReplacementVisitor.cpp (file contents), Revision 1455 by gezelter, Thu Jun 24 20:44:18 2010 UTC vs.
branches/development/src/visitors/ReplacementVisitor.cpp (file contents), Revision 1874 by gezelter, Wed May 15 15:09:35 2013 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  
43   #include <cstring>
# Line 55 | Line 56 | namespace OpenMD {
56      myTypes_.clear();
57    }
58    
59 <  bool ReplacementVisitor::isReplacedAtom(const std::string& atomType) {
59 >  bool ReplacementVisitor::isReplacedAtom(const std::string &atomType) {
60      std::set<std::string>::iterator strIter;
61      strIter = myTypes_.find(atomType);
62      return strIter != myTypes_.end() ? true : false;
63    }
64  
65 <  void ReplacementVisitor::addSite(const std::string& name,
66 <                                   const Vector3d refPos) {
65 >  void ReplacementVisitor::addSite(const std::string &name,
66 >                                   const Vector3d &refPos) {
67      AtomInfo* atomInfo = new AtomInfo();
68      atomInfo->atomTypeName = name;
69      atomInfo->pos = refPos;
70      sites_->addAtomInfo(atomInfo);
71    }
72 <  void ReplacementVisitor::addSite(const std::string& name,
73 <                                   const Vector3d refPos,
74 <                                   const Vector3d refVec) {
72 >  void ReplacementVisitor::addSite(const std::string &name,
73 >                                   const Vector3d &refPos,
74 >                                   const Vector3d &refVec) {
75      AtomInfo* atomInfo = new AtomInfo();
76      atomInfo->atomTypeName = name;
77      atomInfo->pos = refPos;
78 <    atomInfo->dipole = refVec;
78 >    atomInfo->vec = refVec;
79      atomInfo->hasVector = true;
80      sites_->addAtomInfo(atomInfo);
81    }
82    
83    void ReplacementVisitor::visit(DirectionalAtom *datom) {
83 –    std::vector<AtomInfo*>atoms;
84      
85 <    RotMat3x3d   rotMatrix;
86 <    RotMat3x3d   rotTrans;
87 <    AtomInfo *   atomInfo;
85 >    RotMat3x3d   A;
86 >    RotMat3x3d   Atrans;
87 >    Mat3x3d      I;
88      Vector3d     pos;
89 +    Vector3d     vel;
90 +    Vector3d     frc;
91 +    Vector3d     trq;
92 +    Vector3d     j;
93 +    Mat3x3d      skewMat;
94 +
95      Vector3d     newVec;
96 <    Quat4d       q;
96 >    AtomInfo *   atomInfo;
97      AtomData *   atomData;
98      GenericData *data;
99      bool         haveAtomData;
# Line 111 | Line 117 | namespace OpenMD {
117        atomData = new AtomData;
118        haveAtomData = false;
119      }
120 <    
115 <    
120 >        
121      pos = datom->getPos();
122 <    q = datom->getQ();
123 <    rotMatrix = datom->getA();
122 >    vel = datom->getVel();
123 >
124 >    j   = datom->getJ();
125 >    I   = datom->getI();
126 >    A   = datom->getA();
127 >
128 >    skewMat(0, 0) =  0;
129 >    skewMat(0, 1) =  j[2] / I(2, 2);
130 >    skewMat(0, 2) = -j[1] / I(1, 1);    
131 >    skewMat(1, 0) = -j[2] / I(2, 2);
132 >    skewMat(1, 1) =  0;
133 >    skewMat(1, 2) =  j[0] / I(0, 0);    
134 >    skewMat(2, 0) =  j[1] / I(1, 1);
135 >    skewMat(2, 1) = -j[0] / I(0, 0);
136 >    skewMat(2, 2) =  0;
137 >    Mat3x3d mat = (A * skewMat).transpose();
138      
139      // We need A^T to convert from body-fixed to space-fixed:
140 <    rotTrans = rotMatrix.transpose();
140 >    Atrans = A.transpose();
141      
142      AtomInfo* siteInfo;
143      std::vector<AtomInfo*>::iterator iter;
# Line 126 | Line 145 | namespace OpenMD {
145      for( siteInfo = sites_->beginAtomInfo(iter); siteInfo;
146           siteInfo = sites_->nextAtomInfo(iter) ) {
147  
148 <      newVec = rotTrans * siteInfo->pos;    
148 >      newVec = Atrans * siteInfo->pos;    
149        
150        atomInfo = new AtomInfo;
151        atomInfo->atomTypeName = siteInfo->atomTypeName;
152        atomInfo->pos = pos + newVec;
153        
154        if (siteInfo->hasVector) {
155 <        newVec = rotTrans * siteInfo->dipole;
156 <        atomInfo->dipole = newVec;
155 >        newVec = Atrans * siteInfo->vec;
156 >        atomInfo->vec = newVec;
157        } else {
158 <        atomInfo->dipole = V3Zero;
158 >        atomInfo->vec = V3Zero;
159        }
160 <      
160 >
161 >      atomInfo->vel = vel + mat * siteInfo->pos;
162 >      atomInfo->hasVelocity = true;
163 >            
164        atomData->addAtomInfo(atomInfo);
165      }
166      if (!haveAtomData) {

Diff Legend

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