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 978 by tim, Thu May 25 17:02:00 2006 UTC vs.
Revision 1208 by xsun, Wed Jan 16 20:19:28 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 146 | 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
156 >    const RealType convertConstant = 1.439326479e4; //converts Poise angstroms
157 >                                                    // to kcal fs mol^-1 Angstrom^-1
158 >
159      Xiott *= convertConstant;
160      Xiotr *= convertConstant;
161      Xiorr *= convertConstant;
162      
159 –    
160 –    
163      Mat3x3d tmp;
164      Mat3x3d tmpInv;
165      Vector3d tmpVec;
# Line 204 | Line 206 | namespace oopse {
206      Dr6x6.getSubMatrix(0, 3, Drrt);
207      Dr6x6.getSubMatrix(3, 0, Drtr);
208      Dr6x6.getSubMatrix(3, 3, Drrr);
209 <    RealType kt = OOPSEConstant::kB * temperature ;
209 >    RealType kt = OOPSEConstant::kb * temperature ; // in kcal mol^-1
210      Drtt *= kt;
211      Drrt *= kt;
212      Drtr *= kt;
213      Drrr *= kt;
214 <    Xirtt *= OOPSEConstant::kb * temperature;
215 <    Xirtr *= OOPSEConstant::kb * temperature;
216 <    Xirrr *= OOPSEConstant::kb * temperature;
214 >    //Xirtt *= OOPSEConstant::kb * temperature;
215 >    //Xirtr *= OOPSEConstant::kb * temperature;
216 >    //Xirrr *= OOPSEConstant::kb * temperature;
217      
218 +    Mat6x6d Xi, D;
219  
220 <    cr.center = ror;
221 <    cr.Xi.setSubMatrix(0, 0, Xirtt);
222 <    cr.Xi.setSubMatrix(0, 3, Xirtr);
223 <    cr.Xi.setSubMatrix(3, 0, Xirtr);
224 <    cr.Xi.setSubMatrix(3, 3, Xirrr);
225 <    cr.D.setSubMatrix(0, 0, Drtt);
226 <    cr.D.setSubMatrix(0, 3, Drrt);
227 <    cr.D.setSubMatrix(3, 0, Drtr);
228 <    cr.D.setSubMatrix(3, 3, Drrr);    
220 >    cr->setCOR(ror);
221 >
222 >    Xi.setSubMatrix(0, 0, Xirtt);
223 >    Xi.setSubMatrix(0, 3, Xirtr);
224 >    Xi.setSubMatrix(3, 0, Xirtr);
225 >    Xi.setSubMatrix(3, 3, Xirrr);
226 >
227 >    cr->setXi(Xi);
228 >
229 >    D.setSubMatrix(0, 0, Drtt);
230 >    D.setSubMatrix(0, 3, Drrt);
231 >    D.setSubMatrix(3, 0, Drtr);
232 >    D.setSubMatrix(3, 3, Drrr);    
233 >
234 >    cr->setD(D);
235      
236      std::cout << "-----------------------------------------\n";
237      std::cout << "center of resistance :" << std::endl;
# Line 248 | Line 257 | namespace oopse {
257      return true;
258   }
259    
260 <  bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProps& cr) {
260 >  bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProp* cd) {
261      
262      int nbeads = beads.size();
263      DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads);
# Line 312 | Line 321 | namespace oopse {
321              
322          Xitt += Cij;
323          Xitr += U[i] * Cij;
324 <            //Xirr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I;            
324 >        // uncorrected here.  Volume correction is added after we assemble Xiorr
325          Xirr += -U[i] * Cij * U[j];
326        }
327      }
328 +    // add the volume correction here:
329 +    Xirr += (6.0 * viscosity * volume) * I;    
330      
331 <    const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs
331 >    const RealType convertConstant = 1.439326479e4; //converts Poise angstroms
332 >                                                    // to kcal fs mol^-1 Angstrom^-1
333      Xitt *= convertConstant;
334      Xitr *= convertConstant;
335      Xirr *= convertConstant;
336      
337 <    RealType kt = OOPSEConstant::kB * temperature;
337 >    RealType kt = OOPSEConstant::kb * temperature; // in kcal mol^-1
338      
339      Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
340      Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
# Line 398 | Line 410 | namespace oopse {
410      //Xid /= OOPSEConstant::energyConvert;
411      Xid *= OOPSEConstant::kb * temperature;
412  
413 <    cr.center = rod;
402 <    cr.D.setSubMatrix(0, 0, Ddtt);
403 <    cr.D.setSubMatrix(0, 3, Ddtr);
404 <    cr.D.setSubMatrix(3, 0, Ddtr);
405 <    cr.D.setSubMatrix(3, 3, Ddrr);
406 <    cr.Xi = Xid;
413 >    Mat6x6d Xi, D;
414  
415 +    cd->setCOR(rod);
416 +
417 +    cd->setXi(Xid);
418 +
419 +    D.setSubMatrix(0, 0, Ddtt);
420 +    D.setSubMatrix(0, 3, Ddtr);
421 +    D.setSubMatrix(3, 0, Ddtr);
422 +    D.setSubMatrix(3, 3, Ddrr);
423 +
424 +    cd->setD(D);
425 +
426      std::cout << "viscosity = " << viscosity << std::endl;
427      std::cout << "temperature = " << temperature << std::endl;
428      std::cout << "center of diffusion :" << std::endl;
429      std::cout << rod << std::endl;
430      std::cout << "diffusion tensor at center of diffusion " << std::endl;
431 <    std::cout << "translation(A^2/fs) :" << std::endl;
431 >    std::cout << "translation(A^2 / fs) :" << std::endl;
432      std::cout << Ddtt << std::endl;
433 <    std::cout << "translation-rotation(A^3/fs):" << std::endl;
433 >    std::cout << "translation-rotation(A / fs):" << std::endl;
434      std::cout << Ddtr << std::endl;
435 <    std::cout << "rotation(A^4/fs):" << std::endl;
435 >    std::cout << "rotation(fs^-1):" << std::endl;
436      std::cout << Ddrr << std::endl;
437  
438      std::cout << "resistance tensor at center of diffusion " << std::endl;
# Line 430 | Line 448 | namespace oopse {
448      Xid.getSubMatrix(3, 3, Xidrr);
449  
450      std::cout << Xidtt << std::endl;
451 <    std::cout << "rotation-translation (kcal*fs*mol^-1*Ang^-3):" << std::endl;
451 >    std::cout << "rotation-translation (kcal*fs*mol^-1*Ang^-1):" << std::endl;
452      std::cout << Xidrt << std::endl;
453 <    std::cout << "translation-rotation(kcal*fs*mol^-1*Ang^-3):" << std::endl;
453 >    std::cout << "translation-rotation(kcal*fs*mol^-1*Ang^-1):" << std::endl;
454      std::cout << Xidtr << std::endl;
455 <    std::cout << "rotation(kcal*fs*mol^-1*Ang^-4):" << std::endl;
455 >    std::cout << "rotation(kcal*fs*mol^-1):" << std::endl;
456      std::cout << Xidrr << std::endl;
457  
458      return true;

Diff Legend

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