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; |
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); |
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; |
132 |
|
|
133 |
|
//calculate the total volume |
134 |
|
|
135 |
< |
double volume = 0.0; |
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); |
138 |
|
} |
149 |
|
} |
150 |
|
} |
151 |
|
|
152 |
< |
const double convertConstant = 6.023; //convert poise.angstrom to amu/fs |
152 |
> |
const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs |
153 |
|
Xiott *= convertConstant; |
154 |
|
Xiotr *= convertConstant; |
155 |
|
Xiorr *= convertConstant; |
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; |
246 |
|
return true; |
247 |
|
} |
248 |
|
|
249 |
< |
bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, double viscosity, double temperature, HydroProps& cr) { |
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; |
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[j].radius*beads[j].radius)) / rij2; |
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 |
< |
double constant = 8.0 * NumericConstant::PI * viscosity * rij; |
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 |
< |
double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius); |
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; |
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); |
302 |
|
} |
313 |
|
} |
314 |
|
} |
315 |
|
|
316 |
< |
const double convertConstant = 6.023; //convert poise.angstrom to amu/fs |
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; |
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 |
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; |