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 1339 by gezelter, Thu Apr 23 18:31:05 2009 UTC vs.
Revision 1341 by skuang, Fri May 8 19:47:05 2009 UTC

# Line 53 | Line 53
53   #include "math/ParallelRandNumGen.hpp"
54   #endif
55  
56 – /* Remove me after testing*/
57 – /*
58 – #include <cstdio>
59 – #include <iostream>
60 – */
61 – /*End remove me*/
56  
57   namespace oopse {
58    
# Line 74 | Line 68 | namespace oopse {
68      stringToEnumMap_["Unknown"] = rnemdUnknown;
69  
70      rnemdObjectSelection_ = simParams->getRNEMD_objectSelection();
77 –
78 –    std::cerr << "calling  evaluator with " << rnemdObjectSelection_ << "\n";
71      evaluator_.loadScriptString(rnemdObjectSelection_);
72 <    std::cerr << "done\n";
72 >    seleMan_.setSelectionSet(evaluator_.evaluate());
73 >
74 >
75 >    // do some sanity checking
76 >
77 >    int selectionCount = seleMan_.getSelectionCount();
78 >    int nIntegrable = info->getNGlobalIntegrableObjects();
79 >
80 >    if (selectionCount > nIntegrable) {
81 >      sprintf(painCave.errMsg,
82 >              "RNEMD warning: The current RNEMD_objectSelection,\n"
83 >              "\t\t%s\n"
84 >              "\thas resulted in %d selected objects.  However,\n"
85 >              "\tthe total number of integrable objects in the system\n"
86 >              "\tis only %d.  This is almost certainly not what you want\n"
87 >              "\tto do.  A likely cause of this is forgetting the _RB_0\n"
88 >              "\tselector in the selection script!\n",
89 >              rnemdObjectSelection_.c_str(),
90 >              selectionCount, nIntegrable);
91 >      painCave.isFatal = 0;
92 >      simError();
93 >
94 >    }
95      
96      const std::string st = simParams->getRNEMD_swapType();
97  
# Line 89 | Line 103 | namespace oopse {
103      set_RNEMD_nBins(simParams->getRNEMD_nBins());
104      exchangeSum_ = 0.0;
105      counter_ = 0; //added by shenyu
92 –    //profile_.open("profile", std::ios::out);
106  
107   #ifndef IS_MPI
108      if (simParams->haveSeed()) {
# Line 110 | Line 123 | namespace oopse {
123    
124    RNEMD::~RNEMD() {
125      delete randNumGen_;
113 –    //profile_.close();
126    }
127  
128    void RNEMD::doSwap() {
117 –    //std::cerr << "in RNEMD!\n";  
118 –    //std::cerr << "nBins = " << nBins_ << "\n";
129      int midBin = nBins_ / 2;
120 –    //std::cerr << "midBin = " << midBin << "\n";
121 –    //std::cerr << "swapTime = " << swapTime_ << "\n";
122 –    //std::cerr << "swapType = " << rnemdType_ << "\n";
123 –    //std::cerr << "selection = " << rnemdObjectSelection_ << "\n";
130  
131      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
132      Mat3x3d hmat = currentSnap_->getHmat();
133  
128 –    //std::cerr << "hmat = " << hmat << "\n";
129 –
134      seleMan_.setSelectionSet(evaluator_.evaluate());
135  
132 –    //std::cerr << "selectionCount = " << seleMan_.getSelectionCount() << "\n\n";
133 –
136      int selei;
137      StuntDouble* sd;
138      int idx;
# Line 158 | Line 160 | namespace oopse {
160        // which bin is this stuntdouble in?
161        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
162  
163 <      int binNo = int((nBins_-1) * (pos.z() + 0.5*hmat(2,2)) / hmat(2,2));
163 >      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
164  
163 –      //std::cerr << "pos.z() = " << pos.z() << " bin = " << binNo << "\n";
165  
166        // if we're in bin 0 or the middleBin
167        if (binNo == 0 || binNo == midBin) {
# Line 174 | Line 175 | namespace oopse {
175            
176            value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
177                            vel[2]*vel[2]);
177 –          
178            if (sd->isDirectional()) {
179              Vector3d angMom = sd->getJ();
180              Mat3x3d I = sd->getI();
# Line 232 | Line 232 | namespace oopse {
232          }
233        }
234      }
235 <    //std::cerr << "smallest value = " << min_val  << "\n";
236 <    //std::cerr << "largest value = " << max_val << "\n";
237 <    
235 >
236      // missing:  swap information in parallel
237  
238      if (max_found && min_found) {
239        if (min_val< max_val) {
240 +
241          Vector3d min_vel = min_sd->getVel();
242          Vector3d max_vel = max_sd->getVel();
243          RealType temp_vel;
244 +
245          switch(rnemdType_) {
246          case rnemdKinetic :
247            min_sd->setVel(max_vel);
248 <          max_sd->setVel(min_vel);                
248 >          max_sd->setVel(min_vel);
249 >
250            if (min_sd->isDirectional() && max_sd->isDirectional()) {
251              Vector3d min_angMom = min_sd->getJ();
252              Vector3d max_angMom = max_sd->getJ();
253              min_sd->setJ(max_angMom);
254              max_sd->setJ(min_angMom);
255 <          }
255 >          }
256            break;
257          case rnemdPx :
258            temp_vel = min_vel.x();
# Line 285 | Line 286 | namespace oopse {
286      } else {
287        std::cerr << "exchange NOT performed.\none of the two slabs empty.\n";
288      }
288 –    std::cerr << "exchangeSum = " << exchangeSum_ << "\n";
289    }
290  
291    void RNEMD::getStatus() {
292 –    //std::cerr << "in RNEMD!\n";  
293 –    //std::cerr << "nBins = " << nBins_ << "\n";
294 –    int midBin = nBins_ / 2;
295 –    //std::cerr << "midBin = " << midBin << "\n";
296 –    //std::cerr << "swapTime = " << swapTime_ << "\n";
297 –    //std::cerr << "exchangeSum = " << exchangeSum_ << "\n";
298 –    //std::cerr << "swapType = " << rnemdType_ << "\n";
299 –    //std::cerr << "selection = " << rnemdObjectSelection_ << "\n";
292  
293      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
294      Mat3x3d hmat = currentSnap_->getHmat();
295 +    Stats& stat = currentSnap_->statData;
296  
297 <    //std::cerr << "hmat = " << hmat << "\n";
297 >    stat[Stats::RNEMD_SWAP_TOTAL] = exchangeSum_;
298  
299      seleMan_.setSelectionSet(evaluator_.evaluate());
300  
308 –    //std::cerr << "selectionCount = " << seleMan_.getSelectionCount() << "\n\n";
309 –
301      int selei;
302      StuntDouble* sd;
303      int idx;
313 –    /*
314 –    RealType min_val;
315 –    bool min_found = false;  
316 –    StuntDouble* min_sd;
304  
318 –    RealType max_val;
319 –    bool max_found = false;
320 –    StuntDouble* max_sd;
321 –    */
305      std::vector<RealType> valueHist;  // keeps track of what's being averaged
306      std::vector<int> valueCount; // keeps track of the number of degrees of
307                                   // freedom being averaged
# Line 340 | Line 323 | namespace oopse {
323        // which bin is this stuntdouble in?
324        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
325        
326 <      int binNo = int((nBins_-1) * (pos.z()+0.5*hmat(2,2)) / hmat(2,2));
326 >      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
327        
328        //std::cerr << "pos.z() = " << pos.z() << " bin = " << binNo << "\n";
329        
# Line 356 | Line 339 | namespace oopse {
339                          vel[2]*vel[2]);
340          
341          valueCount[binNo] += 3;
342 <
342 >        //std::cerr <<"starting value = " << value << "\n";
343          if (sd->isDirectional()) {
344 +          //std::cerr << "angMom calculated.\n";
345            Vector3d angMom = sd->getJ();
346 +          //std::cerr << "current angMom: " << angMom << "\n";
347            Mat3x3d I = sd->getI();
348            
349            if (sd->isLinear()) {
# Line 370 | Line 355 | namespace oopse {
355  
356              valueCount[binNo] +=2;
357  
358 <          } else {                        
358 >          } else {
359 >            //std::cerr << "non-linear molecule.\n";
360              value += angMom[0]*angMom[0]/I(0, 0)
361                + angMom[1]*angMom[1]/I(1, 1)
362                + angMom[2]*angMom[2]/I(2, 2);
# Line 378 | Line 364 | namespace oopse {
364  
365            }
366          }
367 <        //std::cerr <<"this value = " << value << "\n";
368 <        value *= 0.5 / OOPSEConstant::energyConvert;  // get it in kcal / mol
369 <        value *= 2.0 / OOPSEConstant::kb;             // convert to temperature
370 <        //std::cerr <<"this value = " << value << "\n";
367 >        //std::cerr <<"total value = " << value << "\n";
368 >        //value *= 0.5 / OOPSEConstant::energyConvert;  // get it in kcal / mol
369 >        //value *= 2.0 / OOPSEConstant::kb;            // convert to temperature
370 >        value = value / OOPSEConstant::energyConvert / OOPSEConstant::kb;
371 >        //std::cerr <<"value = " << value << "\n";
372          break;
373        case rnemdPx :
374          value = mass * vel[0];

Diff Legend

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