ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/integrators/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 1396 by gezelter, Sat Dec 5 02:57:05 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 + #ifdef IS_MPI
199 +    if (worldRank == 0) {
200 + #endif
201 +      std::cerr << "total fail trials: " << failTrialCount_ << "\n";
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 c_x = c_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 +      c = 1 - targetFlux_ / Pcz;
682 +      a000 = Kcx;
683 +      a110 = Kcy;
684 +      c0 = Kcz * c * c - Kcx - Kcy - Kcz;
685 +      a001 = px * px * Khx;
686 +      a111 = py * py * Khy;
687 +      b01 = -2.0 * Khx * px * (1.0 + px);
688 +      b11 = -2.0 * Khy * py * (1.0 + py);
689 +      c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py)
690 +        + Khz * (fastpow(c * pz - pz - 1.0, 2) - 1.0);
691 +      break;      
692 +    default :
693 +      break;
694 +    }
695 +
696 +    RealType v1 = a000 * a111 - a001 * a110;
697 +    RealType v2 = a000 * b01;
698 +    RealType v3 = a000 * b11;
699 +    RealType v4 = a000 * c1 - a001 * c0;
700 +    RealType v8 = a110 * b01;
701 +    RealType v10 = - b01 * c0;
702 +
703 +    RealType u0 = v2 * v10 - v4 * v4;
704 +    RealType u1 = -2.0 * v3 * v4;
705 +    RealType u2 = -v2 * v8 - v3 * v3 - 2.0 * v1 * v4;
706 +    RealType u3 = -2.0 * v1 * v3;
707 +    RealType u4 = - v1 * v1;
708 +    //rescale coefficients
709 +    RealType maxAbs = fabs(u0);
710 +    if (maxAbs < fabs(u1)) maxAbs = fabs(u1);
711 +    if (maxAbs < fabs(u2)) maxAbs = fabs(u2);
712 +    if (maxAbs < fabs(u3)) maxAbs = fabs(u3);
713 +    if (maxAbs < fabs(u4)) maxAbs = fabs(u4);
714 +    u0 /= maxAbs;
715 +    u1 /= maxAbs;
716 +    u2 /= maxAbs;
717 +    u3 /= maxAbs;
718 +    u4 /= maxAbs;
719 +    //max_element(start, end) is also available.
720 +    Polynomial<RealType> poly; //same as DoublePolynomial poly;
721 +    poly.setCoefficient(4, u4);
722 +    poly.setCoefficient(3, u3);
723 +    poly.setCoefficient(2, u2);
724 +    poly.setCoefficient(1, u1);
725 +    poly.setCoefficient(0, u0);
726 +    std::vector<RealType> realRoots = poly.FindRealRoots();
727 +
728 +    std::vector<RealType>::iterator ri;
729 +    RealType r1, r2, alpha0;
730 +    std::vector<std::pair<RealType,RealType> > rps;
731 +    for (ri = realRoots.begin(); ri !=realRoots.end(); ri++) {
732 +      r2 = *ri;
733 +      //check if FindRealRoots() give the right answer
734 +      if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) {
735 +        sprintf(painCave.errMsg,
736 +                "RNEMD Warning: polynomial solve seems to have an error!");
737 +        painCave.isFatal = 0;
738 +        simError();
739 +        failRootCount_++;
740 +      }
741 +      //might not be useful w/o rescaling coefficients
742 +      alpha0 = -c0 - a110 * r2 * r2;
743 +      if (alpha0 >= 0.0) {
744 +        r1 = sqrt(alpha0 / a000);
745 +        if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) < 1e-6)
746 +          { rps.push_back(std::make_pair(r1, r2)); }
747 +        if (r1 > 1e-6) { //r1 non-negative
748 +          r1 = -r1;
749 +          if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) <1e-6)
750 +            { rps.push_back(std::make_pair(r1, r2)); }
751 +        }
752 +      }
753 +    }
754 +    // Consider combininig together the solving pair part w/ the searching
755 +    // best solution part so that we don't need the pairs vector
756 +    if (!rps.empty()) {
757 +      RealType smallestDiff = HONKING_LARGE_VALUE;
758 +      RealType diff;
759 +      std::pair<RealType,RealType> bestPair = std::make_pair(1.0, 1.0);
760 +      std::vector<std::pair<RealType,RealType> >::iterator rpi;
761 +      for (rpi = rps.begin(); rpi != rps.end(); rpi++) {
762 +        r1 = (*rpi).first;
763 +        r2 = (*rpi).second;
764 +        switch(rnemdType_) {
765 +        case rnemdKineticScale :
766 +          diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
767 +            + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2)
768 +            + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcy, 2);
769 +          break;
770 +        case rnemdPxScale :
771 +          diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
772 +            + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcy, 2);
773 +          break;
774 +        case rnemdPyScale :
775 +          diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
776 +            + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2);
777 +          break;
778 +        case rnemdPzScale :
779 +        default :
780 +          break;
781 +        }
782 +        if (diff < smallestDiff) {
783 +          smallestDiff = diff;
784 +          bestPair = *rpi;
785 +        }
786 +      }
787 + #ifdef IS_MPI
788 +      if (worldRank == 0) {
789 + #endif
790 +        std::cerr << "we choose r1 = " << bestPair.first
791 +                  << " and r2 = " << bestPair.second << "\n";
792 + #ifdef IS_MPI
793 +      }
794 + #endif
795 +
796 +      RealType x, y, z;
797 +        switch(rnemdType_) {
798 +        case rnemdKineticScale :
799 +          x = bestPair.first;
800 +          y = bestPair.first;
801 +          z = bestPair.second;
802 +          break;
803 +        case rnemdPxScale :
804 +          x = c;
805 +          y = bestPair.first;
806 +          z = bestPair.second;
807 +          break;
808 +        case rnemdPyScale :
809 +          x = bestPair.first;
810 +          y = c;
811 +          z = bestPair.second;
812 +          break;
813 +        case rnemdPzScale :
814 +          x = bestPair.first;
815 +          y = bestPair.second;
816 +          z = c;
817 +          break;          
818 +        default :
819 +          break;
820 +        }
821 +      std::vector<StuntDouble*>::iterator sdi;
822 +      Vector3d vel;
823 +      for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
824 +        vel = (*sdi)->getVel();
825 +        vel.x() *= x;
826 +        vel.y() *= y;
827 +        vel.z() *= z;
828 +        (*sdi)->setVel(vel);
829 +      }
830 +      //convert to hotBin coefficient
831 +      x = 1.0 + px * (1.0 - x);
832 +      y = 1.0 + py * (1.0 - y);
833 +      z = 1.0 + pz * (1.0 - z);
834 +      for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
835 +        vel = (*sdi)->getVel();
836 +        vel.x() *= x;
837 +        vel.y() *= y;
838 +        vel.z() *= z;
839 +        (*sdi)->setVel(vel);
840 +      }
841 +      exchangeSum_ += targetFlux_;
842 +      //we may want to check whether the exchange has been successful
843 +    } else {
844 +      std::cerr << "exchange NOT performed!\n";//MPI incompatible
845 +      failTrialCount_++;
846 +    }
847 +
848 +  }
849 +
850 +  void RNEMD::doRNEMD() {
851 +
852 +    switch(rnemdType_) {
853 +    case rnemdKineticScale :
854 +    case rnemdPxScale :
855 +    case rnemdPyScale :
856 +    case rnemdPzScale :
857 +      doScale();
858 +      break;
859 +    case rnemdKineticSwap :
860 +    case rnemdPx :
861 +    case rnemdPy :
862 +    case rnemdPz :
863 +      doSwap();
864 +      break;
865 +    case rnemdUnknown :
866 +    default :
867 +      break;
868 +    }
869 +  }
870 +
871 +  void RNEMD::collectData() {
872 +
873 +    Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
874 +    Mat3x3d hmat = currentSnap_->getHmat();
875 +
876      seleMan_.setSelectionSet(evaluator_.evaluate());
877  
878      int selei;
879      StuntDouble* sd;
880      int idx;
881  
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?
882      for (sd = seleMan_.beginSelected(selei); sd != NULL;
883           sd = seleMan_.nextSelected(selei)) {
884        
# Line 324 | Line 895 | namespace oopse {
895        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
896        
897        int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
898 <      
899 <      //std::cerr << "pos.z() = " << pos.z() << " bin = " << binNo << "\n";
900 <      
898 >
899 >      if (rnemdLogWidth_ == midBin_ + 1)
900 >        if (binNo > midBin_)
901 >          binNo = nBins_ - binNo;
902 >
903        RealType mass = sd->getMass();
904        Vector3d vel = sd->getVel();
332 –      //std::cerr << "mass = " << mass << " vel = " << vel << "\n";
905        RealType value;
906 +      RealType xVal, yVal, zVal;
907  
908        switch(rnemdType_) {
909 <      case rnemdKinetic :
909 >      case rnemdKineticSwap :
910 >      case rnemdKineticScale :
911          
912          value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
913                          vel[2]*vel[2]);
914          
915 <        valueCount[binNo] += 3;
342 <        //std::cerr <<"starting value = " << value << "\n";
915 >        valueCount_[binNo] += 3;
916          if (sd->isDirectional()) {
344 –          //std::cerr << "angMom calculated.\n";
917            Vector3d angMom = sd->getJ();
346 –          //std::cerr << "current angMom: " << angMom << "\n";
918            Mat3x3d I = sd->getI();
919            
920            if (sd->isLinear()) {
# Line 353 | Line 924 | namespace oopse {
924              value += angMom[j] * angMom[j] / I(j, j) +
925                angMom[k] * angMom[k] / I(k, k);
926  
927 <            valueCount[binNo] +=2;
927 >            valueCount_[binNo] +=2;
928  
929            } else {
359 –            //std::cerr << "non-linear molecule.\n";
930              value += angMom[0]*angMom[0]/I(0, 0)
931                + angMom[1]*angMom[1]/I(1, 1)
932                + angMom[2]*angMom[2]/I(2, 2);
933 <            valueCount[binNo] +=3;
364 <
933 >            valueCount_[binNo] +=3;
934            }
935          }
936 <        //std::cerr <<"total value = " << value << "\n";
937 <        //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";
936 >        value = value / PhysicalConstants::energyConvert / PhysicalConstants::kb;
937 >
938          break;
939        case rnemdPx :
940 +      case rnemdPxScale :
941          value = mass * vel[0];
942 <        valueCount[binNo]++;
942 >        valueCount_[binNo]++;
943 >        xVal = mass * vel.x() * vel.x() / PhysicalConstants::energyConvert
944 >          / PhysicalConstants::kb;
945 >        yVal = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert
946 >          / PhysicalConstants::kb;
947 >        zVal = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert
948 >          / PhysicalConstants::kb;
949 >        xTempHist_[binNo] += xVal;
950 >        yTempHist_[binNo] += yVal;
951 >        zTempHist_[binNo] += zVal;
952          break;
953        case rnemdPy :
954 +      case rnemdPyScale :
955          value = mass * vel[1];
956 <        valueCount[binNo]++;
956 >        valueCount_[binNo]++;
957          break;
958        case rnemdPz :
959 +      case rnemdPzScale :
960          value = mass * vel[2];
961 <        valueCount[binNo]++;
961 >        valueCount_[binNo]++;
962          break;
963        case rnemdUnknown :
964        default :
965          break;
966        }
967 <      //std::cerr << "bin = " << binNo << " value = " << value ;
390 <      valueHist[binNo] += value;
391 <      //std::cerr << " hist = " << valueHist[binNo] << " count = " << valueCount[binNo] << "\n";
967 >      valueHist_[binNo] += value;
968      }
969 <    
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";
969 >
970    }
971 +
972 +  void RNEMD::getStarted() {
973 +    Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
974 +    Stats& stat = currentSnap_->statData;
975 +    stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
976 +  }
977 +
978 +  void RNEMD::getStatus() {
979 +
980 +    Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
981 +    Stats& stat = currentSnap_->statData;
982 +    RealType time = currentSnap_->getTime();
983 +
984 +    stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
985 +    //or to be more meaningful, define another item as exchangeSum_ / time
986 +    int j;
987 +
988 + #ifdef IS_MPI
989 +
990 +    // all processors have the same number of bins, and STL vectors pack their
991 +    // arrays, so in theory, this should be safe:
992 +
993 +    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueHist_[0],
994 +                              rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
995 +    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount_[0],
996 +                              rnemdLogWidth_, MPI::INT, MPI::SUM);
997 +    if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale) {
998 +      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xTempHist_[0],
999 +                                rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1000 +      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &yTempHist_[0],
1001 +                                rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1002 +      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &zTempHist_[0],
1003 +                                rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1004 +    }
1005 +    // If we're the root node, should we print out the results
1006 +    int worldRank = MPI::COMM_WORLD.Get_rank();
1007 +    if (worldRank == 0) {
1008 + #endif
1009 +      rnemdLog_ << time;
1010 +      for (j = 0; j < rnemdLogWidth_; j++) {
1011 +        rnemdLog_ << "\t" << valueHist_[j] / (RealType)valueCount_[j];
1012 +      }
1013 +      rnemdLog_ << "\n";
1014 +      if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale ) {
1015 +        xTempLog_ << time;      
1016 +        for (j = 0; j < rnemdLogWidth_; j++) {
1017 +          xTempLog_ << "\t" << xTempHist_[j] / (RealType)valueCount_[j];
1018 +        }
1019 +        xTempLog_ << "\n";
1020 +        yTempLog_ << time;
1021 +        for (j = 0; j < rnemdLogWidth_; j++) {
1022 +          yTempLog_ << "\t" << yTempHist_[j] / (RealType)valueCount_[j];
1023 +        }
1024 +        yTempLog_ << "\n";
1025 +        zTempLog_ << time;
1026 +        for (j = 0; j < rnemdLogWidth_; j++) {
1027 +          zTempLog_ << "\t" << zTempHist_[j] / (RealType)valueCount_[j];
1028 +        }
1029 +        zTempLog_ << "\n";
1030 +      }
1031 + #ifdef IS_MPI
1032 +    }
1033 + #endif
1034 +    for (j = 0; j < rnemdLogWidth_; j++) {
1035 +      valueCount_[j] = 0;
1036 +      valueHist_[j] = 0.0;
1037 +    }
1038 +    if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale)
1039 +      for (j = 0; j < rnemdLogWidth_; j++) {
1040 +        xTempHist_[j] = 0.0;
1041 +        yTempHist_[j] = 0.0;
1042 +        zTempHist_[j] = 0.0;
1043 +      }
1044 +  }
1045   }

Diff Legend

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