ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/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 908 by tim, Mon Mar 20 19:12:14 2006 UTC

# Line 83 | Line 83 | namespace oopse {
83      while (ifs.getline(buffer, BufferSize)) {
84          StringTokenizer tokenizer(buffer);
85          HydroProp currProp;
86 <        if (tokenizer.countTokens() >= 67) {
86 >        if (tokenizer.countTokens() >= 40) {
87              std::string atomName = tokenizer.nextToken();
88 <            currProp.cod[0] = tokenizer.nextTokenAsDouble();
89 <            currProp.cod[1] = tokenizer.nextTokenAsDouble();
90 <            currProp.cod[2] = tokenizer.nextTokenAsDouble();
88 >            currProp.cor[0] = tokenizer.nextTokenAsDouble();
89 >            currProp.cor[1] = tokenizer.nextTokenAsDouble();
90 >            currProp.cor[2] = tokenizer.nextTokenAsDouble();
91 >            
92 >            currProp.Xirtt(0,0) = tokenizer.nextTokenAsDouble();
93 >            currProp.Xirtt(0,1) = tokenizer.nextTokenAsDouble();
94 >            currProp.Xirtt(0,2) = tokenizer.nextTokenAsDouble();
95 >            currProp.Xirtt(1,0) = tokenizer.nextTokenAsDouble();
96 >            currProp.Xirtt(1,1) = tokenizer.nextTokenAsDouble();
97 >            currProp.Xirtt(1,2) = tokenizer.nextTokenAsDouble();
98 >            currProp.Xirtt(2,0) = tokenizer.nextTokenAsDouble();
99 >            currProp.Xirtt(2,1) = tokenizer.nextTokenAsDouble();
100 >            currProp.Xirtt(2,2) = tokenizer.nextTokenAsDouble();
101  
102 <            currProp.Ddtt(0,0) = tokenizer.nextTokenAsDouble();
103 <            currProp.Ddtt(0,1) = tokenizer.nextTokenAsDouble();
104 <            currProp.Ddtt(0,2) = tokenizer.nextTokenAsDouble();
105 <            currProp.Ddtt(1,0) = tokenizer.nextTokenAsDouble();
106 <            currProp.Ddtt(1,1) = tokenizer.nextTokenAsDouble();
107 <            currProp.Ddtt(1,2) = tokenizer.nextTokenAsDouble();
108 <            currProp.Ddtt(2,0) = tokenizer.nextTokenAsDouble();
109 <            currProp.Ddtt(2,1) = tokenizer.nextTokenAsDouble();
110 <            currProp.Ddtt(2,2) = tokenizer.nextTokenAsDouble();
102 >            currProp.Xirrt(0,0) = tokenizer.nextTokenAsDouble();
103 >            currProp.Xirrt(0,1) = tokenizer.nextTokenAsDouble();
104 >            currProp.Xirrt(0,2) = tokenizer.nextTokenAsDouble();
105 >            currProp.Xirrt(1,0) = tokenizer.nextTokenAsDouble();
106 >            currProp.Xirrt(1,1) = tokenizer.nextTokenAsDouble();
107 >            currProp.Xirrt(1,2) = tokenizer.nextTokenAsDouble();
108 >            currProp.Xirrt(2,0) = tokenizer.nextTokenAsDouble();
109 >            currProp.Xirrt(2,1) = tokenizer.nextTokenAsDouble();
110 >            currProp.Xirrt(2,2) = tokenizer.nextTokenAsDouble();
111 >        
112 >            currProp.Xirtr(0,0) = tokenizer.nextTokenAsDouble();
113 >            currProp.Xirtr(0,1) = tokenizer.nextTokenAsDouble();
114 >            currProp.Xirtr(0,2) = tokenizer.nextTokenAsDouble();
115 >            currProp.Xirtr(1,0) = tokenizer.nextTokenAsDouble();
116 >            currProp.Xirtr(1,1) = tokenizer.nextTokenAsDouble();
117 >            currProp.Xirtr(1,2) = tokenizer.nextTokenAsDouble();
118 >            currProp.Xirtr(2,0) = tokenizer.nextTokenAsDouble();
119 >            currProp.Xirtr(2,1) = tokenizer.nextTokenAsDouble();
120 >            currProp.Xirtr(2,2) = tokenizer.nextTokenAsDouble();
121  
122 <            currProp.Ddtr(0,0) = tokenizer.nextTokenAsDouble();
123 <            currProp.Ddtr(0,1) = tokenizer.nextTokenAsDouble();
124 <            currProp.Ddtr(0,2) = tokenizer.nextTokenAsDouble();
125 <            currProp.Ddtr(1,0) = tokenizer.nextTokenAsDouble();
126 <            currProp.Ddtr(1,1) = tokenizer.nextTokenAsDouble();
127 <            currProp.Ddtr(1,2) = tokenizer.nextTokenAsDouble();
128 <            currProp.Ddtr(2,0) = tokenizer.nextTokenAsDouble();
129 <            currProp.Ddtr(2,1) = tokenizer.nextTokenAsDouble();
130 <            currProp.Ddtr(2,2) = tokenizer.nextTokenAsDouble();
122 >            currProp.Xirrr(0,0) = tokenizer.nextTokenAsDouble();
123 >            currProp.Xirrr(0,1) = tokenizer.nextTokenAsDouble();
124 >            currProp.Xirrr(0,2) = tokenizer.nextTokenAsDouble();
125 >            currProp.Xirrr(1,0) = tokenizer.nextTokenAsDouble();
126 >            currProp.Xirrr(1,1) = tokenizer.nextTokenAsDouble();
127 >            currProp.Xirrr(1,2) = tokenizer.nextTokenAsDouble();
128 >            currProp.Xirrr(2,0) = tokenizer.nextTokenAsDouble();
129 >            currProp.Xirrr(2,1) = tokenizer.nextTokenAsDouble();
130 >            currProp.Xirrr(2,2) = tokenizer.nextTokenAsDouble();
131  
132 <            currProp.Ddrr(0,0) = tokenizer.nextTokenAsDouble();
133 <            currProp.Ddrr(0,1) = tokenizer.nextTokenAsDouble();
134 <            currProp.Ddrr(0,2) = tokenizer.nextTokenAsDouble();
135 <            currProp.Ddrr(1,0) = tokenizer.nextTokenAsDouble();
136 <            currProp.Ddrr(1,1) = tokenizer.nextTokenAsDouble();
137 <            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();
132 >            SquareMatrix<double, 6> Xir;
133 >            Xir.setSubMatrix(0, 0, currProp.Xirtt);
134 >            Xir.setSubMatrix(0, 3, currProp.Xirrt);
135 >            Xir.setSubMatrix(3, 0, currProp.Xirtr);
136 >            Xir.setSubMatrix(3, 3, currProp.Xirrr);
137 >            CholeskyDecomposition(Xir, currProp.S);            
138  
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();
139              props.insert(std::map<std::string, HydroProp>::value_type(atomName, currProp));
140          }
141      }
# Line 203 | Line 181 | namespace oopse {
181                   omega[2] = angMom[2] /I(2, 2);
182               }
183  
184 <             //apply friction force and torque at center of diffusion
184 >             //apply friction force and torque at center of resistance
185               A = integrableObject->getA();
186               Atrans = A.transpose();
187 <             Vector3d rcd = Atrans * hydroProps_[index].cod;  
188 <             Vector3d vcd = vel + cross(omega, rcd);
189 <             vcd = A* vcd;
190 <             Vector3d frictionForce = -(hydroProps_[index].Xidtt * vcd + hydroProps_[index].Xidrt * omega);
191 <             frictionForce = Atrans*frictionForce;
192 <             integrableObject->addFrc(frictionForce);
193 <             Vector3d frictionTorque = - (hydroProps_[index].Xidtr * vcd + hydroProps_[index].Xidrr * omega);
194 <             frictionTorque = Atrans*frictionTorque;
195 <             integrableObject->addTrq(frictionTorque+ cross(rcd, frictionForce));
196 <            
197 <             //apply random force and torque at center of diffustion
198 <             Vector3d randomForce;
199 <             Vector3d randomTorque;
200 <             genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
201 <             randomForce = Atrans*randomForce;
202 <             randomTorque = Atrans* randomTorque;
203 <             integrableObject->addFrc(randomForce);            
204 <             integrableObject->addTrq(randomTorque + cross(rcd, randomForce ));
205 <            
187 >             Vector3d rcr = Atrans * hydroProps_[index].cor;  
188 >             Vector3d vcdLab = vel + cross(omega, rcr);
189 >             Vector3d vcdBody = A* vcdLab;
190 >             Vector3d frictionForceBody = -(hydroProps_[index].Xirtt * vcdBody + hydroProps_[index].Xirrt * omega);
191 >             Vector3d frictionForceLab = Atrans*frictionForceBody;
192 >             integrableObject->addFrc(frictionForceLab);
193 >             Vector3d frictionTorqueBody = - (hydroProps_[index].Xirtr * vcdBody + hydroProps_[index].Xirrr * omega);
194 >             Vector3d frictionTorqueLab = Atrans*frictionTorqueBody;
195 >             integrableObject->addTrq(frictionTorqueLab+ cross(rcr, frictionForceLab));
196 >
197 >             //apply random force and torque at center of resistance
198 >             Vector3d randomForceBody;
199 >             Vector3d randomTorqueBody;
200 >             genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
201 >             Vector3d randomForceLab = Atrans*randomForceBody;
202 >             Vector3d randomTorqueLab = Atrans* randomTorqueBody;
203 >             integrableObject->addFrc(randomForceLab);            
204 >             integrableObject->addTrq(randomTorqueLab + cross(rcr, randomForceLab ));            
205 >
206            } else {
207               //spheric atom
208 <             Vector3d frictionForce = -(hydroProps_[index].Xidtt *vel);    
208 >             Vector3d frictionForce = -(hydroProps_[index].Xirtt *vel);    
209               Vector3d randomForce;
210               Vector3d randomTorque;
211               genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
212  
235 –             //randomForce /= OOPSEConstant::energyConvert;
236 –             //randomTorque /= OOPSEConstant::energyConvert;
213               integrableObject->addFrc(frictionForce+randomForce);            
214            }
215  
# Line 249 | Line 225 | void LDForceManager::genRandomForceAndTorque(Vector3d&
225    }
226  
227   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 –    */
228  
229 <    SquareMatrix<double, 6> Xid;
265 <    SquareMatrix<double, 6> S;
229 >
230      Vector<double, 6> Z;
231      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);
232  
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 –    */
233          
234      Z[0] = randNumGen_.randNorm(0, variance);
235      Z[1] = randNumGen_.randNorm(0, variance);
# Line 289 | Line 239 | void LDForceManager::genRandomForceAndTorque(Vector3d&
239      Z[5] = randNumGen_.randNorm(0, variance);
240      
241  
242 <    generalForce = S*Z;
242 >    generalForce = hydroProps_[index].S*Z;
243      
244      force[0] = generalForce[0];
245      force[1] = generalForce[1];

Diff Legend

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