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

Comparing trunk/src/applications/hydrodynamics/ApproximationModel.cpp (file contents):
Revision 977 by tim, Thu May 25 16:27:27 2006 UTC vs.
Revision 1322 by gezelter, Mon Nov 24 22:01:56 2008 UTC

# Line 72 | Line 72 | namespace oopse {
72    bool ApproximationModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) {
73      
74      bool ret = true;
75 <    HydroProps cr;
76 <    HydroProps cd;
75 >    HydroProp* cr = new HydroProp();
76 >    HydroProp* cd = new HydroProp();
77      calcHydroPropsAtCR(beads_, viscosity, temperature, cr);
78 <    //calcHydroPropsAtCD(beads_, viscosity, temperature, cd);
78 >    calcHydroPropsAtCD(beads_, viscosity, temperature, cd);
79      setCR(cr);
80      setCD(cd);
81 –    
81      return true;    
82    }
83    
84 <  bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProps& cr) {
84 >  bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProp* cr) {
85      
86      int nbeads = beads.size();
87      DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads);
# Line 103 | Line 102 | namespace oopse {
102                Mat3x3d tmpMat;
103                tmpMat = outProduct(Rij, Rij) / rij2;
104                RealType constant = 8.0 * NumericConstant::PI * viscosity * rij;
105 <              Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
105 >              RealType tmp1 = 1.0 + sumSigma2OverRij2/3.0;
106 >              RealType tmp2 = 1.0 - sumSigma2OverRij2;
107 >              Tij = (tmp1 * I + tmp2 * tmpMat ) / constant;
108              }else {
109                RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
110                Tij(0, 0) = constant;
# Line 144 | Line 145 | namespace oopse {
145          
146          Xiott += Cij;
147          Xiotr += U[i] * Cij;
148 <        //Xiorr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I;    
148 >        // uncorrected here.  Volume correction is added after we assemble Xiorr
149          Xiorr += -U[i] * Cij * U[j];
150        }
151      }
152 +
153 +    // add the volume correction
154 +    Xiorr += (6.0 * viscosity * volume) * I;    
155      
156 <    const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs
157 <    Xiott *= convertConstant;
158 <    Xiotr *= convertConstant;
155 <    Xiorr *= convertConstant;
156 >    Xiott *= OOPSEConstant::viscoConvert;
157 >    Xiotr *= OOPSEConstant::viscoConvert;
158 >    Xiorr *= OOPSEConstant::viscoConvert;
159      
157 –    
158 –    
160      Mat3x3d tmp;
161      Mat3x3d tmpInv;
162      Vector3d tmpVec;
# Line 202 | Line 203 | namespace oopse {
203      Dr6x6.getSubMatrix(0, 3, Drrt);
204      Dr6x6.getSubMatrix(3, 0, Drtr);
205      Dr6x6.getSubMatrix(3, 3, Drrr);
206 <    RealType kt = OOPSEConstant::kB * temperature ;
206 >    RealType kt = OOPSEConstant::kb * temperature ; // in kcal mol^-1
207      Drtt *= kt;
208      Drrt *= kt;
209      Drtr *= kt;
210      Drrr *= kt;
211 <    Xirtt *= OOPSEConstant::kb * temperature;
212 <    Xirtr *= OOPSEConstant::kb * temperature;
213 <    Xirrr *= OOPSEConstant::kb * temperature;
211 >    //Xirtt *= OOPSEConstant::kb * temperature;
212 >    //Xirtr *= OOPSEConstant::kb * temperature;
213 >    //Xirrr *= OOPSEConstant::kb * temperature;
214      
215 +    Mat6x6d Xi, D;
216  
217 <    cr.center = ror;
218 <    cr.Xi.setSubMatrix(0, 0, Xirtt);
219 <    cr.Xi.setSubMatrix(0, 3, Xirtr);
220 <    cr.Xi.setSubMatrix(3, 0, Xirtr);
221 <    cr.Xi.setSubMatrix(3, 3, Xirrr);
222 <    cr.D.setSubMatrix(0, 0, Drtt);
223 <    cr.D.setSubMatrix(0, 3, Drrt);
224 <    cr.D.setSubMatrix(3, 0, Drtr);
225 <    cr.D.setSubMatrix(3, 3, Drrr);    
217 >    cr->setCOR(ror);
218 >
219 >    Xi.setSubMatrix(0, 0, Xirtt);
220 >    Xi.setSubMatrix(0, 3, Xirtr);
221 >    Xi.setSubMatrix(3, 0, Xirtr);
222 >    Xi.setSubMatrix(3, 3, Xirrr);
223 >
224 >    cr->setXi(Xi);
225 >
226 >    D.setSubMatrix(0, 0, Drtt);
227 >    D.setSubMatrix(0, 3, Drrt);
228 >    D.setSubMatrix(3, 0, Drtr);
229 >    D.setSubMatrix(3, 3, Drrr);    
230 >
231 >    cr->setD(D);
232      
233      std::cout << "-----------------------------------------\n";
234      std::cout << "center of resistance :" << std::endl;
# Line 246 | Line 254 | namespace oopse {
254      return true;
255   }
256    
257 <  bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProps& cr) {
257 >  bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProp* cd) {
258      
259      int nbeads = beads.size();
260      DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads);
# Line 267 | Line 275 | namespace oopse {
275            Mat3x3d tmpMat;
276            tmpMat = outProduct(Rij, Rij) / rij2;
277            RealType constant = 8.0 * NumericConstant::PI * viscosity * rij;
278 <          Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
278 >          RealType tmp1 = 1.0 + sumSigma2OverRij2/3.0;
279 >          RealType tmp2 = 1.0 - sumSigma2OverRij2;
280 >          Tij = (tmp1 * I + tmp2 * tmpMat ) / constant;
281          }else {
282            RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
283            Tij(0, 0) = constant;
# Line 308 | Line 318 | namespace oopse {
318              
319          Xitt += Cij;
320          Xitr += U[i] * Cij;
321 <            //Xirr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I;            
321 >        // uncorrected here.  Volume correction is added after we assemble Xiorr
322          Xirr += -U[i] * Cij * U[j];
323        }
324      }
325 +    // add the volume correction here:
326 +    Xirr += (6.0 * viscosity * volume) * I;    
327      
328 <    const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs
329 <    Xitt *= convertConstant;
330 <    Xitr *= convertConstant;
319 <    Xirr *= convertConstant;
328 >    Xitt *= OOPSEConstant::viscoConvert;
329 >    Xitr *= OOPSEConstant::viscoConvert;
330 >    Xirr *= OOPSEConstant::viscoConvert;
331      
332 <    RealType kt = OOPSEConstant::kB * temperature;
332 >    RealType kt = OOPSEConstant::kb * temperature; // in kcal mol^-1
333      
334      Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
335      Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
# Line 394 | Line 405 | namespace oopse {
405      //Xid /= OOPSEConstant::energyConvert;
406      Xid *= OOPSEConstant::kb * temperature;
407  
408 <    cr.center = rod;
398 <    cr.D.setSubMatrix(0, 0, Ddtt);
399 <    cr.D.setSubMatrix(0, 3, Ddtr);
400 <    cr.D.setSubMatrix(3, 0, Ddtr);
401 <    cr.D.setSubMatrix(3, 3, Ddrr);
402 <    cr.Xi = Xid;
408 >    Mat6x6d Xi, D;
409  
410 +    cd->setCOR(rod);
411 +
412 +    cd->setXi(Xid);
413 +
414 +    D.setSubMatrix(0, 0, Ddtt);
415 +    D.setSubMatrix(0, 3, Ddtr);
416 +    D.setSubMatrix(3, 0, Ddtr);
417 +    D.setSubMatrix(3, 3, Ddrr);
418 +
419 +    cd->setD(D);
420 +
421      std::cout << "viscosity = " << viscosity << std::endl;
422      std::cout << "temperature = " << temperature << std::endl;
423      std::cout << "center of diffusion :" << std::endl;
424      std::cout << rod << std::endl;
425      std::cout << "diffusion tensor at center of diffusion " << std::endl;
426 <    std::cout << "translation(A^2/fs) :" << std::endl;
426 >    std::cout << "translation(A^2 / fs) :" << std::endl;
427      std::cout << Ddtt << std::endl;
428 <    std::cout << "translation-rotation(A^3/fs):" << std::endl;
428 >    std::cout << "translation-rotation(A / fs):" << std::endl;
429      std::cout << Ddtr << std::endl;
430 <    std::cout << "rotation(A^4/fs):" << std::endl;
430 >    std::cout << "rotation(fs^-1):" << std::endl;
431      std::cout << Ddrr << std::endl;
432  
433      std::cout << "resistance tensor at center of diffusion " << std::endl;
# Line 426 | Line 443 | namespace oopse {
443      Xid.getSubMatrix(3, 3, Xidrr);
444  
445      std::cout << Xidtt << std::endl;
446 <    std::cout << "rotation-translation (kcal*fs*mol^-1*Ang^-3):" << std::endl;
446 >    std::cout << "rotation-translation (kcal*fs*mol^-1*Ang^-1):" << std::endl;
447      std::cout << Xidrt << std::endl;
448 <    std::cout << "translation-rotation(kcal*fs*mol^-1*Ang^-3):" << std::endl;
448 >    std::cout << "translation-rotation(kcal*fs*mol^-1*Ang^-1):" << std::endl;
449      std::cout << Xidtr << std::endl;
450 <    std::cout << "rotation(kcal*fs*mol^-1*Ang^-4):" << std::endl;
450 >    std::cout << "rotation(kcal*fs*mol^-1):" << std::endl;
451      std::cout << Xidrr << std::endl;
452  
453      return true;

Diff Legend

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