# | Line 39 | Line 39 | |
---|---|---|
39 | * such damages. | |
40 | */ | |
41 | ||
42 | + | #include <cmath> |
43 | #include "integrators/RNEMD.hpp" | |
44 | + | #include "math/Vector3.hpp" |
45 | #include "math/SquareMatrix3.hpp" | |
46 | + | #include "math/Polynomial.hpp" |
47 | #include "primitives/Molecule.hpp" | |
48 | #include "primitives/StuntDouble.hpp" | |
49 | + | #include "utils/OOPSEConstant.hpp" |
50 | + | #include "utils/Tuple.hpp" |
51 | ||
52 | #ifndef IS_MPI | |
53 | #include "math/SeqRandNumGen.hpp" | |
# | Line 50 | Line 55 | |
55 | #include "math/ParallelRandNumGen.hpp" | |
56 | #endif | |
57 | ||
58 | < | /* Remove me after testing*/ |
54 | < | /* |
55 | < | #include <cstdio> |
56 | < | #include <iostream> |
57 | < | */ |
58 | < | /*End remove me*/ |
58 | > | #define HONKING_LARGE_VALUE 1.0e10 |
59 | ||
60 | namespace oopse { | |
61 | ||
62 | < | RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info) { |
63 | < | |
62 | > | RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) { |
63 | > | |
64 | > | failTrialCount_ = 0; |
65 | > | failRootCount_ = 0; |
66 | > | |
67 | int seedValue; | |
68 | Globals * simParams = info->getSimParams(); | |
69 | ||
70 | < | stringToEnumMap_["Kinetic"] = rnemdKinetic; |
70 | > | stringToEnumMap_["KineticSwap"] = rnemdKineticSwap; |
71 | > | stringToEnumMap_["KineticScale"] = rnemdKineticScale; |
72 | > | stringToEnumMap_["PxScale"] = rnemdPxScale; |
73 | > | stringToEnumMap_["PyScale"] = rnemdPyScale; |
74 | > | stringToEnumMap_["PzScale"] = rnemdPzScale; |
75 | stringToEnumMap_["Px"] = rnemdPx; | |
76 | stringToEnumMap_["Py"] = rnemdPy; | |
77 | stringToEnumMap_["Pz"] = rnemdPz; | |
78 | stringToEnumMap_["Unknown"] = rnemdUnknown; | |
79 | ||
80 | rnemdObjectSelection_ = simParams->getRNEMD_objectSelection(); | |
74 | – | |
75 | – | std::cerr << "calling evaluator with " << rnemdObjectSelection_ << "\n"; |
81 | evaluator_.loadScriptString(rnemdObjectSelection_); | |
82 | < | std::cerr << "done"; |
82 | > | seleMan_.setSelectionSet(evaluator_.evaluate()); |
83 | > | |
84 | > | // do some sanity checking |
85 | > | |
86 | > | int selectionCount = seleMan_.getSelectionCount(); |
87 | > | int nIntegrable = info->getNGlobalIntegrableObjects(); |
88 | > | |
89 | > | if (selectionCount > nIntegrable) { |
90 | > | sprintf(painCave.errMsg, |
91 | > | "RNEMD warning: The current RNEMD_objectSelection,\n" |
92 | > | "\t\t%s\n" |
93 | > | "\thas resulted in %d selected objects. However,\n" |
94 | > | "\tthe total number of integrable objects in the system\n" |
95 | > | "\tis only %d. This is almost certainly not what you want\n" |
96 | > | "\tto do. A likely cause of this is forgetting the _RB_0\n" |
97 | > | "\tselector in the selection script!\n", |
98 | > | rnemdObjectSelection_.c_str(), |
99 | > | selectionCount, nIntegrable); |
100 | > | painCave.isFatal = 0; |
101 | > | simError(); |
102 | > | |
103 | > | } |
104 | ||
105 | < | const std::string st = simParams->getRNEMD_swapType(); |
105 | > | const std::string st = simParams->getRNEMD_exchangeType(); |
106 | ||
107 | std::map<std::string, RNEMDTypeEnum>::iterator i; | |
108 | i = stringToEnumMap_.find(st); | |
109 | < | rnemdType_ = (i == stringToEnumMap_.end()) ? RNEMD::rnemdUnknown : i->second; |
109 | > | rnemdType_ = (i == stringToEnumMap_.end()) ? RNEMD::rnemdUnknown : i->second; |
110 | > | if (rnemdType_ == rnemdUnknown) { |
111 | > | std::cerr << "WARNING! RNEMD Type Unknown!\n"; |
112 | > | } |
113 | ||
114 | < | set_RNEMD_swapTime(simParams->getRNEMD_swapTime()); |
114 | > | #ifdef IS_MPI |
115 | > | if (worldRank == 0) { |
116 | > | #endif |
117 | > | |
118 | > | std::string rnemdFileName; |
119 | > | std::string xTempFileName; |
120 | > | std::string yTempFileName; |
121 | > | std::string zTempFileName; |
122 | > | switch(rnemdType_) { |
123 | > | case rnemdKineticSwap : |
124 | > | case rnemdKineticScale : |
125 | > | rnemdFileName = "temperature.log"; |
126 | > | break; |
127 | > | case rnemdPx : |
128 | > | case rnemdPxScale : |
129 | > | case rnemdPy : |
130 | > | case rnemdPyScale : |
131 | > | rnemdFileName = "momemtum.log"; |
132 | > | xTempFileName = "temperatureX.log"; |
133 | > | yTempFileName = "temperatureY.log"; |
134 | > | zTempFileName = "temperatureZ.log"; |
135 | > | xTempLog_.open(xTempFileName.c_str()); |
136 | > | yTempLog_.open(yTempFileName.c_str()); |
137 | > | zTempLog_.open(zTempFileName.c_str()); |
138 | > | break; |
139 | > | case rnemdPz : |
140 | > | case rnemdPzScale : |
141 | > | case rnemdUnknown : |
142 | > | default : |
143 | > | rnemdFileName = "rnemd.log"; |
144 | > | break; |
145 | > | } |
146 | > | rnemdLog_.open(rnemdFileName.c_str()); |
147 | > | |
148 | > | #ifdef IS_MPI |
149 | > | } |
150 | > | #endif |
151 | > | |
152 | > | set_RNEMD_exchange_time(simParams->getRNEMD_exchangeTime()); |
153 | set_RNEMD_nBins(simParams->getRNEMD_nBins()); | |
154 | < | exchangeSum_ = 0.0; |
155 | < | |
154 | > | midBin_ = nBins_ / 2; |
155 | > | if (simParams->haveRNEMD_logWidth()) { |
156 | > | rnemdLogWidth_ = simParams->getRNEMD_logWidth(); |
157 | > | if (rnemdLogWidth_ != nBins_ || rnemdLogWidth_ != midBin_ + 1) { |
158 | > | std::cerr << "WARNING! RNEMD_logWidth has abnormal value!\n"; |
159 | > | std::cerr << "Automaically set back to default.\n"; |
160 | > | rnemdLogWidth_ = nBins_; |
161 | > | } |
162 | > | } else { |
163 | > | rnemdLogWidth_ = nBins_; |
164 | > | } |
165 | > | valueHist_.resize(rnemdLogWidth_, 0.0); |
166 | > | valueCount_.resize(rnemdLogWidth_, 0); |
167 | > | xTempHist_.resize(rnemdLogWidth_, 0.0); |
168 | > | yTempHist_.resize(rnemdLogWidth_, 0.0); |
169 | > | zTempHist_.resize(rnemdLogWidth_, 0.0); |
170 | > | |
171 | > | set_RNEMD_exchange_total(0.0); |
172 | > | if (simParams->haveRNEMD_targetFlux()) { |
173 | > | set_RNEMD_target_flux(simParams->getRNEMD_targetFlux()); |
174 | > | } else { |
175 | > | set_RNEMD_target_flux(0.0); |
176 | > | } |
177 | > | |
178 | #ifndef IS_MPI | |
179 | if (simParams->haveSeed()) { | |
180 | seedValue = simParams->getSeed(); | |
# | Line 105 | Line 194 | namespace oopse { | |
194 | ||
195 | RNEMD::~RNEMD() { | |
196 | delete randNumGen_; | |
197 | + | |
198 | + | std::cerr << "total fail trials: " << failTrialCount_ << "\n"; |
199 | + | #ifdef IS_MPI |
200 | + | if (worldRank == 0) { |
201 | + | #endif |
202 | + | rnemdLog_.close(); |
203 | + | if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPyScale) |
204 | + | std::cerr<< "total root-checking warnings: " << failRootCount_ << "\n"; |
205 | + | if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPy || rnemdType_ == rnemdPyScale) { |
206 | + | xTempLog_.close(); |
207 | + | yTempLog_.close(); |
208 | + | zTempLog_.close(); |
209 | + | } |
210 | + | #ifdef IS_MPI |
211 | + | } |
212 | + | #endif |
213 | } | |
214 | ||
215 | void RNEMD::doSwap() { | |
111 | – | std::cerr << "in RNEMD!\n"; |
112 | – | std::cerr << "nBins = " << nBins_ << "\n"; |
113 | – | std::cerr << "swapTime = " << swapTime_ << "\n"; |
114 | – | std::cerr << "exchangeSum = " << exchangeSum_ << "\n"; |
115 | – | std::cerr << "swapType = " << rnemdType_ << "\n"; |
116 | – | std::cerr << "selection = " << rnemdObjectSelection_ << "\n"; |
216 | ||
217 | + | Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); |
218 | + | Mat3x3d hmat = currentSnap_->getHmat(); |
219 | + | |
220 | seleMan_.setSelectionSet(evaluator_.evaluate()); | |
221 | ||
222 | < | std::cerr << "selectionCount = " << seleMan_.getSelectionCount() << "\n\n"; |
121 | < | |
122 | < | int i; |
222 | > | int selei; |
223 | StuntDouble* sd; | |
224 | + | int idx; |
225 | ||
226 | < | for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) { |
226 | > | RealType min_val; |
227 | > | bool min_found = false; |
228 | > | StuntDouble* min_sd; |
229 | > | |
230 | > | RealType max_val; |
231 | > | bool max_found = false; |
232 | > | StuntDouble* max_sd; |
233 | > | |
234 | > | for (sd = seleMan_.beginSelected(selei); sd != NULL; |
235 | > | sd = seleMan_.nextSelected(selei)) { |
236 | > | |
237 | > | idx = sd->getLocalIndex(); |
238 | > | |
239 | Vector3d pos = sd->getPos(); | |
240 | < | //wrap the stuntdoubles into a cell |
240 | > | |
241 | > | // wrap the stuntdouble's position back into the box: |
242 | > | |
243 | if (usePeriodicBoundaryConditions_) | |
244 | < | info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos); |
245 | < | int binNo = int(nBins_ * (pos.z()) / hmat(2,2)); |
246 | < | sliceSDLists_[binNo].push_back(sd); |
244 | > | currentSnap_->wrapVector(pos); |
245 | > | |
246 | > | // which bin is this stuntdouble in? |
247 | > | // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)] |
248 | > | |
249 | > | int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_; |
250 | > | |
251 | > | |
252 | > | // if we're in bin 0 or the middleBin |
253 | > | if (binNo == 0 || binNo == midBin_) { |
254 | > | |
255 | > | RealType mass = sd->getMass(); |
256 | > | Vector3d vel = sd->getVel(); |
257 | > | RealType value; |
258 | > | |
259 | > | switch(rnemdType_) { |
260 | > | case rnemdKineticSwap : |
261 | > | |
262 | > | value = mass * (vel[0]*vel[0] + vel[1]*vel[1] + |
263 | > | vel[2]*vel[2]); |
264 | > | if (sd->isDirectional()) { |
265 | > | Vector3d angMom = sd->getJ(); |
266 | > | Mat3x3d I = sd->getI(); |
267 | > | |
268 | > | if (sd->isLinear()) { |
269 | > | int i = sd->linearAxis(); |
270 | > | int j = (i + 1) % 3; |
271 | > | int k = (i + 2) % 3; |
272 | > | value += angMom[j] * angMom[j] / I(j, j) + |
273 | > | angMom[k] * angMom[k] / I(k, k); |
274 | > | } else { |
275 | > | value += angMom[0]*angMom[0]/I(0, 0) |
276 | > | + angMom[1]*angMom[1]/I(1, 1) |
277 | > | + angMom[2]*angMom[2]/I(2, 2); |
278 | > | } |
279 | > | } |
280 | > | //make exchangeSum_ comparable between swap & scale |
281 | > | //temporarily without using energyConvert |
282 | > | //value = value * 0.5 / OOPSEConstant::energyConvert; |
283 | > | value *= 0.5; |
284 | > | break; |
285 | > | case rnemdPx : |
286 | > | value = mass * vel[0]; |
287 | > | break; |
288 | > | case rnemdPy : |
289 | > | value = mass * vel[1]; |
290 | > | break; |
291 | > | case rnemdPz : |
292 | > | value = mass * vel[2]; |
293 | > | break; |
294 | > | default : |
295 | > | break; |
296 | > | } |
297 | > | |
298 | > | if (binNo == 0) { |
299 | > | if (!min_found) { |
300 | > | min_val = value; |
301 | > | min_sd = sd; |
302 | > | min_found = true; |
303 | > | } else { |
304 | > | if (min_val > value) { |
305 | > | min_val = value; |
306 | > | min_sd = sd; |
307 | > | } |
308 | > | } |
309 | > | } else { //midBin_ |
310 | > | if (!max_found) { |
311 | > | max_val = value; |
312 | > | max_sd = sd; |
313 | > | max_found = true; |
314 | > | } else { |
315 | > | if (max_val < value) { |
316 | > | max_val = value; |
317 | > | max_sd = sd; |
318 | > | } |
319 | > | } |
320 | > | } |
321 | > | } |
322 | } | |
323 | ||
324 | < | } |
324 | > | #ifdef IS_MPI |
325 | > | int nProc, worldRank; |
326 | > | |
327 | > | nProc = MPI::COMM_WORLD.Get_size(); |
328 | > | worldRank = MPI::COMM_WORLD.Get_rank(); |
329 | > | |
330 | > | bool my_min_found = min_found; |
331 | > | bool my_max_found = max_found; |
332 | > | |
333 | > | // Even if we didn't find a minimum, did someone else? |
334 | > | MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, |
335 | > | 1, MPI::BOOL, MPI::LAND); |
336 | > | |
337 | > | // Even if we didn't find a maximum, did someone else? |
338 | > | MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, |
339 | > | 1, MPI::BOOL, MPI::LAND); |
340 | > | |
341 | > | struct { |
342 | > | RealType val; |
343 | > | int rank; |
344 | > | } max_vals, min_vals; |
345 | > | |
346 | > | if (min_found) { |
347 | > | if (my_min_found) |
348 | > | min_vals.val = min_val; |
349 | > | else |
350 | > | min_vals.val = HONKING_LARGE_VALUE; |
351 | > | |
352 | > | min_vals.rank = worldRank; |
353 | > | |
354 | > | // Who had the minimum? |
355 | > | MPI::COMM_WORLD.Allreduce(&min_vals, &min_vals, |
356 | > | 1, MPI::REALTYPE_INT, MPI::MINLOC); |
357 | > | min_val = min_vals.val; |
358 | > | } |
359 | > | |
360 | > | if (max_found) { |
361 | > | if (my_max_found) |
362 | > | max_vals.val = max_val; |
363 | > | else |
364 | > | max_vals.val = -HONKING_LARGE_VALUE; |
365 | > | |
366 | > | max_vals.rank = worldRank; |
367 | > | |
368 | > | // Who had the maximum? |
369 | > | MPI::COMM_WORLD.Allreduce(&max_vals, &max_vals, |
370 | > | 1, MPI::REALTYPE_INT, MPI::MAXLOC); |
371 | > | max_val = max_vals.val; |
372 | > | } |
373 | > | #endif |
374 | > | |
375 | > | if (max_found && min_found) { |
376 | > | if (min_val< max_val) { |
377 | > | |
378 | > | #ifdef IS_MPI |
379 | > | if (max_vals.rank == worldRank && min_vals.rank == worldRank) { |
380 | > | // I have both maximum and minimum, so proceed like a single |
381 | > | // processor version: |
382 | > | #endif |
383 | > | // objects to be swapped: velocity & angular velocity |
384 | > | Vector3d min_vel = min_sd->getVel(); |
385 | > | Vector3d max_vel = max_sd->getVel(); |
386 | > | RealType temp_vel; |
387 | > | |
388 | > | switch(rnemdType_) { |
389 | > | case rnemdKineticSwap : |
390 | > | min_sd->setVel(max_vel); |
391 | > | max_sd->setVel(min_vel); |
392 | > | if (min_sd->isDirectional() && max_sd->isDirectional()) { |
393 | > | Vector3d min_angMom = min_sd->getJ(); |
394 | > | Vector3d max_angMom = max_sd->getJ(); |
395 | > | min_sd->setJ(max_angMom); |
396 | > | max_sd->setJ(min_angMom); |
397 | > | } |
398 | > | break; |
399 | > | case rnemdPx : |
400 | > | temp_vel = min_vel.x(); |
401 | > | min_vel.x() = max_vel.x(); |
402 | > | max_vel.x() = temp_vel; |
403 | > | min_sd->setVel(min_vel); |
404 | > | max_sd->setVel(max_vel); |
405 | > | break; |
406 | > | case rnemdPy : |
407 | > | temp_vel = min_vel.y(); |
408 | > | min_vel.y() = max_vel.y(); |
409 | > | max_vel.y() = temp_vel; |
410 | > | min_sd->setVel(min_vel); |
411 | > | max_sd->setVel(max_vel); |
412 | > | break; |
413 | > | case rnemdPz : |
414 | > | temp_vel = min_vel.z(); |
415 | > | min_vel.z() = max_vel.z(); |
416 | > | max_vel.z() = temp_vel; |
417 | > | min_sd->setVel(min_vel); |
418 | > | max_sd->setVel(max_vel); |
419 | > | break; |
420 | > | default : |
421 | > | break; |
422 | > | } |
423 | > | #ifdef IS_MPI |
424 | > | // the rest of the cases only apply in parallel simulations: |
425 | > | } else if (max_vals.rank == worldRank) { |
426 | > | // I had the max, but not the minimum |
427 | > | |
428 | > | Vector3d min_vel; |
429 | > | Vector3d max_vel = max_sd->getVel(); |
430 | > | MPI::Status status; |
431 | > | |
432 | > | // point-to-point swap of the velocity vector |
433 | > | MPI::COMM_WORLD.Sendrecv(max_vel.getArrayPointer(), 3, MPI::REALTYPE, |
434 | > | min_vals.rank, 0, |
435 | > | min_vel.getArrayPointer(), 3, MPI::REALTYPE, |
436 | > | min_vals.rank, 0, status); |
437 | > | |
438 | > | switch(rnemdType_) { |
439 | > | case rnemdKineticSwap : |
440 | > | max_sd->setVel(min_vel); |
441 | > | |
442 | > | if (max_sd->isDirectional()) { |
443 | > | Vector3d min_angMom; |
444 | > | Vector3d max_angMom = max_sd->getJ(); |
445 | > | |
446 | > | // point-to-point swap of the angular momentum vector |
447 | > | MPI::COMM_WORLD.Sendrecv(max_angMom.getArrayPointer(), 3, |
448 | > | MPI::REALTYPE, min_vals.rank, 1, |
449 | > | min_angMom.getArrayPointer(), 3, |
450 | > | MPI::REALTYPE, min_vals.rank, 1, |
451 | > | status); |
452 | > | |
453 | > | max_sd->setJ(min_angMom); |
454 | > | } |
455 | > | break; |
456 | > | case rnemdPx : |
457 | > | max_vel.x() = min_vel.x(); |
458 | > | max_sd->setVel(max_vel); |
459 | > | break; |
460 | > | case rnemdPy : |
461 | > | max_vel.y() = min_vel.y(); |
462 | > | max_sd->setVel(max_vel); |
463 | > | break; |
464 | > | case rnemdPz : |
465 | > | max_vel.z() = min_vel.z(); |
466 | > | max_sd->setVel(max_vel); |
467 | > | break; |
468 | > | default : |
469 | > | break; |
470 | > | } |
471 | > | } else if (min_vals.rank == worldRank) { |
472 | > | // I had the minimum but not the maximum: |
473 | > | |
474 | > | Vector3d max_vel; |
475 | > | Vector3d min_vel = min_sd->getVel(); |
476 | > | MPI::Status status; |
477 | > | |
478 | > | // point-to-point swap of the velocity vector |
479 | > | MPI::COMM_WORLD.Sendrecv(min_vel.getArrayPointer(), 3, MPI::REALTYPE, |
480 | > | max_vals.rank, 0, |
481 | > | max_vel.getArrayPointer(), 3, MPI::REALTYPE, |
482 | > | max_vals.rank, 0, status); |
483 | > | |
484 | > | switch(rnemdType_) { |
485 | > | case rnemdKineticSwap : |
486 | > | min_sd->setVel(max_vel); |
487 | > | |
488 | > | if (min_sd->isDirectional()) { |
489 | > | Vector3d min_angMom = min_sd->getJ(); |
490 | > | Vector3d max_angMom; |
491 | > | |
492 | > | // point-to-point swap of the angular momentum vector |
493 | > | MPI::COMM_WORLD.Sendrecv(min_angMom.getArrayPointer(), 3, |
494 | > | MPI::REALTYPE, max_vals.rank, 1, |
495 | > | max_angMom.getArrayPointer(), 3, |
496 | > | MPI::REALTYPE, max_vals.rank, 1, |
497 | > | status); |
498 | > | |
499 | > | min_sd->setJ(max_angMom); |
500 | > | } |
501 | > | break; |
502 | > | case rnemdPx : |
503 | > | min_vel.x() = max_vel.x(); |
504 | > | min_sd->setVel(min_vel); |
505 | > | break; |
506 | > | case rnemdPy : |
507 | > | min_vel.y() = max_vel.y(); |
508 | > | min_sd->setVel(min_vel); |
509 | > | break; |
510 | > | case rnemdPz : |
511 | > | min_vel.z() = max_vel.z(); |
512 | > | min_sd->setVel(min_vel); |
513 | > | break; |
514 | > | default : |
515 | > | break; |
516 | > | } |
517 | > | } |
518 | > | #endif |
519 | > | exchangeSum_ += max_val - min_val; |
520 | > | } else { |
521 | > | std::cerr << "exchange NOT performed!\nmin_val > max_val.\n"; |
522 | > | failTrialCount_++; |
523 | > | } |
524 | > | } else { |
525 | > | std::cerr << "exchange NOT performed!\n"; |
526 | > | std::cerr << "at least one of the two slabs empty.\n"; |
527 | > | failTrialCount_++; |
528 | > | } |
529 | > | |
530 | > | } |
531 | > | |
532 | > | void RNEMD::doScale() { |
533 | > | |
534 | > | Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); |
535 | > | Mat3x3d hmat = currentSnap_->getHmat(); |
536 | > | |
537 | > | seleMan_.setSelectionSet(evaluator_.evaluate()); |
538 | > | |
539 | > | int selei; |
540 | > | StuntDouble* sd; |
541 | > | int idx; |
542 | > | |
543 | > | std::vector<StuntDouble*> hotBin, coldBin; |
544 | > | |
545 | > | RealType Phx = 0.0; |
546 | > | RealType Phy = 0.0; |
547 | > | RealType Phz = 0.0; |
548 | > | RealType Khx = 0.0; |
549 | > | RealType Khy = 0.0; |
550 | > | RealType Khz = 0.0; |
551 | > | RealType Pcx = 0.0; |
552 | > | RealType Pcy = 0.0; |
553 | > | RealType Pcz = 0.0; |
554 | > | RealType Kcx = 0.0; |
555 | > | RealType Kcy = 0.0; |
556 | > | RealType Kcz = 0.0; |
557 | > | |
558 | > | for (sd = seleMan_.beginSelected(selei); sd != NULL; |
559 | > | sd = seleMan_.nextSelected(selei)) { |
560 | > | |
561 | > | idx = sd->getLocalIndex(); |
562 | > | |
563 | > | Vector3d pos = sd->getPos(); |
564 | > | |
565 | > | // wrap the stuntdouble's position back into the box: |
566 | > | |
567 | > | if (usePeriodicBoundaryConditions_) |
568 | > | currentSnap_->wrapVector(pos); |
569 | > | |
570 | > | // which bin is this stuntdouble in? |
571 | > | // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)] |
572 | > | |
573 | > | int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_; |
574 | > | |
575 | > | // if we're in bin 0 or the middleBin |
576 | > | if (binNo == 0 || binNo == midBin_) { |
577 | > | |
578 | > | RealType mass = sd->getMass(); |
579 | > | Vector3d vel = sd->getVel(); |
580 | > | |
581 | > | if (binNo == 0) { |
582 | > | hotBin.push_back(sd); |
583 | > | Phx += mass * vel.x(); |
584 | > | Phy += mass * vel.y(); |
585 | > | Phz += mass * vel.z(); |
586 | > | Khx += mass * vel.x() * vel.x(); |
587 | > | Khy += mass * vel.y() * vel.y(); |
588 | > | Khz += mass * vel.z() * vel.z(); |
589 | > | } else { //midBin_ |
590 | > | coldBin.push_back(sd); |
591 | > | Pcx += mass * vel.x(); |
592 | > | Pcy += mass * vel.y(); |
593 | > | Pcz += mass * vel.z(); |
594 | > | Kcx += mass * vel.x() * vel.x(); |
595 | > | Kcy += mass * vel.y() * vel.y(); |
596 | > | Kcz += mass * vel.z() * vel.z(); |
597 | > | } |
598 | > | } |
599 | > | } |
600 | > | |
601 | > | Khx *= 0.5; |
602 | > | Khy *= 0.5; |
603 | > | Khz *= 0.5; |
604 | > | Kcx *= 0.5; |
605 | > | Kcy *= 0.5; |
606 | > | Kcz *= 0.5; |
607 | > | |
608 | > | #ifdef IS_MPI |
609 | > | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phx, 1, MPI::REALTYPE, MPI::SUM); |
610 | > | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phy, 1, MPI::REALTYPE, MPI::SUM); |
611 | > | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phz, 1, MPI::REALTYPE, MPI::SUM); |
612 | > | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcx, 1, MPI::REALTYPE, MPI::SUM); |
613 | > | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcy, 1, MPI::REALTYPE, MPI::SUM); |
614 | > | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcz, 1, MPI::REALTYPE, MPI::SUM); |
615 | > | |
616 | > | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khx, 1, MPI::REALTYPE, MPI::SUM); |
617 | > | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khy, 1, MPI::REALTYPE, MPI::SUM); |
618 | > | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khz, 1, MPI::REALTYPE, MPI::SUM); |
619 | > | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcx, 1, MPI::REALTYPE, MPI::SUM); |
620 | > | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcy, 1, MPI::REALTYPE, MPI::SUM); |
621 | > | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcz, 1, MPI::REALTYPE, MPI::SUM); |
622 | > | #endif |
623 | > | |
624 | > | //use coldBin coeff's |
625 | > | RealType px = Pcx / Phx; |
626 | > | RealType py = Pcy / Phy; |
627 | > | RealType pz = Pcz / Phz; |
628 | > | |
629 | > | RealType a000, a110, c0, a001, a111, b01, b11, c1, c; |
630 | > | switch(rnemdType_) { |
631 | > | case rnemdKineticScale : |
632 | > | /*used hotBin coeff's & only scale x & y dimensions |
633 | > | RealType px = Phx / Pcx; |
634 | > | RealType py = Phy / Pcy; |
635 | > | a110 = Khy; |
636 | > | c0 = - Khx - Khy - targetFlux_; |
637 | > | a000 = Khx; |
638 | > | a111 = Kcy * py * py |
639 | > | b11 = -2.0 * Kcy * py * (1.0 + py); |
640 | > | c1 = Kcy * py * (2.0 + py) + Kcx * px * ( 2.0 + px) + targetFlux_; |
641 | > | b01 = -2.0 * Kcx * px * (1.0 + px); |
642 | > | a001 = Kcx * px * px; |
643 | > | */ |
644 | > | |
645 | > | //scale all three dimensions, let x = y |
646 | > | a000 = Kcx + Kcy; |
647 | > | a110 = Kcz; |
648 | > | c0 = targetFlux_ - Kcx - Kcy - Kcz; |
649 | > | a001 = Khx * px * px + Khy * py * py; |
650 | > | a111 = Khz * pz * pz; |
651 | > | b01 = -2.0 * (Khx * px * (1.0 + px) + Khy * py * (1.0 + py)); |
652 | > | b11 = -2.0 * Khz * pz * (1.0 + pz); |
653 | > | c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py) |
654 | > | + Khz * pz * (2.0 + pz) - targetFlux_; |
655 | > | break; |
656 | > | case rnemdPxScale : |
657 | > | c = 1 - targetFlux_ / Pcx; |
658 | > | a000 = Kcy; |
659 | > | a110 = Kcz; |
660 | > | c0 = Kcx * c * c - Kcx - Kcy - Kcz; |
661 | > | a001 = py * py * Khy; |
662 | > | a111 = pz * pz * Khz; |
663 | > | b01 = -2.0 * Khy * py * (1.0 + py); |
664 | > | b11 = -2.0 * Khz * pz * (1.0 + pz); |
665 | > | c1 = Khy * py * (2.0 + py) + Khz * pz * (2.0 + pz) |
666 | > | + Khx * (fastpow(c * px - px - 1.0, 2) - 1.0); |
667 | > | break; |
668 | > | case rnemdPyScale : |
669 | > | c = 1 - targetFlux_ / Pcy; |
670 | > | a000 = Kcx; |
671 | > | a110 = Kcz; |
672 | > | c0 = Kcy * c * c - Kcx - Kcy - Kcz; |
673 | > | a001 = px * px * Khx; |
674 | > | a111 = pz * pz * Khz; |
675 | > | b01 = -2.0 * Khx * px * (1.0 + px); |
676 | > | b11 = -2.0 * Khz * pz * (1.0 + pz); |
677 | > | c1 = Khx * px * (2.0 + px) + Khz * pz * (2.0 + pz) |
678 | > | + Khy * (fastpow(c * py - py - 1.0, 2) - 1.0); |
679 | > | break; |
680 | > | case rnemdPzScale ://we don't really do this, do we? |
681 | > | default : |
682 | > | break; |
683 | > | } |
684 | > | |
685 | > | RealType v1 = a000 * a111 - a001 * a110; |
686 | > | RealType v2 = a000 * b01; |
687 | > | RealType v3 = a000 * b11; |
688 | > | RealType v4 = a000 * c1 - a001 * c0; |
689 | > | RealType v8 = a110 * b01; |
690 | > | RealType v10 = - b01 * c0; |
691 | > | |
692 | > | RealType u0 = v2 * v10 - v4 * v4; |
693 | > | RealType u1 = -2.0 * v3 * v4; |
694 | > | RealType u2 = -v2 * v8 - v3 * v3 - 2.0 * v1 * v4; |
695 | > | RealType u3 = -2.0 * v1 * v3; |
696 | > | RealType u4 = - v1 * v1; |
697 | > | //rescale coefficients |
698 | > | RealType maxAbs = fabs(u0); |
699 | > | if (maxAbs < fabs(u1)) maxAbs = fabs(u1); |
700 | > | if (maxAbs < fabs(u2)) maxAbs = fabs(u2); |
701 | > | if (maxAbs < fabs(u3)) maxAbs = fabs(u3); |
702 | > | if (maxAbs < fabs(u4)) maxAbs = fabs(u4); |
703 | > | u0 /= maxAbs; |
704 | > | u1 /= maxAbs; |
705 | > | u2 /= maxAbs; |
706 | > | u3 /= maxAbs; |
707 | > | u4 /= maxAbs; |
708 | > | //max_element(start, end) is also available. |
709 | > | Polynomial<RealType> poly; //same as DoublePolynomial poly; |
710 | > | poly.setCoefficient(4, u4); |
711 | > | poly.setCoefficient(3, u3); |
712 | > | poly.setCoefficient(2, u2); |
713 | > | poly.setCoefficient(1, u1); |
714 | > | poly.setCoefficient(0, u0); |
715 | > | std::vector<RealType> realRoots = poly.FindRealRoots(); |
716 | > | |
717 | > | std::vector<RealType>::iterator ri; |
718 | > | RealType r1, r2, alpha0; |
719 | > | std::vector<std::pair<RealType,RealType> > rps; |
720 | > | for (ri = realRoots.begin(); ri !=realRoots.end(); ri++) { |
721 | > | r2 = *ri; |
722 | > | //check if FindRealRoots() give the right answer |
723 | > | if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) { |
724 | > | std::cerr << "WARNING! eq solvers might have mistakes!\n"; |
725 | > | failRootCount_++; |
726 | > | } |
727 | > | //might not be useful w/o rescaling coefficients |
728 | > | alpha0 = -c0 - a110 * r2 * r2; |
729 | > | if (alpha0 >= 0.0) { |
730 | > | r1 = sqrt(alpha0 / a000); |
731 | > | if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) < 1e-6) |
732 | > | { rps.push_back(std::make_pair(r1, r2)); } |
733 | > | if (r1 > 1e-6) { //r1 non-negative |
734 | > | r1 = -r1; |
735 | > | if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) <1e-6) |
736 | > | { rps.push_back(std::make_pair(r1, r2)); } |
737 | > | } |
738 | > | } |
739 | > | } |
740 | > | // Consider combininig together the solving pair part w/ the searching |
741 | > | // best solution part so that we don't need the pairs vector |
742 | > | if (!rps.empty()) { |
743 | > | RealType smallestDiff = HONKING_LARGE_VALUE; |
744 | > | RealType diff; |
745 | > | std::pair<RealType,RealType> bestPair = std::make_pair(1.0, 1.0); |
746 | > | std::vector<std::pair<RealType,RealType> >::iterator rpi; |
747 | > | for (rpi = rps.begin(); rpi != rps.end(); rpi++) { |
748 | > | r1 = (*rpi).first; |
749 | > | r2 = (*rpi).second; |
750 | > | switch(rnemdType_) { |
751 | > | case rnemdKineticScale : |
752 | > | diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2) |
753 | > | + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2) |
754 | > | + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcy, 2); |
755 | > | break; |
756 | > | case rnemdPxScale : |
757 | > | diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2) |
758 | > | + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcy, 2); |
759 | > | break; |
760 | > | case rnemdPyScale : |
761 | > | diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2) |
762 | > | + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2); |
763 | > | break; |
764 | > | case rnemdPzScale : |
765 | > | default : |
766 | > | break; |
767 | > | } |
768 | > | if (diff < smallestDiff) { |
769 | > | smallestDiff = diff; |
770 | > | bestPair = *rpi; |
771 | > | } |
772 | > | } |
773 | > | #ifdef IS_MPI |
774 | > | if (worldRank == 0) { |
775 | > | #endif |
776 | > | std::cerr << "we choose r1 = " << bestPair.first |
777 | > | << " and r2 = " << bestPair.second << "\n"; |
778 | > | #ifdef IS_MPI |
779 | > | } |
780 | > | #endif |
781 | > | |
782 | > | RealType x, y, z; |
783 | > | switch(rnemdType_) { |
784 | > | case rnemdKineticScale : |
785 | > | x = bestPair.first; |
786 | > | y = bestPair.first; |
787 | > | z = bestPair.second; |
788 | > | break; |
789 | > | case rnemdPxScale : |
790 | > | x = c; |
791 | > | y = bestPair.first; |
792 | > | z = bestPair.second; |
793 | > | break; |
794 | > | case rnemdPyScale : |
795 | > | x = bestPair.first; |
796 | > | y = c; |
797 | > | z = bestPair.second; |
798 | > | break; |
799 | > | case rnemdPzScale : |
800 | > | default : |
801 | > | break; |
802 | > | } |
803 | > | std::vector<StuntDouble*>::iterator sdi; |
804 | > | Vector3d vel; |
805 | > | for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) { |
806 | > | vel = (*sdi)->getVel(); |
807 | > | vel.x() *= x; |
808 | > | vel.y() *= y; |
809 | > | vel.z() *= z; |
810 | > | (*sdi)->setVel(vel); |
811 | > | } |
812 | > | //convert to hotBin coefficient |
813 | > | x = 1.0 + px * (1.0 - x); |
814 | > | y = 1.0 + py * (1.0 - y); |
815 | > | z = 1.0 + pz * (1.0 - z); |
816 | > | for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) { |
817 | > | vel = (*sdi)->getVel(); |
818 | > | vel.x() *= x; |
819 | > | vel.y() *= y; |
820 | > | vel.z() *= z; |
821 | > | (*sdi)->setVel(vel); |
822 | > | } |
823 | > | exchangeSum_ += targetFlux_; |
824 | > | //we may want to check whether the exchange has been successful |
825 | > | } else { |
826 | > | std::cerr << "exchange NOT performed!\n"; |
827 | > | failTrialCount_++; |
828 | > | } |
829 | > | |
830 | > | } |
831 | > | |
832 | > | void RNEMD::doRNEMD() { |
833 | > | |
834 | > | switch(rnemdType_) { |
835 | > | case rnemdKineticScale : |
836 | > | case rnemdPxScale : |
837 | > | case rnemdPyScale : |
838 | > | case rnemdPzScale : |
839 | > | doScale(); |
840 | > | break; |
841 | > | case rnemdKineticSwap : |
842 | > | case rnemdPx : |
843 | > | case rnemdPy : |
844 | > | case rnemdPz : |
845 | > | doSwap(); |
846 | > | break; |
847 | > | case rnemdUnknown : |
848 | > | default : |
849 | > | break; |
850 | > | } |
851 | > | } |
852 | > | |
853 | > | void RNEMD::collectData() { |
854 | > | |
855 | > | Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); |
856 | > | Mat3x3d hmat = currentSnap_->getHmat(); |
857 | > | |
858 | > | seleMan_.setSelectionSet(evaluator_.evaluate()); |
859 | > | |
860 | > | int selei; |
861 | > | StuntDouble* sd; |
862 | > | int idx; |
863 | > | |
864 | > | for (sd = seleMan_.beginSelected(selei); sd != NULL; |
865 | > | sd = seleMan_.nextSelected(selei)) { |
866 | > | |
867 | > | idx = sd->getLocalIndex(); |
868 | > | |
869 | > | Vector3d pos = sd->getPos(); |
870 | > | |
871 | > | // wrap the stuntdouble's position back into the box: |
872 | > | |
873 | > | if (usePeriodicBoundaryConditions_) |
874 | > | currentSnap_->wrapVector(pos); |
875 | > | |
876 | > | // which bin is this stuntdouble in? |
877 | > | // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)] |
878 | > | |
879 | > | int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_; |
880 | > | |
881 | > | if (rnemdLogWidth_ == midBin_ + 1) |
882 | > | if (binNo > midBin_) |
883 | > | binNo = nBins_ - binNo; |
884 | > | |
885 | > | RealType mass = sd->getMass(); |
886 | > | Vector3d vel = sd->getVel(); |
887 | > | RealType value; |
888 | > | RealType xVal, yVal, zVal; |
889 | > | |
890 | > | switch(rnemdType_) { |
891 | > | case rnemdKineticSwap : |
892 | > | case rnemdKineticScale : |
893 | > | |
894 | > | value = mass * (vel[0]*vel[0] + vel[1]*vel[1] + |
895 | > | vel[2]*vel[2]); |
896 | > | |
897 | > | valueCount_[binNo] += 3; |
898 | > | if (sd->isDirectional()) { |
899 | > | Vector3d angMom = sd->getJ(); |
900 | > | Mat3x3d I = sd->getI(); |
901 | > | |
902 | > | if (sd->isLinear()) { |
903 | > | int i = sd->linearAxis(); |
904 | > | int j = (i + 1) % 3; |
905 | > | int k = (i + 2) % 3; |
906 | > | value += angMom[j] * angMom[j] / I(j, j) + |
907 | > | angMom[k] * angMom[k] / I(k, k); |
908 | > | |
909 | > | valueCount_[binNo] +=2; |
910 | > | |
911 | > | } else { |
912 | > | value += angMom[0]*angMom[0]/I(0, 0) |
913 | > | + angMom[1]*angMom[1]/I(1, 1) |
914 | > | + angMom[2]*angMom[2]/I(2, 2); |
915 | > | valueCount_[binNo] +=3; |
916 | > | } |
917 | > | } |
918 | > | value = value / OOPSEConstant::energyConvert / OOPSEConstant::kb; |
919 | > | |
920 | > | break; |
921 | > | case rnemdPx : |
922 | > | case rnemdPxScale : |
923 | > | value = mass * vel[0]; |
924 | > | valueCount_[binNo]++; |
925 | > | xVal = mass * vel.x() * vel.x() / OOPSEConstant::energyConvert |
926 | > | / OOPSEConstant::kb; |
927 | > | yVal = mass * vel.y() * vel.y() / OOPSEConstant::energyConvert |
928 | > | / OOPSEConstant::kb; |
929 | > | zVal = mass * vel.z() * vel.z() / OOPSEConstant::energyConvert |
930 | > | / OOPSEConstant::kb; |
931 | > | xTempHist_[binNo] += xVal; |
932 | > | yTempHist_[binNo] += yVal; |
933 | > | zTempHist_[binNo] += zVal; |
934 | > | break; |
935 | > | case rnemdPy : |
936 | > | case rnemdPyScale : |
937 | > | value = mass * vel[1]; |
938 | > | valueCount_[binNo]++; |
939 | > | break; |
940 | > | case rnemdPz : |
941 | > | case rnemdPzScale : |
942 | > | value = mass * vel[2]; |
943 | > | valueCount_[binNo]++; |
944 | > | break; |
945 | > | case rnemdUnknown : |
946 | > | default : |
947 | > | break; |
948 | > | } |
949 | > | valueHist_[binNo] += value; |
950 | > | } |
951 | > | |
952 | > | } |
953 | > | |
954 | > | void RNEMD::getStarted() { |
955 | > | Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); |
956 | > | Stats& stat = currentSnap_->statData; |
957 | > | stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_; |
958 | > | } |
959 | > | |
960 | > | void RNEMD::getStatus() { |
961 | > | |
962 | > | Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); |
963 | > | Stats& stat = currentSnap_->statData; |
964 | > | RealType time = currentSnap_->getTime(); |
965 | > | |
966 | > | stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_; |
967 | > | //or to be more meaningful, define another item as exchangeSum_ / time |
968 | > | |
969 | > | |
970 | > | #ifdef IS_MPI |
971 | > | |
972 | > | // all processors have the same number of bins, and STL vectors pack their |
973 | > | // arrays, so in theory, this should be safe: |
974 | > | |
975 | > | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueHist_[0], |
976 | > | rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); |
977 | > | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount_[0], |
978 | > | rnemdLogWidth_, MPI::INT, MPI::SUM); |
979 | > | |
980 | > | // If we're the root node, should we print out the results |
981 | > | int worldRank = MPI::COMM_WORLD.Get_rank(); |
982 | > | if (worldRank == 0) { |
983 | > | #endif |
984 | > | int j; |
985 | > | rnemdLog_ << time; |
986 | > | for (j = 0; j < rnemdLogWidth_; j++) { |
987 | > | rnemdLog_ << "\t" << valueHist_[j] / (RealType)valueCount_[j]; |
988 | > | valueHist_[j] = 0.0; |
989 | > | } |
990 | > | rnemdLog_ << "\n"; |
991 | > | if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale ) { |
992 | > | xTempLog_ << time; |
993 | > | for (j = 0; j < rnemdLogWidth_; j++) { |
994 | > | xTempLog_ << "\t" << xTempHist_[j] / (RealType)valueCount_[j]; |
995 | > | xTempHist_[j] = 0.0; |
996 | > | } |
997 | > | xTempLog_ << "\n"; |
998 | > | yTempLog_ << time; |
999 | > | for (j = 0; j < rnemdLogWidth_; j++) { |
1000 | > | yTempLog_ << "\t" << yTempHist_[j] / (RealType)valueCount_[j]; |
1001 | > | yTempHist_[j] = 0.0; |
1002 | > | } |
1003 | > | yTempLog_ << "\n"; |
1004 | > | zTempLog_ << time; |
1005 | > | for (j = 0; j < rnemdLogWidth_; j++) { |
1006 | > | zTempLog_ << "\t" << zTempHist_[j] / (RealType)valueCount_[j]; |
1007 | > | zTempHist_[j] = 0.0; |
1008 | > | } |
1009 | > | zTempLog_ << "\n"; |
1010 | > | } |
1011 | > | for (j = 0; j < rnemdLogWidth_; j++) valueCount_[j] = 0; |
1012 | > | #ifdef IS_MPI |
1013 | > | } |
1014 | > | #endif |
1015 | > | } |
1016 | > | |
1017 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |