44 |
|
#include "math/DynamicRectMatrix.hpp" |
45 |
|
#include "math/SquareMatrix3.hpp" |
46 |
|
#include "utils/OOPSEConstant.hpp" |
47 |
< |
#include "applications/hydrodynamics/Spheric.hpp" |
48 |
< |
#include "applications/hydrodynamics/Ellipsoid.hpp" |
47 |
> |
#include "hydrodynamics/Sphere.hpp" |
48 |
> |
#include "hydrodynamics/Ellipsoid.hpp" |
49 |
|
#include "applications/hydrodynamics/CompositeShape.hpp" |
50 |
|
#include "math/LU.hpp" |
51 |
+ |
#include "utils/simError.h" |
52 |
|
namespace oopse { |
53 |
|
/** |
54 |
|
* Reference: |
57 |
|
* Biophysical Journal, 75(6), 3044, 1999 |
58 |
|
*/ |
59 |
|
|
60 |
< |
ApproximationModel::ApproximationModel(StuntDouble* sd, SimInfo* info): HydrodynamicsModel(sd, info){ |
61 |
< |
/* |
62 |
< |
DynamicProperty::const_iterator iter; |
63 |
< |
|
63 |
< |
iter = extraParams.find("Viscosity"); |
64 |
< |
if (iter != extraParams.end()) { |
65 |
< |
boost::any param = iter->second; |
66 |
< |
viscosity = boost::any_cast<double>(param); |
67 |
< |
}else { |
68 |
< |
std::cout << "ApproximationModel Error\n" ; |
69 |
< |
} |
70 |
< |
|
71 |
< |
iter = extraParams.find("Temperature"); |
72 |
< |
if (iter != extraParams.end()) { |
73 |
< |
boost::any param = iter->second; |
74 |
< |
temperature = boost::any_cast<double>(param); |
75 |
< |
}else { |
76 |
< |
std::cout << "ApproximationModel Error\n" ; |
77 |
< |
} |
78 |
< |
*/ |
79 |
< |
} |
80 |
< |
|
81 |
< |
bool ApproximationModel::calcHydroProps(Spheric* spheric, double viscosity, double temperature) { |
82 |
< |
return internalCalcHydroProps(static_cast<Shape*>(spheric), viscosity, temperature); |
83 |
< |
} |
84 |
< |
|
85 |
< |
bool ApproximationModel::calcHydroProps(Ellipsoid* ellipsoid, double viscosity, double temperature) { |
86 |
< |
return internalCalcHydroProps(static_cast<Shape*>(ellipsoid), viscosity, temperature); |
87 |
< |
} |
88 |
< |
bool ApproximationModel::calcHydroProps(CompositeShape* compositeShape, double viscosity, double temperature) { |
89 |
< |
return internalCalcHydroProps(static_cast<Shape*>(compositeShape), viscosity, temperature); |
90 |
< |
} |
91 |
< |
|
92 |
< |
|
93 |
< |
bool ApproximationModel::internalCalcHydroProps(Shape* shape, double viscosity, double temperature) { |
60 |
> |
ApproximationModel::ApproximationModel(StuntDouble* sd, SimInfo* info): HydrodynamicsModel(sd, info){ |
61 |
> |
} |
62 |
> |
|
63 |
> |
void ApproximationModel::init() { |
64 |
|
if (!createBeads(beads_)) { |
65 |
< |
std::cout << "can not create beads" << std::endl; |
66 |
< |
return false; |
65 |
> |
sprintf(painCave.errMsg, "ApproximationModel::init() : Can not create beads\n"); |
66 |
> |
painCave.isFatal = 1; |
67 |
> |
simError(); |
68 |
|
} |
69 |
< |
|
69 |
> |
|
70 |
> |
} |
71 |
> |
|
72 |
> |
bool ApproximationModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) { |
73 |
> |
|
74 |
|
bool ret = true; |
75 |
|
HydroProps cr; |
76 |
|
HydroProps cd; |
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 |
|
|
82 |
|
return true; |
83 |
< |
} |
84 |
< |
|
85 |
< |
bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, double viscosity, double temperature, HydroProps& cr) { |
86 |
< |
|
83 |
> |
} |
84 |
> |
|
85 |
> |
bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProps& cr) { |
86 |
> |
|
87 |
|
int nbeads = beads.size(); |
88 |
< |
DynamicRectMatrix<double> B(3*nbeads, 3*nbeads); |
89 |
< |
DynamicRectMatrix<double> C(3*nbeads, 3*nbeads); |
88 |
> |
DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads); |
89 |
> |
DynamicRectMatrix<RealType> C(3*nbeads, 3*nbeads); |
90 |
|
Mat3x3d I; |
91 |
|
I(0, 0) = 1.0; |
92 |
|
I(1, 1) = 1.0; |
93 |
|
I(2, 2) = 1.0; |
94 |
|
|
95 |
|
for (std::size_t i = 0; i < nbeads; ++i) { |
96 |
< |
for (std::size_t j = 0; j < nbeads; ++j) { |
97 |
< |
Mat3x3d Tij; |
96 |
> |
for (std::size_t j = 0; j < nbeads; ++j) { |
97 |
> |
Mat3x3d Tij; |
98 |
|
if (i != j ) { |
99 |
< |
Vector3d Rij = beads[i].pos - beads[j].pos; |
100 |
< |
double rij = Rij.length(); |
101 |
< |
double rij2 = rij * rij; |
102 |
< |
double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[i].radius*beads[i].radius)) / rij2; |
103 |
< |
Mat3x3d tmpMat; |
104 |
< |
tmpMat = outProduct(Rij, Rij) / rij2; |
105 |
< |
double constant = 8.0 * NumericConstant::PI * viscosity * rij; |
106 |
< |
Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant; |
99 |
> |
Vector3d Rij = beads[i].pos - beads[j].pos; |
100 |
> |
RealType rij = Rij.length(); |
101 |
> |
RealType rij2 = rij * rij; |
102 |
> |
RealType sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2; |
103 |
> |
Mat3x3d tmpMat; |
104 |
> |
tmpMat = outProduct(Rij, Rij) / rij2; |
105 |
> |
RealType constant = 8.0 * NumericConstant::PI * viscosity * rij; |
106 |
> |
Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant; |
107 |
|
}else { |
108 |
< |
double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius); |
109 |
< |
Tij(0, 0) = constant; |
110 |
< |
Tij(1, 1) = constant; |
111 |
< |
Tij(2, 2) = constant; |
108 |
> |
RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius); |
109 |
> |
Tij(0, 0) = constant; |
110 |
> |
Tij(1, 1) = constant; |
111 |
> |
Tij(2, 2) = constant; |
112 |
|
} |
113 |
|
B.setSubMatrix(i*3, j*3, Tij); |
114 |
< |
} |
114 |
> |
} |
115 |
|
} |
116 |
< |
|
116 |
> |
|
117 |
|
//invert B Matrix |
118 |
|
invertMatrix(B, C); |
119 |
< |
|
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); |
123 |
> |
Mat3x3d currU; |
124 |
> |
currU.setupSkewMat(beads[i].pos); |
125 |
> |
U.push_back(currU); |
126 |
|
} |
127 |
|
|
128 |
|
//calculate Xi matrix at arbitrary origin O |
129 |
|
Mat3x3d Xiott; |
130 |
|
Mat3x3d Xiorr; |
131 |
|
Mat3x3d Xiotr; |
132 |
< |
|
132 |
> |
|
133 |
|
//calculate the total volume |
134 |
< |
|
135 |
< |
double volume = 0.0; |
134 |
> |
|
135 |
> |
RealType volume = 0.0; |
136 |
|
for (std::vector<BeadParam>::iterator iter = beads.begin(); iter != beads.end(); ++iter) { |
137 |
< |
volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3); |
137 |
> |
volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3); |
138 |
|
} |
139 |
< |
|
139 |
> |
|
140 |
|
for (std::size_t i = 0; i < nbeads; ++i) { |
141 |
< |
for (std::size_t j = 0; j < nbeads; ++j) { |
142 |
< |
Mat3x3d Cij; |
143 |
< |
C.getSubMatrix(i*3, j*3, Cij); |
144 |
< |
|
145 |
< |
Xiott += Cij; |
146 |
< |
Xiotr += U[i] * Cij; |
147 |
< |
Xiorr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I; |
148 |
< |
} |
141 |
> |
for (std::size_t j = 0; j < nbeads; ++j) { |
142 |
> |
Mat3x3d Cij; |
143 |
> |
C.getSubMatrix(i*3, j*3, Cij); |
144 |
> |
|
145 |
> |
Xiott += Cij; |
146 |
> |
Xiotr += U[i] * Cij; |
147 |
> |
//Xiorr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I; |
148 |
> |
Xiorr += -U[i] * Cij * U[j]; |
149 |
> |
} |
150 |
|
} |
151 |
< |
|
152 |
< |
const double convertConstant = 6.023; //convert poise.angstrom to amu/fs |
151 |
> |
|
152 |
> |
const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs |
153 |
|
Xiott *= convertConstant; |
154 |
|
Xiotr *= convertConstant; |
155 |
|
Xiorr *= convertConstant; |
156 |
|
|
181 |
– |
|
157 |
|
|
158 |
+ |
|
159 |
|
Mat3x3d tmp; |
160 |
|
Mat3x3d tmpInv; |
161 |
|
Vector3d tmpVec; |
185 |
|
Xirrr = Xiorr - Uor * Xiott * Uor + Xiotr * Uor - Uor * Xiotr.transpose(); |
186 |
|
|
187 |
|
|
188 |
< |
SquareMatrix<double,6> Xir6x6; |
189 |
< |
SquareMatrix<double,6> Dr6x6; |
188 |
> |
SquareMatrix<RealType,6> Xir6x6; |
189 |
> |
SquareMatrix<RealType,6> Dr6x6; |
190 |
|
|
191 |
|
Xir6x6.setSubMatrix(0, 0, Xirtt); |
192 |
|
Xir6x6.setSubMatrix(0, 3, Xirtr.transpose()); |
202 |
|
Dr6x6.getSubMatrix(0, 3, Drrt); |
203 |
|
Dr6x6.getSubMatrix(3, 0, Drtr); |
204 |
|
Dr6x6.getSubMatrix(3, 3, Drrr); |
205 |
< |
double kt = OOPSEConstant::kB * temperature ; |
205 |
> |
RealType kt = OOPSEConstant::kB * temperature ; |
206 |
|
Drtt *= kt; |
207 |
|
Drrt *= kt; |
208 |
|
Drtr *= kt; |
245 |
|
|
246 |
|
return true; |
247 |
|
} |
248 |
< |
|
249 |
< |
bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, double viscosity, double temperature, HydroProps& cr) { |
250 |
< |
|
248 |
> |
|
249 |
> |
bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProps& cr) { |
250 |
> |
|
251 |
|
int nbeads = beads.size(); |
252 |
< |
DynamicRectMatrix<double> B(3*nbeads, 3*nbeads); |
253 |
< |
DynamicRectMatrix<double> C(3*nbeads, 3*nbeads); |
252 |
> |
DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads); |
253 |
> |
DynamicRectMatrix<RealType> C(3*nbeads, 3*nbeads); |
254 |
|
Mat3x3d I; |
255 |
|
I(0, 0) = 1.0; |
256 |
|
I(1, 1) = 1.0; |
257 |
|
I(2, 2) = 1.0; |
258 |
|
|
259 |
|
for (std::size_t i = 0; i < nbeads; ++i) { |
260 |
< |
for (std::size_t j = 0; j < nbeads; ++j) { |
261 |
< |
Mat3x3d Tij; |
262 |
< |
if (i != j ) { |
263 |
< |
Vector3d Rij = beads[i].pos - beads[j].pos; |
264 |
< |
double rij = Rij.length(); |
265 |
< |
double rij2 = rij * rij; |
266 |
< |
double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[i].radius*beads[i].radius)) / rij2; |
267 |
< |
Mat3x3d tmpMat; |
268 |
< |
tmpMat = outProduct(Rij, Rij) / rij2; |
269 |
< |
double constant = 8.0 * NumericConstant::PI * viscosity * rij; |
270 |
< |
Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant; |
271 |
< |
}else { |
272 |
< |
double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius); |
273 |
< |
Tij(0, 0) = constant; |
274 |
< |
Tij(1, 1) = constant; |
275 |
< |
Tij(2, 2) = constant; |
300 |
< |
} |
301 |
< |
B.setSubMatrix(i*3, j*3, Tij); |
260 |
> |
for (std::size_t j = 0; j < nbeads; ++j) { |
261 |
> |
Mat3x3d Tij; |
262 |
> |
if (i != j ) { |
263 |
> |
Vector3d Rij = beads[i].pos - beads[j].pos; |
264 |
> |
RealType rij = Rij.length(); |
265 |
> |
RealType rij2 = rij * rij; |
266 |
> |
RealType sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2; |
267 |
> |
Mat3x3d tmpMat; |
268 |
> |
tmpMat = outProduct(Rij, Rij) / rij2; |
269 |
> |
RealType constant = 8.0 * NumericConstant::PI * viscosity * rij; |
270 |
> |
Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant; |
271 |
> |
}else { |
272 |
> |
RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius); |
273 |
> |
Tij(0, 0) = constant; |
274 |
> |
Tij(1, 1) = constant; |
275 |
> |
Tij(2, 2) = constant; |
276 |
|
} |
277 |
+ |
B.setSubMatrix(i*3, j*3, Tij); |
278 |
+ |
} |
279 |
|
} |
280 |
< |
|
280 |
> |
|
281 |
|
//invert B Matrix |
282 |
|
invertMatrix(B, C); |
283 |
< |
|
283 |
> |
|
284 |
|
//prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0) |
285 |
|
std::vector<Mat3x3d> U; |
286 |
|
for (int i = 0; i < nbeads; ++i) { |
287 |
< |
Mat3x3d currU; |
288 |
< |
currU.setupSkewMat(beads[i].pos); |
289 |
< |
U.push_back(currU); |
287 |
> |
Mat3x3d currU; |
288 |
> |
currU.setupSkewMat(beads[i].pos); |
289 |
> |
U.push_back(currU); |
290 |
|
} |
291 |
|
|
292 |
|
//calculate Xi matrix at arbitrary origin O |
296 |
|
|
297 |
|
//calculate the total volume |
298 |
|
|
299 |
< |
double volume = 0.0; |
299 |
> |
RealType volume = 0.0; |
300 |
|
for (std::vector<BeadParam>::iterator iter = beads.begin(); iter != beads.end(); ++iter) { |
301 |
< |
volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3); |
301 |
> |
volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3); |
302 |
|
} |
303 |
< |
|
303 |
> |
|
304 |
|
for (std::size_t i = 0; i < nbeads; ++i) { |
305 |
< |
for (std::size_t j = 0; j < nbeads; ++j) { |
306 |
< |
Mat3x3d Cij; |
307 |
< |
C.getSubMatrix(i*3, j*3, Cij); |
305 |
> |
for (std::size_t j = 0; j < nbeads; ++j) { |
306 |
> |
Mat3x3d Cij; |
307 |
> |
C.getSubMatrix(i*3, j*3, Cij); |
308 |
|
|
309 |
< |
Xitt += Cij; |
310 |
< |
Xitr += U[i] * Cij; |
311 |
< |
Xirr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I; |
312 |
< |
} |
309 |
> |
Xitt += Cij; |
310 |
> |
Xitr += U[i] * Cij; |
311 |
> |
//Xirr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I; |
312 |
> |
Xirr += -U[i] * Cij * U[j]; |
313 |
> |
} |
314 |
|
} |
315 |
< |
|
316 |
< |
const double convertConstant = 6.023; //convert poise.angstrom to amu/fs |
315 |
> |
|
316 |
> |
const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs |
317 |
|
Xitt *= convertConstant; |
318 |
|
Xitr *= convertConstant; |
319 |
|
Xirr *= convertConstant; |
320 |
< |
|
321 |
< |
double kt = OOPSEConstant::kB * temperature; |
322 |
< |
|
320 |
> |
|
321 |
> |
RealType kt = OOPSEConstant::kB * temperature; |
322 |
> |
|
323 |
|
Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O |
324 |
|
Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O |
325 |
|
Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O |
326 |
< |
|
326 |
> |
|
327 |
|
const static Mat3x3d zeroMat(0.0); |
328 |
|
|
329 |
|
Mat3x3d XittInv(0.0); |
377 |
|
Ddrr = Dorr; |
378 |
|
Ddtr = Dotr + Dorr * Uod; |
379 |
|
|
380 |
< |
SquareMatrix<double, 6> Dd; |
380 |
> |
SquareMatrix<RealType, 6> Dd; |
381 |
|
Dd.setSubMatrix(0, 0, Ddtt); |
382 |
|
Dd.setSubMatrix(0, 3, Ddtr.transpose()); |
383 |
|
Dd.setSubMatrix(3, 0, Ddtr); |
384 |
|
Dd.setSubMatrix(3, 3, Ddrr); |
385 |
< |
SquareMatrix<double, 6> Xid; |
385 |
> |
SquareMatrix<RealType, 6> Xid; |
386 |
|
Ddtt *= kt; |
387 |
|
Ddtr *=kt; |
388 |
|
Ddrr *= kt; |
434 |
|
std::cout << Xidrr << std::endl; |
435 |
|
|
436 |
|
return true; |
437 |
< |
|
438 |
< |
} |
437 |
> |
|
438 |
> |
} |
439 |
|
|
440 |
< |
/* |
464 |
< |
void ApproximationModel::writeBeads(std::ostream& os) { |
440 |
> |
void ApproximationModel::writeBeads(std::ostream& os) { |
441 |
|
std::vector<BeadParam>::iterator iter; |
442 |
|
os << beads_.size() << std::endl; |
443 |
|
os << "Generated by Hydro" << std::endl; |
444 |
|
for (iter = beads_.begin(); iter != beads_.end(); ++iter) { |
445 |
< |
os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl; |
445 |
> |
os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl; |
446 |
|
} |
447 |
< |
|
447 |
> |
|
448 |
> |
} |
449 |
|
} |
473 |
– |
*/ |
474 |
– |
|
475 |
– |
|
476 |
– |
} |