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 900 by tim, Wed Mar 15 20:58:47 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 = kt*tmpInv;
179 >    Dotr = -kt*XirrInv * Xitr * tmpInv;
180 >    
181      tmp = Xirr - Xitr * XittInv * Xitr.transpose();    
182      tmpInv = tmp.inverse();
183      
184 <    Dorr = kt * tmpInv*1.0E16;
184 >    Dorr = kt * 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      props_.diffCenter = rod;
219 <    props_.transDiff = Ddtt;
220 <    props_.transRotDiff = Ddtr;
221 <    props_.rotDiff = Ddrr;
219 >    props_.Ddtt = Ddtt;
220 >    props_.Ddtr = Ddtr;
221 >    props_.Ddrr = Ddrr;
222  
223 <    return true;    
223 >    SquareMatrix<double, 6> Dd;
224 >    Dd.setSubMatrix(0, 0, Ddtt);
225 >    Dd.setSubMatrix(0, 3, Ddtr.transpose());
226 >    Dd.setSubMatrix(3, 0, Ddtr);
227 >    Dd.setSubMatrix(3, 3, Ddrr);    
228 >    SquareMatrix<double, 6> Xid;
229 >    invertMatrix(Dd, Xid);
230 >
231 >    //Xidtt in units of kcal*fs*mol^-1*Ang^-2
232 >    Xid *= OOPSEConstant::kb*temperature_;
233 >
234 >    Xid.getSubMatrix(0, 0, props_.Xidtt);
235 >    Xid.getSubMatrix(0, 3, props_.Xidrt);
236 >    Xid.getSubMatrix(3, 0, props_.Xidtr);
237 >    Xid.getSubMatrix(3, 3, props_.Xidrr);
238 >
239 >
240 >    std::cout << "center of diffusion :" << std::endl;
241 >    std::cout << rod << std::endl;
242 >    std::cout << "diffusion tensor at center of diffusion " << std::endl;
243 >    std::cout << "translation(A^2/fs) :" << std::endl;
244 >    std::cout << Ddtt << std::endl;
245 >    std::cout << "translation-rotation(A^3/fs):" << std::endl;
246 >    std::cout << Ddtr << std::endl;
247 >    std::cout << "rotation(A^4/fs):" << std::endl;
248 >    std::cout << Ddrr << std::endl;
249 >
250 >    std::cout << "resistance tensor at center of diffusion " << std::endl;
251 >    std::cout << "translation(kcal*fs*mol^-1*Ang^-2) :" << std::endl;
252 >    std::cout << props_.Xidtt << std::endl;
253 >    std::cout << "rotation-translation (kcal*fs*mol^-1*Ang^-3):" << std::endl;
254 >    std::cout << props_.Xidrt << std::endl;
255 >    std::cout << "translation-rotation(kcal*fs*mol^-1*Ang^-3):" << std::endl;
256 >    std::cout << props_.Xidtr << std::endl;
257 >    std::cout << "rotation(kcal*fs*mol^-1*Ang^-4):" << std::endl;
258 >    std::cout << props_.Xidrr << std::endl;
259 >    
260 >      
261   }
262  
263   void HydrodynamicsModel::writeBeads(std::ostream& os) {
# Line 245 | Line 271 | void HydrodynamicsModel::writeDiffCenterAndDiffTensor(
271   }
272  
273   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";
274  
275 <    os << "//diffusion center" << std::endl;
276 <    os << props_.diffCenter << std::endl;
275 >    os << sd_->getType() << "\t";
276 >    os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\t";
277  
278 <    os << "//translational diffusion tensor" << std::endl;
279 <    os << props_.transDiff << std::endl;
278 >    os << props_.Ddtt(0, 0) << "\t" << props_.Ddtt(0, 1) << "\t" << props_.Ddtt(0, 2) << "\t"
279 >        << props_.Ddtt(1, 0) << "\t" << props_.Ddtt(1, 1) << "\t" << props_.Ddtt(1, 2) << "\t"
280 >        << props_.Ddtt(2, 0) << "\t" << props_.Ddtt(2, 1) << "\t" << props_.Ddtt(2, 2) << "\t";
281 >    
282 >    os << props_.Ddtr(0, 0) << "\t" << props_.Ddtr(0, 1) << "\t" << props_.Ddtr(0, 2) << "\t"
283 >        << props_.Ddtr(1, 0) << "\t" << props_.Ddtr(1, 1) << "\t" << props_.Ddtr(1, 2) << "\t"
284 >        << props_.Ddtr(2, 0) << "\t" << props_.Ddtr(2, 1) << "\t" << props_.Ddtr(2, 2) << "\t";
285  
286 <    os << "//translation-rotation coupling diffusion tensor" << std::endl;
287 <    os << props_.transRotDiff << std::endl;
286 >    os << props_.Ddrr(0, 0) << "\t" << props_.Ddrr(0, 1) << "\t" << props_.Ddrr(0, 2) << "\t"
287 >        << props_.Ddrr(1, 0) << "\t" << props_.Ddrr(1, 1) << "\t" << props_.Ddrr(1, 2) << "\t"
288 >        << props_.Ddrr(2, 0) << "\t" << props_.Ddrr(2, 1) << "\t" << props_.Ddrr(2, 2) <<"\t";        
289  
290 <    os << "//rotational diffusion tensor" << std::endl;
291 <    os << props_.rotDiff << std::endl;
292 <    
265 <    /*
266 <    os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\n"
290 >    os << props_.Xidtt(0, 0) << "\t" << props_.Xidtt(0, 1) << "\t" << props_.Xidtt(0, 2) << "\t"
291 >        << props_.Xidtt(1, 0) << "\t" << props_.Xidtt(1, 1) << "\t" << props_.Xidtt(1, 2) << "\t"
292 >        << props_.Xidtt(2, 0) << "\t" << props_.Xidtt(2, 1) << "\t" << props_.Xidtt(2, 2) << "\t";
293  
294 <    os << props_.transDiff(0, 0) << "\t" << props_.transDiff(0, 1) << "\t" << props_.transDiff(0, 2) << "\t"
295 <        << props_.transDiff(1, 0) << "\t" << props_.transDiff(1, 1) << "\t" << props_.transDiff(1, 2) << "\t"
296 <        << props_.transDiff(2, 0) << "\t" << props_.transDiff(2, 1) << "\t" << props_.transDiff(2, 2) << "\n";
294 >    os << props_.Xidrt(0, 0) << "\t" << props_.Xidrt(0, 1) << "\t" << props_.Xidrt(0, 2) << "\t"
295 >        << props_.Xidrt(1, 0) << "\t" << props_.Xidrt(1, 1) << "\t" << props_.Xidrt(1, 2) << "\t"
296 >        << props_.Xidrt(2, 0) << "\t" << props_.Xidrt(2, 1) << "\t" << props_.Xidrt(2, 2) << "\t";
297      
298 <    os << props_.transRotDiff(0, 0) << "\t" << props_.transRotDiff(0, 1) << "\t" << props_.transRotDiff(0, 2) << "\t"
299 <        << props_.transRotDiff(1, 0) << "\t" << props_.transRotDiff(1, 1) << "\t" << props_.transRotDiff(1, 2) << "\t"
300 <        << props_.transRotDiff(2, 0) << "\t" << props_.transRotDiff(2, 1) << "\t" << props_.transRotDiff(2, 2) << "\t"
298 >    os << props_.Xidtr(0, 0) << "\t" << props_.Xidtr(0, 1) << "\t" << props_.Xidtr(0, 2) << "\t"
299 >        << props_.Xidtr(1, 0) << "\t" << props_.Xidtr(1, 1) << "\t" << props_.Xidtr(1, 2) << "\t"
300 >        << props_.Xidtr(2, 0) << "\t" << props_.Xidtr(2, 1) << "\t" << props_.Xidtr(2, 2) << "\t";
301  
302 <    os << props_.rotDiff(0, 0) << "\t" << props_.rotDiff(0, 1) << "\t" << props_.rotDiff(0, 2) << "\t"
303 <        << props_.rotDiff(1, 0) << "\t" << props_.rotDiff(1, 1) << "\t" << props_.rotDiff(1, 2) << "\t"
304 <        << props_.rotDiff(2, 0) << "\t" << props_.rotDiff(2, 1) << "\t" << props_.rotDiff(2, 2) << ";"
305 <        << std::endl;
280 <    */
302 >    os << props_.Xidrr(0, 0) << "\t" << props_.Xidrr(0, 1) << "\t" << props_.Xidrr(0, 2) << "\t"
303 >        << props_.Xidrr(1, 0) << "\t" << props_.Xidrr(1, 1) << "\t" << props_.Xidrr(1, 2) << "\t"
304 >        << props_.Xidrr(2, 0) << "\t" << props_.Xidrr(2, 1) << "\t" << props_.Xidrr(2, 2) << std::endl;
305 >    
306   }
307  
308   }

Diff Legend

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