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 1350 by gezelter, Thu May 21 18:56:45 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 55 | Line 57
57  
58   #define HONKING_LARGE_VALUE 1.0e10
59  
60 < namespace oopse {
60 > namespace OpenMD {
61    
62    RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
63 <    
63 >
64 >    failTrialCount_ = 0;
65 >    failRootCount_ = 0;
66 >
67      int seedValue;
68      Globals * simParams = info->getSimParams();
69  
70 <    stringToEnumMap_["Kinetic"] = rnemdKinetic;
70 >    stringToEnumMap_["KineticSwap"] = rnemdKineticSwap;
71 >    stringToEnumMap_["KineticScale"] = rnemdKineticScale;
72 >    stringToEnumMap_["PxScale"] = rnemdPxScale;
73 >    stringToEnumMap_["PyScale"] = rnemdPyScale;
74 >    stringToEnumMap_["PzScale"] = rnemdPzScale;
75      stringToEnumMap_["Px"] = rnemdPx;
76      stringToEnumMap_["Py"] = rnemdPy;
77      stringToEnumMap_["Pz"] = rnemdPz;
# Line 72 | Line 81 | namespace oopse {
81      evaluator_.loadScriptString(rnemdObjectSelection_);
82      seleMan_.setSelectionSet(evaluator_.evaluate());
83  
75 –
84      // do some sanity checking
85  
86      int selectionCount = seleMan_.getSelectionCount();
# Line 94 | Line 102 | namespace oopse {
102  
103      }
104      
105 <    const std::string st = simParams->getRNEMD_swapType();
105 >    const std::string st = simParams->getRNEMD_exchangeType();
106  
107      std::map<std::string, RNEMDTypeEnum>::iterator i;
108      i = stringToEnumMap_.find(st);
109 <    rnemdType_  = (i == stringToEnumMap_.end()) ? RNEMD::rnemdUnknown : i->second;
109 >    rnemdType_ = (i == stringToEnumMap_.end()) ? RNEMD::rnemdUnknown : i->second;
110 >    if (rnemdType_ == rnemdUnknown) {
111 >      std::cerr << "WARNING! RNEMD Type Unknown!\n";
112 >    }
113  
114 <    set_RNEMD_swapTime(simParams->getRNEMD_swapTime());
114 > #ifdef IS_MPI
115 >    if (worldRank == 0) {
116 > #endif
117 >
118 >      std::string rnemdFileName;
119 >      std::string xTempFileName;
120 >      std::string yTempFileName;
121 >      std::string zTempFileName;
122 >      switch(rnemdType_) {
123 >      case rnemdKineticSwap :
124 >      case rnemdKineticScale :
125 >        rnemdFileName = "temperature.log";
126 >        break;
127 >      case rnemdPx :
128 >      case rnemdPxScale :
129 >      case rnemdPy :
130 >      case rnemdPyScale :
131 >        rnemdFileName = "momemtum.log";
132 >        xTempFileName = "temperatureX.log";
133 >        yTempFileName = "temperatureY.log";
134 >        zTempFileName = "temperatureZ.log";
135 >        xTempLog_.open(xTempFileName.c_str());
136 >        yTempLog_.open(yTempFileName.c_str());
137 >        zTempLog_.open(zTempFileName.c_str());
138 >        break;
139 >      case rnemdPz :
140 >      case rnemdPzScale :
141 >      case rnemdUnknown :
142 >      default :
143 >        rnemdFileName = "rnemd.log";
144 >        break;
145 >      }
146 >      rnemdLog_.open(rnemdFileName.c_str());
147 >
148 > #ifdef IS_MPI
149 >    }
150 > #endif
151 >
152 >    set_RNEMD_exchange_time(simParams->getRNEMD_exchangeTime());
153      set_RNEMD_nBins(simParams->getRNEMD_nBins());
154 <    exchangeSum_ = 0.0;
154 >    midBin_ = nBins_ / 2;
155 >    if (simParams->haveRNEMD_logWidth()) {
156 >      rnemdLogWidth_ = simParams->getRNEMD_logWidth();
157 >      if (rnemdLogWidth_ != nBins_ || rnemdLogWidth_ != midBin_ + 1) {
158 >        std::cerr << "WARNING! RNEMD_logWidth has abnormal value!\n";
159 >        std::cerr << "Automaically set back to default.\n";
160 >        rnemdLogWidth_ = nBins_;
161 >      }
162 >    } else {
163 >      rnemdLogWidth_ = nBins_;
164 >    }
165 >    valueHist_.resize(rnemdLogWidth_, 0.0);
166 >    valueCount_.resize(rnemdLogWidth_, 0);
167 >    xTempHist_.resize(rnemdLogWidth_, 0.0);
168 >    yTempHist_.resize(rnemdLogWidth_, 0.0);
169 >    zTempHist_.resize(rnemdLogWidth_, 0.0);
170  
171 +    set_RNEMD_exchange_total(0.0);
172 +    if (simParams->haveRNEMD_targetFlux()) {
173 +      set_RNEMD_target_flux(simParams->getRNEMD_targetFlux());
174 +    } else {
175 +      set_RNEMD_target_flux(0.0);
176 +    }
177 +
178   #ifndef IS_MPI
179      if (simParams->haveSeed()) {
180        seedValue = simParams->getSeed();
# Line 123 | Line 194 | namespace oopse {
194    
195    RNEMD::~RNEMD() {
196      delete randNumGen_;
197 +
198 +    std::cerr << "total fail trials: " << failTrialCount_ << "\n";
199 + #ifdef IS_MPI
200 +    if (worldRank == 0) {
201 + #endif
202 +      rnemdLog_.close();
203 +      if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPyScale)
204 +        std::cerr<< "total root-checking warnings: " << failRootCount_ << "\n";
205 +      if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPy || rnemdType_ == rnemdPyScale) {
206 +        xTempLog_.close();
207 +        yTempLog_.close();
208 +        zTempLog_.close();
209 +      }
210 + #ifdef IS_MPI
211 +    }
212 + #endif
213    }
214  
215    void RNEMD::doSwap() {
129 –    int midBin = nBins_ / 2;
216  
217      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
218      Mat3x3d hmat = currentSnap_->getHmat();
# Line 164 | Line 250 | namespace oopse {
250  
251  
252        // if we're in bin 0 or the middleBin
253 <      if (binNo == 0 || binNo == midBin) {
253 >      if (binNo == 0 || binNo == midBin_) {
254          
255          RealType mass = sd->getMass();
256          Vector3d vel = sd->getVel();
257          RealType value;
258  
259          switch(rnemdType_) {
260 <        case rnemdKinetic :
260 >        case rnemdKineticSwap :
261            
262            value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
263                            vel[2]*vel[2]);
# Line 191 | Line 277 | namespace oopse {
277                  + angMom[2]*angMom[2]/I(2, 2);
278              }
279            }
280 <          value = value * 0.5 / OOPSEConstant::energyConvert;
280 >          //make exchangeSum_ comparable between swap & scale
281 >          //temporarily without using energyConvert
282 >          //value = value * 0.5 / PhysicalConstants::energyConvert;
283 >          value *= 0.5;
284            break;
285          case rnemdPx :
286            value = mass * vel[0];
# Line 202 | Line 291 | namespace oopse {
291          case rnemdPz :
292            value = mass * vel[2];
293            break;
205 –        case rnemdUnknown :
294          default :
295            break;
296          }
# Line 218 | Line 306 | namespace oopse {
306                min_sd = sd;
307              }
308            }
309 <        } else {
309 >        } else { //midBin_
310            if (!max_found) {
311              max_val = value;
312              max_sd = sd;
# Line 298 | Line 386 | namespace oopse {
386            RealType temp_vel;
387            
388            switch(rnemdType_) {
389 <          case rnemdKinetic :
389 >          case rnemdKineticSwap :
390              min_sd->setVel(max_vel);
391              max_sd->setVel(min_vel);
392              if (min_sd->isDirectional() && max_sd->isDirectional()) {
# Line 329 | Line 417 | namespace oopse {
417              min_sd->setVel(min_vel);
418              max_sd->setVel(max_vel);
419              break;
332 –          case rnemdUnknown :
420            default :
421              break;
422            }
# Line 349 | Line 436 | namespace oopse {
436                                     min_vals.rank, 0, status);
437            
438            switch(rnemdType_) {
439 <          case rnemdKinetic :
439 >          case rnemdKineticSwap :
440              max_sd->setVel(min_vel);
441              
442              if (max_sd->isDirectional()) {
# Line 378 | Line 465 | namespace oopse {
465              max_vel.z() = min_vel.z();
466              max_sd->setVel(max_vel);
467              break;
381 –          case rnemdUnknown :
468            default :
469              break;
470            }
# Line 396 | Line 482 | namespace oopse {
482                                     max_vals.rank, 0, status);
483            
484            switch(rnemdType_) {
485 <          case rnemdKinetic :
485 >          case rnemdKineticSwap :
486              min_sd->setVel(max_vel);
487              
488              if (min_sd->isDirectional()) {
# Line 425 | Line 511 | namespace oopse {
511              min_vel.z() = max_vel.z();
512              min_sd->setVel(min_vel);
513              break;
428 –          case rnemdUnknown :
514            default :
515              break;
516            }
# Line 433 | Line 518 | namespace oopse {
518   #endif
519          exchangeSum_ += max_val - min_val;
520        } else {
521 <        std::cerr << "exchange NOT performed.\nmin_val > max_val.\n";
521 >        std::cerr << "exchange NOT performed!\nmin_val > max_val.\n";
522 >        failTrialCount_++;
523        }
524      } else {
525 <      std::cerr << "exchange NOT performed.\none of the two slabs empty.\n";
525 >      std::cerr << "exchange NOT performed!\n";
526 >      std::cerr << "at least one of the two slabs empty.\n";
527 >      failTrialCount_++;
528      }
529      
530    }
531    
532 <  void RNEMD::getStatus() {
532 >  void RNEMD::doScale() {
533  
534      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
535      Mat3x3d hmat = currentSnap_->getHmat();
448 –    Stats& stat = currentSnap_->statData;
449 –    RealType time = currentSnap_->getTime();
536  
451 –    stat[Stats::RNEMD_SWAP_TOTAL] = exchangeSum_;
452 –
537      seleMan_.setSelectionSet(evaluator_.evaluate());
538  
539      int selei;
540      StuntDouble* sd;
541      int idx;
542  
543 <    std::vector<RealType> valueHist(nBins_, 0.0); // keeps track of what's
544 <                                                  // being averaged
545 <    std::vector<int> valueCount(nBins_, 0);       // keeps track of the
546 <                                                  // number of degrees of
547 <                                                  // freedom being averaged
543 >    std::vector<StuntDouble*> hotBin, coldBin;
544 >
545 >    RealType Phx = 0.0;
546 >    RealType Phy = 0.0;
547 >    RealType Phz = 0.0;
548 >    RealType Khx = 0.0;
549 >    RealType Khy = 0.0;
550 >    RealType Khz = 0.0;
551 >    RealType Pcx = 0.0;
552 >    RealType Pcy = 0.0;
553 >    RealType Pcz = 0.0;
554 >    RealType Kcx = 0.0;
555 >    RealType Kcy = 0.0;
556 >    RealType Kcz = 0.0;
557 >
558 >    for (sd = seleMan_.beginSelected(selei); sd != NULL;
559 >         sd = seleMan_.nextSelected(selei)) {
560 >
561 >      idx = sd->getLocalIndex();
562 >
563 >      Vector3d pos = sd->getPos();
564 >
565 >      // wrap the stuntdouble's position back into the box:
566 >
567 >      if (usePeriodicBoundaryConditions_)
568 >        currentSnap_->wrapVector(pos);
569 >
570 >      // which bin is this stuntdouble in?
571 >      // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
572 >
573 >      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
574 >
575 >      // if we're in bin 0 or the middleBin
576 >      if (binNo == 0 || binNo == midBin_) {
577 >        
578 >        RealType mass = sd->getMass();
579 >        Vector3d vel = sd->getVel();
580 >      
581 >        if (binNo == 0) {
582 >          hotBin.push_back(sd);
583 >          Phx += mass * vel.x();
584 >          Phy += mass * vel.y();
585 >          Phz += mass * vel.z();
586 >          Khx += mass * vel.x() * vel.x();
587 >          Khy += mass * vel.y() * vel.y();
588 >          Khz += mass * vel.z() * vel.z();
589 >        } else { //midBin_
590 >          coldBin.push_back(sd);
591 >          Pcx += mass * vel.x();
592 >          Pcy += mass * vel.y();
593 >          Pcz += mass * vel.z();
594 >          Kcx += mass * vel.x() * vel.x();
595 >          Kcy += mass * vel.y() * vel.y();
596 >          Kcz += mass * vel.z() * vel.z();
597 >        }
598 >      }
599 >    }
600 >
601 >    Khx *= 0.5;
602 >    Khy *= 0.5;
603 >    Khz *= 0.5;
604 >    Kcx *= 0.5;
605 >    Kcy *= 0.5;
606 >    Kcz *= 0.5;
607 >
608 > #ifdef IS_MPI
609 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phx, 1, MPI::REALTYPE, MPI::SUM);
610 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phy, 1, MPI::REALTYPE, MPI::SUM);
611 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phz, 1, MPI::REALTYPE, MPI::SUM);
612 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcx, 1, MPI::REALTYPE, MPI::SUM);
613 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcy, 1, MPI::REALTYPE, MPI::SUM);
614 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcz, 1, MPI::REALTYPE, MPI::SUM);
615 >
616 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khx, 1, MPI::REALTYPE, MPI::SUM);
617 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khy, 1, MPI::REALTYPE, MPI::SUM);
618 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khz, 1, MPI::REALTYPE, MPI::SUM);
619 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcx, 1, MPI::REALTYPE, MPI::SUM);
620 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcy, 1, MPI::REALTYPE, MPI::SUM);
621 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcz, 1, MPI::REALTYPE, MPI::SUM);
622 > #endif
623 >
624 >    //use coldBin coeff's
625 >    RealType px = Pcx / Phx;
626 >    RealType py = Pcy / Phy;
627 >    RealType pz = Pcz / Phz;
628 >
629 >    RealType a000, a110, c0, a001, a111, b01, b11, c1, c;
630 >    switch(rnemdType_) {
631 >    case rnemdKineticScale :
632 >    /*used hotBin coeff's & only scale x & y dimensions
633 >      RealType px = Phx / Pcx;
634 >      RealType py = Phy / Pcy;
635 >      a110 = Khy;
636 >      c0 = - Khx - Khy - targetFlux_;
637 >      a000 = Khx;
638 >      a111 = Kcy * py * py
639 >      b11 = -2.0 * Kcy * py * (1.0 + py);
640 >      c1 = Kcy * py * (2.0 + py) + Kcx * px * ( 2.0 + px) + targetFlux_;
641 >      b01 = -2.0 * Kcx * px * (1.0 + px);
642 >      a001 = Kcx * px * px;
643 >    */
644 >
645 >      //scale all three dimensions, let x = y
646 >      a000 = Kcx + Kcy;
647 >      a110 = Kcz;
648 >      c0 = targetFlux_ - Kcx - Kcy - Kcz;
649 >      a001 = Khx * px * px + Khy * py * py;
650 >      a111 = Khz * pz * pz;
651 >      b01 = -2.0 * (Khx * px * (1.0 + px) + Khy * py * (1.0 + py));
652 >      b11 = -2.0 * Khz * pz * (1.0 + pz);
653 >      c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py)
654 >         + Khz * pz * (2.0 + pz) - targetFlux_;
655 >      break;
656 >    case rnemdPxScale :
657 >      c = 1 - targetFlux_ / Pcx;
658 >      a000 = Kcy;
659 >      a110 = Kcz;
660 >      c0 = Kcx * c * c - Kcx - Kcy - Kcz;
661 >      a001 = py * py * Khy;
662 >      a111 = pz * pz * Khz;
663 >      b01 = -2.0 * Khy * py * (1.0 + py);
664 >      b11 = -2.0 * Khz * pz * (1.0 + pz);
665 >      c1 = Khy * py * (2.0 + py) + Khz * pz * (2.0 + pz)
666 >         + Khx * (fastpow(c * px - px - 1.0, 2) - 1.0);
667 >      break;
668 >    case rnemdPyScale :
669 >      c = 1 - targetFlux_ / Pcy;
670 >      a000 = Kcx;
671 >      a110 = Kcz;
672 >      c0 = Kcy * c * c - Kcx - Kcy - Kcz;
673 >      a001 = px * px * Khx;
674 >      a111 = pz * pz * Khz;
675 >      b01 = -2.0 * Khx * px * (1.0 + px);
676 >      b11 = -2.0 * Khz * pz * (1.0 + pz);
677 >      c1 = Khx * px * (2.0 + px) + Khz * pz * (2.0 + pz)
678 >         + Khy * (fastpow(c * py - py - 1.0, 2) - 1.0);
679 >      break;
680 >    case rnemdPzScale ://we don't really do this, do we?
681 >    default :
682 >      break;
683 >    }
684  
685 +    RealType v1 = a000 * a111 - a001 * a110;
686 +    RealType v2 = a000 * b01;
687 +    RealType v3 = a000 * b11;
688 +    RealType v4 = a000 * c1 - a001 * c0;
689 +    RealType v8 = a110 * b01;
690 +    RealType v10 = - b01 * c0;
691 +
692 +    RealType u0 = v2 * v10 - v4 * v4;
693 +    RealType u1 = -2.0 * v3 * v4;
694 +    RealType u2 = -v2 * v8 - v3 * v3 - 2.0 * v1 * v4;
695 +    RealType u3 = -2.0 * v1 * v3;
696 +    RealType u4 = - v1 * v1;
697 +    //rescale coefficients
698 +    RealType maxAbs = fabs(u0);
699 +    if (maxAbs < fabs(u1)) maxAbs = fabs(u1);
700 +    if (maxAbs < fabs(u2)) maxAbs = fabs(u2);
701 +    if (maxAbs < fabs(u3)) maxAbs = fabs(u3);
702 +    if (maxAbs < fabs(u4)) maxAbs = fabs(u4);
703 +    u0 /= maxAbs;
704 +    u1 /= maxAbs;
705 +    u2 /= maxAbs;
706 +    u3 /= maxAbs;
707 +    u4 /= maxAbs;
708 +    //max_element(start, end) is also available.
709 +    Polynomial<RealType> poly; //same as DoublePolynomial poly;
710 +    poly.setCoefficient(4, u4);
711 +    poly.setCoefficient(3, u3);
712 +    poly.setCoefficient(2, u2);
713 +    poly.setCoefficient(1, u1);
714 +    poly.setCoefficient(0, u0);
715 +    std::vector<RealType> realRoots = poly.FindRealRoots();
716 +
717 +    std::vector<RealType>::iterator ri;
718 +    RealType r1, r2, alpha0;
719 +    std::vector<std::pair<RealType,RealType> > rps;
720 +    for (ri = realRoots.begin(); ri !=realRoots.end(); ri++) {
721 +      r2 = *ri;
722 +      //check if FindRealRoots() give the right answer
723 +      if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) {
724 +        std::cerr << "WARNING! eq solvers might have mistakes!\n";
725 +        failRootCount_++;
726 +      }
727 +      //might not be useful w/o rescaling coefficients
728 +      alpha0 = -c0 - a110 * r2 * r2;
729 +      if (alpha0 >= 0.0) {
730 +        r1 = sqrt(alpha0 / a000);
731 +        if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) < 1e-6)
732 +          { rps.push_back(std::make_pair(r1, r2)); }
733 +        if (r1 > 1e-6) { //r1 non-negative
734 +          r1 = -r1;
735 +          if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) <1e-6)
736 +            { rps.push_back(std::make_pair(r1, r2)); }
737 +        }
738 +      }
739 +    }
740 +    // Consider combininig together the solving pair part w/ the searching
741 +    // best solution part so that we don't need the pairs vector
742 +    if (!rps.empty()) {
743 +      RealType smallestDiff = HONKING_LARGE_VALUE;
744 +      RealType diff;
745 +      std::pair<RealType,RealType> bestPair = std::make_pair(1.0, 1.0);
746 +      std::vector<std::pair<RealType,RealType> >::iterator rpi;
747 +      for (rpi = rps.begin(); rpi != rps.end(); rpi++) {
748 +        r1 = (*rpi).first;
749 +        r2 = (*rpi).second;
750 +        switch(rnemdType_) {
751 +        case rnemdKineticScale :
752 +          diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
753 +            + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2)
754 +            + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcy, 2);
755 +          break;
756 +        case rnemdPxScale :
757 +          diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
758 +            + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcy, 2);
759 +          break;
760 +        case rnemdPyScale :
761 +          diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
762 +            + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2);
763 +          break;
764 +        case rnemdPzScale :
765 +        default :
766 +          break;
767 +        }
768 +        if (diff < smallestDiff) {
769 +          smallestDiff = diff;
770 +          bestPair = *rpi;
771 +        }
772 +      }
773 + #ifdef IS_MPI
774 +      if (worldRank == 0) {
775 + #endif
776 +        std::cerr << "we choose r1 = " << bestPair.first
777 +                  << " and r2 = " << bestPair.second << "\n";
778 + #ifdef IS_MPI
779 +      }
780 + #endif
781 +
782 +      RealType x, y, z;
783 +        switch(rnemdType_) {
784 +        case rnemdKineticScale :
785 +          x = bestPair.first;
786 +          y = bestPair.first;
787 +          z = bestPair.second;
788 +          break;
789 +        case rnemdPxScale :
790 +          x = c;
791 +          y = bestPair.first;
792 +          z = bestPair.second;
793 +          break;
794 +        case rnemdPyScale :
795 +          x = bestPair.first;
796 +          y = c;
797 +          z = bestPair.second;
798 +          break;
799 +        case rnemdPzScale :
800 +        default :
801 +          break;
802 +        }
803 +      std::vector<StuntDouble*>::iterator sdi;
804 +      Vector3d vel;
805 +      for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
806 +        vel = (*sdi)->getVel();
807 +        vel.x() *= x;
808 +        vel.y() *= y;
809 +        vel.z() *= z;
810 +        (*sdi)->setVel(vel);
811 +      }
812 +      //convert to hotBin coefficient
813 +      x = 1.0 + px * (1.0 - x);
814 +      y = 1.0 + py * (1.0 - y);
815 +      z = 1.0 + pz * (1.0 - z);
816 +      for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
817 +        vel = (*sdi)->getVel();
818 +        vel.x() *= x;
819 +        vel.y() *= y;
820 +        vel.z() *= z;
821 +        (*sdi)->setVel(vel);
822 +      }
823 +      exchangeSum_ += targetFlux_;
824 +      //we may want to check whether the exchange has been successful
825 +    } else {
826 +      std::cerr << "exchange NOT performed!\n";
827 +      failTrialCount_++;
828 +    }
829 +
830 +  }
831 +
832 +  void RNEMD::doRNEMD() {
833 +
834 +    switch(rnemdType_) {
835 +    case rnemdKineticScale :
836 +    case rnemdPxScale :
837 +    case rnemdPyScale :
838 +    case rnemdPzScale :
839 +      doScale();
840 +      break;
841 +    case rnemdKineticSwap :
842 +    case rnemdPx :
843 +    case rnemdPy :
844 +    case rnemdPz :
845 +      doSwap();
846 +      break;
847 +    case rnemdUnknown :
848 +    default :
849 +      break;
850 +    }
851 +  }
852 +
853 +  void RNEMD::collectData() {
854 +
855 +    Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
856 +    Mat3x3d hmat = currentSnap_->getHmat();
857 +
858 +    seleMan_.setSelectionSet(evaluator_.evaluate());
859 +
860 +    int selei;
861 +    StuntDouble* sd;
862 +    int idx;
863 +
864      for (sd = seleMan_.beginSelected(selei); sd != NULL;
865           sd = seleMan_.nextSelected(selei)) {
866        
# Line 477 | Line 876 | namespace oopse {
876        // which bin is this stuntdouble in?
877        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
878        
879 <      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;    
880 <      
879 >      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
880 >
881 >      if (rnemdLogWidth_ == midBin_ + 1)
882 >        if (binNo > midBin_)
883 >          binNo = nBins_ - binNo;
884 >
885        RealType mass = sd->getMass();
886        Vector3d vel = sd->getVel();
887        RealType value;
888 +      RealType xVal, yVal, zVal;
889  
890        switch(rnemdType_) {
891 <      case rnemdKinetic :
891 >      case rnemdKineticSwap :
892 >      case rnemdKineticScale :
893          
894          value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
895                          vel[2]*vel[2]);
896          
897 <        valueCount[binNo] += 3;
897 >        valueCount_[binNo] += 3;
898          if (sd->isDirectional()) {
899            Vector3d angMom = sd->getJ();
900            Mat3x3d I = sd->getI();
# Line 501 | Line 906 | namespace oopse {
906              value += angMom[j] * angMom[j] / I(j, j) +
907                angMom[k] * angMom[k] / I(k, k);
908  
909 <            valueCount[binNo] +=2;
909 >            valueCount_[binNo] +=2;
910  
911            } else {
912              value += angMom[0]*angMom[0]/I(0, 0)
913                + angMom[1]*angMom[1]/I(1, 1)
914                + angMom[2]*angMom[2]/I(2, 2);
915 <            valueCount[binNo] +=3;
915 >            valueCount_[binNo] +=3;
916            }
917          }
918 <        value = value / OOPSEConstant::energyConvert / OOPSEConstant::kb;
918 >        value = value / PhysicalConstants::energyConvert / PhysicalConstants::kb;
919  
920          break;
921        case rnemdPx :
922 +      case rnemdPxScale :
923          value = mass * vel[0];
924 <        valueCount[binNo]++;
924 >        valueCount_[binNo]++;
925 >        xVal = mass * vel.x() * vel.x() / PhysicalConstants::energyConvert
926 >          / PhysicalConstants::kb;
927 >        yVal = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert
928 >          / PhysicalConstants::kb;
929 >        zVal = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert
930 >          / PhysicalConstants::kb;
931 >        xTempHist_[binNo] += xVal;
932 >        yTempHist_[binNo] += yVal;
933 >        zTempHist_[binNo] += zVal;
934          break;
935        case rnemdPy :
936 +      case rnemdPyScale :
937          value = mass * vel[1];
938 <        valueCount[binNo]++;
938 >        valueCount_[binNo]++;
939          break;
940        case rnemdPz :
941 +      case rnemdPzScale :
942          value = mass * vel[2];
943 <        valueCount[binNo]++;
943 >        valueCount_[binNo]++;
944          break;
945        case rnemdUnknown :
946        default :
947          break;
948        }
949 <      valueHist[binNo] += value;
949 >      valueHist_[binNo] += value;
950      }
951  
952 +  }
953 +
954 +  void RNEMD::getStarted() {
955 +    Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
956 +    Stats& stat = currentSnap_->statData;
957 +    stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
958 +  }
959 +
960 +  void RNEMD::getStatus() {
961 +
962 +    Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
963 +    Stats& stat = currentSnap_->statData;
964 +    RealType time = currentSnap_->getTime();
965 +
966 +    stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
967 +    //or to be more meaningful, define another item as exchangeSum_ / time
968 +
969 +
970   #ifdef IS_MPI
971  
972      // all processors have the same number of bins, and STL vectors pack their
973      // arrays, so in theory, this should be safe:
974  
975 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueHist[0],
976 <                              nBins_, MPI::REALTYPE, MPI::SUM);
977 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount[0],
978 <                              nBins_, MPI::INT, MPI::SUM);
975 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueHist_[0],
976 >                              rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
977 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount_[0],
978 >                              rnemdLogWidth_, MPI::INT, MPI::SUM);
979  
980      // If we're the root node, should we print out the results
981      int worldRank = MPI::COMM_WORLD.Get_rank();
982      if (worldRank == 0) {
983   #endif
984 <      
985 <      std::cout << time;
986 <      for (int j = 0; j < nBins_; j++)
987 <        std::cout << "\t" << valueHist[j] / (RealType)valueCount[j];
988 <      std::cout << "\n";
989 <      
984 >      int j;
985 >      rnemdLog_ << time;
986 >      for (j = 0; j < rnemdLogWidth_; j++) {
987 >        rnemdLog_ << "\t" << valueHist_[j] / (RealType)valueCount_[j];
988 >        valueHist_[j] = 0.0;
989 >      }
990 >      rnemdLog_ << "\n";
991 >      if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale ) {
992 >        xTempLog_ << time;      
993 >        for (j = 0; j < rnemdLogWidth_; j++) {
994 >          xTempLog_ << "\t" << xTempHist_[j] / (RealType)valueCount_[j];
995 >          xTempHist_[j] = 0.0;
996 >        }
997 >        xTempLog_ << "\n";
998 >        yTempLog_ << time;
999 >        for (j = 0; j < rnemdLogWidth_; j++) {
1000 >          yTempLog_ << "\t" << yTempHist_[j] / (RealType)valueCount_[j];
1001 >          yTempHist_[j] = 0.0;
1002 >        }
1003 >        yTempLog_ << "\n";
1004 >        zTempLog_ << time;
1005 >        for (j = 0; j < rnemdLogWidth_; j++) {
1006 >          zTempLog_ << "\t" << zTempHist_[j] / (RealType)valueCount_[j];
1007 >          zTempHist_[j] = 0.0;
1008 >        }
1009 >        zTempLog_ << "\n";
1010 >      }
1011 >      for (j = 0; j < rnemdLogWidth_; j++) valueCount_[j] = 0;
1012   #ifdef IS_MPI
1013 <    }
1013 >    }    
1014   #endif
1015 +
1016 +      
1017    }
1018 +
1019   }

Diff Legend

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