# | Line 6 | Line 6 | |
---|---|---|
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 |
9 | > | * 1. Redistributions of source code must retain the above copyright |
10 | * notice, this list of conditions and the following disclaimer. | |
11 | * | |
12 | < | * 3. Redistributions in binary form must reproduce the above copyright |
12 | > | * 2. Redistributions in binary form must reproduce the above copyright |
13 | * notice, this list of conditions and the following disclaimer in the | |
14 | * documentation and/or other materials provided with the | |
15 | * distribution. | |
# | Line 37 | Line 28 | |
28 | * arising out of the use of or inability to use software, even if the | |
29 | * University of Notre Dame has been advised of the possibility of | |
30 | * such damages. | |
31 | + | * |
32 | + | * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your |
33 | + | * research, please cite the appropriate papers when you publish your |
34 | + | * work. Good starting points are: |
35 | + | * |
36 | + | * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). |
37 | + | * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). |
38 | + | * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). |
39 | + | * [4] Vardeman & Gezelter, in progress (2009). |
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" |
49 | > | #include "utils/PhysicalConstants.hpp" |
50 | #include "utils/Tuple.hpp" | |
51 | ||
52 | #ifndef IS_MPI | |
# | Line 55 | Line 57 | |
57 | ||
58 | #define HONKING_LARGE_VALUE 1.0e10 | |
59 | ||
60 | < | namespace oopse { |
60 | > | namespace OpenMD { |
61 | ||
62 | RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) { | |
63 | < | |
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; | |
# | Line 72 | Line 81 | namespace oopse { | |
81 | evaluator_.loadScriptString(rnemdObjectSelection_); | |
82 | seleMan_.setSelectionSet(evaluator_.evaluate()); | |
83 | ||
75 | – | |
84 | // do some sanity checking | |
85 | ||
86 | int selectionCount = seleMan_.getSelectionCount(); | |
# | Line 94 | Line 102 | namespace oopse { | |
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; |
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 123 | 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() { | |
129 | – | int midBin = nBins_ / 2; |
216 | ||
217 | Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); | |
218 | Mat3x3d hmat = currentSnap_->getHmat(); | |
# | Line 164 | Line 250 | namespace oopse { | |
250 | ||
251 | ||
252 | // if we're in bin 0 or the middleBin | |
253 | < | if (binNo == 0 || binNo == midBin) { |
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 rnemdKinetic : |
260 | > | case rnemdKineticSwap : |
261 | ||
262 | value = mass * (vel[0]*vel[0] + vel[1]*vel[1] + | |
263 | vel[2]*vel[2]); | |
# | Line 191 | Line 277 | namespace oopse { | |
277 | + angMom[2]*angMom[2]/I(2, 2); | |
278 | } | |
279 | } | |
280 | < | value = value * 0.5 / OOPSEConstant::energyConvert; |
280 | > | //make exchangeSum_ comparable between swap & scale |
281 | > | //temporarily without using energyConvert |
282 | > | //value = value * 0.5 / PhysicalConstants::energyConvert; |
283 | > | value *= 0.5; |
284 | break; | |
285 | case rnemdPx : | |
286 | value = mass * vel[0]; | |
# | Line 202 | Line 291 | namespace oopse { | |
291 | case rnemdPz : | |
292 | value = mass * vel[2]; | |
293 | break; | |
205 | – | case rnemdUnknown : |
294 | default : | |
295 | break; | |
296 | } | |
# | Line 218 | Line 306 | namespace oopse { | |
306 | min_sd = sd; | |
307 | } | |
308 | } | |
309 | < | } else { |
309 | > | } else { //midBin_ |
310 | if (!max_found) { | |
311 | max_val = value; | |
312 | max_sd = sd; | |
# | Line 298 | Line 386 | namespace oopse { | |
386 | RealType temp_vel; | |
387 | ||
388 | switch(rnemdType_) { | |
389 | < | case rnemdKinetic : |
389 | > | case rnemdKineticSwap : |
390 | min_sd->setVel(max_vel); | |
391 | max_sd->setVel(min_vel); | |
392 | if (min_sd->isDirectional() && max_sd->isDirectional()) { | |
# | Line 329 | Line 417 | namespace oopse { | |
417 | min_sd->setVel(min_vel); | |
418 | max_sd->setVel(max_vel); | |
419 | break; | |
332 | – | case rnemdUnknown : |
420 | default : | |
421 | break; | |
422 | } | |
# | Line 349 | Line 436 | namespace oopse { | |
436 | min_vals.rank, 0, status); | |
437 | ||
438 | switch(rnemdType_) { | |
439 | < | case rnemdKinetic : |
439 | > | case rnemdKineticSwap : |
440 | max_sd->setVel(min_vel); | |
441 | ||
442 | if (max_sd->isDirectional()) { | |
# | Line 378 | Line 465 | namespace oopse { | |
465 | max_vel.z() = min_vel.z(); | |
466 | max_sd->setVel(max_vel); | |
467 | break; | |
381 | – | case rnemdUnknown : |
468 | default : | |
469 | break; | |
470 | } | |
# | Line 396 | Line 482 | namespace oopse { | |
482 | max_vals.rank, 0, status); | |
483 | ||
484 | switch(rnemdType_) { | |
485 | < | case rnemdKinetic : |
485 | > | case rnemdKineticSwap : |
486 | min_sd->setVel(max_vel); | |
487 | ||
488 | if (min_sd->isDirectional()) { | |
# | Line 425 | Line 511 | namespace oopse { | |
511 | min_vel.z() = max_vel.z(); | |
512 | min_sd->setVel(min_vel); | |
513 | break; | |
428 | – | case rnemdUnknown : |
514 | default : | |
515 | break; | |
516 | } | |
# | Line 433 | Line 518 | namespace oopse { | |
518 | #endif | |
519 | exchangeSum_ += max_val - min_val; | |
520 | } else { | |
521 | < | std::cerr << "exchange NOT performed.\nmin_val > max_val.\n"; |
521 | > | std::cerr << "exchange NOT performed!\nmin_val > max_val.\n"; |
522 | > | failTrialCount_++; |
523 | } | |
524 | } else { | |
525 | < | std::cerr << "exchange NOT performed.\none of the two slabs empty.\n"; |
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::getStatus() { |
532 | > | void RNEMD::doScale() { |
533 | ||
534 | Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); | |
535 | Mat3x3d hmat = currentSnap_->getHmat(); | |
448 | – | Stats& stat = currentSnap_->statData; |
449 | – | RealType time = currentSnap_->getTime(); |
536 | ||
451 | – | stat[Stats::RNEMD_SWAP_TOTAL] = exchangeSum_; |
452 | – | |
537 | seleMan_.setSelectionSet(evaluator_.evaluate()); | |
538 | ||
539 | int selei; | |
540 | StuntDouble* sd; | |
541 | int idx; | |
542 | ||
543 | < | std::vector<RealType> valueHist(nBins_, 0.0); // keeps track of what's |
544 | < | // being averaged |
545 | < | std::vector<int> valueCount(nBins_, 0); // keeps track of the |
546 | < | // number of degrees of |
547 | < | // freedom being averaged |
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 | ||
# | Line 477 | Line 876 | namespace oopse { | |
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 | < | |
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 rnemdKinetic : |
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; |
897 | > | valueCount_[binNo] += 3; |
898 | if (sd->isDirectional()) { | |
899 | Vector3d angMom = sd->getJ(); | |
900 | Mat3x3d I = sd->getI(); | |
# | Line 501 | Line 906 | namespace oopse { | |
906 | value += angMom[j] * angMom[j] / I(j, j) + | |
907 | angMom[k] * angMom[k] / I(k, k); | |
908 | ||
909 | < | valueCount[binNo] +=2; |
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; |
915 | > | valueCount_[binNo] +=3; |
916 | } | |
917 | } | |
918 | < | value = value / OOPSEConstant::energyConvert / OOPSEConstant::kb; |
918 | > | value = value / PhysicalConstants::energyConvert / PhysicalConstants::kb; |
919 | ||
920 | break; | |
921 | case rnemdPx : | |
922 | + | case rnemdPxScale : |
923 | value = mass * vel[0]; | |
924 | < | valueCount[binNo]++; |
924 | > | valueCount_[binNo]++; |
925 | > | xVal = mass * vel.x() * vel.x() / PhysicalConstants::energyConvert |
926 | > | / PhysicalConstants::kb; |
927 | > | yVal = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert |
928 | > | / PhysicalConstants::kb; |
929 | > | zVal = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert |
930 | > | / PhysicalConstants::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]++; |
938 | > | valueCount_[binNo]++; |
939 | break; | |
940 | case rnemdPz : | |
941 | + | case rnemdPzScale : |
942 | value = mass * vel[2]; | |
943 | < | valueCount[binNo]++; |
943 | > | valueCount_[binNo]++; |
944 | break; | |
945 | case rnemdUnknown : | |
946 | default : | |
947 | break; | |
948 | } | |
949 | < | valueHist[binNo] += value; |
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 | < | nBins_, MPI::REALTYPE, MPI::SUM); |
977 | < | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount[0], |
978 | < | nBins_, MPI::INT, MPI::SUM); |
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 | < | |
985 | < | std::cout << time; |
986 | < | for (int j = 0; j < nBins_; j++) |
987 | < | std::cout << "\t" << valueHist[j] / (RealType)valueCount[j]; |
988 | < | std::cout << "\n"; |
989 | < | |
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 | < | } |
1013 | > | } |
1014 | #endif | |
1015 | + | |
1016 | + | |
1017 | } | |
1018 | + | |
1019 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |