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

Comparing trunk/src/integrators/LDForceManager.cpp (file contents):
Revision 904 by tim, Thu Mar 16 22:50:48 2006 UTC vs.
Revision 909 by tim, Tue Mar 21 00:26:55 2006 UTC

# Line 50 | Line 50 | namespace oopse {
50      if (simParams->haveHydroPropFile()) {
51          hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
52      } else {
53 <        //error
53 >                sprintf( painCave.errMsg,
54 >                       "HydroPropFile keyword must be set if Lagevin Dynamics is used\n");
55 >                painCave.severity = OOPSE_ERROR;
56 >                painCave.isFatal = 1;
57 >                simError();  
58      }
59  
60      SimInfo::MoleculeIterator i;
# Line 64 | Line 68 | namespace oopse {
68              if (iter != hydroPropMap.end()) {
69                  hydroProps_.push_back(iter->second);
70              } else {
71 <                //error
71 >                sprintf( painCave.errMsg,
72 >                       "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
73 >                painCave.severity = OOPSE_ERROR;
74 >                painCave.isFatal = 1;
75 >                simError();  
76              }
77              
78             }
# Line 83 | Line 91 | namespace oopse {
91      while (ifs.getline(buffer, BufferSize)) {
92          StringTokenizer tokenizer(buffer);
93          HydroProp currProp;
94 <        if (tokenizer.countTokens() >= 67) {
94 >        if (tokenizer.countTokens() >= 40) {
95              std::string atomName = tokenizer.nextToken();
96 <            currProp.cod[0] = tokenizer.nextTokenAsDouble();
97 <            currProp.cod[1] = tokenizer.nextTokenAsDouble();
98 <            currProp.cod[2] = tokenizer.nextTokenAsDouble();
96 >            currProp.cor[0] = tokenizer.nextTokenAsDouble();
97 >            currProp.cor[1] = tokenizer.nextTokenAsDouble();
98 >            currProp.cor[2] = tokenizer.nextTokenAsDouble();
99 >            
100 >            currProp.Xirtt(0,0) = tokenizer.nextTokenAsDouble();
101 >            currProp.Xirtt(0,1) = tokenizer.nextTokenAsDouble();
102 >            currProp.Xirtt(0,2) = tokenizer.nextTokenAsDouble();
103 >            currProp.Xirtt(1,0) = tokenizer.nextTokenAsDouble();
104 >            currProp.Xirtt(1,1) = tokenizer.nextTokenAsDouble();
105 >            currProp.Xirtt(1,2) = tokenizer.nextTokenAsDouble();
106 >            currProp.Xirtt(2,0) = tokenizer.nextTokenAsDouble();
107 >            currProp.Xirtt(2,1) = tokenizer.nextTokenAsDouble();
108 >            currProp.Xirtt(2,2) = tokenizer.nextTokenAsDouble();
109  
110 <            currProp.Ddtt(0,0) = tokenizer.nextTokenAsDouble();
111 <            currProp.Ddtt(0,1) = tokenizer.nextTokenAsDouble();
112 <            currProp.Ddtt(0,2) = tokenizer.nextTokenAsDouble();
113 <            currProp.Ddtt(1,0) = tokenizer.nextTokenAsDouble();
114 <            currProp.Ddtt(1,1) = tokenizer.nextTokenAsDouble();
115 <            currProp.Ddtt(1,2) = tokenizer.nextTokenAsDouble();
116 <            currProp.Ddtt(2,0) = tokenizer.nextTokenAsDouble();
117 <            currProp.Ddtt(2,1) = tokenizer.nextTokenAsDouble();
118 <            currProp.Ddtt(2,2) = tokenizer.nextTokenAsDouble();
110 >            currProp.Xirrt(0,0) = tokenizer.nextTokenAsDouble();
111 >            currProp.Xirrt(0,1) = tokenizer.nextTokenAsDouble();
112 >            currProp.Xirrt(0,2) = tokenizer.nextTokenAsDouble();
113 >            currProp.Xirrt(1,0) = tokenizer.nextTokenAsDouble();
114 >            currProp.Xirrt(1,1) = tokenizer.nextTokenAsDouble();
115 >            currProp.Xirrt(1,2) = tokenizer.nextTokenAsDouble();
116 >            currProp.Xirrt(2,0) = tokenizer.nextTokenAsDouble();
117 >            currProp.Xirrt(2,1) = tokenizer.nextTokenAsDouble();
118 >            currProp.Xirrt(2,2) = tokenizer.nextTokenAsDouble();
119 >        
120 >            currProp.Xirtr(0,0) = tokenizer.nextTokenAsDouble();
121 >            currProp.Xirtr(0,1) = tokenizer.nextTokenAsDouble();
122 >            currProp.Xirtr(0,2) = tokenizer.nextTokenAsDouble();
123 >            currProp.Xirtr(1,0) = tokenizer.nextTokenAsDouble();
124 >            currProp.Xirtr(1,1) = tokenizer.nextTokenAsDouble();
125 >            currProp.Xirtr(1,2) = tokenizer.nextTokenAsDouble();
126 >            currProp.Xirtr(2,0) = tokenizer.nextTokenAsDouble();
127 >            currProp.Xirtr(2,1) = tokenizer.nextTokenAsDouble();
128 >            currProp.Xirtr(2,2) = tokenizer.nextTokenAsDouble();
129  
130 <            currProp.Ddtr(0,0) = tokenizer.nextTokenAsDouble();
131 <            currProp.Ddtr(0,1) = tokenizer.nextTokenAsDouble();
132 <            currProp.Ddtr(0,2) = tokenizer.nextTokenAsDouble();
133 <            currProp.Ddtr(1,0) = tokenizer.nextTokenAsDouble();
134 <            currProp.Ddtr(1,1) = tokenizer.nextTokenAsDouble();
135 <            currProp.Ddtr(1,2) = tokenizer.nextTokenAsDouble();
136 <            currProp.Ddtr(2,0) = tokenizer.nextTokenAsDouble();
137 <            currProp.Ddtr(2,1) = tokenizer.nextTokenAsDouble();
138 <            currProp.Ddtr(2,2) = tokenizer.nextTokenAsDouble();
130 >            currProp.Xirrr(0,0) = tokenizer.nextTokenAsDouble();
131 >            currProp.Xirrr(0,1) = tokenizer.nextTokenAsDouble();
132 >            currProp.Xirrr(0,2) = tokenizer.nextTokenAsDouble();
133 >            currProp.Xirrr(1,0) = tokenizer.nextTokenAsDouble();
134 >            currProp.Xirrr(1,1) = tokenizer.nextTokenAsDouble();
135 >            currProp.Xirrr(1,2) = tokenizer.nextTokenAsDouble();
136 >            currProp.Xirrr(2,0) = tokenizer.nextTokenAsDouble();
137 >            currProp.Xirrr(2,1) = tokenizer.nextTokenAsDouble();
138 >            currProp.Xirrr(2,2) = tokenizer.nextTokenAsDouble();
139  
140 <            currProp.Ddrr(0,0) = tokenizer.nextTokenAsDouble();
141 <            currProp.Ddrr(0,1) = tokenizer.nextTokenAsDouble();
142 <            currProp.Ddrr(0,2) = tokenizer.nextTokenAsDouble();
143 <            currProp.Ddrr(1,0) = tokenizer.nextTokenAsDouble();
144 <            currProp.Ddrr(1,1) = tokenizer.nextTokenAsDouble();
145 <            currProp.Ddrr(1,2) = tokenizer.nextTokenAsDouble();
118 <            currProp.Ddrr(2,0) = tokenizer.nextTokenAsDouble();
119 <            currProp.Ddrr(2,1) = tokenizer.nextTokenAsDouble();
120 <            currProp.Ddrr(2,2) = tokenizer.nextTokenAsDouble();                
121 <
122 <            currProp.Xidtt(0,0) = tokenizer.nextTokenAsDouble();
123 <            currProp.Xidtt(0,1) = tokenizer.nextTokenAsDouble();
124 <            currProp.Xidtt(0,2) = tokenizer.nextTokenAsDouble();
125 <            currProp.Xidtt(1,0) = tokenizer.nextTokenAsDouble();
126 <            currProp.Xidtt(1,1) = tokenizer.nextTokenAsDouble();
127 <            currProp.Xidtt(1,2) = tokenizer.nextTokenAsDouble();
128 <            currProp.Xidtt(2,0) = tokenizer.nextTokenAsDouble();
129 <            currProp.Xidtt(2,1) = tokenizer.nextTokenAsDouble();
130 <            currProp.Xidtt(2,2) = tokenizer.nextTokenAsDouble();
131 <
132 <            currProp.Xidrt(0,0) = tokenizer.nextTokenAsDouble();
133 <            currProp.Xidrt(0,1) = tokenizer.nextTokenAsDouble();
134 <            currProp.Xidrt(0,2) = tokenizer.nextTokenAsDouble();
135 <            currProp.Xidrt(1,0) = tokenizer.nextTokenAsDouble();
136 <            currProp.Xidrt(1,1) = tokenizer.nextTokenAsDouble();
137 <            currProp.Xidrt(1,2) = tokenizer.nextTokenAsDouble();
138 <            currProp.Xidrt(2,0) = tokenizer.nextTokenAsDouble();
139 <            currProp.Xidrt(2,1) = tokenizer.nextTokenAsDouble();
140 <            currProp.Xidrt(2,2) = tokenizer.nextTokenAsDouble();
141 <            
142 <            currProp.Xidtr(0,0) = tokenizer.nextTokenAsDouble();
143 <            currProp.Xidtr(0,1) = tokenizer.nextTokenAsDouble();
144 <            currProp.Xidtr(0,2) = tokenizer.nextTokenAsDouble();
145 <            currProp.Xidtr(1,0) = tokenizer.nextTokenAsDouble();
146 <            currProp.Xidtr(1,1) = tokenizer.nextTokenAsDouble();
147 <            currProp.Xidtr(1,2) = tokenizer.nextTokenAsDouble();
148 <            currProp.Xidtr(2,0) = tokenizer.nextTokenAsDouble();
149 <            currProp.Xidtr(2,1) = tokenizer.nextTokenAsDouble();
150 <            currProp.Xidtr(2,2) = tokenizer.nextTokenAsDouble();
140 >            SquareMatrix<double, 6> Xir;
141 >            Xir.setSubMatrix(0, 0, currProp.Xirtt);
142 >            Xir.setSubMatrix(0, 3, currProp.Xirrt);
143 >            Xir.setSubMatrix(3, 0, currProp.Xirtr);
144 >            Xir.setSubMatrix(3, 3, currProp.Xirrr);
145 >            CholeskyDecomposition(Xir, currProp.S);            
146  
152 –            currProp.Xidrr(0,0) = tokenizer.nextTokenAsDouble();
153 –            currProp.Xidrr(0,1) = tokenizer.nextTokenAsDouble();
154 –            currProp.Xidrr(0,2) = tokenizer.nextTokenAsDouble();
155 –            currProp.Xidrr(1,0) = tokenizer.nextTokenAsDouble();
156 –            currProp.Xidrr(1,1) = tokenizer.nextTokenAsDouble();
157 –            currProp.Xidrr(1,2) = tokenizer.nextTokenAsDouble();
158 –            currProp.Xidrr(2,0) = tokenizer.nextTokenAsDouble();
159 –            currProp.Xidrr(2,1) = tokenizer.nextTokenAsDouble();
160 –            currProp.Xidrr(2,2) = tokenizer.nextTokenAsDouble();
147              props.insert(std::map<std::string, HydroProp>::value_type(atomName, currProp));
148          }
149      }
# Line 203 | Line 189 | namespace oopse {
189                   omega[2] = angMom[2] /I(2, 2);
190               }
191  
192 <             //apply friction force and torque at center of diffusion
192 >             //apply friction force and torque at center of resistance
193               A = integrableObject->getA();
194               Atrans = A.transpose();
195 <             Vector3d rcd = Atrans * hydroProps_[index].cod;  
196 <             Vector3d vcd = vel + cross(omega, rcd);
197 <             vcd = A* vcd;
198 <             Vector3d frictionForce = -(hydroProps_[index].Xidtt * vcd + hydroProps_[index].Xidrt * omega);
199 <             frictionForce = Atrans*frictionForce;
200 <             integrableObject->addFrc(frictionForce);
201 <             Vector3d frictionTorque = - (hydroProps_[index].Xidtr * vcd + hydroProps_[index].Xidrr * omega);
202 <             frictionTorque = Atrans*frictionTorque;
203 <             integrableObject->addTrq(frictionTorque+ cross(rcd, frictionForce));
204 <            
205 <             //apply random force and torque at center of diffustion
206 <             Vector3d randomForce;
207 <             Vector3d randomTorque;
208 <             genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
209 <             randomForce = Atrans*randomForce;
210 <             randomTorque = Atrans* randomTorque;
211 <             integrableObject->addFrc(randomForce);            
212 <             integrableObject->addTrq(randomTorque + cross(rcd, randomForce ));
213 <            
195 >             Vector3d rcr = Atrans * hydroProps_[index].cor;  
196 >             Vector3d vcdLab = vel + cross(omega, rcr);
197 >             Vector3d vcdBody = A* vcdLab;
198 >             Vector3d frictionForceBody = -(hydroProps_[index].Xirtt * vcdBody + hydroProps_[index].Xirrt * omega);
199 >             Vector3d frictionForceLab = Atrans*frictionForceBody;
200 >             integrableObject->addFrc(frictionForceLab);
201 >             Vector3d frictionTorqueBody = - (hydroProps_[index].Xirtr * vcdBody + hydroProps_[index].Xirrr * omega);
202 >             Vector3d frictionTorqueLab = Atrans*frictionTorqueBody;
203 >             integrableObject->addTrq(frictionTorqueLab+ cross(rcr, frictionForceLab));
204 >
205 >             //apply random force and torque at center of resistance
206 >             Vector3d randomForceBody;
207 >             Vector3d randomTorqueBody;
208 >             genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
209 >             Vector3d randomForceLab = Atrans*randomForceBody;
210 >             Vector3d randomTorqueLab = Atrans* randomTorqueBody;
211 >             integrableObject->addFrc(randomForceLab);            
212 >             integrableObject->addTrq(randomTorqueLab + cross(rcr, randomForceLab ));            
213 >
214            } else {
215               //spheric atom
216 <             Vector3d frictionForce = -(hydroProps_[index].Xidtt *vel);    
216 >             Vector3d frictionForce = -(hydroProps_[index].Xirtt *vel);    
217               Vector3d randomForce;
218               Vector3d randomTorque;
219               genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
220  
235 –             //randomForce /= OOPSEConstant::energyConvert;
236 –             //randomTorque /= OOPSEConstant::energyConvert;
221               integrableObject->addFrc(frictionForce+randomForce);            
222            }
223  
# Line 249 | Line 233 | void LDForceManager::genRandomForceAndTorque(Vector3d&
233    }
234  
235   void LDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, double variance) {
252 –    /*
253 –    SquareMatrix<double, 6> Dd;
254 –    SquareMatrix<double, 6> S;
255 –    Vector<double, 6> Z;
256 –    Vector<double, 6> generalForce;
257 –    Dd.setSubMatrix(0, 0, hydroProps_[index].Ddtt);
258 –    Dd.setSubMatrix(0, 3, hydroProps_[index].Ddtr.transpose());
259 –    Dd.setSubMatrix(3, 0, hydroProps_[index].Ddtr);
260 –    Dd.setSubMatrix(3, 3, hydroProps_[index].Ddrr);
261 –    CholeskyDecomposition(Dd, S);
262 –    */
236  
237 <    SquareMatrix<double, 6> Xid;
265 <    SquareMatrix<double, 6> S;
237 >
238      Vector<double, 6> Z;
239      Vector<double, 6> generalForce;
268 –    Xid.setSubMatrix(0, 0, hydroProps_[index].Xidtt);
269 –    Xid.setSubMatrix(0, 3, hydroProps_[index].Xidrt);
270 –    Xid.setSubMatrix(3, 0, hydroProps_[index].Xidtr);
271 –    Xid.setSubMatrix(3, 3, hydroProps_[index].Xidrr);
272 –    CholeskyDecomposition(Xid, S);
240  
274 –    /*
275 –    Xid *= variance;
276 –    Z[0] = randNumGen_.randNorm(0, 1.0);
277 –    Z[1] = randNumGen_.randNorm(0, 1.0);
278 –    Z[2] = randNumGen_.randNorm(0, 1.0);
279 –    Z[3] = randNumGen_.randNorm(0, 1.0);
280 –    Z[4] = randNumGen_.randNorm(0, 1.0);
281 –    Z[5] = randNumGen_.randNorm(0, 1.0);
282 –    */
241          
242      Z[0] = randNumGen_.randNorm(0, variance);
243      Z[1] = randNumGen_.randNorm(0, variance);
# Line 289 | Line 247 | void LDForceManager::genRandomForceAndTorque(Vector3d&
247      Z[5] = randNumGen_.randNorm(0, variance);
248      
249  
250 <    generalForce = S*Z;
250 >    generalForce = hydroProps_[index].S*Z;
251      
252      force[0] = generalForce[0];
253      force[1] = generalForce[1];

Diff Legend

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