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 904 by tim, Thu Mar 16 22:50:48 2006 UTC vs.
Revision 981 by gezelter, Mon Jun 5 18:24:45 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/HydrodynamicsModel.hpp"
42 + #include "hydrodynamics/Shape.hpp"
43 + #include "hydrodynamics/Sphere.hpp"
44 + #include "hydrodynamics/Ellipsoid.hpp"
45 + #include "applications/hydrodynamics/CompositeShape.hpp"
46  
42 – #include "applications/hydrodynamics/HydrodynamicsModel.hpp"
43 – #include "math/LU.hpp"
44 – #include "math/DynamicRectMatrix.hpp"
45 – #include "math/SquareMatrix3.hpp"
46 – #include "utils/OOPSEConstant.hpp"
47   namespace oopse {
48 < /**
49 < * Reference:
50 < * Beatriz Carrasco and Jose Gracia de la Torre, Hydrodynamic Properties of Rigid Particles:
51 < * Comparison of Different Modeling and Computational Procedures.
52 < * Biophysical Journal, 75(6), 3044, 1999
53 < */
48 >  
49 >  bool HydrodynamicsModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) {
50 >    return false;
51 >  }
52 >  
53 >  void HydrodynamicsModel::writeHydroProps(std::ostream& os) {
54  
55 < HydrodynamicsModel::HydrodynamicsModel(StuntDouble* sd, const DynamicProperty& extraParams) : sd_(sd){
56 <    DynamicProperty::const_iterator iter;
57 <
58 <    iter = extraParams.find("Viscosity");
59 <    if (iter != extraParams.end()) {
60 <        boost::any param = iter->second;
61 <        viscosity_ = boost::any_cast<double>(param);
62 <    }else {
63 <        std::cout << "HydrodynamicsModel Error\n" ;
64 <    }
65 <
66 <    iter = extraParams.find("Temperature");
67 <    if (iter != extraParams.end()) {
68 <        boost::any param = iter->second;
69 <        temperature_ = boost::any_cast<double>(param);
70 <    }else {
71 <        std::cout << "HydrodynamicsModel Error\n" ;
72 <    }    
73 < }
74 <
75 < bool HydrodynamicsModel::calcHydrodyanmicsProps() {
76 <    if (!createBeads(beads_)) {
77 <        std::cout << "can not create beads" << std::endl;
78 <        return false;
79 <    }
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);
93 <    Mat3x3d I;
94 <    I(0, 0) = 1.0;
95 <    I(1, 1) = 1.0;
96 <    I(2, 2) = 1.0;
55 >    Vector3d center;
56 >    Mat6x6d Xi, D;
57      
58 <    for (std::size_t i = 0; i < nbeads; ++i) {
59 <        for (std::size_t j = 0; j < nbeads; ++j) {
60 <            Mat3x3d Tij;
101 <            if (i != j ) {
102 <                Vector3d Rij = beads_[i].pos - beads_[j].pos;
103 <                double rij = Rij.length();
104 <                double rij2 = rij * rij;
105 <                double sumSigma2OverRij2 = ((beads_[i].radius*beads_[i].radius) + (beads_[i].radius*beads_[i].radius)) / rij2;                
106 <                Mat3x3d tmpMat;
107 <                tmpMat = outProduct(Rij, Rij) / rij2;
108 <                double constant = 8.0 * NumericConstant::PI * viscosity_ * rij;
109 <                Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
110 <            }else {
111 <                double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity_ * beads_[i].radius);
112 <                Tij(0, 0) = constant;
113 <                Tij(1, 1) = constant;
114 <                Tij(2, 2) = constant;
115 <            }
116 <            B.setSubMatrix(i*3, j*3, Tij);
117 <        }
118 <    }
58 >    os << sd_->getType() << "\t";
59 >    
60 >    //center of resistance
61  
62 <    //invert B Matrix
121 <    invertMatrix(B, C);
62 >    center = cr_->getCOR();
63  
64 <    //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) {
126 <        Mat3x3d currU;
127 <        currU.setupSkewMat(beads_[i].pos);
128 <        U.push_back(currU);
129 <    }
64 >    os << center[0] <<  "\t" << center[1] <<  "\t" << center[2] <<  "\t";
65      
66 <    //calculate Xi matrix at arbitrary origin O
67 <    Mat3x3d Xitt;
133 <    Mat3x3d Xirr;
134 <    Mat3x3d Xitr;
66 >    //resistance tensor at center of resistance
67 >    //translation
68  
69 <    //calculate the total volume
69 >    Xi = cr_->getXi();
70  
71 <    double volume = 0.0;
72 <    for (std::vector<BeadParam>::iterator iter = beads_.begin(); iter != beads_.end(); ++iter) {
73 <        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
141 <    }
142 <        
143 <    for (std::size_t i = 0; i < nbeads; ++i) {
144 <        for (std::size_t j = 0; j < nbeads; ++j) {
145 <            Mat3x3d Cij;
146 <            C.getSubMatrix(i*3, j*3, Cij);
147 <            
148 <            Xitt += Cij;
149 <            Xitr += U[i] * Cij;
150 <            Xirr += -U[i] * Cij * U[j] + (6 * viscosity_ * volume) * I;            
151 <        }
152 <    }
153 <
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
163 <    Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O
164 <
165 <    const static Mat3x3d zeroMat(0.0);
71 >    os << Xi(0, 0) <<  "\t" << Xi(0, 1) <<  "\t" << Xi(0, 2) <<  "\t"
72 >       << Xi(1, 0) <<  "\t" << Xi(1, 1) <<  "\t" << Xi(1, 2) <<  "\t"
73 >       << Xi(2, 0) <<  "\t" << Xi(2, 1) <<  "\t" << Xi(2, 2) <<  "\t";
74      
75 <    Mat3x3d XittInv(0.0);
76 <    XittInv = Xitt.inverse();
75 >    //rotation-translation
76 >    os << Xi(0, 3) <<  "\t" << Xi(0, 4) <<  "\t" << Xi(0, 5) <<  "\t"
77 >       << Xi(1, 3) <<  "\t" << Xi(1, 4) <<  "\t" << Xi(1, 5) <<  "\t"
78 >       << Xi(2, 3) <<  "\t" << Xi(2, 4) <<  "\t" << Xi(2, 5) <<  "\t";
79      
80 <    Mat3x3d XirrInv;
81 <    XirrInv = Xirr.inverse();
82 <
83 <    Mat3x3d tmp;
174 <    Mat3x3d tmpInv;
175 <    tmp = Xitt - Xitr.transpose() * XirrInv * Xitr;
176 <    tmpInv = tmp.inverse();
177 <
178 <    Dott = tmpInv;
179 <    Dotr = -XirrInv * Xitr * tmpInv;
80 >    //translation-rotation
81 >    os << Xi(3, 0) <<  "\t" << Xi(3, 1) <<  "\t" << Xi(3, 2) <<  "\t"
82 >       << Xi(4, 0) <<  "\t" << Xi(4, 1) <<  "\t" << Xi(4, 2) <<  "\t"
83 >       << Xi(5, 0) <<  "\t" << Xi(5, 1) <<  "\t" << Xi(5, 2) <<  "\t";
84      
85 <    tmp = Xirr - Xitr * XittInv * Xitr.transpose();    
86 <    tmpInv = tmp.inverse();
85 >    //rotation
86 >    os << Xi(3, 3) <<  "\t" << Xi(3, 4) <<  "\t" << Xi(3, 5) <<  "\t"
87 >       << Xi(4, 3) <<  "\t" << Xi(4, 4) <<  "\t" << Xi(4, 5) <<  "\t"
88 >       << Xi(5, 3) <<  "\t" << Xi(5, 4) <<  "\t" << Xi(5, 5) <<  "\t";
89      
90 <    Dorr = tmpInv;
90 >    
91 >    //diffusion tensor at center of resistance
92 >    //translation
93  
94 <    //calculate center of diffusion
187 <    tmp(0, 0) = Dorr(1, 1) + Dorr(2, 2);
188 <    tmp(0, 1) = - Dorr(0, 1);
189 <    tmp(0, 2) = -Dorr(0, 2);
190 <    tmp(1, 0) = -Dorr(0, 1);
191 <    tmp(1, 1) = Dorr(0, 0)  + Dorr(2, 2);
192 <    tmp(1, 2) = -Dorr(1, 2);
193 <    tmp(2, 0) = -Dorr(0, 2);
194 <    tmp(2, 1) = -Dorr(1, 2);
195 <    tmp(2, 2) = Dorr(1, 1) + Dorr(0, 0);
94 >    D = cr_->getD();
95  
96 <    Vector3d tmpVec;
97 <    tmpVec[0] = Dotr(1, 2) - Dotr(2, 1);
98 <    tmpVec[1] = Dotr(2, 0) - Dotr(0, 2);
200 <    tmpVec[2] = Dotr(0, 1) - Dotr(1, 0);
201 <
202 <    tmpInv = tmp.inverse();
96 >    os << D(0, 0) <<  "\t" << D(0, 1) <<  "\t" << D(0, 2) <<  "\t"
97 >       << D(1, 0) <<  "\t" << D(1, 1) <<  "\t" << D(1, 2) <<  "\t"
98 >       << D(2, 0) <<  "\t" << D(2, 1) <<  "\t" << D(2, 2) <<  "\t";
99      
100 <    Vector3d rod = tmpInv * tmpVec;
101 <
102 <    //calculate Diffusion Tensor at center of diffusion
103 <    Mat3x3d Uod;
208 <    Uod.setupSkewMat(rod);
100 >    //rotation-translation
101 >    os << D(0, 3) <<  "\t" << D(0, 4) <<  "\t" << D(0, 5) <<  "\t"
102 >       << D(1, 3) <<  "\t" << D(1, 4) <<  "\t" << D(1, 5) <<  "\t"
103 >       << D(2, 3) <<  "\t" << D(2, 4) <<  "\t" << D(2, 5) <<  "\t";
104      
105 <    Mat3x3d Ddtt; //translational diffusion tensor at diffusion center
106 <    Mat3x3d Ddtr; //rotational diffusion tensor at diffusion center
107 <    Mat3x3d Ddrr; //translation-rotation couplingl diffusion tensor at diffusion tensor
105 >    //translation-rotation
106 >    os << D(3, 0) <<  "\t" << D(3, 1) <<  "\t" << D(3, 2) <<  "\t"
107 >       << D(4, 0) <<  "\t" << D(4, 1) <<  "\t" << D(4, 2) <<  "\t"
108 >       << D(5, 0) <<  "\t" << D(5, 1) <<  "\t" << D(5, 2) <<  "\t";
109      
110 <    Ddtt = Dott - Uod * Dorr * Uod + Dotr.transpose() * Uod - Uod * Dotr;
111 <    Ddrr = Dorr;
112 <    Ddtr = Dotr + Dorr * Uod;
110 >    //rotation
111 >    os << D(3, 3) <<  "\t" << D(3, 4) <<  "\t" << D(3, 5) <<  "\t"
112 >       << D(4, 3) <<  "\t" << D(4, 4) <<  "\t" << D(4, 5) <<  "\t"
113 >       << D(5, 3) <<  "\t" << D(5, 4) <<  "\t" << D(5, 5) <<  "\t";
114 >    
115 >    //---------------------------------------------------------------------
116 >    
117 >    //center of diffusion
118  
119 <    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);
119 >    center = cd_->getCOR();
120  
121 <
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_.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 <
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;
121 >    os << center[0] <<  "\t" << center[1] <<  "\t" << center[2] <<  "\t";
122      
123 <      
124 < }
123 >    //resistance tensor at center of diffusion
124 >    //translation
125  
126 < void HydrodynamicsModel::writeBeads(std::ostream& os) {
270 <    std::vector<BeadParam>::iterator iter;
271 <    os << beads_.size() << std::endl;
272 <    os << "Generated by Hydro" << std::endl;
273 <    for (iter = beads_.begin(); iter != beads_.end(); ++iter) {
274 <        os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
275 <    }
126 >    Xi = cd_->getXi();
127  
128 < }
129 <
130 < void HydrodynamicsModel::writeDiffCenterAndDiffTensor(std::ostream& os) {
280 <
281 <    os << sd_->getType() << "\t";
282 <    os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\t";
283 <
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";
128 >    os << Xi(0, 0) <<  "\t" << Xi(0, 1) <<  "\t" << Xi(0, 2) <<  "\t"
129 >       << Xi(1, 0) <<  "\t" << Xi(1, 1) <<  "\t" << Xi(1, 2) <<  "\t"
130 >       << Xi(2, 0) <<  "\t" << Xi(2, 1) <<  "\t" << Xi(2, 2) <<  "\t";
131      
132 <    os << props_.Ddtr(0, 0) << "\t" << props_.Ddtr(0, 1) << "\t" << props_.Ddtr(0, 2) << "\t"
133 <        << props_.Ddtr(1, 0) << "\t" << props_.Ddtr(1, 1) << "\t" << props_.Ddtr(1, 2) << "\t"
134 <        << props_.Ddtr(2, 0) << "\t" << props_.Ddtr(2, 1) << "\t" << props_.Ddtr(2, 2) << "\t";
132 >    //rotation-translation
133 >    os << Xi(0, 3) <<  "\t" << Xi(0, 4) <<  "\t" << Xi(0, 5) <<  "\t"
134 >       << Xi(1, 3) <<  "\t" << Xi(1, 4) <<  "\t" << Xi(1, 5) <<  "\t"
135 >       << Xi(2, 3) <<  "\t" << Xi(2, 4) <<  "\t" << Xi(2, 5) <<  "\t";
136 >    
137 >    //translation-rotation
138 >    os << Xi(3, 0) <<  "\t" << Xi(3, 1) <<  "\t" << Xi(3, 2) <<  "\t"
139 >       << Xi(4, 0) <<  "\t" << Xi(4, 1) <<  "\t" << Xi(4, 2) <<  "\t"
140 >       << Xi(5, 0) <<  "\t" << Xi(5, 1) <<  "\t" << Xi(5, 2) <<  "\t";
141 >    
142 >    //rotation
143 >    os << Xi(3, 3) <<  "\t" << Xi(3, 4) <<  "\t" << Xi(3, 5) <<  "\t"
144 >       << Xi(4, 3) <<  "\t" << Xi(4, 4) <<  "\t" << Xi(4, 5) <<  "\t"
145 >       << Xi(5, 3) <<  "\t" << Xi(5, 4) <<  "\t" << Xi(5, 5) <<  "\t";
146 >    
147 >    
148 >    //diffusion tensor at center of diffusion
149 >    //translation
150  
151 <    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";        
151 >    D = cd_->getD();
152  
153 <    os << props_.Xidtt(0, 0) << "\t" << props_.Xidtt(0, 1) << "\t" << props_.Xidtt(0, 2) << "\t"
154 <        << props_.Xidtt(1, 0) << "\t" << props_.Xidtt(1, 1) << "\t" << props_.Xidtt(1, 2) << "\t"
155 <        << props_.Xidtt(2, 0) << "\t" << props_.Xidtt(2, 1) << "\t" << props_.Xidtt(2, 2) << "\t";
299 <
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";
153 >    os << D(0, 0) <<  "\t" << D(0, 1) <<  "\t" << D(0, 2) <<  "\t"
154 >       << D(1, 0) <<  "\t" << D(1, 1) <<  "\t" << D(1, 2) <<  "\t"
155 >       << D(2, 0) <<  "\t" << D(2, 1) <<  "\t" << D(2, 2) <<  "\t";
156      
157 <    os << props_.Xidtr(0, 0) << "\t" << props_.Xidtr(0, 1) << "\t" << props_.Xidtr(0, 2) << "\t"
158 <        << props_.Xidtr(1, 0) << "\t" << props_.Xidtr(1, 1) << "\t" << props_.Xidtr(1, 2) << "\t"
159 <        << props_.Xidtr(2, 0) << "\t" << props_.Xidtr(2, 1) << "\t" << props_.Xidtr(2, 2) << "\t";
160 <
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;
157 >    //rotation-translation
158 >    os << D(0, 3) <<  "\t" << D(0, 4) <<  "\t" << D(0, 5) <<  "\t"
159 >       << D(1, 3) <<  "\t" << D(1, 4) <<  "\t" << D(1, 5) <<  "\t"
160 >       << D(2, 3) <<  "\t" << D(2, 4) <<  "\t" << D(2, 5) <<  "\t";
161      
162 +    //translation-rotation
163 +    os << D(3, 0) <<  "\t" << D(3, 1) <<  "\t" << D(3, 2) <<  "\t"
164 +       << D(4, 0) <<  "\t" << D(4, 1) <<  "\t" << D(4, 2) <<  "\t"
165 +       << D(5, 0) <<  "\t" << D(5, 1) <<  "\t" << D(5, 2) <<  "\t";
166 +    
167 +    //rotation
168 +    os << D(3, 3) <<  "\t" << D(3, 4) <<  "\t" << D(3, 5) <<  "\t"
169 +       << D(4, 3) <<  "\t" << D(4, 4) <<  "\t" << D(4, 5) <<  "\t"
170 +       << D(5, 3) <<  "\t" << D(5, 4) <<  "\t" << D(5, 5) <<  "\n";    
171 +    
172 +  }
173 +  
174   }
313 –
314 – }

Diff Legend

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