ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/hydrodynamics/HydrodynamicsModel.cpp
(Generate patch)

Comparing trunk/src/applications/hydrodynamics/HydrodynamicsModel.cpp (file contents):
Revision 893 by tim, Fri Feb 24 21:17:05 2006 UTC vs.
Revision 904 by tim, Thu Mar 16 22:50:48 2006 UTC

# Line 77 | Line 77 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
77          std::cout << "can not create beads" << std::endl;
78          return false;
79      }
80 <    
80 >
81 >    //calcResistanceTensor();
82 >    calcDiffusionTensor();    
83 >    return true;    
84 > }
85 >
86 > void HydrodynamicsModel::calcResistanceTensor() {
87 > }
88 >
89 > void HydrodynamicsModel::calcDiffusionTensor() {
90      int nbeads = beads_.size();
91      DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
92      DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
# Line 105 | Line 114 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
114                  Tij(2, 2) = constant;
115              }
116              B.setSubMatrix(i*3, j*3, Tij);
108 –            std::cout << Tij << std::endl;
117          }
118      }
119  
112 –    std::cout << "B=\n"
113 –                  << B << std::endl;
120      //invert B Matrix
121      invertMatrix(B, C);
122  
117 –    std::cout << "C=\n"
118 –                  << C << std::endl;
119 –
123      //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
124      std::vector<Mat3x3d> U;
125      for (int i = 0; i < nbeads; ++i) {
# Line 134 | Line 137 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
137  
138      double volume = 0.0;
139      for (std::vector<BeadParam>::iterator iter = beads_.begin(); iter != beads_.end(); ++iter) {
140 <        volume = 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
140 >        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
141      }
142          
143      for (std::size_t i = 0; i < nbeads; ++i) {
# Line 144 | Line 147 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
147              
148              Xitt += Cij;
149              Xitr += U[i] * Cij;
150 <            Xirr += -U[i] * Cij * U[j];            
148 <            //Xirr += -U[i] * Cij * U[j] + (0.166*6 * viscosity_ * volume) * I;            
150 >            Xirr += -U[i] * Cij * U[j] + (6 * viscosity_ * volume) * I;            
151          }
152      }
153  
154 <    //invert Xi to get Diffusion Tensor at arbitrary origin O
155 <    RectMatrix<double, 6, 6> Xi;    
156 <    RectMatrix<double, 6, 6> Do;
157 <    Xi.setSubMatrix(0, 0, Xitt);
156 <    Xi.setSubMatrix(0, 3, Xitr.transpose());
157 <    Xi.setSubMatrix(3, 0, Xitr);
158 <    Xi.setSubMatrix(3, 3, Xirr);
159 <    //invertMatrix(Xi, Do);
160 <    double kt = OOPSEConstant::kB * temperature_ * 1.66E-2;
161 <    //Do *= kt;    
154 >    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
155 >    Xitt *= convertConstant;
156 >    Xitr *= convertConstant;
157 >    Xirr *= convertConstant;
158  
159 +    double kt = OOPSEConstant::kB * temperature_;
160  
161      Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
162      Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
# Line 170 | Line 167 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
167      Mat3x3d XittInv(0.0);
168      XittInv = Xitt.inverse();
169      
173 –    //Xirr may not be inverted,if it one of the diagonal element is zero, for example
174 –    //( a11 a12 0)
175 –    //( a21 a22 0)
176 –    //( 0    0    0)
170      Mat3x3d XirrInv;
171      XirrInv = Xirr.inverse();
172  
# Line 182 | Line 175 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
175      tmp = Xitt - Xitr.transpose() * XirrInv * Xitr;
176      tmpInv = tmp.inverse();
177  
178 <    Dott = kt * tmpInv;
179 <    Dotr = -kt*XirrInv * Xitr * tmpInv* 1.0E8;
180 <
178 >    Dott = tmpInv;
179 >    Dotr = -XirrInv * Xitr * tmpInv;
180 >    
181      tmp = Xirr - Xitr * XittInv * Xitr.transpose();    
182      tmpInv = tmp.inverse();
183      
184 <    Dorr = kt * tmpInv*1.0E16;
184 >    Dorr = tmpInv;
185  
193 –    //Do.getSubMatrix(0, 0 , Dott);
194 –    //Do.getSubMatrix(3, 0, Dotr);
195 –    //Do.getSubMatrix(3, 3, Dorr);
196 –
186      //calculate center of diffusion
187      tmp(0, 0) = Dorr(1, 1) + Dorr(2, 2);
188      tmp(0, 1) = - Dorr(0, 1);
# Line 225 | Line 214 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
214      Ddtt = Dott - Uod * Dorr * Uod + Dotr.transpose() * Uod - Uod * Dotr;
215      Ddrr = Dorr;
216      Ddtr = Dotr + Dorr * Uod;
217 <    
217 >
218 >    SquareMatrix<double, 6> Dd;
219 >    Dd.setSubMatrix(0, 0, Ddtt);
220 >    Dd.setSubMatrix(0, 3, Ddtr.transpose());
221 >    Dd.setSubMatrix(3, 0, Ddtr);
222 >    Dd.setSubMatrix(3, 3, Ddrr);    
223 >    SquareMatrix<double, 6> Xid;
224 >    Ddtt *= kt;
225 >    Ddtr *=kt;
226 >    Ddrr *= kt;
227 >    invertMatrix(Dd, Xid);
228 >
229 >
230 >
231 >    //Xidtt in units of kcal*fs*mol^-1*Ang^-2
232 >    //Xid /= OOPSEConstant::energyConvert;
233 >    Xid *= OOPSEConstant::kb * temperature_;
234      props_.diffCenter = rod;
235 <    props_.transDiff = Ddtt;
236 <    props_.transRotDiff = Ddtr;
237 <    props_.rotDiff = Ddrr;
235 >    props_.Ddtt = Ddtt;
236 >    props_.Ddtr = Ddtr;
237 >    props_.Ddrr = Ddrr;
238 >    Xid.getSubMatrix(0, 0, props_.Xidtt);
239 >    Xid.getSubMatrix(0, 3, props_.Xidrt);
240 >    Xid.getSubMatrix(3, 0, props_.Xidtr);
241 >    Xid.getSubMatrix(3, 3, props_.Xidrr);
242  
243 <    return true;    
243 >
244 >    std::cout << "viscosity = " << viscosity_ << std::endl;
245 >    std::cout << "temperature = " << temperature_ << std::endl;
246 >    std::cout << "center of diffusion :" << std::endl;
247 >    std::cout << rod << std::endl;
248 >    std::cout << "diffusion tensor at center of diffusion " << std::endl;
249 >    std::cout << "translation(A^2/fs) :" << std::endl;
250 >    std::cout << Ddtt << std::endl;
251 >    std::cout << "translation-rotation(A^3/fs):" << std::endl;
252 >    std::cout << Ddtr << std::endl;
253 >    std::cout << "rotation(A^4/fs):" << std::endl;
254 >    std::cout << Ddrr << std::endl;
255 >
256 >    std::cout << "resistance tensor at center of diffusion " << std::endl;
257 >    std::cout << "translation(kcal*fs*mol^-1*Ang^-2) :" << std::endl;
258 >    std::cout << props_.Xidtt << std::endl;
259 >    std::cout << "rotation-translation (kcal*fs*mol^-1*Ang^-3):" << std::endl;
260 >    std::cout << props_.Xidrt << std::endl;
261 >    std::cout << "translation-rotation(kcal*fs*mol^-1*Ang^-3):" << std::endl;
262 >    std::cout << props_.Xidtr << std::endl;
263 >    std::cout << "rotation(kcal*fs*mol^-1*Ang^-4):" << std::endl;
264 >    std::cout << props_.Xidrr << std::endl;
265 >    
266 >      
267   }
268  
269   void HydrodynamicsModel::writeBeads(std::ostream& os) {
# Line 245 | Line 277 | void HydrodynamicsModel::writeDiffCenterAndDiffTensor(
277   }
278  
279   void HydrodynamicsModel::writeDiffCenterAndDiffTensor(std::ostream& os) {
248 –    os << "//viscosity = " << viscosity_ << std::endl;
249 –    os << "//temperature = " << temperature_<< std::endl;
250 –    std::vector<BeadParam>::iterator iter;
251 –    os << sd_->getType() << "\n";
280  
281 <    os << "//diffusion center" << std::endl;
282 <    os << props_.diffCenter << std::endl;
281 >    os << sd_->getType() << "\t";
282 >    os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\t";
283  
284 <    os << "//translational diffusion tensor" << std::endl;
285 <    os << props_.transDiff << std::endl;
284 >    os << props_.Ddtt(0, 0) << "\t" << props_.Ddtt(0, 1) << "\t" << props_.Ddtt(0, 2) << "\t"
285 >        << props_.Ddtt(1, 0) << "\t" << props_.Ddtt(1, 1) << "\t" << props_.Ddtt(1, 2) << "\t"
286 >        << props_.Ddtt(2, 0) << "\t" << props_.Ddtt(2, 1) << "\t" << props_.Ddtt(2, 2) << "\t";
287 >    
288 >    os << props_.Ddtr(0, 0) << "\t" << props_.Ddtr(0, 1) << "\t" << props_.Ddtr(0, 2) << "\t"
289 >        << props_.Ddtr(1, 0) << "\t" << props_.Ddtr(1, 1) << "\t" << props_.Ddtr(1, 2) << "\t"
290 >        << props_.Ddtr(2, 0) << "\t" << props_.Ddtr(2, 1) << "\t" << props_.Ddtr(2, 2) << "\t";
291  
292 <    os << "//translation-rotation coupling diffusion tensor" << std::endl;
293 <    os << props_.transRotDiff << std::endl;
292 >    os << props_.Ddrr(0, 0) << "\t" << props_.Ddrr(0, 1) << "\t" << props_.Ddrr(0, 2) << "\t"
293 >        << props_.Ddrr(1, 0) << "\t" << props_.Ddrr(1, 1) << "\t" << props_.Ddrr(1, 2) << "\t"
294 >        << props_.Ddrr(2, 0) << "\t" << props_.Ddrr(2, 1) << "\t" << props_.Ddrr(2, 2) <<"\t";        
295  
296 <    os << "//rotational diffusion tensor" << std::endl;
297 <    os << props_.rotDiff << std::endl;
298 <    
265 <    /*
266 <    os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\n"
296 >    os << props_.Xidtt(0, 0) << "\t" << props_.Xidtt(0, 1) << "\t" << props_.Xidtt(0, 2) << "\t"
297 >        << props_.Xidtt(1, 0) << "\t" << props_.Xidtt(1, 1) << "\t" << props_.Xidtt(1, 2) << "\t"
298 >        << props_.Xidtt(2, 0) << "\t" << props_.Xidtt(2, 1) << "\t" << props_.Xidtt(2, 2) << "\t";
299  
300 <    os << props_.transDiff(0, 0) << "\t" << props_.transDiff(0, 1) << "\t" << props_.transDiff(0, 2) << "\t"
301 <        << props_.transDiff(1, 0) << "\t" << props_.transDiff(1, 1) << "\t" << props_.transDiff(1, 2) << "\t"
302 <        << props_.transDiff(2, 0) << "\t" << props_.transDiff(2, 1) << "\t" << props_.transDiff(2, 2) << "\n";
300 >    os << props_.Xidrt(0, 0) << "\t" << props_.Xidrt(0, 1) << "\t" << props_.Xidrt(0, 2) << "\t"
301 >        << props_.Xidrt(1, 0) << "\t" << props_.Xidrt(1, 1) << "\t" << props_.Xidrt(1, 2) << "\t"
302 >        << props_.Xidrt(2, 0) << "\t" << props_.Xidrt(2, 1) << "\t" << props_.Xidrt(2, 2) << "\t";
303      
304 <    os << props_.transRotDiff(0, 0) << "\t" << props_.transRotDiff(0, 1) << "\t" << props_.transRotDiff(0, 2) << "\t"
305 <        << props_.transRotDiff(1, 0) << "\t" << props_.transRotDiff(1, 1) << "\t" << props_.transRotDiff(1, 2) << "\t"
306 <        << props_.transRotDiff(2, 0) << "\t" << props_.transRotDiff(2, 1) << "\t" << props_.transRotDiff(2, 2) << "\t"
304 >    os << props_.Xidtr(0, 0) << "\t" << props_.Xidtr(0, 1) << "\t" << props_.Xidtr(0, 2) << "\t"
305 >        << props_.Xidtr(1, 0) << "\t" << props_.Xidtr(1, 1) << "\t" << props_.Xidtr(1, 2) << "\t"
306 >        << props_.Xidtr(2, 0) << "\t" << props_.Xidtr(2, 1) << "\t" << props_.Xidtr(2, 2) << "\t";
307  
308 <    os << props_.rotDiff(0, 0) << "\t" << props_.rotDiff(0, 1) << "\t" << props_.rotDiff(0, 2) << "\t"
309 <        << props_.rotDiff(1, 0) << "\t" << props_.rotDiff(1, 1) << "\t" << props_.rotDiff(1, 2) << "\t"
310 <        << props_.rotDiff(2, 0) << "\t" << props_.rotDiff(2, 1) << "\t" << props_.rotDiff(2, 2) << ";"
311 <        << std::endl;
280 <    */
308 >    os << props_.Xidrr(0, 0) << "\t" << props_.Xidrr(0, 1) << "\t" << props_.Xidrr(0, 2) << "\t"
309 >        << props_.Xidrr(1, 0) << "\t" << props_.Xidrr(1, 1) << "\t" << props_.Xidrr(1, 2) << "\t"
310 >        << props_.Xidrr(2, 0) << "\t" << props_.Xidrr(2, 1) << "\t" << props_.Xidrr(2, 2) << std::endl;
311 >    
312   }
313  
314   }

Diff Legend

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