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 1341 by skuang, Fri May 8 19:47:05 2009 UTC vs.
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC

# 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 53 | Line 55
55   #include "math/ParallelRandNumGen.hpp"
56   #endif
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 71 | Line 81 | namespace oopse {
81      evaluator_.loadScriptString(rnemdObjectSelection_);
82      seleMan_.setSelectionSet(evaluator_.evaluate());
83  
74 –
84      // do some sanity checking
85  
86      int selectionCount = seleMan_.getSelectionCount();
# Line 93 | 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;
155 <    counter_ = 0; //added by shenyu
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 233 | Line 321 | namespace oopse {
321        }
322      }
323  
324 <    // missing:  swap information in parallel
324 > #ifdef IS_MPI
325 >    int nProc, worldRank;
326  
327 <    if (max_found && min_found) {
328 <      if (min_val< max_val) {
327 >    nProc = MPI::COMM_WORLD.Get_size();
328 >    worldRank = MPI::COMM_WORLD.Get_rank();
329  
330 <        Vector3d min_vel = min_sd->getVel();
331 <        Vector3d max_vel = max_sd->getVel();
243 <        RealType temp_vel;
330 >    bool my_min_found = min_found;
331 >    bool my_max_found = max_found;
332  
333 <        switch(rnemdType_) {
334 <        case rnemdKinetic :
335 <          min_sd->setVel(max_vel);
336 <          max_sd->setVel(min_vel);
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 (min_sd->isDirectional() && max_sd->isDirectional()) {
376 <            Vector3d min_angMom = min_sd->getJ();
377 <            Vector3d max_angMom = max_sd->getJ();
378 <            min_sd->setJ(max_angMom);
379 <            max_sd->setJ(min_angMom);
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 <          break;
424 <        case rnemdPx :
425 <          temp_vel = min_vel.x();
426 <          min_vel.x() = max_vel.x();
427 <          max_vel.x() = temp_vel;
428 <          min_sd->setVel(min_vel);
429 <          max_sd->setVel(max_vel);
430 <          break;
431 <        case rnemdPy :
432 <          temp_vel = min_vel.y();
433 <          min_vel.y() = max_vel.y();
434 <          max_vel.y() = temp_vel;
435 <          min_sd->setVel(min_vel);
436 <          max_sd->setVel(max_vel);
437 <          break;
438 <        case rnemdPz :
439 <          temp_vel = min_vel.z();
440 <          min_vel.z() = max_vel.z();
441 <          max_vel.z() = temp_vel;
442 <          min_sd->setVel(min_vel);
443 <          max_sd->setVel(max_vel);
444 <          break;
445 <        case rnemdUnknown :
446 <        default :
447 <          break;
448 <        }
449 <      exchangeSum_ += max_val - min_val;
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";
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::doScale() {
533  
291 –  void RNEMD::getStatus() {
292 –
534      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
535      Mat3x3d hmat = currentSnap_->getHmat();
295 –    Stats& stat = currentSnap_->statData;
536  
537 <    stat[Stats::RNEMD_SWAP_TOTAL] = exchangeSum_;
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  
305 –    std::vector<RealType> valueHist;  // keeps track of what's being averaged
306 –    std::vector<int> valueCount; // keeps track of the number of degrees of
307 –                                 // freedom being averaged
308 –    valueHist.resize(nBins_);
309 –    valueCount.resize(nBins_);
310 –    //do they initialize themselves to zero automatically?
864      for (sd = seleMan_.beginSelected(selei); sd != NULL;
865           sd = seleMan_.nextSelected(selei)) {
866        
# Line 324 | Line 877 | namespace oopse {
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 <      //std::cerr << "pos.z() = " << pos.z() << " bin = " << binNo << "\n";
882 <      
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();
332 –      //std::cerr << "mass = " << mass << " vel = " << vel << "\n";
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;
342 <        //std::cerr <<"starting value = " << value << "\n";
897 >        valueCount_[binNo] += 3;
898          if (sd->isDirectional()) {
344 –          //std::cerr << "angMom calculated.\n";
899            Vector3d angMom = sd->getJ();
346 –          //std::cerr << "current angMom: " << angMom << "\n";
900            Mat3x3d I = sd->getI();
901            
902            if (sd->isLinear()) {
# Line 353 | 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 {
359 –            //std::cerr << "non-linear molecule.\n";
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;
364 <
915 >            valueCount_[binNo] +=3;
916            }
917          }
918 <        //std::cerr <<"total value = " << value << "\n";
919 <        //value *= 0.5 / OOPSEConstant::energyConvert;  // get it in kcal / mol
369 <        //value *= 2.0 / OOPSEConstant::kb;            // convert to temperature
370 <        value = value / OOPSEConstant::energyConvert / OOPSEConstant::kb;
371 <        //std::cerr <<"value = " << value << "\n";
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 <      //std::cerr << "bin = " << binNo << " value = " << value ;
390 <      valueHist[binNo] += value;
391 <      //std::cerr << " hist = " << valueHist[binNo] << " count = " << valueCount[binNo] << "\n";
949 >      valueHist_[binNo] += value;
950      }
951 <    
394 <    std::cout << counter_++;
395 <    for (int j = 0; j < nBins_; j++)
396 <      std::cout << "\t" << valueHist[j] / (RealType)valueCount[j];
397 <    std::cout << "\n";
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 +  }
1018 +
1019   }

Diff Legend

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