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 905 by tim, Thu Mar 16 22:50:48 2006 UTC vs.
Revision 906 by tim, Fri Mar 17 23:20:35 2006 UTC

# Line 80 | 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  
92 –            currProp.Ddtt(0,0) = tokenizer.nextTokenAsDouble();
93 –            currProp.Ddtt(0,1) = tokenizer.nextTokenAsDouble();
94 –            currProp.Ddtt(0,2) = tokenizer.nextTokenAsDouble();
95 –            currProp.Ddtt(1,0) = tokenizer.nextTokenAsDouble();
96 –            currProp.Ddtt(1,1) = tokenizer.nextTokenAsDouble();
97 –            currProp.Ddtt(1,2) = tokenizer.nextTokenAsDouble();
98 –            currProp.Ddtt(2,0) = tokenizer.nextTokenAsDouble();
99 –            currProp.Ddtt(2,1) = tokenizer.nextTokenAsDouble();
100 –            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();
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();
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 203 | 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               Atrans = A.transpose();
221 <             Vector3d rcd = Atrans * hydroProps_[index].cod;  
222 <             Vector3d vcd = vel + cross(omega, rcd);
223 <             vcd = A* vcd;
224 <             Vector3d frictionForce = -(hydroProps_[index].Xidtt * vcd + hydroProps_[index].Xidrt * omega);
225 <             frictionForce = Atrans*frictionForce;
226 <             integrableObject->addFrc(frictionForce);
227 <             Vector3d frictionTorque = - (hydroProps_[index].Xidtr * vcd + hydroProps_[index].Xidrr * omega);
228 <             frictionTorque = Atrans*frictionTorque;
229 <             integrableObject->addTrq(frictionTorque+ cross(rcd, frictionForce));
230 <            
231 <             //apply random force and torque at center of diffustion
232 <             Vector3d randomForce;
233 <             Vector3d randomTorque;
234 <             genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
235 <             randomForce = Atrans*randomForce;
236 <             randomTorque = Atrans* randomTorque;
237 <             integrableObject->addFrc(randomForce);            
238 <             integrableObject->addTrq(randomTorque + cross(rcd, randomForce ));
239 <            
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  
235 –             //randomForce /= OOPSEConstant::energyConvert;
236 –             //randomTorque /= OOPSEConstant::energyConvert;
247               integrableObject->addFrc(frictionForce+randomForce);            
248            }
249  
# Line 249 | Line 259 | void LDForceManager::genRandomForceAndTorque(Vector3d&
259    }
260  
261   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 –    */
262  
263 <    SquareMatrix<double, 6> Xid;
265 <    SquareMatrix<double, 6> S;
263 >
264      Vector<double, 6> Z;
265      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);
266  
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 –    */
267          
268      Z[0] = randNumGen_.randNorm(0, variance);
269      Z[1] = randNumGen_.randNorm(0, variance);
# Line 289 | Line 273 | void LDForceManager::genRandomForceAndTorque(Vector3d&
273      Z[5] = randNumGen_.randNorm(0, variance);
274      
275  
276 <    generalForce = S*Z;
276 >    generalForce = hydroProps_[index].S*Z;
277      
278      force[0] = generalForce[0];
279      force[1] = generalForce[1];

Diff Legend

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