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 1442 by gezelter, Mon May 10 17:28:26 2010 UTC vs.
Revision 1560 by skuang, Wed May 11 17:55:32 2011 UTC

# Line 111 | Line 111 | namespace OpenMD {
111        std::cerr << "WARNING! RNEMD Type Unknown!\n";
112      }
113  
114 +    output3DTemp_ = false;
115 +    if (simParams->haveRNEMD_outputDimensionalTemperature()) {
116 +      output3DTemp_ = simParams->getRNEMD_outputDimensionalTemperature();
117 +    }
118 +
119   #ifdef IS_MPI
120      if (worldRank == 0) {
121   #endif
122  
123        std::string rnemdFileName;
119 –      std::string xTempFileName;
120 –      std::string yTempFileName;
121 –      std::string zTempFileName;
124        switch(rnemdType_) {
125        case rnemdKineticSwap :
126        case rnemdKineticScale :
# Line 129 | Line 131 | namespace OpenMD {
131        case rnemdPy :
132        case rnemdPyScale :
133          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());
134          break;
135        case rnemdPz :
136        case rnemdPzScale :
# Line 145 | Line 141 | namespace OpenMD {
141        }
142        rnemdLog_.open(rnemdFileName.c_str());
143  
144 +      std::string xTempFileName;
145 +      std::string yTempFileName;
146 +      std::string zTempFileName;
147 +      if (output3DTemp_) {
148 +        xTempFileName = "temperatureX.log";
149 +        yTempFileName = "temperatureY.log";
150 +        zTempFileName = "temperatureZ.log";
151 +        xTempLog_.open(xTempFileName.c_str());
152 +        yTempLog_.open(yTempFileName.c_str());
153 +        zTempLog_.open(zTempFileName.c_str());
154 +      }
155 +
156   #ifdef IS_MPI
157      }
158   #endif
# Line 152 | Line 160 | namespace OpenMD {
160      set_RNEMD_exchange_time(simParams->getRNEMD_exchangeTime());
161      set_RNEMD_nBins(simParams->getRNEMD_nBins());
162      midBin_ = nBins_ / 2;
163 +    if (simParams->haveRNEMD_binShift()) {
164 +        if (simParams->getRNEMD_binShift()) {
165 +          zShift_ = 0.5 / (RealType)(nBins_);
166 +        } else {
167 +          zShift_ = 0.0;
168 +        }
169 +    } else {
170 +      zShift_ = 0.0;
171 +    }
172 +    //std::cerr << "we have zShift_ = " << zShift_ << "\n";
173 +    //shift slabs by half slab width, might be useful in heterogeneous systems
174 +    //set to 0.0 if not using it; can NOT be used in status output yet
175      if (simParams->haveRNEMD_logWidth()) {
176 <      rnemdLogWidth_ = simParams->getRNEMD_logWidth();
176 >      set_RNEMD_logWidth(simParams->getRNEMD_logWidth());
177 >      /*arbitary rnemdLogWidth_ no checking
178        if (rnemdLogWidth_ != nBins_ && rnemdLogWidth_ != midBin_ + 1) {
179          std::cerr << "WARNING! RNEMD_logWidth has abnormal value!\n";
180          std::cerr << "Automaically set back to default.\n";
181          rnemdLogWidth_ = nBins_;
182 <      }
182 >        }*/
183      } else {
184 <      rnemdLogWidth_ = nBins_;
184 >      set_RNEMD_logWidth(nBins_);
185      }
186      valueHist_.resize(rnemdLogWidth_, 0.0);
187      valueCount_.resize(rnemdLogWidth_, 0);
188      xTempHist_.resize(rnemdLogWidth_, 0.0);
189      yTempHist_.resize(rnemdLogWidth_, 0.0);
190      zTempHist_.resize(rnemdLogWidth_, 0.0);
191 +    xyzTempCount_.resize(rnemdLogWidth_, 0);
192  
193      set_RNEMD_exchange_total(0.0);
194      if (simParams->haveRNEMD_targetFlux()) {
# Line 202 | Line 224 | namespace OpenMD {
224        rnemdLog_.close();
225        if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPyScale)
226          std::cerr<< "total root-checking warnings: " << failRootCount_ << "\n";
227 <      if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPy || rnemdType_ == rnemdPyScale) {
227 >      if (output3DTemp_) {
228          xTempLog_.close();
229          yTempLog_.close();
230          zTempLog_.close();
# Line 246 | Line 268 | namespace OpenMD {
268        // which bin is this stuntdouble in?
269        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
270  
271 <      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
271 >      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_;
272  
273  
274        // if we're in bin 0 or the middleBin
# Line 261 | Line 283 | namespace OpenMD {
283            
284            value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
285                            vel[2]*vel[2]);
286 +          /*
287            if (sd->isDirectional()) {
288              Vector3d angMom = sd->getJ();
289              Mat3x3d I = sd->getI();
# Line 276 | Line 299 | namespace OpenMD {
299                  + angMom[1]*angMom[1]/I(1, 1)
300                  + angMom[2]*angMom[2]/I(2, 2);
301              }
302 <          }
302 >          } no exchange of angular momenta
303 >          */
304            //make exchangeSum_ comparable between swap & scale
305            //temporarily without using energyConvert
306            //value = value * 0.5 / PhysicalConstants::energyConvert;
# Line 330 | Line 354 | namespace OpenMD {
354      bool my_min_found = min_found;
355      bool my_max_found = max_found;
356  
357 <    // Even if we didn't find a minimum, did someone else?
358 <    MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found,
359 <                              1, MPI::BOOL, MPI::LAND);
360 <    
357 >    // Even if we didn't find a minimum, did someone else? debugging...
358 >    //MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found,
359 >    //                          1, MPI::BOOL, MPI::LAND);
360 >    MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, 1, MPI::BOOL, MPI::LOR);
361      // Even if we didn't find a maximum, did someone else?
362 <    MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found,
363 <                              1, MPI::BOOL, MPI::LAND);
364 <    
362 >    //MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found,
363 >    //                          1, MPI::BOOL, MPI::LAND);
364 >    MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, 1, MPI::BOOL, MPI::LOR);
365      struct {
366        RealType val;
367        int rank;
# Line 373 | Line 397 | namespace OpenMD {
397   #endif
398  
399      if (max_found && min_found) {
400 <      if (min_val< max_val) {
400 >      if (min_val < max_val) {
401  
402   #ifdef IS_MPI      
403          if (max_vals.rank == worldRank && min_vals.rank == worldRank) {
404            // I have both maximum and minimum, so proceed like a single
405            // processor version:
406   #endif
407 <          // objects to be swapped: velocity & angular velocity
407 >          // objects to be swapped: velocity ONLY
408            Vector3d min_vel = min_sd->getVel();
409            Vector3d max_vel = max_sd->getVel();
410            RealType temp_vel;
# Line 389 | Line 413 | namespace OpenMD {
413            case rnemdKineticSwap :
414              min_sd->setVel(max_vel);
415              max_sd->setVel(min_vel);
416 +            /*
417              if (min_sd->isDirectional() && max_sd->isDirectional()) {
418                Vector3d min_angMom = min_sd->getJ();
419                Vector3d max_angMom = max_sd->getJ();
420                min_sd->setJ(max_angMom);
421                max_sd->setJ(min_angMom);
422 <            }
422 >            } no angular momentum exchange
423 >            */
424              break;
425            case rnemdPx :
426              temp_vel = min_vel.x();
# Line 438 | Line 464 | namespace OpenMD {
464            switch(rnemdType_) {
465            case rnemdKineticSwap :
466              max_sd->setVel(min_vel);
467 <            
467 >            /*            
468              if (max_sd->isDirectional()) {
469                Vector3d min_angMom;
470                Vector3d max_angMom = max_sd->getJ();
# Line 451 | Line 477 | namespace OpenMD {
477                                         status);
478  
479                max_sd->setJ(min_angMom);
480 <            }
480 >            } no angular momentum exchange
481 >            */
482              break;
483            case rnemdPx :
484              max_vel.x() = min_vel.x();
# Line 484 | Line 511 | namespace OpenMD {
511            switch(rnemdType_) {
512            case rnemdKineticSwap :
513              min_sd->setVel(max_vel);
514 <            
514 >            /*            
515              if (min_sd->isDirectional()) {
516                Vector3d min_angMom = min_sd->getJ();
517                Vector3d max_angMom;
# Line 497 | Line 524 | namespace OpenMD {
524                                         status);
525  
526                min_sd->setJ(max_angMom);
527 <            }
527 >            } no angular momentum exchange
528 >            */
529              break;
530            case rnemdPx :
531              min_vel.x() = max_vel.x();
# Line 570 | Line 598 | namespace OpenMD {
598        // which bin is this stuntdouble in?
599        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
600  
601 <      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
601 >      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_;
602  
603        // if we're in bin 0 or the middleBin
604        if (binNo == 0 || binNo == midBin_) {
# Line 776 | Line 804 | namespace OpenMD {
804              + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2);
805            break;
806          case rnemdPzScale :
807 +          diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
808 +            + fastpow(r1 * r1 / r2 / r2 - Kcy/Kcx, 2);
809          default :
810            break;
811          }
# Line 879 | Line 909 | namespace OpenMD {
909      StuntDouble* sd;
910      int idx;
911  
912 +    // alternative approach, track all molecules instead of only those selected for scaling/swapping
913 +    //SimInfo::MoleculeIterator miter;
914 +    //std::vector<StuntDouble*>::iterator iiter;
915 +    //Molecule* mol;
916 +    //StuntDouble* integrableObject;
917 +    //for (mol = info_->beginMolecule(miter); mol != NULL;
918 +    //     mol = info_->nextMolecule(miter))
919 +    // integrableObject is essentially sd
920 +    //for (integrableObject = mol->beginIntegrableObject(iiter);
921 +    //     integrableObject != NULL;
922 +    //     integrableObject = mol->nextIntegrableObject(iiter))
923      for (sd = seleMan_.beginSelected(selei); sd != NULL;
924           sd = seleMan_.nextSelected(selei)) {
925        
# Line 894 | Line 935 | namespace OpenMD {
935        // which bin is this stuntdouble in?
936        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
937        
938 <      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
939 <
938 >      int binNo = int(rnemdLogWidth_ * (pos.z() / hmat(2,2) + 0.5)) %
939 >        rnemdLogWidth_;
940 >      /* no symmetrization allowed due to arbitary rnemdLogWidth_ value
941        if (rnemdLogWidth_ == midBin_ + 1)
942          if (binNo > midBin_)
943            binNo = nBins_ - binNo;
944 <
944 >      */
945        RealType mass = sd->getMass();
946        Vector3d vel = sd->getVel();
947        RealType value;
# Line 940 | Line 982 | namespace OpenMD {
982        case rnemdPxScale :
983          value = mass * vel[0];
984          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;
985          break;
986        case rnemdPy :
987        case rnemdPyScale :
# Line 957 | Line 990 | namespace OpenMD {
990          break;
991        case rnemdPz :
992        case rnemdPzScale :
993 <        value = mass * vel[2];
993 >        value = pos.z(); //temporarily for homogeneous systems ONLY
994          valueCount_[binNo]++;
995          break;
996        case rnemdUnknown :
997        default :
998 +        value = 1.0;
999 +        valueCount_[binNo]++;
1000          break;
1001        }
1002        valueHist_[binNo] += value;
968 –    }
1003  
1004 +      if (output3DTemp_) {
1005 +        xVal = mass * vel.x() * vel.x() / PhysicalConstants::energyConvert
1006 +          / PhysicalConstants::kb;
1007 +        yVal = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert
1008 +          / PhysicalConstants::kb;
1009 +        zVal = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert
1010 +          / PhysicalConstants::kb;
1011 +        xTempHist_[binNo] += xVal;
1012 +        yTempHist_[binNo] += yVal;
1013 +        zTempHist_[binNo] += zVal;
1014 +        xyzTempCount_[binNo]++;
1015 +      }
1016 +    }
1017    }
1018  
1019    void RNEMD::getStarted() {
1020 +    collectData();
1021 +    /* now should be able to output profile in step 0, but might not be useful
1022      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1023      Stats& stat = currentSnap_->statData;
1024      stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
1025 +    */
1026 +    getStatus();
1027    }
1028  
1029    void RNEMD::getStatus() {
# Line 994 | Line 1045 | namespace OpenMD {
1045                                rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1046      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount_[0],
1047                                rnemdLogWidth_, MPI::INT, MPI::SUM);
1048 <    if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale) {
1048 >    if (output3DTemp_) {
1049        MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xTempHist_[0],
1050                                  rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1051        MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &yTempHist_[0],
1052                                  rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1053        MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &zTempHist_[0],
1054                                  rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1055 +      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xyzTempCount_[0],
1056 +                                rnemdLogWidth_, MPI::INT, MPI::SUM);
1057      }
1058      // If we're the root node, should we print out the results
1059      int worldRank = MPI::COMM_WORLD.Get_rank();
# Line 1011 | Line 1064 | namespace OpenMD {
1064          rnemdLog_ << "\t" << valueHist_[j] / (RealType)valueCount_[j];
1065        }
1066        rnemdLog_ << "\n";
1067 <      if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale ) {
1067 >      if (output3DTemp_) {
1068          xTempLog_ << time;      
1069          for (j = 0; j < rnemdLogWidth_; j++) {
1070 <          xTempLog_ << "\t" << xTempHist_[j] / (RealType)valueCount_[j];
1070 >          xTempLog_ << "\t" << xTempHist_[j] / (RealType)xyzTempCount_[j];
1071          }
1072          xTempLog_ << "\n";
1073          yTempLog_ << time;
1074          for (j = 0; j < rnemdLogWidth_; j++) {
1075 <          yTempLog_ << "\t" << yTempHist_[j] / (RealType)valueCount_[j];
1075 >          yTempLog_ << "\t" << yTempHist_[j] / (RealType)xyzTempCount_[j];
1076          }
1077          yTempLog_ << "\n";
1078          zTempLog_ << time;
1079          for (j = 0; j < rnemdLogWidth_; j++) {
1080 <          zTempLog_ << "\t" << zTempHist_[j] / (RealType)valueCount_[j];
1080 >          zTempLog_ << "\t" << zTempHist_[j] / (RealType)xyzTempCount_[j];
1081          }
1082          zTempLog_ << "\n";
1083        }
# Line 1035 | Line 1088 | namespace OpenMD {
1088        valueCount_[j] = 0;
1089        valueHist_[j] = 0.0;
1090      }
1091 <    if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale)
1091 >    if (output3DTemp_)
1092        for (j = 0; j < rnemdLogWidth_; j++) {
1093          xTempHist_[j] = 0.0;
1094          yTempHist_[j] = 0.0;
1095          zTempHist_[j] = 0.0;
1096 +        xyzTempCount_[j] = 0;
1097        }
1098    }
1099   }

Diff Legend

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