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 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC vs.
Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC

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

Diff Legend

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