1 |
/* |
2 |
* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. |
3 |
* |
4 |
* The University of Notre Dame grants you ("Licensee") a |
5 |
* non-exclusive, royalty free, license to use, modify and |
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 |
19 |
* notice, this list of conditions and the following disclaimer. |
20 |
* |
21 |
* 3. Redistributions in binary form must reproduce the above copyright |
22 |
* notice, this list of conditions and the following disclaimer in the |
23 |
* documentation and/or other materials provided with the |
24 |
* distribution. |
25 |
* |
26 |
* This software is provided "AS IS," without a warranty of any |
27 |
* kind. All express or implied conditions, representations and |
28 |
* warranties, including any implied warranty of merchantability, |
29 |
* fitness for a particular purpose or non-infringement, are hereby |
30 |
* excluded. The University of Notre Dame and its licensors shall not |
31 |
* be liable for any damages suffered by licensee as a result of |
32 |
* using, modifying or distributing the software or its |
33 |
* derivatives. In no event will the University of Notre Dame or its |
34 |
* licensors be liable for any lost revenue, profit or data, or for |
35 |
* direct, indirect, special, consequential, incidental or punitive |
36 |
* damages, however caused and regardless of the theory of liability, |
37 |
* arising out of the use of or inability to use software, even if the |
38 |
* University of Notre Dame has been advised of the possibility of |
39 |
* such damages. |
40 |
*/ |
41 |
|
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 |
*/ |
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; |
97 |
|
98 |
for (std::size_t i = 0; i < nbeads; ++i) { |
99 |
for (std::size_t j = 0; j < nbeads; ++j) { |
100 |
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 |
} |
119 |
|
120 |
//invert B Matrix |
121 |
invertMatrix(B, C); |
122 |
|
123 |
//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 |
} |
130 |
|
131 |
//calculate Xi matrix at arbitrary origin O |
132 |
Mat3x3d Xitt; |
133 |
Mat3x3d Xirr; |
134 |
Mat3x3d Xitr; |
135 |
|
136 |
//calculate the total volume |
137 |
|
138 |
double volume = 0.0; |
139 |
for (std::vector<BeadParam>::iterator iter = beads_.begin(); iter != beads_.end(); ++iter) { |
140 |
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); |
166 |
|
167 |
Mat3x3d XittInv(0.0); |
168 |
XittInv = Xitt.inverse(); |
169 |
|
170 |
Mat3x3d XirrInv; |
171 |
XirrInv = Xirr.inverse(); |
172 |
|
173 |
Mat3x3d tmp; |
174 |
Mat3x3d tmpInv; |
175 |
tmp = Xitt - Xitr.transpose() * XirrInv * Xitr; |
176 |
tmpInv = tmp.inverse(); |
177 |
|
178 |
Dott = kt*tmpInv; |
179 |
Dotr = -kt*XirrInv * Xitr * tmpInv; |
180 |
|
181 |
tmp = Xirr - Xitr * XittInv * Xitr.transpose(); |
182 |
tmpInv = tmp.inverse(); |
183 |
|
184 |
Dorr = kt * tmpInv; |
185 |
|
186 |
//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); |
196 |
|
197 |
Vector3d tmpVec; |
198 |
tmpVec[0] = Dotr(1, 2) - Dotr(2, 1); |
199 |
tmpVec[1] = Dotr(2, 0) - Dotr(0, 2); |
200 |
tmpVec[2] = Dotr(0, 1) - Dotr(1, 0); |
201 |
|
202 |
tmpInv = tmp.inverse(); |
203 |
|
204 |
Vector3d rod = tmpInv * tmpVec; |
205 |
|
206 |
//calculate Diffusion Tensor at center of diffusion |
207 |
Mat3x3d Uod; |
208 |
Uod.setupSkewMat(rod); |
209 |
|
210 |
Mat3x3d Ddtt; //translational diffusion tensor at diffusion center |
211 |
Mat3x3d Ddtr; //rotational diffusion tensor at diffusion center |
212 |
Mat3x3d Ddrr; //translation-rotation couplingl diffusion tensor at diffusion tensor |
213 |
|
214 |
Ddtt = Dott - Uod * Dorr * Uod + Dotr.transpose() * Uod - Uod * Dotr; |
215 |
Ddrr = Dorr; |
216 |
Ddtr = Dotr + Dorr * Uod; |
217 |
|
218 |
props_.diffCenter = rod; |
219 |
props_.Ddtt = Ddtt; |
220 |
props_.Ddtr = Ddtr; |
221 |
props_.Ddrr = Ddrr; |
222 |
|
223 |
SquareMatrix<double, 6> Dd; |
224 |
Dd.setSubMatrix(0, 0, Ddtt); |
225 |
Dd.setSubMatrix(0, 3, Ddtr.transpose()); |
226 |
Dd.setSubMatrix(3, 0, Ddtr); |
227 |
Dd.setSubMatrix(3, 3, Ddrr); |
228 |
SquareMatrix<double, 6> Xid; |
229 |
invertMatrix(Dd, Xid); |
230 |
|
231 |
//Xidtt in units of kcal*fs*mol^-1*Ang^-2 |
232 |
Xid *= OOPSEConstant::kb*temperature_; |
233 |
|
234 |
Xid.getSubMatrix(0, 0, props_.Xidtt); |
235 |
Xid.getSubMatrix(0, 3, props_.Xidrt); |
236 |
Xid.getSubMatrix(3, 0, props_.Xidtr); |
237 |
Xid.getSubMatrix(3, 3, props_.Xidrr); |
238 |
|
239 |
|
240 |
std::cout << "center of diffusion :" << std::endl; |
241 |
std::cout << rod << std::endl; |
242 |
std::cout << "diffusion tensor at center of diffusion " << std::endl; |
243 |
std::cout << "translation(A^2/fs) :" << std::endl; |
244 |
std::cout << Ddtt << std::endl; |
245 |
std::cout << "translation-rotation(A^3/fs):" << std::endl; |
246 |
std::cout << Ddtr << std::endl; |
247 |
std::cout << "rotation(A^4/fs):" << std::endl; |
248 |
std::cout << Ddrr << std::endl; |
249 |
|
250 |
std::cout << "resistance tensor at center of diffusion " << std::endl; |
251 |
std::cout << "translation(kcal*fs*mol^-1*Ang^-2) :" << std::endl; |
252 |
std::cout << props_.Xidtt << std::endl; |
253 |
std::cout << "rotation-translation (kcal*fs*mol^-1*Ang^-3):" << std::endl; |
254 |
std::cout << props_.Xidrt << std::endl; |
255 |
std::cout << "translation-rotation(kcal*fs*mol^-1*Ang^-3):" << std::endl; |
256 |
std::cout << props_.Xidtr << std::endl; |
257 |
std::cout << "rotation(kcal*fs*mol^-1*Ang^-4):" << std::endl; |
258 |
std::cout << props_.Xidrr << std::endl; |
259 |
|
260 |
|
261 |
} |
262 |
|
263 |
void HydrodynamicsModel::writeBeads(std::ostream& os) { |
264 |
std::vector<BeadParam>::iterator iter; |
265 |
os << beads_.size() << std::endl; |
266 |
os << "Generated by Hydro" << std::endl; |
267 |
for (iter = beads_.begin(); iter != beads_.end(); ++iter) { |
268 |
os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl; |
269 |
} |
270 |
|
271 |
} |
272 |
|
273 |
void HydrodynamicsModel::writeDiffCenterAndDiffTensor(std::ostream& os) { |
274 |
|
275 |
os << sd_->getType() << "\t"; |
276 |
os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\t"; |
277 |
|
278 |
os << props_.Ddtt(0, 0) << "\t" << props_.Ddtt(0, 1) << "\t" << props_.Ddtt(0, 2) << "\t" |
279 |
<< props_.Ddtt(1, 0) << "\t" << props_.Ddtt(1, 1) << "\t" << props_.Ddtt(1, 2) << "\t" |
280 |
<< props_.Ddtt(2, 0) << "\t" << props_.Ddtt(2, 1) << "\t" << props_.Ddtt(2, 2) << "\t"; |
281 |
|
282 |
os << props_.Ddtr(0, 0) << "\t" << props_.Ddtr(0, 1) << "\t" << props_.Ddtr(0, 2) << "\t" |
283 |
<< props_.Ddtr(1, 0) << "\t" << props_.Ddtr(1, 1) << "\t" << props_.Ddtr(1, 2) << "\t" |
284 |
<< props_.Ddtr(2, 0) << "\t" << props_.Ddtr(2, 1) << "\t" << props_.Ddtr(2, 2) << "\t"; |
285 |
|
286 |
os << props_.Ddrr(0, 0) << "\t" << props_.Ddrr(0, 1) << "\t" << props_.Ddrr(0, 2) << "\t" |
287 |
<< props_.Ddrr(1, 0) << "\t" << props_.Ddrr(1, 1) << "\t" << props_.Ddrr(1, 2) << "\t" |
288 |
<< props_.Ddrr(2, 0) << "\t" << props_.Ddrr(2, 1) << "\t" << props_.Ddrr(2, 2) <<"\t"; |
289 |
|
290 |
os << props_.Xidtt(0, 0) << "\t" << props_.Xidtt(0, 1) << "\t" << props_.Xidtt(0, 2) << "\t" |
291 |
<< props_.Xidtt(1, 0) << "\t" << props_.Xidtt(1, 1) << "\t" << props_.Xidtt(1, 2) << "\t" |
292 |
<< props_.Xidtt(2, 0) << "\t" << props_.Xidtt(2, 1) << "\t" << props_.Xidtt(2, 2) << "\t"; |
293 |
|
294 |
os << props_.Xidrt(0, 0) << "\t" << props_.Xidrt(0, 1) << "\t" << props_.Xidrt(0, 2) << "\t" |
295 |
<< props_.Xidrt(1, 0) << "\t" << props_.Xidrt(1, 1) << "\t" << props_.Xidrt(1, 2) << "\t" |
296 |
<< props_.Xidrt(2, 0) << "\t" << props_.Xidrt(2, 1) << "\t" << props_.Xidrt(2, 2) << "\t"; |
297 |
|
298 |
os << props_.Xidtr(0, 0) << "\t" << props_.Xidtr(0, 1) << "\t" << props_.Xidtr(0, 2) << "\t" |
299 |
<< props_.Xidtr(1, 0) << "\t" << props_.Xidtr(1, 1) << "\t" << props_.Xidtr(1, 2) << "\t" |
300 |
<< props_.Xidtr(2, 0) << "\t" << props_.Xidtr(2, 1) << "\t" << props_.Xidtr(2, 2) << "\t"; |
301 |
|
302 |
os << props_.Xidrr(0, 0) << "\t" << props_.Xidrr(0, 1) << "\t" << props_.Xidrr(0, 2) << "\t" |
303 |
<< props_.Xidrr(1, 0) << "\t" << props_.Xidrr(1, 1) << "\t" << props_.Xidrr(1, 2) << "\t" |
304 |
<< props_.Xidrr(2, 0) << "\t" << props_.Xidrr(2, 1) << "\t" << props_.Xidrr(2, 2) << std::endl; |
305 |
|
306 |
} |
307 |
|
308 |
} |