ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Shake.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/Shake.cpp (file contents):
Revision 1253 by tim, Thu Jun 3 21:51:55 2004 UTC vs.
Revision 1254 by tim, Wed Jun 9 16:16:33 2004 UTC

# Line 2 | Line 2
2   #include "SimInfo.hpp"
3   #include <cmath>
4   #include "simError.h"
5 + #include "MatVec3.h"
6   ////////////////////////////////////////////////////////////////////////////////
7   //Implementation of DCShakeFunctor
8   ////////////////////////////////////////////////////////////////////////////////
# Line 71 | Line 72 | int DCShakeFunctor::operator()(ConstraintAtom* consAto
72      dy = rab[1] * gab;
73      dz = rab[2] * gab;
74  
75 +    //set atom1's position
76      posA[0] += rma * dx;
77      posA[1] += rma * dy;
78      posA[2] += rma * dz;
77 –
79      consAtom1->setPos(posA);
80  
81 +    //set atom2's position
82      posB[0] -= rmb * dx;
83      posB[1] -= rmb * dy;
84      posB[2] -= rmb * dz;
83 –
85      consAtom2->setPos(posB);
86  
87 <    dx = dx / dt;
88 <    dy = dy / dt;
89 <    dz = dz / dt;
87 >    //dx = dx / dt;
88 >    //dy = dy / dt;
89 >    //dz = dz / dt;
90  
91 <    consAtom1->getVel(velA);
91 >    ////set atom1's velocity
92 >    //consAtom1->getVel(velA);
93 >    //velA[0] += rma * dx;
94 >    //velA[1] += rma * dy;
95 >    //velA[2] += rma * dz;
96 >    //consAtom1->setVel(velA);
97  
98 <    velA[0] += rma * dx;
99 <    velA[1] += rma * dy;
100 <    velA[2] += rma * dz;
98 >    ////set atom2's velocity
99 >    //consAtom2->getVel(velB);
100 >    //velB[0] -= rmb * dx;
101 >    //velB[1] -= rmb * dy;
102 >    //velB[2] -= rmb * dz;
103 >    //consAtom2->setVel(velB);
104  
96 –    consAtom1->setVel(velA);
97 –
98 –    consAtom2->getVel(velB);
99 –
100 –    velB[0] -= rmb * dx;
101 –    velB[1] -= rmb * dy;
102 –    velB[2] -= rmb * dz;
103 –
104 –    consAtom2->setVel(velB);
105 –
105      return consSuccess;
106    }
107    else
# Line 115 | Line 114 | int DCShakeFunctor::operator()(ConstraintAtom* consAto
114    return consElemHandlerFail;
115   }
116  
117 + /**
118 +  * QSHAKE Algorithm
119 +  * Reference
120 +  * T.R. Forester and W. Smith, SHAKE, Rattle and Roll: Efficient Constraint Algorithms for Linked
121 +  * Rigid Bodies, J. Comp. Chem., 19(1), 102 -111 (1998)
122 +  */
123   int DCShakeFunctor::operator()(ConstraintRigidBody* consRB1, ConstraintRigidBody* consRB2){
124 <  return consElemHandlerFail;
124 >  double posA[3];
125 >  double posB[3];
126 >  double oldPosA[3];
127 >  double oldPosB[3];
128 >  double velA[3];
129 >  double velB[3];
130 >  double pab[3];
131 >  double tempPab[3];
132 >  double rab[3];
133 >  double rma, rmb;
134 >  double dx, dy, dz;
135 >  double rpab;
136 >  double rabsq, pabsq, rpabsq;
137 >  double diffsq;
138 >  double gab;
139 >  double dt;
140 >  double dt2;
141 >  double consForce[3];
142 >
143 >  const int conRBMaxIter = 10;
144 >  
145 >  dt = info->dt;
146 >  dt2 = dt * dt;
147 >
148 >  consRB1->getOldAtomPos(oldPosA);
149 >  consRB2->getOldAtomPos(oldPosB);      
150 >
151 >
152 >  for(int i=0 ; i < conRBMaxIter; i++){  
153 >  consRB1->getCurAtomPos(posA);
154 >  consRB2->getCurAtomPos(posB);
155 >
156 >  
157 >  pab[0] = posA[0] - posB[0];
158 >  pab[1] = posA[1] - posB[1];
159 >  pab[2] = posA[2] - posB[2];
160 >
161 >  //periodic boundary condition
162 >
163 >  info->wrapVector(pab);
164 >
165 >  pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
166 >
167 >  rabsq = curPair->getBondLength2();
168 >  diffsq = rabsq - pabsq;
169 >
170 >  // the original rattle code from alan tidesley
171 >  if (fabs(diffsq) > (consTolerance * rabsq * 2)){
172 >      
173 >    rab[0] = oldPosA[0] - oldPosB[0];
174 >    rab[1] = oldPosA[1] - oldPosB[1];
175 >    rab[2] = oldPosA[2] - oldPosB[2];
176 >
177 >    info->wrapVector(rab);
178 >
179 >    rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
180 >
181 >    rpabsq = rpab * rpab;
182 >
183 >
184 >    if (rpabsq < (rabsq * -diffsq)){
185 >      return consFail;
186 >    }
187 >
188 >    //rma = 1.0 / consRB1->getMass();
189 >    //rmb = 1.0 / consRB2->getMass();
190 >    
191 >    tempPab[0] = pab[0] / sqrt(pabsq);
192 >    tempPab[1] = pab[1] / sqrt(pabsq);
193 >    tempPab[2] = pab[2] / sqrt(pabsq);
194 >    
195 >    rma = getEffInvMass(consRB1, tempPab);
196 >
197 >    tempPab[0] = -tempPab[0];
198 >    tempPab[1] = -tempPab[1];
199 >    tempPab[2] = -tempPab[2];    
200 >    rmb = getEffInvMass(consRB2, tempPab);
201 >
202 >    gab = diffsq / (2.0* dt * dt * (rma + rmb) * rpab) ;
203 >    consForce[0] = rab[0] * gab;
204 >    consForce[1] = rab[1] * gab;
205 >    consForce[2] = rab[2] * gab;
206 >
207 >    //integrate consRB1 using constraint force;
208 >    integrate(consRB1,consForce);
209 >
210 >    //integrate consRB2 using constraint force;
211 >    consForce[0] = -consForce[0];
212 >    consForce[1] = -consForce[1];
213 >    consForce[2] = -consForce[2];
214 >    integrate(consRB2,consForce);    
215 >    
216 >  }
217 >  else{
218 >    if (i ==0)
219 >      return consAlready;
220 >    else
221 >      return consSuccess;
222 >  }
223 >  }
224 >
225 >  return consExceedMaxIter;
226   }
227  
228 + double DCShakeFunctor::getEffInvMass(ConstraintRigidBody* consRB, double bondDir[3]){
229 +  double effInvMass;  //effective inversse mass
230 +  double effInvMassCorr;  //correction for effective inverse mass
231 +  double aTrans[3][3];
232 +  double a[3][3];
233 +
234 +  double IFrame[3][3];
235 +  double IBody[3][3];
236 +  double invI[3][3];
237 +  double refCoor[3];
238 +  double refCrossBond[3];
239 +  double tempVec1[3];
240 +  double tempVec2[3];
241 +
242 +  effInvMass = 1.0 / consRB ->getMass();
243 +  
244 +  consRB->getRefCoor(refCoor);
245 +  consRB->getA(a);
246 +  consRB->getI(IBody);
247 +
248 +  crossProduct3(refCoor, bondDir, refCrossBond);
249 +
250 +  matMul3(a, IBody, IFrame);
251 +
252 +  invertMat3(IFrame, invI);
253 +
254 +  matVecMul3(invI, refCrossBond, tempVec1);
255 +  
256 +  crossProduct3(tempVec1, refCoor, tempVec2);
257 +
258 +  effInvMassCorr = dotProduct3(tempVec1, bondDir);
259 +    
260 +  effInvMass += effInvMassCorr;
261 +
262 +  return effInvMass;
263 + }
264 +
265 + void DCShakeFunctor::integrate(ConstraintRigidBody* consRB, double force[3]){
266 +
267 + }
268   ////////////////////////////////////////////////////////////////////////////////
269   //Implementation of JCShakeFunctor
270   ////////////////////////////////////////////////////////////////////////////////
# Line 135 | Line 281 | int JCShakeFunctor::operator()(ConstraintRigidBody* co
281  
282   int JCShakeFunctor::operator()(ConstraintRigidBody* consRB1, ConstraintRigidBody* consRB2){
283    return consElemHandlerFail;
284 < }
284 > }

Diff Legend

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