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

Comparing trunk/src/applications/hydrodynamics/AnalyticalModel.cpp (file contents):
Revision 906 by tim, Fri Mar 17 23:20:35 2006 UTC vs.
Revision 956 by gezelter, Tue May 16 02:06:37 2006 UTC

# Line 38 | Line 38
38   * University of Notre Dame has been advised of the possibility of
39   * such damages.
40   */
41 < #include "applications/hydrodynamics/AnalyticalModel.hpp"
42 < #include "applications/hydrodynamics/Spheric.hpp"
43 < #include "applications/hydrodynamics/Ellipsoid.hpp"
44 < #include "applications/hydrodynamics/CompositeShape.hpp"
41 > #include "hydrodynamics/AnalyticalModel.hpp"
42 > #include "hydrodynamics/Sphere.hpp"
43 > #include "hydrodynamics/Ellipsoid.hpp"
44 > #include "hydrodynamics/CompositeShape.hpp"
45   #include "math/LU.hpp"
46   namespace oopse {
47 < bool AnalyticalModel::calcHydroProps(Spheric* spheric, double viscosity, double temperature) {
48 <
49 <    double radius = spheric->getRadius();
47 >  
48 >  bool AnalyticalModel::calcHydroProps(Sphere* sphere, double viscosity, double temperature) {
49 >    
50 >    double radius = sphere->getRadius();
51      HydroProps props;
52      props.center =V3Zero;
53      double Xitt  = 6.0 * NumericConstant::PI * viscosity * radius;
# Line 65 | Line 66 | bool AnalyticalModel::calcHydroProps(Spheric* spheric,
66      double kt = OOPSEConstant::kB * temperature;
67      props.D *= kt;
68      props.Xi *= OOPSEConstant::kb * temperature;
69 <
69 >    
70      setCR(props);
71      setCD(props);
71 –
72 –    return true;
72      
73 +    return true;    
74   }
75  
76 < /**
77 < * calculate the ratio of friction coeffiction constant between ellipsoid and spheric
78 < * with same volume.
79 < * @param m
80 < * @param n
81 < * @note
82 < * Reference:
83 < *
84 < * (1) Victor A. Bloomfield, On-Line Biophysics Textbook, Volume: Separations and Hydrodynamics
85 < * Chapter 1,Survey of Biomolecular Hydrodynamics
86 < * http://www.biophysics.org/education/vbloomfield.pdf
87 < * (2) F. Perrin , J. Phys. Radium, [7] 5, 497-511, 1934
88 < * (3) F. Perrin, J. Phys. Radium, [7] 7, 1-11, 1936
89 < */        
90 < bool AnalyticalModel::calcHydroProps(Ellipsoid* ellipsoid, double viscosity, double temperature) {
91 <    double ft;
92 <    double fra;
93 <    double frb;
94 <    double a = ellipsoid->getA();
95 <    double b = ellipsoid->getB();
96 <    double q = a/b; //?
97 <    if (q > 1.0) {//prolate
98 <        ft = sqrt(1-q*q)/(pow(q, 2.0/3.0)*log((1 + sqrt(1-q*q))/q));
99 <        fra = 4*(1-q*q)/(3*(2 - 2*pow(q, 4.0/3.0)/ft)); //not sure
100 <        frb = 4*(1-q*q*q*q) /(3*q*q*(2*pow(q, -2.0/3.0)*(2-q*q)/ft-2));
101 <    } else {//oblate
102 <        ft = sqrt(1-q*q)/(pow(q, 2.0/3.0)*atan(sqrt(q*q-1)));
103 <        fra = 4*(1-q*q)/(3*(2 - 2*pow(q, 4.0/3.0)/ft)); //not sure
104 <        frb = 4*(1-q*q*q*q) /(3*q*q*(2*pow(q, -2.0/3.0)*(2-q*q)/ft-2));
76 >  /**
77 >   * Reference:
78 >   * (2) F. Perrin , J. Phys. Radium, [7] 5, 497-511, 1934
79 >   * (3) F. Perrin, J. Phys. Radium, [7] 7, 1-11, 1936
80 >   */        
81 >  bool AnalyticalModel::calcHydroProps(Ellipsoid* ellipsoid, double viscosity, double temperature) {
82 >    
83 >    double rMajor = ellipsoid->getRMajor();
84 >    double rMinor = ellipsoid->getRMinor();
85 >    
86 >    double a = rMinor;
87 >    double b = rMajor;
88 >    double a2 = a * a;
89 >    double b2 = b* b;
90 >    
91 >    double p = a /b;
92 >    double S;
93 >    if (p > 1.0) { //prolate
94 >      S = 2.0/sqrt(a2 - b2) * log((a + sqrt(a2-b2))/b);
95 >    } else { //oblate
96 >      S = 2.0/sqrt(b2 - a2) * atan(sqrt(b2-a2)/a);
97      }
98 <                    
99 <    double radius = pow(a*a*b, 1.0/3.0);
98 >    
99 >    double P = 1.0/(a2 - b2) * (S - 2.0/a);
100 >    double Q = 0.5/(a2-b2) * (2.0*a/b2 - S);
101 >    
102 >    double transMinor = 16.0 * NumericConstant::PI * viscosity * (a2 - b2) /((2.0*a2-b2)*S -2.0*a);
103 >    double transMajor = 32.0 * NumericConstant::PI * viscosity * (a2 - b2) /((2.0*a2-3.0*b2)*S +2.0*a);
104 >    double rotMinor = 32.0/3.0 * NumericConstant::PI * viscosity *(a2 - b2) * b2 /(2.0*a -b2*S);
105 >    double rotMajor = 32.0/3.0 * NumericConstant::PI * viscosity *(a2*a2 - b2*b2)/((2.0*a2-b2)*S-2.0*a);
106 >    
107 >    
108      HydroProps props;
109 –    double Xitt  = 6.0 * NumericConstant::PI * viscosity * radius;
110 –    double Xirr = 8.0 * NumericConstant::PI * viscosity * radius * radius * radius;
111 –    props.Xi(0, 0) = Xitt;
112 –    props.Xi(1, 1) = Xitt;
113 –    props.Xi(2, 2) = Xitt;
114 –    props.Xi(3, 3) = Xirr;
115 –    props.Xi(4, 4) = Xirr;
116 –    props.Xi(5, 5) = Xirr;
109      
110 +    props.Xi(0,0) = transMajor;
111 +    props.Xi(1,1) = transMajor;
112 +    props.Xi(2,2) = transMinor;
113 +    props.Xi(3,3) = rotMajor;
114 +    props.Xi(4,4) = rotMajor;
115 +    props.Xi(5,5) = rotMinor;
116 +    
117      const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
118      props.Xi *= convertConstant;    
120 –    props.Xi(0,0) *= ft;
121 –    props.Xi(1,1) *= ft;
122 –    props.Xi(2,2) *= ft;
123 –    props.Xi(3,3) *= fra;
124 –    props.Xi(4,4) *= fra;
125 –    props.Xi(5,5) *= frb;
119      
120      Mat6x6d XiCopy = props.Xi;
128 –    XiCopy /= OOPSEConstant::kb * temperature;
121      invertMatrix(XiCopy, props.D);
122      double kt = OOPSEConstant::kB * temperature;
123      props.D *= kt;
124 <
124 >    props.Xi *= OOPSEConstant::kb * temperature;
125 >    
126      setCR(props);
127      setCD(props);
128 <
128 >    
129      return true;
130 < }
131 <
132 < bool AnalyticalModel::calcHydroProps(CompositeShape* compositexShape, double viscosity, double temperature) {
130 >  }
131 >  
132 >  bool AnalyticalModel::calcHydroProps(CompositeShape* compositeShape, double viscosity, double temperature) {
133      return false;
134 +  }
135 +  
136 +  void AnalyticalModel::writeBeads(std::ostream& os) {
137 +    os << "1\n";
138 +    os << "Generated by Hydro\n";
139 +    Vector3d pos = sd_->getPos();
140 +    os << sd_->getType() << "\t" << pos[0] << "\t" << pos[1] << "\t" << pos[2] << std::endl;
141 +  }  
142   }
142 –        
143 –
144 –
145 – }

Diff Legend

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