ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/rnemd/RNEMD.cpp
(Generate patch)

Comparing trunk/src/integrators/RNEMD.cpp (file contents):
Revision 1330 by skuang, Thu Mar 19 21:03:36 2009 UTC vs.
Revision 1368 by skuang, Mon Oct 19 13:39:04 2009 UTC

# 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) {
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 <    const std::string st = simParams->getRNEMD_swapType();
80 >    rnemdObjectSelection_ = simParams->getRNEMD_objectSelection();
81 >    evaluator_.loadScriptString(rnemdObjectSelection_);
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_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 + #ifdef IS_MPI
115 +    if (worldRank == 0) {
116 + #endif
117  
118 <    set_RNEMD_swapTime(simParams->getRNEMD_swapTime());
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 100 | 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() {
216 <    std::cerr << "in RNEMD!\n";  
217 <    std::cerr << "nBins = " << nBins_ << "\n";
218 <    std::cerr << "swapTime = " << swapTime_ << "\n";
219 <    std::cerr << "exchangeSum = " << exchangeSum_ << "\n";
220 <    std::cerr << "swapType = " << rnemdType_ << "\n";
221 <  }  
216 >
217 >    Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
218 >    Mat3x3d hmat = currentSnap_->getHmat();
219 >
220 >    seleMan_.setSelectionSet(evaluator_.evaluate());
221 >
222 >    int selei;
223 >    StuntDouble* sd;
224 >    int idx;
225 >
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 >
241 >      // wrap the stuntdouble's position back into the box:
242 >
243 >      if (usePeriodicBoundaryConditions_)
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 > #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   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines