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 1331 by gezelter, Thu Apr 2 16:04:52 2009 UTC vs.
Revision 1350 by gezelter, Thu May 21 18:56:45 2009 UTC

# Line 40 | Line 40
40   */
41  
42   #include "integrators/RNEMD.hpp"
43 + #include "math/Vector3.hpp"
44   #include "math/SquareMatrix3.hpp"
45   #include "primitives/Molecule.hpp"
46   #include "primitives/StuntDouble.hpp"
47 + #include "utils/OOPSEConstant.hpp"
48 + #include "utils/Tuple.hpp"
49  
50   #ifndef IS_MPI
51   #include "math/SeqRandNumGen.hpp"
# Line 50 | Line 53
53   #include "math/ParallelRandNumGen.hpp"
54   #endif
55  
56 < /* Remove me after testing*/
54 < /*
55 < #include <cstdio>
56 < #include <iostream>
57 < */
58 < /*End remove me*/
56 > #define HONKING_LARGE_VALUE 1.0e10
57  
58   namespace oopse {
59    
60 <  RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info) {
60 >  RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
61      
62      int seedValue;
63      Globals * simParams = info->getSimParams();
# Line 71 | Line 69 | namespace oopse {
69      stringToEnumMap_["Unknown"] = rnemdUnknown;
70  
71      rnemdObjectSelection_ = simParams->getRNEMD_objectSelection();
74 –
75 –    std::cerr << "calling  evaluator with " << rnemdObjectSelection_ << "\n";
72      evaluator_.loadScriptString(rnemdObjectSelection_);
73 <    std::cerr << "done";
73 >    seleMan_.setSelectionSet(evaluator_.evaluate());
74 >
75 >
76 >    // do some sanity checking
77 >
78 >    int selectionCount = seleMan_.getSelectionCount();
79 >    int nIntegrable = info->getNGlobalIntegrableObjects();
80 >
81 >    if (selectionCount > nIntegrable) {
82 >      sprintf(painCave.errMsg,
83 >              "RNEMD warning: The current RNEMD_objectSelection,\n"
84 >              "\t\t%s\n"
85 >              "\thas resulted in %d selected objects.  However,\n"
86 >              "\tthe total number of integrable objects in the system\n"
87 >              "\tis only %d.  This is almost certainly not what you want\n"
88 >              "\tto do.  A likely cause of this is forgetting the _RB_0\n"
89 >              "\tselector in the selection script!\n",
90 >              rnemdObjectSelection_.c_str(),
91 >              selectionCount, nIntegrable);
92 >      painCave.isFatal = 0;
93 >      simError();
94 >
95 >    }
96      
97      const std::string st = simParams->getRNEMD_swapType();
98  
# Line 85 | Line 103 | namespace oopse {
103      set_RNEMD_swapTime(simParams->getRNEMD_swapTime());
104      set_RNEMD_nBins(simParams->getRNEMD_nBins());
105      exchangeSum_ = 0.0;
106 <    
106 >
107   #ifndef IS_MPI
108      if (simParams->haveSeed()) {
109        seedValue = simParams->getSeed();
# Line 108 | Line 126 | namespace oopse {
126    }
127  
128    void RNEMD::doSwap() {
129 <    std::cerr << "in RNEMD!\n";  
112 <    std::cerr << "nBins = " << nBins_ << "\n";
113 <    std::cerr << "swapTime = " << swapTime_ << "\n";
114 <    std::cerr << "exchangeSum = " << exchangeSum_ << "\n";
115 <    std::cerr << "swapType = " << rnemdType_ << "\n";
116 <    std::cerr << "selection = " << rnemdObjectSelection_ << "\n";
129 >    int midBin = nBins_ / 2;
130  
131 +    Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
132 +    Mat3x3d hmat = currentSnap_->getHmat();
133 +
134      seleMan_.setSelectionSet(evaluator_.evaluate());
135  
136 <    std::cerr << "selectionCount = " << seleMan_.getSelectionCount() << "\n\n";
121 <
122 <    int i;
136 >    int selei;
137      StuntDouble* sd;
138 +    int idx;
139  
140 <    for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
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 >
151 >      idx = sd->getLocalIndex();
152 >
153        Vector3d pos = sd->getPos();
154 <      //wrap the stuntdoubles into a cell      
154 >
155 >      // wrap the stuntdouble's position back into the box:
156 >
157        if (usePeriodicBoundaryConditions_)
158 <        info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
159 <      int binNo = int(nBins_ * (pos.z()) / hmat(2,2));
160 <      sliceSDLists_[binNo].push_back(sd);
158 >        currentSnap_->wrapVector(pos);
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) + 0.5)) % nBins_;
164 >
165 >
166 >      // if we're in bin 0 or the middleBin
167 >      if (binNo == 0 || binNo == midBin) {
168 >        
169 >        RealType mass = sd->getMass();
170 >        Vector3d vel = sd->getVel();
171 >        RealType value;
172 >
173 >        switch(rnemdType_) {
174 >        case rnemdKinetic :
175 >          
176 >          value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
177 >                          vel[2]*vel[2]);
178 >          if (sd->isDirectional()) {
179 >            Vector3d angMom = sd->getJ();
180 >            Mat3x3d I = sd->getI();
181 >            
182 >            if (sd->isLinear()) {
183 >              int i = sd->linearAxis();
184 >              int j = (i + 1) % 3;
185 >              int k = (i + 2) % 3;
186 >              value += angMom[j] * angMom[j] / I(j, j) +
187 >                angMom[k] * angMom[k] / I(k, k);
188 >            } else {                        
189 >              value += angMom[0]*angMom[0]/I(0, 0)
190 >                + angMom[1]*angMom[1]/I(1, 1)
191 >                + angMom[2]*angMom[2]/I(2, 2);
192 >            }
193 >          }
194 >          value = value * 0.5 / OOPSEConstant::energyConvert;
195 >          break;
196 >        case rnemdPx :
197 >          value = mass * vel[0];
198 >          break;
199 >        case rnemdPy :
200 >          value = mass * vel[1];
201 >          break;
202 >        case rnemdPz :
203 >          value = mass * vel[2];
204 >          break;
205 >        case rnemdUnknown :
206 >        default :
207 >          break;
208 >        }
209 >        
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 <  }  
236 > #ifdef IS_MPI
237 >    int nProc, worldRank;
238 >
239 >    nProc = MPI::COMM_WORLD.Get_size();
240 >    worldRank = MPI::COMM_WORLD.Get_rank();
241 >
242 >    bool my_min_found = min_found;
243 >    bool my_max_found = max_found;
244 >
245 >    // Even if we didn't find a minimum, did someone else?
246 >    MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found,
247 >                              1, MPI::BOOL, MPI::LAND);
248 >    
249 >    // Even if we didn't find a maximum, did someone else?
250 >    MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found,
251 >                              1, MPI::BOOL, MPI::LAND);
252 >    
253 >    struct {
254 >      RealType val;
255 >      int rank;
256 >    } max_vals, min_vals;
257 >    
258 >    if (min_found) {
259 >      if (my_min_found)
260 >        min_vals.val = min_val;
261 >      else
262 >        min_vals.val = HONKING_LARGE_VALUE;
263 >      
264 >      min_vals.rank = worldRank;    
265 >      
266 >      // Who had the minimum?
267 >      MPI::COMM_WORLD.Allreduce(&min_vals, &min_vals,
268 >                                1, MPI::REALTYPE_INT, MPI::MINLOC);
269 >      min_val = min_vals.val;
270 >    }
271 >      
272 >    if (max_found) {
273 >      if (my_max_found)
274 >        max_vals.val = max_val;
275 >      else
276 >        max_vals.val = -HONKING_LARGE_VALUE;
277 >      
278 >      max_vals.rank = worldRank;    
279 >      
280 >      // Who had the maximum?
281 >      MPI::COMM_WORLD.Allreduce(&max_vals, &max_vals,
282 >                                1, MPI::REALTYPE_INT, MPI::MAXLOC);
283 >      max_val = max_vals.val;
284 >    }
285 > #endif
286 >
287 >    if (max_found && min_found) {
288 >      if (min_val< max_val) {
289 >
290 > #ifdef IS_MPI      
291 >        if (max_vals.rank == worldRank && min_vals.rank == worldRank) {
292 >          // I have both maximum and minimum, so proceed like a single
293 >          // processor version:
294 > #endif
295 >          // objects to be swapped: velocity & angular velocity
296 >          Vector3d min_vel = min_sd->getVel();
297 >          Vector3d max_vel = max_sd->getVel();
298 >          RealType temp_vel;
299 >          
300 >          switch(rnemdType_) {
301 >          case rnemdKinetic :
302 >            min_sd->setVel(max_vel);
303 >            max_sd->setVel(min_vel);
304 >            if (min_sd->isDirectional() && max_sd->isDirectional()) {
305 >              Vector3d min_angMom = min_sd->getJ();
306 >              Vector3d max_angMom = max_sd->getJ();
307 >              min_sd->setJ(max_angMom);
308 >              max_sd->setJ(min_angMom);
309 >            }
310 >            break;
311 >          case rnemdPx :
312 >            temp_vel = min_vel.x();
313 >            min_vel.x() = max_vel.x();
314 >            max_vel.x() = temp_vel;
315 >            min_sd->setVel(min_vel);
316 >            max_sd->setVel(max_vel);
317 >            break;
318 >          case rnemdPy :
319 >            temp_vel = min_vel.y();
320 >            min_vel.y() = max_vel.y();
321 >            max_vel.y() = temp_vel;
322 >            min_sd->setVel(min_vel);
323 >            max_sd->setVel(max_vel);
324 >            break;
325 >          case rnemdPz :
326 >            temp_vel = min_vel.z();
327 >            min_vel.z() = max_vel.z();
328 >            max_vel.z() = temp_vel;
329 >            min_sd->setVel(min_vel);
330 >            max_sd->setVel(max_vel);
331 >            break;
332 >          case rnemdUnknown :
333 >          default :
334 >            break;
335 >          }
336 > #ifdef IS_MPI
337 >          // the rest of the cases only apply in parallel simulations:
338 >        } else if (max_vals.rank == worldRank) {
339 >          // I had the max, but not the minimum
340 >          
341 >          Vector3d min_vel;
342 >          Vector3d max_vel = max_sd->getVel();
343 >          MPI::Status status;
344 >
345 >          // point-to-point swap of the velocity vector
346 >          MPI::COMM_WORLD.Sendrecv(max_vel.getArrayPointer(), 3, MPI::REALTYPE,
347 >                                   min_vals.rank, 0,
348 >                                   min_vel.getArrayPointer(), 3, MPI::REALTYPE,
349 >                                   min_vals.rank, 0, status);
350 >          
351 >          switch(rnemdType_) {
352 >          case rnemdKinetic :
353 >            max_sd->setVel(min_vel);
354 >            
355 >            if (max_sd->isDirectional()) {
356 >              Vector3d min_angMom;
357 >              Vector3d max_angMom = max_sd->getJ();
358 >
359 >              // point-to-point swap of the angular momentum vector
360 >              MPI::COMM_WORLD.Sendrecv(max_angMom.getArrayPointer(), 3,
361 >                                       MPI::REALTYPE, min_vals.rank, 1,
362 >                                       min_angMom.getArrayPointer(), 3,
363 >                                       MPI::REALTYPE, min_vals.rank, 1,
364 >                                       status);
365 >
366 >              max_sd->setJ(min_angMom);
367 >            }
368 >            break;
369 >          case rnemdPx :
370 >            max_vel.x() = min_vel.x();
371 >            max_sd->setVel(max_vel);
372 >            break;
373 >          case rnemdPy :
374 >            max_vel.y() = min_vel.y();
375 >            max_sd->setVel(max_vel);
376 >            break;
377 >          case rnemdPz :
378 >            max_vel.z() = min_vel.z();
379 >            max_sd->setVel(max_vel);
380 >            break;
381 >          case rnemdUnknown :
382 >          default :
383 >            break;
384 >          }
385 >        } else if (min_vals.rank == worldRank) {
386 >          // I had the minimum but not the maximum:
387 >          
388 >          Vector3d max_vel;
389 >          Vector3d min_vel = min_sd->getVel();
390 >          MPI::Status status;
391 >          
392 >          // point-to-point swap of the velocity vector
393 >          MPI::COMM_WORLD.Sendrecv(min_vel.getArrayPointer(), 3, MPI::REALTYPE,
394 >                                   max_vals.rank, 0,
395 >                                   max_vel.getArrayPointer(), 3, MPI::REALTYPE,
396 >                                   max_vals.rank, 0, status);
397 >          
398 >          switch(rnemdType_) {
399 >          case rnemdKinetic :
400 >            min_sd->setVel(max_vel);
401 >            
402 >            if (min_sd->isDirectional()) {
403 >              Vector3d min_angMom = min_sd->getJ();
404 >              Vector3d max_angMom;
405 >
406 >              // point-to-point swap of the angular momentum vector
407 >              MPI::COMM_WORLD.Sendrecv(min_angMom.getArrayPointer(), 3,
408 >                                       MPI::REALTYPE, max_vals.rank, 1,
409 >                                       max_angMom.getArrayPointer(), 3,
410 >                                       MPI::REALTYPE, max_vals.rank, 1,
411 >                                       status);
412 >
413 >              min_sd->setJ(max_angMom);
414 >            }
415 >            break;
416 >          case rnemdPx :
417 >            min_vel.x() = max_vel.x();
418 >            min_sd->setVel(min_vel);
419 >            break;
420 >          case rnemdPy :
421 >            min_vel.y() = max_vel.y();
422 >            min_sd->setVel(min_vel);
423 >            break;
424 >          case rnemdPz :
425 >            min_vel.z() = max_vel.z();
426 >            min_sd->setVel(min_vel);
427 >            break;
428 >          case rnemdUnknown :
429 >          default :
430 >            break;
431 >          }
432 >        }
433 > #endif
434 >        exchangeSum_ += max_val - min_val;
435 >      } else {
436 >        std::cerr << "exchange NOT performed.\nmin_val > max_val.\n";
437 >      }
438 >    } else {
439 >      std::cerr << "exchange NOT performed.\none of the two slabs empty.\n";
440 >    }
441 >    
442 >  }
443 >  
444 >  void RNEMD::getStatus() {
445 >
446 >    Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
447 >    Mat3x3d hmat = currentSnap_->getHmat();
448 >    Stats& stat = currentSnap_->statData;
449 >    RealType time = currentSnap_->getTime();
450 >
451 >    stat[Stats::RNEMD_SWAP_TOTAL] = exchangeSum_;
452 >
453 >    seleMan_.setSelectionSet(evaluator_.evaluate());
454 >
455 >    int selei;
456 >    StuntDouble* sd;
457 >    int idx;
458 >
459 >    std::vector<RealType> valueHist(nBins_, 0.0); // keeps track of what's
460 >                                                  // being averaged
461 >    std::vector<int> valueCount(nBins_, 0);       // keeps track of the
462 >                                                  // number of degrees of
463 >                                                  // freedom being averaged
464 >
465 >    for (sd = seleMan_.beginSelected(selei); sd != NULL;
466 >         sd = seleMan_.nextSelected(selei)) {
467 >      
468 >      idx = sd->getLocalIndex();
469 >      
470 >      Vector3d pos = sd->getPos();
471 >
472 >      // wrap the stuntdouble's position back into the box:
473 >      
474 >      if (usePeriodicBoundaryConditions_)
475 >        currentSnap_->wrapVector(pos);
476 >      
477 >      // which bin is this stuntdouble in?
478 >      // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
479 >      
480 >      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;    
481 >      
482 >      RealType mass = sd->getMass();
483 >      Vector3d vel = sd->getVel();
484 >      RealType value;
485 >
486 >      switch(rnemdType_) {
487 >      case rnemdKinetic :
488 >        
489 >        value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
490 >                        vel[2]*vel[2]);
491 >        
492 >        valueCount[binNo] += 3;
493 >        if (sd->isDirectional()) {
494 >          Vector3d angMom = sd->getJ();
495 >          Mat3x3d I = sd->getI();
496 >          
497 >          if (sd->isLinear()) {
498 >            int i = sd->linearAxis();
499 >            int j = (i + 1) % 3;
500 >            int k = (i + 2) % 3;
501 >            value += angMom[j] * angMom[j] / I(j, j) +
502 >              angMom[k] * angMom[k] / I(k, k);
503 >
504 >            valueCount[binNo] +=2;
505 >
506 >          } else {
507 >            value += angMom[0]*angMom[0]/I(0, 0)
508 >              + angMom[1]*angMom[1]/I(1, 1)
509 >              + angMom[2]*angMom[2]/I(2, 2);
510 >            valueCount[binNo] +=3;
511 >          }
512 >        }
513 >        value = value / OOPSEConstant::energyConvert / OOPSEConstant::kb;
514 >
515 >        break;
516 >      case rnemdPx :
517 >        value = mass * vel[0];
518 >        valueCount[binNo]++;
519 >        break;
520 >      case rnemdPy :
521 >        value = mass * vel[1];
522 >        valueCount[binNo]++;
523 >        break;
524 >      case rnemdPz :
525 >        value = mass * vel[2];
526 >        valueCount[binNo]++;
527 >        break;
528 >      case rnemdUnknown :
529 >      default :
530 >        break;
531 >      }
532 >      valueHist[binNo] += value;
533 >    }
534 >
535 > #ifdef IS_MPI
536 >
537 >    // all processors have the same number of bins, and STL vectors pack their
538 >    // arrays, so in theory, this should be safe:
539 >
540 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueHist[0],
541 >                              nBins_, MPI::REALTYPE, MPI::SUM);
542 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount[0],
543 >                              nBins_, MPI::INT, MPI::SUM);
544 >
545 >    // If we're the root node, should we print out the results
546 >    int worldRank = MPI::COMM_WORLD.Get_rank();
547 >    if (worldRank == 0) {
548 > #endif
549 >      
550 >      std::cout << time;
551 >      for (int j = 0; j < nBins_; j++)
552 >        std::cout << "\t" << valueHist[j] / (RealType)valueCount[j];
553 >      std::cout << "\n";
554 >      
555 > #ifdef IS_MPI
556 >    }
557 > #endif
558 >  }
559   }

Diff Legend

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