69 |
|
|
70 |
|
} |
71 |
|
|
72 |
< |
bool ApproximationModel::calcHydroProps(Shape* shape, double viscosity, double temperature) { |
72 |
> |
bool ApproximationModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) { |
73 |
|
|
74 |
|
bool ret = true; |
75 |
|
HydroProps cr; |
82 |
|
return true; |
83 |
|
} |
84 |
|
|
85 |
< |
bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, double viscosity, double temperature, HydroProps& cr) { |
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; |
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[j].radius*beads[j].radius)) / rij2; |
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 |
< |
double constant = 8.0 * NumericConstant::PI * viscosity * rij; |
106 |
< |
Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant; |
105 |
> |
RealType constant = 8.0 * NumericConstant::PI * viscosity * rij; |
106 |
> |
RealType tmp1 = 1.0 + sumSigma2OverRij2/3.0; |
107 |
> |
RealType tmp2 = 1.0 - sumSigma2OverRij2; |
108 |
> |
Tij = (tmp1 * I + tmp2 * tmpMat ) / constant; |
109 |
|
}else { |
110 |
< |
double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius); |
110 |
> |
RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius); |
111 |
|
Tij(0, 0) = constant; |
112 |
|
Tij(1, 1) = constant; |
113 |
|
Tij(2, 2) = constant; |
134 |
|
|
135 |
|
//calculate the total volume |
136 |
|
|
137 |
< |
double volume = 0.0; |
137 |
> |
RealType volume = 0.0; |
138 |
|
for (std::vector<BeadParam>::iterator iter = beads.begin(); iter != beads.end(); ++iter) { |
139 |
|
volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3); |
140 |
|
} |
151 |
|
} |
152 |
|
} |
153 |
|
|
154 |
< |
const double convertConstant = 6.023; //convert poise.angstrom to amu/fs |
154 |
> |
const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs |
155 |
|
Xiott *= convertConstant; |
156 |
|
Xiotr *= convertConstant; |
157 |
|
Xiorr *= convertConstant; |
187 |
|
Xirrr = Xiorr - Uor * Xiott * Uor + Xiotr * Uor - Uor * Xiotr.transpose(); |
188 |
|
|
189 |
|
|
190 |
< |
SquareMatrix<double,6> Xir6x6; |
191 |
< |
SquareMatrix<double,6> Dr6x6; |
190 |
> |
SquareMatrix<RealType,6> Xir6x6; |
191 |
> |
SquareMatrix<RealType,6> Dr6x6; |
192 |
|
|
193 |
|
Xir6x6.setSubMatrix(0, 0, Xirtt); |
194 |
|
Xir6x6.setSubMatrix(0, 3, Xirtr.transpose()); |
204 |
|
Dr6x6.getSubMatrix(0, 3, Drrt); |
205 |
|
Dr6x6.getSubMatrix(3, 0, Drtr); |
206 |
|
Dr6x6.getSubMatrix(3, 3, Drrr); |
207 |
< |
double kt = OOPSEConstant::kB * temperature ; |
207 |
> |
RealType kt = OOPSEConstant::kB * temperature ; |
208 |
|
Drtt *= kt; |
209 |
|
Drrt *= kt; |
210 |
|
Drtr *= kt; |
248 |
|
return true; |
249 |
|
} |
250 |
|
|
251 |
< |
bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, double viscosity, double temperature, HydroProps& cr) { |
251 |
> |
bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProps& cr) { |
252 |
|
|
253 |
|
int nbeads = beads.size(); |
254 |
< |
DynamicRectMatrix<double> B(3*nbeads, 3*nbeads); |
255 |
< |
DynamicRectMatrix<double> C(3*nbeads, 3*nbeads); |
254 |
> |
DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads); |
255 |
> |
DynamicRectMatrix<RealType> C(3*nbeads, 3*nbeads); |
256 |
|
Mat3x3d I; |
257 |
|
I(0, 0) = 1.0; |
258 |
|
I(1, 1) = 1.0; |
263 |
|
Mat3x3d Tij; |
264 |
|
if (i != j ) { |
265 |
|
Vector3d Rij = beads[i].pos - beads[j].pos; |
266 |
< |
double rij = Rij.length(); |
267 |
< |
double rij2 = rij * rij; |
268 |
< |
double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2; |
266 |
> |
RealType rij = Rij.length(); |
267 |
> |
RealType rij2 = rij * rij; |
268 |
> |
RealType sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2; |
269 |
|
Mat3x3d tmpMat; |
270 |
|
tmpMat = outProduct(Rij, Rij) / rij2; |
271 |
< |
double constant = 8.0 * NumericConstant::PI * viscosity * rij; |
272 |
< |
Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant; |
271 |
> |
RealType constant = 8.0 * NumericConstant::PI * viscosity * rij; |
272 |
> |
RealType tmp1 = 1.0 + sumSigma2OverRij2/3.0; |
273 |
> |
RealType tmp2 = 1.0 - sumSigma2OverRij2; |
274 |
> |
Tij = (tmp1 * I + tmp2 * tmpMat ) / constant; |
275 |
|
}else { |
276 |
< |
double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius); |
276 |
> |
RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius); |
277 |
|
Tij(0, 0) = constant; |
278 |
|
Tij(1, 1) = constant; |
279 |
|
Tij(2, 2) = constant; |
300 |
|
|
301 |
|
//calculate the total volume |
302 |
|
|
303 |
< |
double volume = 0.0; |
303 |
> |
RealType volume = 0.0; |
304 |
|
for (std::vector<BeadParam>::iterator iter = beads.begin(); iter != beads.end(); ++iter) { |
305 |
|
volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3); |
306 |
|
} |
317 |
|
} |
318 |
|
} |
319 |
|
|
320 |
< |
const double convertConstant = 6.023; //convert poise.angstrom to amu/fs |
320 |
> |
const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs |
321 |
|
Xitt *= convertConstant; |
322 |
|
Xitr *= convertConstant; |
323 |
|
Xirr *= convertConstant; |
324 |
|
|
325 |
< |
double kt = OOPSEConstant::kB * temperature; |
325 |
> |
RealType kt = OOPSEConstant::kB * temperature; |
326 |
|
|
327 |
|
Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O |
328 |
|
Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O |
381 |
|
Ddrr = Dorr; |
382 |
|
Ddtr = Dotr + Dorr * Uod; |
383 |
|
|
384 |
< |
SquareMatrix<double, 6> Dd; |
384 |
> |
SquareMatrix<RealType, 6> Dd; |
385 |
|
Dd.setSubMatrix(0, 0, Ddtt); |
386 |
|
Dd.setSubMatrix(0, 3, Ddtr.transpose()); |
387 |
|
Dd.setSubMatrix(3, 0, Ddtr); |
388 |
|
Dd.setSubMatrix(3, 3, Ddrr); |
389 |
< |
SquareMatrix<double, 6> Xid; |
389 |
> |
SquareMatrix<RealType, 6> Xid; |
390 |
|
Ddtt *= kt; |
391 |
|
Ddtr *=kt; |
392 |
|
Ddrr *= kt; |