ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/rnemd/RNEMD.cpp
(Generate patch)

Comparing trunk/src/rnemd/RNEMD.cpp (file contents):
Revision 2026 by gezelter, Wed Oct 22 12:23:59 2014 UTC vs.
Revision 2068 by gezelter, Thu Mar 5 15:40:58 2015 UTC

# Line 327 | Line 327 | namespace OpenMD {
327            break;
328          }
329        }
330 <      if (hasAngularMomentumFluxVector) {
331 <        std::vector<RealType> amf = rnemdParams->getAngularMomentumFluxVector();
332 <        if (amf.size() != 3) {
333 <          sprintf(painCave.errMsg,
334 <                  "RNEMD: Incorrect number of parameters specified for angularMomentumFluxVector.\n"
335 <                  "\tthere should be 3 parameters, but %lu were specified.\n",
336 <                  amf.size());
337 <          painCave.isFatal = 1;
338 <          simError();      
339 <        }
340 <        angularMomentumFluxVector_.x() = amf[0];
341 <        angularMomentumFluxVector_.y() = amf[1];
342 <        angularMomentumFluxVector_.z() = amf[2];
343 <      } else {
344 <        angularMomentumFluxVector_ = V3Zero;
345 <        if (hasAngularMomentumFlux) {
346 <          RealType angularMomentumFlux = rnemdParams->getAngularMomentumFlux();
347 <          switch (rnemdFluxType_) {
348 <          case rnemdLx:
349 <            angularMomentumFluxVector_.x() = angularMomentumFlux;
350 <            break;
351 <          case rnemdLy:
352 <            angularMomentumFluxVector_.y() = angularMomentumFlux;
353 <            break;
354 <          case rnemdLz:
355 <            angularMomentumFluxVector_.z() = angularMomentumFlux;
356 <            break;
357 <          case rnemdKeLx:
358 <            angularMomentumFluxVector_.x() = angularMomentumFlux;
359 <            break;
360 <          case rnemdKeLy:
361 <            angularMomentumFluxVector_.y() = angularMomentumFlux;
362 <            break;
363 <          case rnemdKeLz:
364 <            angularMomentumFluxVector_.z() = angularMomentumFlux;
365 <            break;
366 <          default:
367 <            break;
368 <          }
369 <        }        
370 <      }
371 <
372 <      if (hasCoordinateOrigin) {
373 <        std::vector<RealType> co = rnemdParams->getCoordinateOrigin();
374 <        if (co.size() != 3) {
375 <          sprintf(painCave.errMsg,
376 <                  "RNEMD: Incorrect number of parameters specified for coordinateOrigin.\n"
377 <                  "\tthere should be 3 parameters, but %lu were specified.\n",
378 <                  co.size());
379 <          painCave.isFatal = 1;
380 <          simError();      
381 <        }
382 <        coordinateOrigin_.x() = co[0];
383 <        coordinateOrigin_.y() = co[1];
384 <        coordinateOrigin_.z() = co[2];
385 <      } else {
386 <        coordinateOrigin_ = V3Zero;
387 <      }
388 <
389 <      // do some sanity checking
390 <
391 <      int selectionCount = seleMan_.getSelectionCount();
392 <
393 <      int nIntegrable = info->getNGlobalIntegrableObjects();
394 <
395 <      if (selectionCount > nIntegrable) {
396 <        sprintf(painCave.errMsg,
397 <                "RNEMD: The current objectSelection,\n"
398 <                "\t\t%s\n"
399 <                "\thas resulted in %d selected objects.  However,\n"
400 <                "\tthe total number of integrable objects in the system\n"
401 <                "\tis only %d.  This is almost certainly not what you want\n"
402 <                "\tto do.  A likely cause of this is forgetting the _RB_0\n"
403 <                "\tselector in the selection script!\n",
404 <                rnemdObjectSelection_.c_str(),
405 <                selectionCount, nIntegrable);
406 <        painCave.isFatal = 0;
407 <        painCave.severity = OPENMD_WARNING;
408 <        simError();
330 >    }
331 >    if (hasAngularMomentumFluxVector) {
332 >      std::vector<RealType> amf = rnemdParams->getAngularMomentumFluxVector();
333 >      if (amf.size() != 3) {
334 >        sprintf(painCave.errMsg,
335 >                "RNEMD: Incorrect number of parameters specified for angularMomentumFluxVector.\n"
336 >                "\tthere should be 3 parameters, but %lu were specified.\n",
337 >                amf.size());
338 >        painCave.isFatal = 1;
339 >        simError();      
340        }
341 <
342 <      areaAccumulator_ = new Accumulator();
343 <
344 <      nBins_ = rnemdParams->getOutputBins();
345 <      binWidth_ = rnemdParams->getOutputBinWidth();
346 <
347 <      data_.resize(RNEMD::ENDINDEX);
417 <      OutputData z;
418 <      z.units =  "Angstroms";
419 <      z.title =  "Z";
420 <      z.dataType = "RealType";
421 <      z.accumulator.reserve(nBins_);
422 <      for (int i = 0; i < nBins_; i++)
423 <        z.accumulator.push_back( new Accumulator() );
424 <      data_[Z] = z;
425 <      outputMap_["Z"] =  Z;
426 <
427 <      OutputData r;
428 <      r.units =  "Angstroms";
429 <      r.title =  "R";
430 <      r.dataType = "RealType";
431 <      r.accumulator.reserve(nBins_);
432 <      for (int i = 0; i < nBins_; i++)
433 <        r.accumulator.push_back( new Accumulator() );
434 <      data_[R] = r;
435 <      outputMap_["R"] =  R;
436 <
437 <      OutputData temperature;
438 <      temperature.units =  "K";
439 <      temperature.title =  "Temperature";
440 <      temperature.dataType = "RealType";
441 <      temperature.accumulator.reserve(nBins_);
442 <      for (int i = 0; i < nBins_; i++)
443 <        temperature.accumulator.push_back( new Accumulator() );
444 <      data_[TEMPERATURE] = temperature;
445 <      outputMap_["TEMPERATURE"] =  TEMPERATURE;
446 <
447 <      OutputData velocity;
448 <      velocity.units = "angstroms/fs";
449 <      velocity.title =  "Velocity";  
450 <      velocity.dataType = "Vector3d";
451 <      velocity.accumulator.reserve(nBins_);
452 <      for (int i = 0; i < nBins_; i++)
453 <        velocity.accumulator.push_back( new VectorAccumulator() );
454 <      data_[VELOCITY] = velocity;
455 <      outputMap_["VELOCITY"] = VELOCITY;
456 <
457 <      OutputData angularVelocity;
458 <      angularVelocity.units = "angstroms^2/fs";
459 <      angularVelocity.title =  "AngularVelocity";  
460 <      angularVelocity.dataType = "Vector3d";
461 <      angularVelocity.accumulator.reserve(nBins_);
462 <      for (int i = 0; i < nBins_; i++)
463 <        angularVelocity.accumulator.push_back( new VectorAccumulator() );
464 <      data_[ANGULARVELOCITY] = angularVelocity;
465 <      outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY;
466 <
467 <      OutputData density;
468 <      density.units =  "g cm^-3";
469 <      density.title =  "Density";
470 <      density.dataType = "RealType";
471 <      density.accumulator.reserve(nBins_);
472 <      for (int i = 0; i < nBins_; i++)
473 <        density.accumulator.push_back( new Accumulator() );
474 <      data_[DENSITY] = density;
475 <      outputMap_["DENSITY"] =  DENSITY;
476 <
477 <      if (hasOutputFields) {
478 <        parseOutputFileFormat(rnemdParams->getOutputFields());
479 <      } else {
480 <        if (usePeriodicBoundaryConditions_)
481 <          outputMask_.set(Z);
482 <        else
483 <          outputMask_.set(R);
341 >      angularMomentumFluxVector_.x() = amf[0];
342 >      angularMomentumFluxVector_.y() = amf[1];
343 >      angularMomentumFluxVector_.z() = amf[2];
344 >    } else {
345 >      angularMomentumFluxVector_ = V3Zero;
346 >      if (hasAngularMomentumFlux) {
347 >        RealType angularMomentumFlux = rnemdParams->getAngularMomentumFlux();
348          switch (rnemdFluxType_) {
485        case rnemdKE:
486        case rnemdRotKE:
487        case rnemdFullKE:
488          outputMask_.set(TEMPERATURE);
489          break;
490        case rnemdPx:
491        case rnemdPy:
492          outputMask_.set(VELOCITY);
493          break;
494        case rnemdPz:        
495        case rnemdPvector:
496          outputMask_.set(VELOCITY);
497          outputMask_.set(DENSITY);
498          break;
349          case rnemdLx:
350 +          angularMomentumFluxVector_.x() = angularMomentumFlux;
351 +          break;
352          case rnemdLy:
353 +          angularMomentumFluxVector_.y() = angularMomentumFlux;
354 +          break;
355          case rnemdLz:
356 <        case rnemdLvector:
503 <          outputMask_.set(ANGULARVELOCITY);
356 >          angularMomentumFluxVector_.z() = angularMomentumFlux;
357            break;
358          case rnemdKeLx:
359 +          angularMomentumFluxVector_.x() = angularMomentumFlux;
360 +          break;
361          case rnemdKeLy:
362 +          angularMomentumFluxVector_.y() = angularMomentumFlux;
363 +          break;
364          case rnemdKeLz:
365 <        case rnemdKeLvector:
509 <          outputMask_.set(TEMPERATURE);
510 <          outputMask_.set(ANGULARVELOCITY);
365 >          angularMomentumFluxVector_.z() = angularMomentumFlux;
366            break;
512        case rnemdKePx:
513        case rnemdKePy:
514          outputMask_.set(TEMPERATURE);
515          outputMask_.set(VELOCITY);
516          break;
517        case rnemdKePvector:
518          outputMask_.set(TEMPERATURE);
519          outputMask_.set(VELOCITY);
520          outputMask_.set(DENSITY);        
521          break;
367          default:
368            break;
369          }
370 <      }
371 <      
527 <      if (hasOutputFileName) {
528 <        rnemdFileName_ = rnemdParams->getOutputFileName();
529 <      } else {
530 <        rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd";
531 <      }          
370 >      }        
371 >    }
372  
373 <      exchangeTime_ = rnemdParams->getExchangeTime();
374 <
375 <      Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot();
376 <      // total exchange sums are zeroed out at the beginning:
377 <
378 <      kineticExchange_ = 0.0;
379 <      momentumExchange_ = V3Zero;
380 <      angularMomentumExchange_ = V3Zero;
381 <
382 <      std::ostringstream selectionAstream;
383 <      std::ostringstream selectionBstream;
373 >    if (hasCoordinateOrigin) {
374 >      std::vector<RealType> co = rnemdParams->getCoordinateOrigin();
375 >      if (co.size() != 3) {
376 >        sprintf(painCave.errMsg,
377 >                "RNEMD: Incorrect number of parameters specified for coordinateOrigin.\n"
378 >                "\tthere should be 3 parameters, but %lu were specified.\n",
379 >                co.size());
380 >        painCave.isFatal = 1;
381 >        simError();      
382 >      }
383 >      coordinateOrigin_.x() = co[0];
384 >      coordinateOrigin_.y() = co[1];
385 >      coordinateOrigin_.z() = co[2];
386 >    } else {
387 >      coordinateOrigin_ = V3Zero;
388 >    }
389      
390 <      if (hasSelectionA_) {
391 <        selectionA_ = rnemdParams->getSelectionA();
392 <      } else {
393 <        if (usePeriodicBoundaryConditions_) {    
394 <          Mat3x3d hmat = currentSnap_->getHmat();
395 <        
396 <          if (hasSlabWidth)
397 <            slabWidth_ = rnemdParams->getSlabWidth();
398 <          else
399 <            slabWidth_ = hmat(2,2) / 10.0;
400 <        
401 <          if (hasSlabACenter)
402 <            slabACenter_ = rnemdParams->getSlabACenter();
403 <          else
404 <            slabACenter_ = 0.0;
405 <        
406 <          selectionAstream << "select wrappedz > "
407 <                           << slabACenter_ - 0.5*slabWidth_
408 <                           <<  " && wrappedz < "
409 <                           << slabACenter_ + 0.5*slabWidth_;
410 <          selectionA_ = selectionAstream.str();
411 <        } else {
412 <          if (hasSphereARadius)
413 <            sphereARadius_ = rnemdParams->getSphereARadius();
414 <          else {
415 <            // use an initial guess to the size of the inner slab to be 1/10 the
416 <            // radius of an approximately spherical hull:
417 <            Thermo thermo(info);
418 <            RealType hVol = thermo.getHullVolume();
419 <            sphereARadius_ = 0.1 * pow((3.0 * hVol / (4.0 * M_PI)), 1.0/3.0);
420 <          }
421 <          selectionAstream << "select r < " << sphereARadius_;
422 <          selectionA_ = selectionAstream.str();
423 <        }
390 >    // do some sanity checking
391 >    
392 >    int selectionCount = seleMan_.getSelectionCount();    
393 >    int nIntegrable = info->getNGlobalIntegrableObjects();
394 >    if (selectionCount > nIntegrable) {
395 >      sprintf(painCave.errMsg,
396 >              "RNEMD: The current objectSelection,\n"
397 >              "\t\t%s\n"
398 >              "\thas resulted in %d selected objects.  However,\n"
399 >              "\tthe total number of integrable objects in the system\n"
400 >              "\tis only %d.  This is almost certainly not what you want\n"
401 >              "\tto do.  A likely cause of this is forgetting the _RB_0\n"
402 >              "\tselector in the selection script!\n",
403 >              rnemdObjectSelection_.c_str(),
404 >              selectionCount, nIntegrable);
405 >      painCave.isFatal = 0;
406 >      painCave.severity = OPENMD_WARNING;
407 >      simError();
408 >    }
409 >    
410 >    areaAccumulator_ = new Accumulator();
411 >    
412 >    nBins_ = rnemdParams->getOutputBins();
413 >    binWidth_ = rnemdParams->getOutputBinWidth();
414 >    
415 >    data_.resize(RNEMD::ENDINDEX);
416 >    OutputData z;
417 >    z.units =  "Angstroms";
418 >    z.title =  "Z";
419 >    z.dataType = "RealType";
420 >    z.accumulator.reserve(nBins_);
421 >    for (int i = 0; i < nBins_; i++)
422 >      z.accumulator.push_back( new Accumulator() );
423 >    data_[Z] = z;
424 >    outputMap_["Z"] =  Z;
425 >    
426 >    OutputData r;
427 >    r.units =  "Angstroms";
428 >    r.title =  "R";
429 >    r.dataType = "RealType";
430 >    r.accumulator.reserve(nBins_);
431 >    for (int i = 0; i < nBins_; i++)
432 >      r.accumulator.push_back( new Accumulator() );
433 >    data_[R] = r;
434 >    outputMap_["R"] =  R;
435 >    
436 >    OutputData temperature;
437 >    temperature.units =  "K";
438 >    temperature.title =  "Temperature";
439 >    temperature.dataType = "RealType";
440 >    temperature.accumulator.reserve(nBins_);
441 >    for (int i = 0; i < nBins_; i++)
442 >      temperature.accumulator.push_back( new Accumulator() );
443 >    data_[TEMPERATURE] = temperature;
444 >    outputMap_["TEMPERATURE"] =  TEMPERATURE;
445 >    
446 >    OutputData velocity;
447 >    velocity.units = "angstroms/fs";
448 >    velocity.title =  "Velocity";  
449 >    velocity.dataType = "Vector3d";
450 >    velocity.accumulator.reserve(nBins_);
451 >    for (int i = 0; i < nBins_; i++)
452 >      velocity.accumulator.push_back( new VectorAccumulator() );
453 >    data_[VELOCITY] = velocity;
454 >    outputMap_["VELOCITY"] = VELOCITY;
455 >    
456 >    OutputData angularVelocity;
457 >    angularVelocity.units = "angstroms^2/fs";
458 >    angularVelocity.title =  "AngularVelocity";  
459 >    angularVelocity.dataType = "Vector3d";
460 >    angularVelocity.accumulator.reserve(nBins_);
461 >    for (int i = 0; i < nBins_; i++)
462 >      angularVelocity.accumulator.push_back( new VectorAccumulator() );
463 >    data_[ANGULARVELOCITY] = angularVelocity;
464 >    outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY;
465 >    
466 >    OutputData density;
467 >    density.units =  "g cm^-3";
468 >    density.title =  "Density";
469 >    density.dataType = "RealType";
470 >    density.accumulator.reserve(nBins_);
471 >    for (int i = 0; i < nBins_; i++)
472 >      density.accumulator.push_back( new Accumulator() );
473 >    data_[DENSITY] = density;
474 >    outputMap_["DENSITY"] =  DENSITY;
475 >    
476 >    if (hasOutputFields) {
477 >      parseOutputFileFormat(rnemdParams->getOutputFields());
478 >    } else {
479 >      if (usePeriodicBoundaryConditions_)
480 >        outputMask_.set(Z);
481 >      else
482 >        outputMask_.set(R);
483 >      switch (rnemdFluxType_) {
484 >      case rnemdKE:
485 >      case rnemdRotKE:
486 >      case rnemdFullKE:
487 >        outputMask_.set(TEMPERATURE);
488 >        break;
489 >      case rnemdPx:
490 >      case rnemdPy:
491 >        outputMask_.set(VELOCITY);
492 >        break;
493 >      case rnemdPz:        
494 >      case rnemdPvector:
495 >        outputMask_.set(VELOCITY);
496 >        outputMask_.set(DENSITY);
497 >        break;
498 >      case rnemdLx:
499 >      case rnemdLy:
500 >      case rnemdLz:
501 >      case rnemdLvector:
502 >        outputMask_.set(ANGULARVELOCITY);
503 >        break;
504 >      case rnemdKeLx:
505 >      case rnemdKeLy:
506 >      case rnemdKeLz:
507 >      case rnemdKeLvector:
508 >        outputMask_.set(TEMPERATURE);
509 >        outputMask_.set(ANGULARVELOCITY);
510 >        break;
511 >      case rnemdKePx:
512 >      case rnemdKePy:
513 >        outputMask_.set(TEMPERATURE);
514 >        outputMask_.set(VELOCITY);
515 >        break;
516 >      case rnemdKePvector:
517 >        outputMask_.set(TEMPERATURE);
518 >        outputMask_.set(VELOCITY);
519 >        outputMask_.set(DENSITY);        
520 >        break;
521 >      default:
522 >        break;
523        }
524 +    }
525      
526 <      if (hasSelectionB_) {
527 <        selectionB_ = rnemdParams->getSelectionB();
528 <
529 <      } else {
530 <        if (usePeriodicBoundaryConditions_) {    
531 <          Mat3x3d hmat = currentSnap_->getHmat();
526 >    if (hasOutputFileName) {
527 >      rnemdFileName_ = rnemdParams->getOutputFileName();
528 >    } else {
529 >      rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd";
530 >    }          
531 >    
532 >    exchangeTime_ = rnemdParams->getExchangeTime();
533 >    
534 >    Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot();
535 >    // total exchange sums are zeroed out at the beginning:
536 >    
537 >    kineticExchange_ = 0.0;
538 >    momentumExchange_ = V3Zero;
539 >    angularMomentumExchange_ = V3Zero;
540 >    
541 >    std::ostringstream selectionAstream;
542 >    std::ostringstream selectionBstream;
543 >    
544 >    if (hasSelectionA_) {
545 >      selectionA_ = rnemdParams->getSelectionA();
546 >    } else {
547 >      if (usePeriodicBoundaryConditions_) {    
548 >        Mat3x3d hmat = currentSnap_->getHmat();
549          
550 <          if (hasSlabWidth)
551 <            slabWidth_ = rnemdParams->getSlabWidth();
552 <          else
553 <            slabWidth_ = hmat(2,2) / 10.0;
550 >        if (hasSlabWidth)
551 >          slabWidth_ = rnemdParams->getSlabWidth();
552 >        else
553 >          slabWidth_ = hmat(2,2) / 10.0;
554          
555 <          if (hasSlabBCenter)
556 <            slabBCenter_ = rnemdParams->getSlabBCenter();
557 <          else
558 <            slabBCenter_ = hmat(2,2) / 2.0;
555 >        if (hasSlabACenter)
556 >          slabACenter_ = rnemdParams->getSlabACenter();
557 >        else
558 >          slabACenter_ = 0.0;
559          
560 <          selectionBstream << "select wrappedz > "
561 <                           << slabBCenter_ - 0.5*slabWidth_
562 <                           <<  " && wrappedz < "
563 <                           << slabBCenter_ + 0.5*slabWidth_;
560 >        selectionAstream << "select wrappedz > "
561 >                         << slabACenter_ - 0.5*slabWidth_
562 >                         <<  " && wrappedz < "
563 >                         << slabACenter_ + 0.5*slabWidth_;
564 >        selectionA_ = selectionAstream.str();
565 >      } else {
566 >        if (hasSphereARadius)
567 >          sphereARadius_ = rnemdParams->getSphereARadius();
568 >        else {
569 >          // use an initial guess to the size of the inner slab to be 1/10 the
570 >          // radius of an approximately spherical hull:
571 >          Thermo thermo(info);
572 >          RealType hVol = thermo.getHullVolume();
573 >          sphereARadius_ = 0.1 * pow((3.0 * hVol / (4.0 * M_PI)), 1.0/3.0);
574 >        }
575 >        selectionAstream << "select r < " << sphereARadius_;
576 >        selectionA_ = selectionAstream.str();
577 >      }
578 >    }
579 >    
580 >    if (hasSelectionB_) {
581 >      selectionB_ = rnemdParams->getSelectionB();
582 >      
583 >    } else {
584 >      if (usePeriodicBoundaryConditions_) {    
585 >        Mat3x3d hmat = currentSnap_->getHmat();
586 >        
587 >        if (hasSlabWidth)
588 >          slabWidth_ = rnemdParams->getSlabWidth();
589 >        else
590 >          slabWidth_ = hmat(2,2) / 10.0;
591 >        
592 >        if (hasSlabBCenter)
593 >          slabBCenter_ = rnemdParams->getSlabBCenter();
594 >        else
595 >          slabBCenter_ = hmat(2,2) / 2.0;
596 >        
597 >        selectionBstream << "select wrappedz > "
598 >                         << slabBCenter_ - 0.5*slabWidth_
599 >                         <<  " && wrappedz < "
600 >                         << slabBCenter_ + 0.5*slabWidth_;
601 >        selectionB_ = selectionBstream.str();
602 >      } else {
603 >        if (hasSphereBRadius_) {
604 >          sphereBRadius_ = rnemdParams->getSphereBRadius();
605 >          selectionBstream << "select r > " << sphereBRadius_;
606            selectionB_ = selectionBstream.str();
607          } else {
608 <          if (hasSphereBRadius_) {
609 <            sphereBRadius_ = rnemdParams->getSphereBRadius();
610 <            selectionBstream << "select r > " << sphereBRadius_;
607 <            selectionB_ = selectionBstream.str();
608 <          } else {
609 <            selectionB_ = "select hull";
610 <            BisHull_ = true;
611 <            hasSelectionB_ = true;
612 <          }
608 >          selectionB_ = "select hull";
609 >          BisHull_ = true;
610 >          hasSelectionB_ = true;
611          }
612        }
613      }
614 <
614 >  
615 >  
616      // object evaluator:
617      evaluator_.loadScriptString(rnemdObjectSelection_);
618      seleMan_.setSelectionSet(evaluator_.evaluate());
# Line 655 | Line 654 | namespace OpenMD {
654  
655      StuntDouble* sd;
656  
657 <    RealType min_val;
657 >    RealType min_val(0.0);
658      int min_found = 0;  
659 <    StuntDouble* min_sd;
659 >    StuntDouble* min_sd = NULL;
660  
661 <    RealType max_val;
661 >    RealType max_val(0.0);
662      int max_found = 0;
663 <    StuntDouble* max_sd;
663 >    StuntDouble* max_sd = NULL;
664  
665      for (sd = seleManA_.beginSelected(selei); sd != NULL;
666           sd = seleManA_.nextSelected(selei)) {
# Line 1021 | Line 1020 | namespace OpenMD {
1020      int selej;
1021  
1022      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1024    RealType time = currentSnap_->getTime();    
1023      Mat3x3d hmat = currentSnap_->getHmat();
1024  
1025      StuntDouble* sd;
# Line 1469 | Line 1467 | namespace OpenMD {
1467      int selej;
1468  
1469      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1472    RealType time = currentSnap_->getTime();    
1470      Mat3x3d hmat = currentSnap_->getHmat();
1471  
1472      StuntDouble* sd;
# Line 1776 | Line 1773 | namespace OpenMD {
1773          areaA = evaluatorA_.getSurfaceArea();
1774        else {
1775          
1779        cerr << "selection A did not have surface area, recomputing\n";
1776          int isd;
1777          StuntDouble* sd;
1778          vector<StuntDouble*> aSites;
# Line 1788 | Line 1784 | namespace OpenMD {
1784   #if defined(HAVE_QHULL)
1785          ConvexHull* surfaceMeshA = new ConvexHull();
1786          surfaceMeshA->computeHull(aSites);
1791        cerr << "flag1\n";
1787          areaA = surfaceMeshA->getArea();
1793        cerr << "Flag2 " << areaA << "\n";
1788          delete surfaceMeshA;
1789   #else
1790          sprintf( painCave.errMsg,
# Line 1816 | Line 1810 | namespace OpenMD {
1810      }
1811  
1812      if (hasSelectionB_) {
1813 <      if (evaluatorB_.hasSurfaceArea())
1813 >      if (evaluatorB_.hasSurfaceArea()) {
1814          areaB = evaluatorB_.getSurfaceArea();
1815 <      else {
1822 <        cerr << "selection B did not have surface area, recomputing\n";
1815 >      } else {
1816  
1817          int isd;
1818          StuntDouble* sd;
# Line 1910 | Line 1903 | namespace OpenMD {
1903    void RNEMD::collectData() {
1904      if (!doRNEMD_) return;
1905      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1906 <    
1906 >
1907      // collectData can be called more frequently than the doRNEMD, so use the
1908      // computed area from the last exchange time:
1909      RealType area = getDividingArea();
# Line 2025 | Line 2018 | namespace OpenMD {
2018          }
2019        }
2020      }
2021 <    
2021 >
2022   #ifdef IS_MPI
2023  
2024      for (int i = 0; i < nBins_; i++) {
# Line 2132 | Line 2125 | namespace OpenMD {
2125          painCave.severity = OPENMD_ERROR;
2126          simError();            
2127        }
2128 <    }  
2128 >    }
2129    }
2130    
2131    void RNEMD::writeOutputFile() {

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines