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 1333 by gezelter, Thu Apr 2 18:11:59 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    
59 <  RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info) {
59 >  RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
60      
61      int seedValue;
62      Globals * simParams = info->getSimParams();
# 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";
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 88 | Line 102 | namespace oopse {
102      set_RNEMD_swapTime(simParams->getRNEMD_swapTime());
103      set_RNEMD_nBins(simParams->getRNEMD_nBins());
104      exchangeSum_ = 0.0;
105 <    
105 >    counter_ = 0; //added by shenyu
106 >
107   #ifndef IS_MPI
108      if (simParams->haveSeed()) {
109        seedValue = simParams->getSeed();
# Line 111 | Line 126 | namespace oopse {
126    }
127  
128    void RNEMD::doSwap() {
114 –    std::cerr << "in RNEMD!\n";  
115 –    std::cerr << "nBins = " << nBins_ << "\n";
129      int midBin = nBins_ / 2;
117 –    std::cerr << "midBin = " << midBin << "\n";
118 –    std::cerr << "swapTime = " << swapTime_ << "\n";
119 –    std::cerr << "exchangeSum = " << exchangeSum_ << "\n";
120 –    std::cerr << "swapType = " << rnemdType_ << "\n";
121 –    std::cerr << "selection = " << rnemdObjectSelection_ << "\n";
130  
131      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
132      Mat3x3d hmat = currentSnap_->getHmat();
133  
126 –    std::cerr << "hmat = " << hmat << "\n";
127 –
134      seleMan_.setSelectionSet(evaluator_.evaluate());
135  
130 –    std::cerr << "selectionCount = " << seleMan_.getSelectionCount() << "\n\n";
131 –
136      int selei;
137      StuntDouble* sd;
138      int idx;
139  
140 <    std::vector<tuple3<RealType, int, StuntDouble* > > endSlice;
141 <    std::vector<tuple3<RealType, int, StuntDouble* > > midSlice;
142 <    
140 >    RealType min_val;
141 >    bool min_found = false;  
142 >    StuntDouble* min_sd;
143 >
144 >    RealType max_val;
145 >    bool max_found = false;
146 >    StuntDouble* max_sd;
147 >
148      for (sd = seleMan_.beginSelected(selei); sd != NULL;
149           sd = seleMan_.nextSelected(selei)) {
150  
# Line 143 | Line 152 | namespace oopse {
152  
153        Vector3d pos = sd->getPos();
154  
146 –      std::cerr << "idx = " << idx << "pos = " << pos << "\n";
155        // wrap the stuntdouble's position back into the box:
156  
157        if (usePeriodicBoundaryConditions_)
158          currentSnap_->wrapVector(pos);
151 –      std::cerr << "new pos.z = " << pos.z() << "\n";
159  
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_ * (pos.z()) / hmat(2,2));
163 >      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
164  
157 –      std::cerr << "bin = " << binNo << "\n";
165  
166        // if we're in bin 0 or the middleBin
167        if (binNo == 0 || binNo == midBin) {
# Line 168 | Line 175 | namespace oopse {
175            
176            value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
177                            vel[2]*vel[2]);
171 –          
178            if (sd->isDirectional()) {
179              Vector3d angMom = sd->getJ();
180              Mat3x3d I = sd->getI();
# Line 201 | Line 207 | namespace oopse {
207            break;
208          }
209          
210 <        if (binNo == 0)
211 <          endSlice.push_back( make_tuple3(value, idx, sd));
212 <        
213 <        if (binNo == midBin)
214 <          midSlice.push_back( make_tuple3(value, idx, sd));              
210 >        if (binNo == 0) {
211 >          if (!min_found) {
212 >            min_val = value;
213 >            min_sd = sd;
214 >            min_found = true;
215 >          } else {
216 >            if (min_val > value) {
217 >              min_val = value;
218 >              min_sd = sd;
219 >            }
220 >          }
221 >        } else {
222 >          if (!max_found) {
223 >            max_val = value;
224 >            max_sd = sd;
225 >            max_found = true;
226 >          } else {
227 >            if (max_val < value) {
228 >              max_val = value;
229 >              max_sd = sd;
230 >            }
231 >          }      
232 >        }
233        }
234 +    }
235  
236 <      // find smallest value in endSlice:
212 <      std::sort(endSlice.begin(), endSlice.end());    
213 <      RealType min_val = endSlice[0].first;
214 <      int min_idx = endSlice[0].second;
215 <      StuntDouble* min_sd = endSlice[0].third;
236 >    // missing:  swap information in parallel
237  
238 <      std::cerr << "smallest value = " << min_val << " idx = " << min_idx << "\n";
238 >    if (max_found && min_found) {
239 >      if (min_val< max_val) {
240  
241 <      // find largest value in midSlice:
242 <      std::sort(midSlice.rbegin(), midSlice.rend());
243 <      RealType max_val = midSlice[0].first;
222 <      int max_idx = midSlice[0].second;
223 <      StuntDouble* max_sd = midSlice[0].third;
241 >        Vector3d min_vel = min_sd->getVel();
242 >        Vector3d max_vel = max_sd->getVel();
243 >        RealType temp_vel;
244  
245 <      std::cerr << "largest value = " << max_val << " idx = " << max_idx << "\n";
245 >        switch(rnemdType_) {
246 >        case rnemdKinetic :
247 >          min_sd->setVel(max_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 >          }
256 >          break;
257 >        case rnemdPx :
258 >          temp_vel = min_vel.x();
259 >          min_vel.x() = max_vel.x();
260 >          max_vel.x() = temp_vel;
261 >          min_sd->setVel(min_vel);
262 >          max_sd->setVel(max_vel);
263 >          break;
264 >        case rnemdPy :
265 >          temp_vel = min_vel.y();
266 >          min_vel.y() = max_vel.y();
267 >          max_vel.y() = temp_vel;
268 >          min_sd->setVel(min_vel);
269 >          max_sd->setVel(max_vel);
270 >          break;
271 >        case rnemdPz :
272 >          temp_vel = min_vel.z();
273 >          min_vel.z() = max_vel.z();
274 >          max_vel.z() = temp_vel;
275 >          min_sd->setVel(min_vel);
276 >          max_sd->setVel(max_vel);
277 >          break;
278 >        case rnemdUnknown :
279 >        default :
280 >          break;
281 >        }
282 >      exchangeSum_ += max_val - min_val;
283 >      } else {
284 >        std::cerr << "exchange NOT performed.\nmin_val > max_val.\n";
285 >      }
286 >    } else {
287 >      std::cerr << "exchange NOT performed.\none of the two slabs empty.\n";
288 >    }
289 >  }
290  
291 +  void RNEMD::getStatus() {
292 +
293 +    Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
294 +    Mat3x3d hmat = currentSnap_->getHmat();
295 +    Stats& stat = currentSnap_->statData;
296 +
297 +    stat[Stats::RNEMD_SWAP_TOTAL] = exchangeSum_;
298 +
299 +    seleMan_.setSelectionSet(evaluator_.evaluate());
300 +
301 +    int selei;
302 +    StuntDouble* sd;
303 +    int idx;
304 +
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
308 +    valueHist.resize(nBins_);
309 +    valueCount.resize(nBins_);
310 +    //do they initialize themselves to zero automatically?
311 +    for (sd = seleMan_.beginSelected(selei); sd != NULL;
312 +         sd = seleMan_.nextSelected(selei)) {
313 +      
314 +      idx = sd->getLocalIndex();
315 +      
316 +      Vector3d pos = sd->getPos();
317 +
318 +      // wrap the stuntdouble's position back into the box:
319 +      
320 +      if (usePeriodicBoundaryConditions_)
321 +        currentSnap_->wrapVector(pos);
322 +      
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_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
327 +      
328 +      //std::cerr << "pos.z() = " << pos.z() << " bin = " << binNo << "\n";
329 +      
330 +      RealType mass = sd->getMass();
331 +      Vector3d vel = sd->getVel();
332 +      //std::cerr << "mass = " << mass << " vel = " << vel << "\n";
333 +      RealType value;
334 +
335 +      switch(rnemdType_) {
336 +      case rnemdKinetic :
337 +        
338 +        value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
339 +                        vel[2]*vel[2]);
340 +        
341 +        valueCount[binNo] += 3;
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()) {
350 +            int i = sd->linearAxis();
351 +            int j = (i + 1) % 3;
352 +            int k = (i + 2) % 3;
353 +            value += angMom[j] * angMom[j] / I(j, j) +
354 +              angMom[k] * angMom[k] / I(k, k);
355 +
356 +            valueCount[binNo] +=2;
357 +
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);
363 +            valueCount[binNo] +=3;
364 +
365 +          }
366 +        }
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];
375 +        valueCount[binNo]++;
376 +        break;
377 +      case rnemdPy :
378 +        value = mass * vel[1];
379 +        valueCount[binNo]++;
380 +        break;
381 +      case rnemdPz :
382 +        value = mass * vel[2];
383 +        valueCount[binNo]++;
384 +        break;
385 +      case rnemdUnknown :
386 +      default :
387 +        break;
388 +      }
389 +      //std::cerr << "bin = " << binNo << " value = " << value ;
390 +      valueHist[binNo] += value;
391 +      //std::cerr << " hist = " << valueHist[binNo] << " count = " << valueCount[binNo] << "\n";
392      }
393 <  }  
393 >    
394 >    std::cout << counter_++;
395 >    for (int j = 0; j < nBins_; j++)
396 >      std::cout << "\t" << valueHist[j] / (RealType)valueCount[j];
397 >    std::cout << "\n";
398 >  }
399   }

Diff Legend

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