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 895 by tim, Mon Mar 13 22:42:40 2006 UTC vs.
Revision 904 by tim, Thu Mar 16 22:50:48 2006 UTC

# Line 41 | Line 41
41   #include <fstream>
42   #include "integrators/LDForceManager.hpp"
43   #include "math/CholeskyDecomposition.hpp"
44 + #include "utils/OOPSEConstant.hpp"
45   namespace oopse {
46  
47    LDForceManager::LDForceManager(SimInfo* info) : ForceManager(info){
# Line 68 | Line 69 | namespace oopse {
69              
70             }
71      }
72 <    variance_ = 2.0*simParams->getDt();
72 >    variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
73    }
74    std::map<std::string, HydroProp> LDForceManager::parseFrictionFile(const std::string& filename) {
75      std::map<std::string, HydroProp> props;
# Line 173 | Line 174 | namespace oopse {
174      Vector3d pos;
175      Vector3d frc;
176      Mat3x3d A;
177 +    Mat3x3d Atrans;
178      Vector3d Tb;
179      Vector3d ji;
180      double mass;
# Line 203 | Line 205 | namespace oopse {
205  
206               //apply friction force and torque at center of diffusion
207               A = integrableObject->getA();
208 <             Vector3d rcd = A.transpose() * hydroProps_[index].cod;  
208 >             Atrans = A.transpose();
209 >             Vector3d rcd = Atrans * hydroProps_[index].cod;  
210               Vector3d vcd = vel + cross(omega, rcd);
211 +             vcd = A* vcd;
212               Vector3d frictionForce = -(hydroProps_[index].Xidtt * vcd + hydroProps_[index].Xidrt * omega);
213 +             frictionForce = Atrans*frictionForce;
214               integrableObject->addFrc(frictionForce);
215               Vector3d frictionTorque = - (hydroProps_[index].Xidtr * vcd + hydroProps_[index].Xidrr * omega);
216 <             integrableObject->addTrq(frictionTorque);
216 >             frictionTorque = Atrans*frictionTorque;
217 >             integrableObject->addTrq(frictionTorque+ cross(rcd, frictionForce));
218              
219               //apply random force and torque at center of diffustion
220               Vector3d randomForce;
221               Vector3d randomTorque;
222               genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
223 <             integrableObject->addFrc(randomForce);
223 >             randomForce = Atrans*randomForce;
224 >             randomTorque = Atrans* randomTorque;
225 >             integrableObject->addFrc(randomForce);            
226               integrableObject->addTrq(randomTorque + cross(rcd, randomForce ));
227              
228            } else {
# Line 223 | Line 231 | namespace oopse {
231               Vector3d randomForce;
232               Vector3d randomTorque;
233               genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
234 +
235 +             //randomForce /= OOPSEConstant::energyConvert;
236 +             //randomTorque /= OOPSEConstant::energyConvert;
237               integrableObject->addFrc(frictionForce+randomForce);            
238            }
239  
# Line 238 | Line 249 | void LDForceManager::genRandomForceAndTorque(Vector3d&
249    }
250  
251   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;
# Line 247 | Line 259 | void LDForceManager::genRandomForceAndTorque(Vector3d&
259      Dd.setSubMatrix(3, 0, hydroProps_[index].Ddtr);
260      Dd.setSubMatrix(3, 3, hydroProps_[index].Ddrr);
261      CholeskyDecomposition(Dd, S);
262 <    
262 >    */
263 >
264 >    SquareMatrix<double, 6> Xid;
265 >    SquareMatrix<double, 6> S;
266 >    Vector<double, 6> Z;
267 >    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);
273 >
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 >    */
283 >        
284      Z[0] = randNumGen_.randNorm(0, variance);
285      Z[1] = randNumGen_.randNorm(0, variance);
286      Z[2] = randNumGen_.randNorm(0, variance);
287      Z[3] = randNumGen_.randNorm(0, variance);
288      Z[4] = randNumGen_.randNorm(0, variance);
289      Z[5] = randNumGen_.randNorm(0, variance);
290 <        
290 >    
291 >
292      generalForce = S*Z;
293 +    
294      force[0] = generalForce[0];
295      force[1] = generalForce[1];
296      force[2] = generalForce[2];

Diff Legend

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