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 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 + * [4]  Vardeman & Gezelter, in progress (2009).                        
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  
47 < #include "applications/hydrodynamics/HydrodynamicsModel.hpp"
48 < #include "math/LU.hpp"
49 < #include "math/DynamicRectMatrix.hpp"
50 < #include "math/SquareMatrix3.hpp"
51 < #include "utils/OOPSEConstant.hpp"
52 < namespace oopse {
53 < /**
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 < */
47 > namespace OpenMD {
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 <    }
55 >    Vector3d center;
56 >    Mat6x6d Xi, D;
57      
58 <    int nbeads = beads_.size();
82 <    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
83 <    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
84 <    Mat3x3d I;
85 <    I(0, 0) = 1.0;
86 <    I(1, 1) = 1.0;
87 <    I(2, 2) = 1.0;
58 >    os << sd_->getType() << "\t";
59      
60 <    for (std::size_t i = 0; i < nbeads; ++i) {
90 <        for (std::size_t j = 0; j < nbeads; ++j) {
91 <            Mat3x3d Tij;
92 <            if (i != j ) {
93 <                Vector3d Rij = beads_[i].pos - beads_[j].pos;
94 <                double rij = Rij.length();
95 <                double rij2 = rij * rij;
96 <                double sumSigma2OverRij2 = ((beads_[i].radius*beads_[i].radius) + (beads_[i].radius*beads_[i].radius)) / rij2;                
97 <                Mat3x3d tmpMat;
98 <                tmpMat = outProduct(Rij, Rij) / rij2;
99 <                double constant = 8.0 * NumericConstant::PI * viscosity_ * rij;
100 <                Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
101 <            }else {
102 <                double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity_ * beads_[i].radius);
103 <                Tij(0, 0) = constant;
104 <                Tij(1, 1) = constant;
105 <                Tij(2, 2) = constant;
106 <            }
107 <            B.setSubMatrix(i*3, j*3, Tij);
108 <            std::cout << Tij << std::endl;
109 <        }
110 <    }
60 >    //center of resistance
61  
62 <    std::cout << "B=\n"
113 <                  << B << std::endl;
114 <    //invert B Matrix
115 <    invertMatrix(B, C);
62 >    center = cr_->getCOR();
63  
64 <    std::cout << "C=\n"
118 <                  << C << std::endl;
119 <
120 <    //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
121 <    std::vector<Mat3x3d> U;
122 <    for (int i = 0; i < nbeads; ++i) {
123 <        Mat3x3d currU;
124 <        currU.setupSkewMat(beads_[i].pos);
125 <        U.push_back(currU);
126 <    }
64 >    os << center[0] <<  "\t" << center[1] <<  "\t" << center[2] <<  "\t";
65      
66 <    //calculate Xi matrix at arbitrary origin O
67 <    Mat3x3d Xitt;
130 <    Mat3x3d Xirr;
131 <    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);
74 <    }
75 <        
76 <    for (std::size_t i = 0; i < nbeads; ++i) {
77 <        for (std::size_t j = 0; j < nbeads; ++j) {
78 <            Mat3x3d Cij;
79 <            C.getSubMatrix(i*3, j*3, Cij);
80 <            
81 <            Xitt += Cij;
82 <            Xitr += U[i] * Cij;
83 <            Xirr += -U[i] * Cij * U[j];            
84 <            //Xirr += -U[i] * Cij * U[j] + (0.166*6 * viscosity_ * volume) * I;            
85 <        }
86 <    }
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 >    //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 >    //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 >    //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 >    
91 >    //diffusion tensor at center of resistance
92 >    //translation
93  
94 <    //invert Xi to get Diffusion Tensor at arbitrary origin O
153 <    RectMatrix<double, 6, 6> Xi;    
154 <    RectMatrix<double, 6, 6> Do;
155 <    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;    
94 >    D = cr_->getD();
95  
96 <
97 <    Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
98 <    Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
166 <    Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O
167 <
168 <    const static Mat3x3d zeroMat(0.0);
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 <    Mat3x3d XittInv(0.0);
101 <    XittInv = Xitt.inverse();
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 <    //Xirr may not be inverted,if it one of the diagonal element is zero, for example
106 <    //( a11 a12 0)
107 <    //( a21 a22 0)
108 <    //( 0    0    0)
109 <    Mat3x3d XirrInv;
110 <    XirrInv = Xirr.inverse();
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 >    //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 <    Mat3x3d tmp;
181 <    Mat3x3d tmpInv;
182 <    tmp = Xitt - Xitr.transpose() * XirrInv * Xitr;
183 <    tmpInv = tmp.inverse();
119 >    center = cd_->getCOR();
120  
121 <    Dott = kt * tmpInv;
186 <    Dotr = -kt*XirrInv * Xitr * tmpInv* 1.0E8;
187 <
188 <    tmp = Xirr - Xitr * XittInv * Xitr.transpose();    
189 <    tmpInv = tmp.inverse();
121 >    os << center[0] <<  "\t" << center[1] <<  "\t" << center[2] <<  "\t";
122      
123 <    Dorr = kt * tmpInv*1.0E16;
123 >    //resistance tensor at center of diffusion
124 >    //translation
125  
126 <    //Do.getSubMatrix(0, 0 , Dott);
194 <    //Do.getSubMatrix(3, 0, Dotr);
195 <    //Do.getSubMatrix(3, 3, Dorr);
126 >    Xi = cd_->getXi();
127  
128 <    //calculate center of diffusion
129 <    tmp(0, 0) = Dorr(1, 1) + Dorr(2, 2);
130 <    tmp(0, 1) = - Dorr(0, 1);
200 <    tmp(0, 2) = -Dorr(0, 2);
201 <    tmp(1, 0) = -Dorr(0, 1);
202 <    tmp(1, 1) = Dorr(0, 0)  + Dorr(2, 2);
203 <    tmp(1, 2) = -Dorr(1, 2);
204 <    tmp(2, 0) = -Dorr(0, 2);
205 <    tmp(2, 1) = -Dorr(1, 2);
206 <    tmp(2, 2) = Dorr(1, 1) + Dorr(0, 0);
207 <
208 <    Vector3d tmpVec;
209 <    tmpVec[0] = Dotr(1, 2) - Dotr(2, 1);
210 <    tmpVec[1] = Dotr(2, 0) - Dotr(0, 2);
211 <    tmpVec[2] = Dotr(0, 1) - Dotr(1, 0);
212 <
213 <    tmpInv = tmp.inverse();
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 <    Vector3d rod = tmpInv * tmpVec;
133 <
134 <    //calculate Diffusion Tensor at center of diffusion
135 <    Mat3x3d Uod;
219 <    Uod.setupSkewMat(rod);
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 <    Mat3x3d Ddtt; //translational diffusion tensor at diffusion center
138 <    Mat3x3d Ddtr; //rotational diffusion tensor at diffusion center
139 <    Mat3x3d Ddrr; //translation-rotation couplingl diffusion tensor at diffusion tensor
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 <    Ddtt = Dott - Uod * Dorr * Uod + Dotr.transpose() * Uod - Uod * Dotr;
143 <    Ddrr = Dorr;
144 <    Ddtr = Dotr + Dorr * Uod;
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 <    props_.diffCenter = rod;
148 <    props_.transDiff = Ddtt;
149 <    props_.transRotDiff = Ddtr;
232 <    props_.rotDiff = Ddrr;
147 >    
148 >    //diffusion tensor at center of diffusion
149 >    //translation
150  
151 <    return true;    
235 < }
151 >    D = cd_->getD();
152  
153 < void HydrodynamicsModel::writeBeads(std::ostream& os) {
154 <    std::vector<BeadParam>::iterator iter;
155 <    os << beads_.size() << std::endl;
240 <    os << "Generated by Hydro" << std::endl;
241 <    for (iter = beads_.begin(); iter != beads_.end(); ++iter) {
242 <        os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
243 <    }
244 <
245 < }
246 <
247 < 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";
252 <
253 <    os << "//diffusion center" << std::endl;
254 <    os << props_.diffCenter << std::endl;
255 <
256 <    os << "//translational diffusion tensor" << std::endl;
257 <    os << props_.transDiff << std::endl;
258 <
259 <    os << "//translation-rotation coupling diffusion tensor" << std::endl;
260 <    os << props_.transRotDiff << std::endl;
261 <
262 <    os << "//rotational diffusion tensor" << std::endl;
263 <    os << props_.rotDiff << std::endl;
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 <    /*
158 <    os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\n"
159 <
160 <    os << props_.transDiff(0, 0) << "\t" << props_.transDiff(0, 1) << "\t" << props_.transDiff(0, 2) << "\t"
269 <        << props_.transDiff(1, 0) << "\t" << props_.transDiff(1, 1) << "\t" << props_.transDiff(1, 2) << "\t"
270 <        << props_.transDiff(2, 0) << "\t" << props_.transDiff(2, 1) << "\t" << props_.transDiff(2, 2) << "\n";
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 <    os << props_.transRotDiff(0, 0) << "\t" << props_.transRotDiff(0, 1) << "\t" << props_.transRotDiff(0, 2) << "\t"
163 <        << props_.transRotDiff(1, 0) << "\t" << props_.transRotDiff(1, 1) << "\t" << props_.transRotDiff(1, 2) << "\t"
164 <        << props_.transRotDiff(2, 0) << "\t" << props_.transRotDiff(2, 1) << "\t" << props_.transRotDiff(2, 2) << "\t"
165 <
166 <    os << props_.rotDiff(0, 0) << "\t" << props_.rotDiff(0, 1) << "\t" << props_.rotDiff(0, 2) << "\t"
167 <        << props_.rotDiff(1, 0) << "\t" << props_.rotDiff(1, 1) << "\t" << props_.rotDiff(1, 2) << "\t"
168 <        << props_.rotDiff(2, 0) << "\t" << props_.rotDiff(2, 1) << "\t" << props_.rotDiff(2, 2) << ";"
169 <        << std::endl;
170 <    */
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   }
282
283 }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines