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 895 by tim, Mon Mar 13 22:42:40 2006 UTC vs.
Revision 908 by tim, Mon Mar 20 19:12:14 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 82 | 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();
117 <            currProp.Ddrr(2,0) = tokenizer.nextTokenAsDouble();
118 <            currProp.Ddrr(2,1) = tokenizer.nextTokenAsDouble();
119 <            currProp.Ddrr(2,2) = tokenizer.nextTokenAsDouble();                
120 <
121 <            currProp.Xidtt(0,0) = tokenizer.nextTokenAsDouble();
122 <            currProp.Xidtt(0,1) = tokenizer.nextTokenAsDouble();
123 <            currProp.Xidtt(0,2) = tokenizer.nextTokenAsDouble();
124 <            currProp.Xidtt(1,0) = tokenizer.nextTokenAsDouble();
125 <            currProp.Xidtt(1,1) = tokenizer.nextTokenAsDouble();
126 <            currProp.Xidtt(1,2) = tokenizer.nextTokenAsDouble();
127 <            currProp.Xidtt(2,0) = tokenizer.nextTokenAsDouble();
128 <            currProp.Xidtt(2,1) = tokenizer.nextTokenAsDouble();
129 <            currProp.Xidtt(2,2) = tokenizer.nextTokenAsDouble();
130 <
131 <            currProp.Xidrt(0,0) = tokenizer.nextTokenAsDouble();
132 <            currProp.Xidrt(0,1) = tokenizer.nextTokenAsDouble();
133 <            currProp.Xidrt(0,2) = tokenizer.nextTokenAsDouble();
134 <            currProp.Xidrt(1,0) = tokenizer.nextTokenAsDouble();
135 <            currProp.Xidrt(1,1) = tokenizer.nextTokenAsDouble();
136 <            currProp.Xidrt(1,2) = tokenizer.nextTokenAsDouble();
137 <            currProp.Xidrt(2,0) = tokenizer.nextTokenAsDouble();
138 <            currProp.Xidrt(2,1) = tokenizer.nextTokenAsDouble();
139 <            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();
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  
151 –            currProp.Xidrr(0,0) = tokenizer.nextTokenAsDouble();
152 –            currProp.Xidrr(0,1) = tokenizer.nextTokenAsDouble();
153 –            currProp.Xidrr(0,2) = tokenizer.nextTokenAsDouble();
154 –            currProp.Xidrr(1,0) = tokenizer.nextTokenAsDouble();
155 –            currProp.Xidrr(1,1) = tokenizer.nextTokenAsDouble();
156 –            currProp.Xidrr(1,2) = tokenizer.nextTokenAsDouble();
157 –            currProp.Xidrr(2,0) = tokenizer.nextTokenAsDouble();
158 –            currProp.Xidrr(2,1) = tokenizer.nextTokenAsDouble();
159 –            currProp.Xidrr(2,2) = tokenizer.nextTokenAsDouble();
139              props.insert(std::map<std::string, HydroProp>::value_type(atomName, currProp));
140          }
141      }
# Line 173 | Line 152 | namespace oopse {
152      Vector3d pos;
153      Vector3d frc;
154      Mat3x3d A;
155 +    Mat3x3d Atrans;
156      Vector3d Tb;
157      Vector3d ji;
158      double mass;
# Line 201 | 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 <             Vector3d rcd = A.transpose() * hydroProps_[index].cod;  
187 <             Vector3d vcd = vel + cross(omega, rcd);
188 <             Vector3d frictionForce = -(hydroProps_[index].Xidtt * vcd + hydroProps_[index].Xidrt * omega);
189 <             integrableObject->addFrc(frictionForce);
190 <             Vector3d frictionTorque = - (hydroProps_[index].Xidtr * vcd + hydroProps_[index].Xidrr * omega);
191 <             integrableObject->addTrq(frictionTorque);
192 <            
193 <             //apply random force and torque at center of diffustion
194 <             Vector3d randomForce;
195 <             Vector3d randomTorque;
196 <             genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
197 <             integrableObject->addFrc(randomForce);
198 <             integrableObject->addTrq(randomTorque + cross(rcd, randomForce ));
199 <            
186 >             Atrans = A.transpose();
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 +
213               integrableObject->addFrc(frictionForce+randomForce);            
214            }
215  
# Line 238 | Line 225 | void LDForceManager::genRandomForceAndTorque(Vector3d&
225    }
226  
227   void LDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, double variance) {
228 <    SquareMatrix<double, 6> Dd;
229 <    SquareMatrix<double, 6> S;
228 >
229 >
230      Vector<double, 6> Z;
231      Vector<double, 6> generalForce;
232 <    Dd.setSubMatrix(0, 0, hydroProps_[index].Ddtt);
233 <    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 <    
232 >
233 >        
234      Z[0] = randNumGen_.randNorm(0, variance);
235      Z[1] = randNumGen_.randNorm(0, variance);
236      Z[2] = randNumGen_.randNorm(0, variance);
237      Z[3] = randNumGen_.randNorm(0, variance);
238      Z[4] = randNumGen_.randNorm(0, variance);
239      Z[5] = randNumGen_.randNorm(0, variance);
240 <        
241 <    generalForce = S*Z;
240 >    
241 >
242 >    generalForce = hydroProps_[index].S*Z;
243 >    
244      force[0] = generalForce[0];
245      force[1] = generalForce[1];
246      force[2] = generalForce[2];

Diff Legend

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