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 1341 by skuang, Fri May 8 19:47:05 2009 UTC vs.
Revision 1350 by gezelter, Thu May 21 18:56:45 2009 UTC

# Line 53 | Line 53
53   #include "math/ParallelRandNumGen.hpp"
54   #endif
55  
56 + #define HONKING_LARGE_VALUE 1.0e10
57  
58   namespace oopse {
59    
# Line 102 | Line 103 | namespace oopse {
103      set_RNEMD_swapTime(simParams->getRNEMD_swapTime());
104      set_RNEMD_nBins(simParams->getRNEMD_nBins());
105      exchangeSum_ = 0.0;
105 –    counter_ = 0; //added by shenyu
106  
107   #ifndef IS_MPI
108      if (simParams->haveSeed()) {
# Line 233 | Line 233 | namespace oopse {
233        }
234      }
235  
236 <    // missing:  swap information in parallel
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 <        Vector3d min_vel = min_sd->getVel();
291 <        Vector3d max_vel = max_sd->getVel();
292 <        RealType temp_vel;
293 <
294 <        switch(rnemdType_) {
295 <        case rnemdKinetic :
296 <          min_sd->setVel(max_vel);
297 <          max_sd->setVel(min_vel);
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 <          if (min_sd->isDirectional() && max_sd->isDirectional()) {
346 <            Vector3d min_angMom = min_sd->getJ();
347 <            Vector3d max_angMom = max_sd->getJ();
348 <            min_sd->setJ(max_angMom);
349 <            max_sd->setJ(min_angMom);
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 <          break;
386 <        case rnemdPx :
387 <          temp_vel = min_vel.x();
388 <          min_vel.x() = max_vel.x();
389 <          max_vel.x() = temp_vel;
390 <          min_sd->setVel(min_vel);
391 <          max_sd->setVel(max_vel);
392 <          break;
393 <        case rnemdPy :
394 <          temp_vel = min_vel.y();
395 <          min_vel.y() = max_vel.y();
396 <          max_vel.y() = temp_vel;
397 <          min_sd->setVel(min_vel);
398 <          max_sd->setVel(max_vel);
399 <          break;
400 <        case rnemdPz :
401 <          temp_vel = min_vel.z();
402 <          min_vel.z() = max_vel.z();
403 <          max_vel.z() = temp_vel;
404 <          min_sd->setVel(min_vel);
405 <          max_sd->setVel(max_vel);
406 <          break;
407 <        case rnemdUnknown :
408 <        default :
409 <          break;
410 <        }
411 <      exchangeSum_ += max_val - min_val;
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";
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 <
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  
# Line 302 | Line 456 | namespace oopse {
456      StuntDouble* sd;
457      int idx;
458  
459 <    std::vector<RealType> valueHist;  // keeps track of what's being averaged
460 <    std::vector<int> valueCount; // keeps track of the number of degrees of
461 <                                 // freedom being averaged
462 <    valueHist.resize(nBins_);
463 <    valueCount.resize(nBins_);
464 <    //do they initialize themselves to zero automatically?
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        
# Line 323 | Line 477 | namespace oopse {
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_;
480 >      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;    
481        
328 –      //std::cerr << "pos.z() = " << pos.z() << " bin = " << binNo << "\n";
329 –      
482        RealType mass = sd->getMass();
483        Vector3d vel = sd->getVel();
332 –      //std::cerr << "mass = " << mass << " vel = " << vel << "\n";
484        RealType value;
485  
486        switch(rnemdType_) {
# Line 339 | Line 490 | namespace oopse {
490                          vel[2]*vel[2]);
491          
492          valueCount[binNo] += 3;
342 –        //std::cerr <<"starting value = " << value << "\n";
493          if (sd->isDirectional()) {
344 –          //std::cerr << "angMom calculated.\n";
494            Vector3d angMom = sd->getJ();
346 –          //std::cerr << "current angMom: " << angMom << "\n";
495            Mat3x3d I = sd->getI();
496            
497            if (sd->isLinear()) {
# Line 356 | Line 504 | namespace oopse {
504              valueCount[binNo] +=2;
505  
506            } else {
359 –            //std::cerr << "non-linear molecule.\n";
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;
364 –
511            }
512          }
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
513          value = value / OOPSEConstant::energyConvert / OOPSEConstant::kb;
514 <        //std::cerr <<"value = " << value << "\n";
514 >
515          break;
516        case rnemdPx :
517          value = mass * vel[0];
# Line 386 | Line 529 | namespace oopse {
529        default :
530          break;
531        }
389 –      //std::cerr << "bin = " << binNo << " value = " << value ;
532        valueHist[binNo] += value;
391 –      //std::cerr << " hist = " << valueHist[binNo] << " count = " << valueCount[binNo] << "\n";
533      }
534 <    
535 <    std::cout << counter_++;
536 <    for (int j = 0; j < nBins_; j++)
537 <      std::cout << "\t" << valueHist[j] / (RealType)valueCount[j];
538 <    std::cout << "\n";
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