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 906 by tim, Fri Mar 17 23:20:35 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 79 | Line 80 | namespace oopse {
80  
81      const unsigned int BufferSize = 65535;
82      char buffer[BufferSize];  
83 +    Mat3x3d Ddtt;
84 +    Mat3x3d Ddtr;
85 +    Mat3x3d Ddrr;
86      while (ifs.getline(buffer, BufferSize)) {
87          StringTokenizer tokenizer(buffer);
88          HydroProp currProp;
89 <        if (tokenizer.countTokens() >= 67) {
89 >        if (tokenizer.countTokens() >= 40) {
90              std::string atomName = tokenizer.nextToken();
91 <            currProp.cod[0] = tokenizer.nextTokenAsDouble();
92 <            currProp.cod[1] = tokenizer.nextTokenAsDouble();
93 <            currProp.cod[2] = tokenizer.nextTokenAsDouble();
91 >            currProp.cor[0] = tokenizer.nextTokenAsDouble();
92 >            currProp.cor[1] = tokenizer.nextTokenAsDouble();
93 >            currProp.cor[2] = tokenizer.nextTokenAsDouble();
94  
91 –            currProp.Ddtt(0,0) = tokenizer.nextTokenAsDouble();
92 –            currProp.Ddtt(0,1) = tokenizer.nextTokenAsDouble();
93 –            currProp.Ddtt(0,2) = tokenizer.nextTokenAsDouble();
94 –            currProp.Ddtt(1,0) = tokenizer.nextTokenAsDouble();
95 –            currProp.Ddtt(1,1) = tokenizer.nextTokenAsDouble();
96 –            currProp.Ddtt(1,2) = tokenizer.nextTokenAsDouble();
97 –            currProp.Ddtt(2,0) = tokenizer.nextTokenAsDouble();
98 –            currProp.Ddtt(2,1) = tokenizer.nextTokenAsDouble();
99 –            currProp.Ddtt(2,2) = tokenizer.nextTokenAsDouble();
95  
96 <            currProp.Ddtr(0,0) = tokenizer.nextTokenAsDouble();
97 <            currProp.Ddtr(0,1) = tokenizer.nextTokenAsDouble();
98 <            currProp.Ddtr(0,2) = tokenizer.nextTokenAsDouble();
99 <            currProp.Ddtr(1,0) = tokenizer.nextTokenAsDouble();
100 <            currProp.Ddtr(1,1) = tokenizer.nextTokenAsDouble();
101 <            currProp.Ddtr(1,2) = tokenizer.nextTokenAsDouble();
102 <            currProp.Ddtr(2,0) = tokenizer.nextTokenAsDouble();
103 <            currProp.Ddtr(2,1) = tokenizer.nextTokenAsDouble();
104 <            currProp.Ddtr(2,2) = tokenizer.nextTokenAsDouble();
96 >            Ddtt(0,0) = tokenizer.nextTokenAsDouble();
97 >            Ddtt(0,1) = tokenizer.nextTokenAsDouble();
98 >            Ddtt(0,2) = tokenizer.nextTokenAsDouble();
99 >            Ddtt(1,0) = tokenizer.nextTokenAsDouble();
100 >            Ddtt(1,1) = tokenizer.nextTokenAsDouble();
101 >            Ddtt(1,2) = tokenizer.nextTokenAsDouble();
102 >            Ddtt(2,0) = tokenizer.nextTokenAsDouble();
103 >            Ddtt(2,1) = tokenizer.nextTokenAsDouble();
104 >            Ddtt(2,2) = tokenizer.nextTokenAsDouble();
105  
106 <            currProp.Ddrr(0,0) = tokenizer.nextTokenAsDouble();
107 <            currProp.Ddrr(0,1) = tokenizer.nextTokenAsDouble();
108 <            currProp.Ddrr(0,2) = tokenizer.nextTokenAsDouble();
109 <            currProp.Ddrr(1,0) = tokenizer.nextTokenAsDouble();
110 <            currProp.Ddrr(1,1) = tokenizer.nextTokenAsDouble();
111 <            currProp.Ddrr(1,2) = tokenizer.nextTokenAsDouble();
112 <            currProp.Ddrr(2,0) = tokenizer.nextTokenAsDouble();
113 <            currProp.Ddrr(2,1) = tokenizer.nextTokenAsDouble();
114 <            currProp.Ddrr(2,2) = tokenizer.nextTokenAsDouble();                
106 >            Ddtr(0,0) = tokenizer.nextTokenAsDouble();
107 >            Ddtr(0,1) = tokenizer.nextTokenAsDouble();
108 >            Ddtr(0,2) = tokenizer.nextTokenAsDouble();
109 >            Ddtr(1,0) = tokenizer.nextTokenAsDouble();
110 >            Ddtr(1,1) = tokenizer.nextTokenAsDouble();
111 >            Ddtr(1,2) = tokenizer.nextTokenAsDouble();
112 >            Ddtr(2,0) = tokenizer.nextTokenAsDouble();
113 >            Ddtr(2,1) = tokenizer.nextTokenAsDouble();
114 >            Ddtr(2,2) = tokenizer.nextTokenAsDouble();
115  
116 <            currProp.Xidtt(0,0) = tokenizer.nextTokenAsDouble();
117 <            currProp.Xidtt(0,1) = tokenizer.nextTokenAsDouble();
118 <            currProp.Xidtt(0,2) = tokenizer.nextTokenAsDouble();
119 <            currProp.Xidtt(1,0) = tokenizer.nextTokenAsDouble();
120 <            currProp.Xidtt(1,1) = tokenizer.nextTokenAsDouble();
121 <            currProp.Xidtt(1,2) = tokenizer.nextTokenAsDouble();
122 <            currProp.Xidtt(2,0) = tokenizer.nextTokenAsDouble();
123 <            currProp.Xidtt(2,1) = tokenizer.nextTokenAsDouble();
124 <            currProp.Xidtt(2,2) = tokenizer.nextTokenAsDouble();
116 >            Ddrr(0,0) = tokenizer.nextTokenAsDouble();
117 >            Ddrr(0,1) = tokenizer.nextTokenAsDouble();
118 >            Ddrr(0,2) = tokenizer.nextTokenAsDouble();
119 >            Ddrr(1,0) = tokenizer.nextTokenAsDouble();
120 >            Ddrr(1,1) = tokenizer.nextTokenAsDouble();
121 >            Ddrr(1,2) = tokenizer.nextTokenAsDouble();
122 >            Ddrr(2,0) = tokenizer.nextTokenAsDouble();
123 >            Ddrr(2,1) = tokenizer.nextTokenAsDouble();
124 >            Ddrr(2,2) = tokenizer.nextTokenAsDouble();                
125  
126 <            currProp.Xidrt(0,0) = tokenizer.nextTokenAsDouble();
127 <            currProp.Xidrt(0,1) = tokenizer.nextTokenAsDouble();
128 <            currProp.Xidrt(0,2) = tokenizer.nextTokenAsDouble();
129 <            currProp.Xidrt(1,0) = tokenizer.nextTokenAsDouble();
130 <            currProp.Xidrt(1,1) = tokenizer.nextTokenAsDouble();
131 <            currProp.Xidrt(1,2) = tokenizer.nextTokenAsDouble();
132 <            currProp.Xidrt(2,0) = tokenizer.nextTokenAsDouble();
133 <            currProp.Xidrt(2,1) = tokenizer.nextTokenAsDouble();
134 <            currProp.Xidrt(2,2) = tokenizer.nextTokenAsDouble();
140 <            
141 <            currProp.Xidtr(0,0) = tokenizer.nextTokenAsDouble();
142 <            currProp.Xidtr(0,1) = tokenizer.nextTokenAsDouble();
143 <            currProp.Xidtr(0,2) = tokenizer.nextTokenAsDouble();
144 <            currProp.Xidtr(1,0) = tokenizer.nextTokenAsDouble();
145 <            currProp.Xidtr(1,1) = tokenizer.nextTokenAsDouble();
146 <            currProp.Xidtr(1,2) = tokenizer.nextTokenAsDouble();
147 <            currProp.Xidtr(2,0) = tokenizer.nextTokenAsDouble();
148 <            currProp.Xidtr(2,1) = tokenizer.nextTokenAsDouble();
149 <            currProp.Xidtr(2,2) = tokenizer.nextTokenAsDouble();
126 >            currProp.Xirtt(0,0) = tokenizer.nextTokenAsDouble();
127 >            currProp.Xirtt(0,1) = tokenizer.nextTokenAsDouble();
128 >            currProp.Xirtt(0,2) = tokenizer.nextTokenAsDouble();
129 >            currProp.Xirtt(1,0) = tokenizer.nextTokenAsDouble();
130 >            currProp.Xirtt(1,1) = tokenizer.nextTokenAsDouble();
131 >            currProp.Xirtt(1,2) = tokenizer.nextTokenAsDouble();
132 >            currProp.Xirtt(2,0) = tokenizer.nextTokenAsDouble();
133 >            currProp.Xirtt(2,1) = tokenizer.nextTokenAsDouble();
134 >            currProp.Xirtt(2,2) = tokenizer.nextTokenAsDouble();
135  
136 <            currProp.Xidrr(0,0) = tokenizer.nextTokenAsDouble();
137 <            currProp.Xidrr(0,1) = tokenizer.nextTokenAsDouble();
138 <            currProp.Xidrr(0,2) = tokenizer.nextTokenAsDouble();
139 <            currProp.Xidrr(1,0) = tokenizer.nextTokenAsDouble();
140 <            currProp.Xidrr(1,1) = tokenizer.nextTokenAsDouble();
141 <            currProp.Xidrr(1,2) = tokenizer.nextTokenAsDouble();
142 <            currProp.Xidrr(2,0) = tokenizer.nextTokenAsDouble();
143 <            currProp.Xidrr(2,1) = tokenizer.nextTokenAsDouble();
144 <            currProp.Xidrr(2,2) = tokenizer.nextTokenAsDouble();
136 >            currProp.Xirrt(0,0) = tokenizer.nextTokenAsDouble();
137 >            currProp.Xirrt(0,1) = tokenizer.nextTokenAsDouble();
138 >            currProp.Xirrt(0,2) = tokenizer.nextTokenAsDouble();
139 >            currProp.Xirrt(1,0) = tokenizer.nextTokenAsDouble();
140 >            currProp.Xirrt(1,1) = tokenizer.nextTokenAsDouble();
141 >            currProp.Xirrt(1,2) = tokenizer.nextTokenAsDouble();
142 >            currProp.Xirrt(2,0) = tokenizer.nextTokenAsDouble();
143 >            currProp.Xirrt(2,1) = tokenizer.nextTokenAsDouble();
144 >            currProp.Xirrt(2,2) = tokenizer.nextTokenAsDouble();
145 >        
146 >            currProp.Xirtr(0,0) = tokenizer.nextTokenAsDouble();
147 >            currProp.Xirtr(0,1) = tokenizer.nextTokenAsDouble();
148 >            currProp.Xirtr(0,2) = tokenizer.nextTokenAsDouble();
149 >            currProp.Xirtr(1,0) = tokenizer.nextTokenAsDouble();
150 >            currProp.Xirtr(1,1) = tokenizer.nextTokenAsDouble();
151 >            currProp.Xirtr(1,2) = tokenizer.nextTokenAsDouble();
152 >            currProp.Xirtr(2,0) = tokenizer.nextTokenAsDouble();
153 >            currProp.Xirtr(2,1) = tokenizer.nextTokenAsDouble();
154 >            currProp.Xirtr(2,2) = tokenizer.nextTokenAsDouble();
155 >
156 >            currProp.Xirrr(0,0) = tokenizer.nextTokenAsDouble();
157 >            currProp.Xirrr(0,1) = tokenizer.nextTokenAsDouble();
158 >            currProp.Xirrr(0,2) = tokenizer.nextTokenAsDouble();
159 >            currProp.Xirrr(1,0) = tokenizer.nextTokenAsDouble();
160 >            currProp.Xirrr(1,1) = tokenizer.nextTokenAsDouble();
161 >            currProp.Xirrr(1,2) = tokenizer.nextTokenAsDouble();
162 >            currProp.Xirrr(2,0) = tokenizer.nextTokenAsDouble();
163 >            currProp.Xirrr(2,1) = tokenizer.nextTokenAsDouble();
164 >            currProp.Xirrr(2,2) = tokenizer.nextTokenAsDouble();
165 >
166 >            SquareMatrix<double, 6> Xir;
167 >            Xir.setSubMatrix(0, 0, currProp.Xirtt);
168 >            Xir.setSubMatrix(0, 3, currProp.Xirrt);
169 >            Xir.setSubMatrix(3, 0, currProp.Xirtr);
170 >            Xir.setSubMatrix(3, 3, currProp.Xirrr);
171 >            CholeskyDecomposition(Xir, currProp.S);            
172 >
173              props.insert(std::map<std::string, HydroProp>::value_type(atomName, currProp));
174          }
175      }
# Line 173 | Line 186 | namespace oopse {
186      Vector3d pos;
187      Vector3d frc;
188      Mat3x3d A;
189 +    Mat3x3d Atrans;
190      Vector3d Tb;
191      Vector3d ji;
192      double mass;
# Line 201 | Line 215 | namespace oopse {
215                   omega[2] = angMom[2] /I(2, 2);
216               }
217  
218 <             //apply friction force and torque at center of diffusion
218 >             //apply friction force and torque at center of resistance
219               A = integrableObject->getA();
220 <             Vector3d rcd = A.transpose() * hydroProps_[index].cod;  
221 <             Vector3d vcd = vel + cross(omega, rcd);
222 <             Vector3d frictionForce = -(hydroProps_[index].Xidtt * vcd + hydroProps_[index].Xidrt * omega);
223 <             integrableObject->addFrc(frictionForce);
224 <             Vector3d frictionTorque = - (hydroProps_[index].Xidtr * vcd + hydroProps_[index].Xidrr * omega);
225 <             integrableObject->addTrq(frictionTorque);
226 <            
227 <             //apply random force and torque at center of diffustion
228 <             Vector3d randomForce;
229 <             Vector3d randomTorque;
230 <             genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
231 <             integrableObject->addFrc(randomForce);
232 <             integrableObject->addTrq(randomTorque + cross(rcd, randomForce ));
233 <            
220 >             Atrans = A.transpose();
221 >             Vector3d rcr = Atrans * hydroProps_[index].cor;  
222 >             Vector3d vcdLab = vel + cross(omega, rcr);
223 >             Vector3d vcdBody = A* vcdLab;
224 >             Vector3d frictionForceBody = -(hydroProps_[index].Xirtt * vcdBody + hydroProps_[index].Xirrt * omega);
225 >             Vector3d frictionForceLab = Atrans*frictionForceBody;
226 >             integrableObject->addFrc(frictionForceLab);
227 >             Vector3d frictionTorqueBody = - (hydroProps_[index].Xirtr * vcdBody + hydroProps_[index].Xirrr * omega);
228 >             Vector3d frictionTorqueLab = Atrans*frictionTorqueBody;
229 >             integrableObject->addTrq(frictionTorqueLab+ cross(rcr, frictionForceLab));
230 >
231 >             //apply random force and torque at center of resistance
232 >             Vector3d randomForceBody;
233 >             Vector3d randomTorqueBody;
234 >             genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
235 >             Vector3d randomForceLab = Atrans*randomForceBody;
236 >             Vector3d randomTorqueLab = Atrans* randomTorqueBody;
237 >             integrableObject->addFrc(randomForceLab);            
238 >             integrableObject->addTrq(randomTorqueLab + cross(rcr, randomForceLab ));            
239 >
240            } else {
241               //spheric atom
242 <             Vector3d frictionForce = -(hydroProps_[index].Xidtt *vel);    
242 >             Vector3d frictionForce = -(hydroProps_[index].Xirtt *vel);    
243               Vector3d randomForce;
244               Vector3d randomTorque;
245               genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
246 +
247               integrableObject->addFrc(frictionForce+randomForce);            
248            }
249  
# Line 238 | Line 259 | void LDForceManager::genRandomForceAndTorque(Vector3d&
259    }
260  
261   void LDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, double variance) {
262 <    SquareMatrix<double, 6> Dd;
263 <    SquareMatrix<double, 6> S;
262 >
263 >
264      Vector<double, 6> Z;
265      Vector<double, 6> generalForce;
266 <    Dd.setSubMatrix(0, 0, hydroProps_[index].Ddtt);
267 <    Dd.setSubMatrix(0, 3, hydroProps_[index].Ddtr.transpose());
247 <    Dd.setSubMatrix(3, 0, hydroProps_[index].Ddtr);
248 <    Dd.setSubMatrix(3, 3, hydroProps_[index].Ddrr);
249 <    CholeskyDecomposition(Dd, S);
250 <    
266 >
267 >        
268      Z[0] = randNumGen_.randNorm(0, variance);
269      Z[1] = randNumGen_.randNorm(0, variance);
270      Z[2] = randNumGen_.randNorm(0, variance);
271      Z[3] = randNumGen_.randNorm(0, variance);
272      Z[4] = randNumGen_.randNorm(0, variance);
273      Z[5] = randNumGen_.randNorm(0, variance);
274 <        
275 <    generalForce = S*Z;
274 >    
275 >
276 >    generalForce = hydroProps_[index].S*Z;
277 >    
278      force[0] = generalForce[0];
279      force[1] = generalForce[1];
280      force[2] = generalForce[2];

Diff Legend

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