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

Comparing branches/development/src/integrators/RNEMD.cpp (file contents):
Revision 1628 by gezelter, Tue Sep 13 22:05:04 2011 UTC vs.
Revision 1629 by gezelter, Wed Sep 14 21:15:17 2011 UTC

# Line 58 | Line 58
58  
59   #define HONKING_LARGE_VALUE 1.0e10
60  
61 + using namespace std;
62   namespace OpenMD {
63    
64 <  RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
64 >  RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info),
65 >                                usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
66  
67      failTrialCount_ = 0;
68      failRootCount_ = 0;
# Line 89 | Line 91 | namespace OpenMD {
91  
92      if (selectionCount > nIntegrable) {
93        sprintf(painCave.errMsg,
94 <              "RNEMD warning: The current RNEMD_objectSelection,\n"
94 >              "RNEMD: The current RNEMD_objectSelection,\n"
95                "\t\t%s\n"
96                "\thas resulted in %d selected objects.  However,\n"
97                "\tthe total number of integrable objects in the system\n"
# Line 99 | Line 101 | namespace OpenMD {
101                rnemdObjectSelection_.c_str(),
102                selectionCount, nIntegrable);
103        painCave.isFatal = 0;
104 +      painCave.severity = OPENMD_WARNING;
105        simError();
103 –
106      }
107      
108 <    const std::string st = simParams->getRNEMD_exchangeType();
108 >    const string st = simParams->getRNEMD_exchangeType();
109  
110 <    std::map<std::string, RNEMDTypeEnum>::iterator i;
110 >    map<string, RNEMDTypeEnum>::iterator i;
111      i = stringToEnumMap_.find(st);
112      rnemdType_ = (i == stringToEnumMap_.end()) ? RNEMD::rnemdUnknown : i->second;
113      if (rnemdType_ == rnemdUnknown) {
114 <      std::cerr << "WARNING! RNEMD Type Unknown!\n";
114 >      sprintf(painCave.errMsg,
115 >              "RNEMD: The current RNEMD_exchangeType,\n"
116 >              "\t\t%s\n"
117 >              "\tis not one of the recognized exchange types.\n",
118 >              st.c_str());
119 >      painCave.isFatal = 1;
120 >      painCave.severity = OPENMD_ERROR;
121 >      simError();
122      }
123 +    
124 +    output3DTemp_ = false;
125 +    if (simParams->haveRNEMD_outputDimensionalTemperature()) {
126 +      output3DTemp_ = simParams->getRNEMD_outputDimensionalTemperature();
127 +    }
128  
129   #ifdef IS_MPI
130      if (worldRank == 0) {
131   #endif
132  
133 <      std::string rnemdFileName;
120 <      std::string xTempFileName;
121 <      std::string yTempFileName;
122 <      std::string zTempFileName;
133 >      string rnemdFileName;
134        switch(rnemdType_) {
135        case rnemdKineticSwap :
136        case rnemdKineticScale :
# Line 130 | Line 141 | namespace OpenMD {
141        case rnemdPy :
142        case rnemdPyScale :
143          rnemdFileName = "momemtum.log";
133 –        xTempFileName = "temperatureX.log";
134 –        yTempFileName = "temperatureY.log";
135 –        zTempFileName = "temperatureZ.log";
136 –        xTempLog_.open(xTempFileName.c_str());
137 –        yTempLog_.open(yTempFileName.c_str());
138 –        zTempLog_.open(zTempFileName.c_str());
144          break;
145        case rnemdPz :
146        case rnemdPzScale :
# Line 146 | Line 151 | namespace OpenMD {
151        }
152        rnemdLog_.open(rnemdFileName.c_str());
153  
154 +      string xTempFileName;
155 +      string yTempFileName;
156 +      string zTempFileName;
157 +      if (output3DTemp_) {
158 +        xTempFileName = "temperatureX.log";
159 +        yTempFileName = "temperatureY.log";
160 +        zTempFileName = "temperatureZ.log";
161 +        xTempLog_.open(xTempFileName.c_str());
162 +        yTempLog_.open(yTempFileName.c_str());
163 +        zTempLog_.open(zTempFileName.c_str());
164 +      }
165 +
166   #ifdef IS_MPI
167      }
168   #endif
# Line 153 | Line 170 | namespace OpenMD {
170      set_RNEMD_exchange_time(simParams->getRNEMD_exchangeTime());
171      set_RNEMD_nBins(simParams->getRNEMD_nBins());
172      midBin_ = nBins_ / 2;
173 +    if (simParams->haveRNEMD_binShift()) {
174 +      if (simParams->getRNEMD_binShift()) {
175 +        zShift_ = 0.5 / (RealType)(nBins_);
176 +      } else {
177 +        zShift_ = 0.0;
178 +      }
179 +    } else {
180 +      zShift_ = 0.0;
181 +    }
182 +    //cerr << "we have zShift_ = " << zShift_ << "\n";
183 +    //shift slabs by half slab width, might be useful in heterogeneous systems
184 +    //set to 0.0 if not using it; can NOT be used in status output yet
185      if (simParams->haveRNEMD_logWidth()) {
186 <      rnemdLogWidth_ = simParams->getRNEMD_logWidth();
187 <      if (rnemdLogWidth_ != nBins_ && rnemdLogWidth_ != midBin_ + 1) {
188 <        std::cerr << "WARNING! RNEMD_logWidth has abnormal value!\n";
189 <        std::cerr << "Automaically set back to default.\n";
186 >      set_RNEMD_logWidth(simParams->getRNEMD_logWidth());
187 >      /*arbitary rnemdLogWidth_ no checking
188 >        if (rnemdLogWidth_ != nBins_ && rnemdLogWidth_ != midBin_ + 1) {
189 >        cerr << "WARNING! RNEMD_logWidth has abnormal value!\n";
190 >        cerr << "Automaically set back to default.\n";
191          rnemdLogWidth_ = nBins_;
192 <      }
192 >        }*/
193      } else {
194 <      rnemdLogWidth_ = nBins_;
194 >      set_RNEMD_logWidth(nBins_);
195      }
196      valueHist_.resize(rnemdLogWidth_, 0.0);
197      valueCount_.resize(rnemdLogWidth_, 0);
198      xTempHist_.resize(rnemdLogWidth_, 0.0);
199      yTempHist_.resize(rnemdLogWidth_, 0.0);
200      zTempHist_.resize(rnemdLogWidth_, 0.0);
201 +    xyzTempCount_.resize(rnemdLogWidth_, 0);
202  
203      set_RNEMD_exchange_total(0.0);
204      if (simParams->haveRNEMD_targetFlux()) {
# Line 199 | Line 230 | namespace OpenMD {
230   #ifdef IS_MPI
231      if (worldRank == 0) {
232   #endif
233 <      std::cerr << "total fail trials: " << failTrialCount_ << "\n";
233 >      
234 >      sprintf(painCave.errMsg,
235 >              "RNEMD: total failed trials: %d\n",
236 >              failTrialCount_);
237 >      painCave.isFatal = 0;
238 >      painCave.severity = OPENMD_INFO;
239 >      simError();
240 >
241        rnemdLog_.close();
242 <      if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPyScale)
243 <        std::cerr<< "total root-checking warnings: " << failRootCount_ << "\n";
244 <      if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPy || rnemdType_ == rnemdPyScale) {
242 >      if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPyScale) {
243 >        sprintf(painCave.errMsg,
244 >                "RNEMD: total root-checking warnings: %d\n",
245 >                failRootCount_);
246 >        painCave.isFatal = 0;
247 >        painCave.severity = OPENMD_INFO;
248 >        simError();
249 >      }
250 >      if (output3DTemp_) {
251          xTempLog_.close();
252          yTempLog_.close();
253          zTempLog_.close();
# Line 247 | Line 291 | namespace OpenMD {
291        // which bin is this stuntdouble in?
292        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
293  
294 <      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
294 >      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_;
295  
296  
297        // if we're in bin 0 or the middleBin
# Line 262 | Line 306 | namespace OpenMD {
306            
307            value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
308                            vel[2]*vel[2]);
309 <          if (sd->isDirectional()) {
309 >          /*
310 >            if (sd->isDirectional()) {
311              Vector3d angMom = sd->getJ();
312              Mat3x3d I = sd->getI();
313              
314              if (sd->isLinear()) {
315 <              int i = sd->linearAxis();
316 <              int j = (i + 1) % 3;
317 <              int k = (i + 2) % 3;
318 <              value += angMom[j] * angMom[j] / I(j, j) +
319 <                angMom[k] * angMom[k] / I(k, k);
315 >            int i = sd->linearAxis();
316 >            int j = (i + 1) % 3;
317 >            int k = (i + 2) % 3;
318 >            value += angMom[j] * angMom[j] / I(j, j) +
319 >            angMom[k] * angMom[k] / I(k, k);
320              } else {                        
321 <              value += angMom[0]*angMom[0]/I(0, 0)
322 <                + angMom[1]*angMom[1]/I(1, 1)
323 <                + angMom[2]*angMom[2]/I(2, 2);
321 >            value += angMom[0]*angMom[0]/I(0, 0)
322 >            + angMom[1]*angMom[1]/I(1, 1)
323 >            + angMom[2]*angMom[2]/I(2, 2);
324              }
325 <          }
325 >            } no exchange of angular momenta
326 >          */
327            //make exchangeSum_ comparable between swap & scale
328            //temporarily without using energyConvert
329            //value = value * 0.5 / PhysicalConstants::energyConvert;
# Line 332 | Line 378 | namespace OpenMD {
378      bool my_max_found = max_found;
379  
380      // Even if we didn't find a minimum, did someone else?
381 <    MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found,
336 <                              1, MPI::BOOL, MPI::LAND);
337 <    
381 >    MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, 1, MPI::BOOL, MPI::LOR);
382      // Even if we didn't find a maximum, did someone else?
383 <    MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found,
340 <                              1, MPI::BOOL, MPI::LAND);
341 <    
383 >    MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, 1, MPI::BOOL, MPI::LOR);
384      struct {
385        RealType val;
386        int rank;
# Line 374 | Line 416 | namespace OpenMD {
416   #endif
417  
418      if (max_found && min_found) {
419 <      if (min_val< max_val) {
419 >      if (min_val < max_val) {
420  
421   #ifdef IS_MPI      
422          if (max_vals.rank == worldRank && min_vals.rank == worldRank) {
423            // I have both maximum and minimum, so proceed like a single
424            // processor version:
425   #endif
426 <          // objects to be swapped: velocity & angular velocity
426 >          // objects to be swapped: velocity ONLY
427            Vector3d min_vel = min_sd->getVel();
428            Vector3d max_vel = max_sd->getVel();
429            RealType temp_vel;
# Line 390 | Line 432 | namespace OpenMD {
432            case rnemdKineticSwap :
433              min_sd->setVel(max_vel);
434              max_sd->setVel(min_vel);
435 <            if (min_sd->isDirectional() && max_sd->isDirectional()) {
435 >            /*
436 >              if (min_sd->isDirectional() && max_sd->isDirectional()) {
437                Vector3d min_angMom = min_sd->getJ();
438                Vector3d max_angMom = max_sd->getJ();
439                min_sd->setJ(max_angMom);
440                max_sd->setJ(min_angMom);
441 <            }
441 >              } no angular momentum exchange
442 >            */
443              break;
444            case rnemdPx :
445              temp_vel = min_vel.x();
# Line 439 | Line 483 | namespace OpenMD {
483            switch(rnemdType_) {
484            case rnemdKineticSwap :
485              max_sd->setVel(min_vel);
486 <            
486 >            //no angular momentum exchange for now
487 >            /*
488              if (max_sd->isDirectional()) {
489                Vector3d min_angMom;
490                Vector3d max_angMom = max_sd->getJ();
491 <
491 >              
492                // point-to-point swap of the angular momentum vector
493                MPI::COMM_WORLD.Sendrecv(max_angMom.getArrayPointer(), 3,
494                                         MPI::REALTYPE, min_vals.rank, 1,
495                                         min_angMom.getArrayPointer(), 3,
496                                         MPI::REALTYPE, min_vals.rank, 1,
497                                         status);
498 <
498 >              
499                max_sd->setJ(min_angMom);
500 <            }
500 >             }
501 >             */            
502              break;
503            case rnemdPx :
504              max_vel.x() = min_vel.x();
# Line 485 | Line 531 | namespace OpenMD {
531            switch(rnemdType_) {
532            case rnemdKineticSwap :
533              min_sd->setVel(max_vel);
534 <            
534 >            // no angular momentum exchange for now
535 >            /*
536              if (min_sd->isDirectional()) {
537                Vector3d min_angMom = min_sd->getJ();
538                Vector3d max_angMom;
539 <
539 >              
540                // point-to-point swap of the angular momentum vector
541                MPI::COMM_WORLD.Sendrecv(min_angMom.getArrayPointer(), 3,
542                                         MPI::REALTYPE, max_vals.rank, 1,
543                                         max_angMom.getArrayPointer(), 3,
544                                         MPI::REALTYPE, max_vals.rank, 1,
545                                         status);
546 <
546 >              
547                min_sd->setJ(max_angMom);
548              }
549 +            */
550              break;
551            case rnemdPx :
552              min_vel.x() = max_vel.x();
# Line 518 | Line 566 | namespace OpenMD {
566          }
567   #endif
568          exchangeSum_ += max_val - min_val;
569 <      } else {
570 <        std::cerr << "exchange NOT performed!\nmin_val > max_val.\n";
569 >      } else {        
570 >        sprintf(painCave.errMsg,
571 >                "RNEMD: exchange NOT performed because min_val > max_val\n");
572 >        painCave.isFatal = 0;
573 >        painCave.severity = OPENMD_INFO;
574 >        simError();        
575          failTrialCount_++;
576        }
577      } else {
578 <      std::cerr << "exchange NOT performed!\n";
579 <      std::cerr << "at least one of the two slabs empty.\n";
578 >      sprintf(painCave.errMsg,
579 >              "RNEMD: exchange NOT performed because at least one\n"
580 >              "\tof the two slabs is empty\n");
581 >      painCave.isFatal = 0;
582 >      painCave.severity = OPENMD_INFO;
583 >      simError();        
584        failTrialCount_++;
585      }
586      
# Line 541 | Line 597 | namespace OpenMD {
597      StuntDouble* sd;
598      int idx;
599  
600 <    std::vector<StuntDouble*> hotBin, coldBin;
600 >    vector<StuntDouble*> hotBin, coldBin;
601  
602      RealType Phx = 0.0;
603      RealType Phy = 0.0;
# Line 571 | Line 627 | namespace OpenMD {
627        // which bin is this stuntdouble in?
628        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
629  
630 <      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
630 >      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_;
631  
632        // if we're in bin 0 or the middleBin
633        if (binNo == 0 || binNo == midBin_) {
# Line 630 | Line 686 | namespace OpenMD {
686      RealType a000, a110, c0, a001, a111, b01, b11, c1, c;
687      switch(rnemdType_) {
688      case rnemdKineticScale :
689 <    /*used hotBin coeff's & only scale x & y dimensions
689 >      // used hotBin coeff's & only scale x & y dimensions
690 >      /*
691        RealType px = Phx / Pcx;
692        RealType py = Phy / Pcy;
693        a110 = Khy;
694        c0 = - Khx - Khy - targetFlux_;
695        a000 = Khx;
696 <      a111 = Kcy * py * py
696 >      a111 = Kcy * py * py;
697        b11 = -2.0 * Kcy * py * (1.0 + py);
698        c1 = Kcy * py * (2.0 + py) + Kcx * px * ( 2.0 + px) + targetFlux_;
699        b01 = -2.0 * Kcx * px * (1.0 + px);
700        a001 = Kcx * px * px;
701 <    */
645 <
701 >      */
702        //scale all three dimensions, let c_x = c_y
703        a000 = Kcx + Kcy;
704        a110 = Kcz;
# Line 652 | Line 708 | namespace OpenMD {
708        b01 = -2.0 * (Khx * px * (1.0 + px) + Khy * py * (1.0 + py));
709        b11 = -2.0 * Khz * pz * (1.0 + pz);
710        c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py)
711 <         + Khz * pz * (2.0 + pz) - targetFlux_;
711 >        + Khz * pz * (2.0 + pz) - targetFlux_;
712        break;
713      case rnemdPxScale :
714        c = 1 - targetFlux_ / Pcx;
# Line 664 | Line 720 | namespace OpenMD {
720        b01 = -2.0 * Khy * py * (1.0 + py);
721        b11 = -2.0 * Khz * pz * (1.0 + pz);
722        c1 = Khy * py * (2.0 + py) + Khz * pz * (2.0 + pz)
723 <         + Khx * (fastpow(c * px - px - 1.0, 2) - 1.0);
723 >        + Khx * (fastpow(c * px - px - 1.0, 2) - 1.0);
724        break;
725      case rnemdPyScale :
726        c = 1 - targetFlux_ / Pcy;
# Line 676 | Line 732 | namespace OpenMD {
732        b01 = -2.0 * Khx * px * (1.0 + px);
733        b11 = -2.0 * Khz * pz * (1.0 + pz);
734        c1 = Khx * px * (2.0 + px) + Khz * pz * (2.0 + pz)
735 <         + Khy * (fastpow(c * py - py - 1.0, 2) - 1.0);
735 >        + Khy * (fastpow(c * py - py - 1.0, 2) - 1.0);
736        break;
737      case rnemdPzScale ://we don't really do this, do we?
738        c = 1 - targetFlux_ / Pcz;
# Line 724 | Line 780 | namespace OpenMD {
780      poly.setCoefficient(2, u2);
781      poly.setCoefficient(1, u1);
782      poly.setCoefficient(0, u0);
783 <    std::vector<RealType> realRoots = poly.FindRealRoots();
783 >    vector<RealType> realRoots = poly.FindRealRoots();
784  
785 <    std::vector<RealType>::iterator ri;
785 >    vector<RealType>::iterator ri;
786      RealType r1, r2, alpha0;
787 <    std::vector<std::pair<RealType,RealType> > rps;
787 >    vector<pair<RealType,RealType> > rps;
788      for (ri = realRoots.begin(); ri !=realRoots.end(); ri++) {
789        r2 = *ri;
790        //check if FindRealRoots() give the right answer
# Line 744 | Line 800 | namespace OpenMD {
800        if (alpha0 >= 0.0) {
801          r1 = sqrt(alpha0 / a000);
802          if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) < 1e-6)
803 <          { rps.push_back(std::make_pair(r1, r2)); }
803 >          { rps.push_back(make_pair(r1, r2)); }
804          if (r1 > 1e-6) { //r1 non-negative
805            r1 = -r1;
806            if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) <1e-6)
807 <            { rps.push_back(std::make_pair(r1, r2)); }
807 >            { rps.push_back(make_pair(r1, r2)); }
808          }
809        }
810      }
811 <    // Consider combininig together the solving pair part w/ the searching
811 >    // Consider combining together the solving pair part w/ the searching
812      // best solution part so that we don't need the pairs vector
813      if (!rps.empty()) {
814        RealType smallestDiff = HONKING_LARGE_VALUE;
815        RealType diff;
816 <      std::pair<RealType,RealType> bestPair = std::make_pair(1.0, 1.0);
817 <      std::vector<std::pair<RealType,RealType> >::iterator rpi;
816 >      pair<RealType,RealType> bestPair = make_pair(1.0, 1.0);
817 >      vector<pair<RealType,RealType> >::iterator rpi;
818        for (rpi = rps.begin(); rpi != rps.end(); rpi++) {
819          r1 = (*rpi).first;
820          r2 = (*rpi).second;
# Line 777 | Line 833 | namespace OpenMD {
833              + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2);
834            break;
835          case rnemdPzScale :
836 +          diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
837 +            + fastpow(r1 * r1 / r2 / r2 - Kcy/Kcx, 2);
838          default :
839            break;
840          }
# Line 788 | Line 846 | namespace OpenMD {
846   #ifdef IS_MPI
847        if (worldRank == 0) {
848   #endif
849 <        std::cerr << "we choose r1 = " << bestPair.first
850 <                  << " and r2 = " << bestPair.second << "\n";
849 >        sprintf(painCave.errMsg,
850 >                "RNEMD: roots r1= %lf\tr2 = %lf\n",
851 >                bestPair.first, bestPair.second);
852 >        painCave.isFatal = 0;
853 >        painCave.severity = OPENMD_INFO;
854 >        simError();
855   #ifdef IS_MPI
856        }
857   #endif
858 <
858 >      
859        RealType x, y, z;
860 <        switch(rnemdType_) {
861 <        case rnemdKineticScale :
862 <          x = bestPair.first;
863 <          y = bestPair.first;
864 <          z = bestPair.second;
865 <          break;
866 <        case rnemdPxScale :
867 <          x = c;
868 <          y = bestPair.first;
869 <          z = bestPair.second;
870 <          break;
871 <        case rnemdPyScale :
872 <          x = bestPair.first;
873 <          y = c;
874 <          z = bestPair.second;
875 <          break;
876 <        case rnemdPzScale :
877 <          x = bestPair.first;
878 <          y = bestPair.second;
879 <          z = c;
880 <          break;          
881 <        default :
882 <          break;
883 <        }
884 <      std::vector<StuntDouble*>::iterator sdi;
860 >      switch(rnemdType_) {
861 >      case rnemdKineticScale :
862 >        x = bestPair.first;
863 >        y = bestPair.first;
864 >        z = bestPair.second;
865 >        break;
866 >      case rnemdPxScale :
867 >        x = c;
868 >        y = bestPair.first;
869 >        z = bestPair.second;
870 >        break;
871 >      case rnemdPyScale :
872 >        x = bestPair.first;
873 >        y = c;
874 >        z = bestPair.second;
875 >        break;
876 >      case rnemdPzScale :
877 >        x = bestPair.first;
878 >        y = bestPair.second;
879 >        z = c;
880 >        break;          
881 >      default :
882 >        break;
883 >      }
884 >      vector<StuntDouble*>::iterator sdi;
885        Vector3d vel;
886        for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
887          vel = (*sdi)->getVel();
# Line 842 | Line 904 | namespace OpenMD {
904        exchangeSum_ += targetFlux_;
905        //we may want to check whether the exchange has been successful
906      } else {
907 <      std::cerr << "exchange NOT performed!\n";//MPI incompatible
907 >      sprintf(painCave.errMsg,
908 >              "RNEMD: exchange NOT performed!\n");
909 >      painCave.isFatal = 0;
910 >      painCave.severity = OPENMD_INFO;
911 >      simError();        
912        failTrialCount_++;
913      }
914  
# Line 880 | Line 946 | namespace OpenMD {
946      StuntDouble* sd;
947      int idx;
948  
949 +    // alternative approach, track all molecules instead of only those
950 +    // selected for scaling/swapping:
951 +    /*
952 +    SimInfo::MoleculeIterator miter;
953 +    vector<StuntDouble*>::iterator iiter;
954 +    Molecule* mol;
955 +    StuntDouble* integrableObject;
956 +    for (mol = info_->beginMolecule(miter); mol != NULL;
957 +         mol = info_->nextMolecule(miter))
958 +      integrableObject is essentially sd
959 +        for (integrableObject = mol->beginIntegrableObject(iiter);
960 +             integrableObject != NULL;
961 +             integrableObject = mol->nextIntegrableObject(iiter))
962 +    */
963      for (sd = seleMan_.beginSelected(selei); sd != NULL;
964           sd = seleMan_.nextSelected(selei)) {
965        
# Line 895 | Line 975 | namespace OpenMD {
975        // which bin is this stuntdouble in?
976        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
977        
978 <      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
979 <
978 >      int binNo = int(rnemdLogWidth_ * (pos.z() / hmat(2,2) + 0.5)) %
979 >        rnemdLogWidth_;
980 >      // no symmetrization allowed due to arbitary rnemdLogWidth_ value
981 >      /*
982        if (rnemdLogWidth_ == midBin_ + 1)
983          if (binNo > midBin_)
984            binNo = nBins_ - binNo;
985 <
985 >      */
986        RealType mass = sd->getMass();
987        Vector3d vel = sd->getVel();
988        RealType value;
# Line 910 | Line 992 | namespace OpenMD {
992        case rnemdKineticSwap :
993        case rnemdKineticScale :
994          
995 <        value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
914 <                        vel[2]*vel[2]);
995 >        value = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
996          
997          valueCount_[binNo] += 3;
998          if (sd->isDirectional()) {
# Line 924 | Line 1005 | namespace OpenMD {
1005              int k = (i + 2) % 3;
1006              value += angMom[j] * angMom[j] / I(j, j) +
1007                angMom[k] * angMom[k] / I(k, k);
1008 <
1008 >            
1009              valueCount_[binNo] +=2;
1010 <
1010 >            
1011            } else {
1012              value += angMom[0]*angMom[0]/I(0, 0)
1013                + angMom[1]*angMom[1]/I(1, 1)
# Line 935 | Line 1016 | namespace OpenMD {
1016            }
1017          }
1018          value = value / PhysicalConstants::energyConvert / PhysicalConstants::kb;
1019 <
1019 >        
1020          break;
1021        case rnemdPx :
1022        case rnemdPxScale :
1023          value = mass * vel[0];
1024          valueCount_[binNo]++;
944 –        xVal = mass * vel.x() * vel.x() / PhysicalConstants::energyConvert
945 –          / PhysicalConstants::kb;
946 –        yVal = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert
947 –          / PhysicalConstants::kb;
948 –        zVal = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert
949 –          / PhysicalConstants::kb;
950 –        xTempHist_[binNo] += xVal;
951 –        yTempHist_[binNo] += yVal;
952 –        zTempHist_[binNo] += zVal;
1025          break;
1026        case rnemdPy :
1027        case rnemdPyScale :
# Line 958 | Line 1030 | namespace OpenMD {
1030          break;
1031        case rnemdPz :
1032        case rnemdPzScale :
1033 <        value = mass * vel[2];
1033 >        value = pos.z(); //temporarily for homogeneous systems ONLY
1034          valueCount_[binNo]++;
1035          break;
1036        case rnemdUnknown :
1037        default :
1038 +        value = 1.0;
1039 +        valueCount_[binNo]++;
1040          break;
1041        }
1042        valueHist_[binNo] += value;
969 –    }
1043  
1044 +      if (output3DTemp_) {
1045 +        xVal = mass * vel.x() * vel.x() / PhysicalConstants::energyConvert
1046 +          / PhysicalConstants::kb;
1047 +        yVal = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert
1048 +          / PhysicalConstants::kb;
1049 +        zVal = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert
1050 +          / PhysicalConstants::kb;
1051 +        xTempHist_[binNo] += xVal;
1052 +        yTempHist_[binNo] += yVal;
1053 +        zTempHist_[binNo] += zVal;
1054 +        xyzTempCount_[binNo]++;
1055 +      }
1056 +    }
1057    }
1058  
1059    void RNEMD::getStarted() {
1060 <    Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1061 <    Stats& stat = currentSnap_->statData;
1062 <    stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
1060 >    collectData();
1061 >    /* now should be able to output profile in step 0, but might not be useful
1062 >       Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1063 >       Stats& stat = currentSnap_->statData;
1064 >       stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
1065 >    */
1066 >    getStatus();
1067    }
1068  
1069    void RNEMD::getStatus() {
# Line 995 | Line 1085 | namespace OpenMD {
1085                                rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1086      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount_[0],
1087                                rnemdLogWidth_, MPI::INT, MPI::SUM);
1088 <    if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale) {
1088 >    if (output3DTemp_) {
1089        MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xTempHist_[0],
1090                                  rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1091        MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &yTempHist_[0],
1092                                  rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1093        MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &zTempHist_[0],
1094                                  rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1095 +      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xyzTempCount_[0],
1096 +                                rnemdLogWidth_, MPI::INT, MPI::SUM);
1097      }
1098      // If we're the root node, should we print out the results
1099      int worldRank = MPI::COMM_WORLD.Get_rank();
# Line 1012 | Line 1104 | namespace OpenMD {
1104          rnemdLog_ << "\t" << valueHist_[j] / (RealType)valueCount_[j];
1105        }
1106        rnemdLog_ << "\n";
1107 <      if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale ) {
1107 >      if (output3DTemp_) {
1108          xTempLog_ << time;      
1109          for (j = 0; j < rnemdLogWidth_; j++) {
1110 <          xTempLog_ << "\t" << xTempHist_[j] / (RealType)valueCount_[j];
1110 >          xTempLog_ << "\t" << xTempHist_[j] / (RealType)xyzTempCount_[j];
1111          }
1112          xTempLog_ << "\n";
1113          yTempLog_ << time;
1114          for (j = 0; j < rnemdLogWidth_; j++) {
1115 <          yTempLog_ << "\t" << yTempHist_[j] / (RealType)valueCount_[j];
1115 >          yTempLog_ << "\t" << yTempHist_[j] / (RealType)xyzTempCount_[j];
1116          }
1117          yTempLog_ << "\n";
1118          zTempLog_ << time;
1119          for (j = 0; j < rnemdLogWidth_; j++) {
1120 <          zTempLog_ << "\t" << zTempHist_[j] / (RealType)valueCount_[j];
1120 >          zTempLog_ << "\t" << zTempHist_[j] / (RealType)xyzTempCount_[j];
1121          }
1122          zTempLog_ << "\n";
1123        }
# Line 1036 | Line 1128 | namespace OpenMD {
1128        valueCount_[j] = 0;
1129        valueHist_[j] = 0.0;
1130      }
1131 <    if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale)
1131 >    if (output3DTemp_)
1132        for (j = 0; j < rnemdLogWidth_; j++) {
1133          xTempHist_[j] = 0.0;
1134          yTempHist_[j] = 0.0;
1135          zTempHist_[j] = 0.0;
1136 +        xyzTempCount_[j] = 0;
1137        }
1138    }
1139   }

Diff Legend

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