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 1350 by gezelter, Thu May 21 18:56:45 2009 UTC vs.
branches/development/src/integrators/RNEMD.cpp (file contents), Revision 1629 by gezelter, Wed Sep 14 21:15:17 2011 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 + * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41  
42 + #include <cmath>
43   #include "integrators/RNEMD.hpp"
44   #include "math/Vector3.hpp"
45   #include "math/SquareMatrix3.hpp"
46 + #include "math/Polynomial.hpp"
47   #include "primitives/Molecule.hpp"
48   #include "primitives/StuntDouble.hpp"
49 < #include "utils/OOPSEConstant.hpp"
49 > #include "utils/PhysicalConstants.hpp"
50   #include "utils/Tuple.hpp"
51  
52   #ifndef IS_MPI
53   #include "math/SeqRandNumGen.hpp"
54   #else
55 + #include <mpi.h>
56   #include "math/ParallelRandNumGen.hpp"
57   #endif
58  
59   #define HONKING_LARGE_VALUE 1.0e10
60  
61 < namespace oopse {
61 > using namespace std;
62 > namespace OpenMD {
63    
64 <  RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
65 <    
64 >  RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info),
65 >                                usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
66 >
67 >    failTrialCount_ = 0;
68 >    failRootCount_ = 0;
69 >
70      int seedValue;
71      Globals * simParams = info->getSimParams();
72  
73 <    stringToEnumMap_["Kinetic"] = rnemdKinetic;
73 >    stringToEnumMap_["KineticSwap"] = rnemdKineticSwap;
74 >    stringToEnumMap_["KineticScale"] = rnemdKineticScale;
75 >    stringToEnumMap_["PxScale"] = rnemdPxScale;
76 >    stringToEnumMap_["PyScale"] = rnemdPyScale;
77 >    stringToEnumMap_["PzScale"] = rnemdPzScale;
78      stringToEnumMap_["Px"] = rnemdPx;
79      stringToEnumMap_["Py"] = rnemdPy;
80      stringToEnumMap_["Pz"] = rnemdPz;
# Line 72 | Line 84 | namespace oopse {
84      evaluator_.loadScriptString(rnemdObjectSelection_);
85      seleMan_.setSelectionSet(evaluator_.evaluate());
86  
75
87      // do some sanity checking
88  
89      int selectionCount = seleMan_.getSelectionCount();
# Line 80 | Line 91 | namespace oopse {
91  
92      if (selectionCount > nIntegrable) {
93        sprintf(painCave.errMsg,
94 <              "RNEMD warning: The current RNEMD_objectSelection,\n"
94 >              "RNEMD: The current RNEMD_objectSelection,\n"
95                "\t\t%s\n"
96                "\thas resulted in %d selected objects.  However,\n"
97                "\tthe total number of integrable objects in the system\n"
# Line 90 | Line 101 | namespace oopse {
101                rnemdObjectSelection_.c_str(),
102                selectionCount, nIntegrable);
103        painCave.isFatal = 0;
104 +      painCave.severity = OPENMD_WARNING;
105        simError();
94
106      }
107      
108 <    const std::string st = simParams->getRNEMD_swapType();
108 >    const string st = simParams->getRNEMD_exchangeType();
109  
110 <    std::map<std::string, RNEMDTypeEnum>::iterator i;
110 >    map<string, RNEMDTypeEnum>::iterator i;
111      i = stringToEnumMap_.find(st);
112 <    rnemdType_  = (i == stringToEnumMap_.end()) ? RNEMD::rnemdUnknown : i->second;
112 >    rnemdType_ = (i == stringToEnumMap_.end()) ? RNEMD::rnemdUnknown : i->second;
113 >    if (rnemdType_ == rnemdUnknown) {
114 >      sprintf(painCave.errMsg,
115 >              "RNEMD: The current RNEMD_exchangeType,\n"
116 >              "\t\t%s\n"
117 >              "\tis not one of the recognized exchange types.\n",
118 >              st.c_str());
119 >      painCave.isFatal = 1;
120 >      painCave.severity = OPENMD_ERROR;
121 >      simError();
122 >    }
123 >    
124 >    output3DTemp_ = false;
125 >    if (simParams->haveRNEMD_outputDimensionalTemperature()) {
126 >      output3DTemp_ = simParams->getRNEMD_outputDimensionalTemperature();
127 >    }
128  
129 <    set_RNEMD_swapTime(simParams->getRNEMD_swapTime());
129 > #ifdef IS_MPI
130 >    if (worldRank == 0) {
131 > #endif
132 >
133 >      string rnemdFileName;
134 >      switch(rnemdType_) {
135 >      case rnemdKineticSwap :
136 >      case rnemdKineticScale :
137 >        rnemdFileName = "temperature.log";
138 >        break;
139 >      case rnemdPx :
140 >      case rnemdPxScale :
141 >      case rnemdPy :
142 >      case rnemdPyScale :
143 >        rnemdFileName = "momemtum.log";
144 >        break;
145 >      case rnemdPz :
146 >      case rnemdPzScale :
147 >      case rnemdUnknown :
148 >      default :
149 >        rnemdFileName = "rnemd.log";
150 >        break;
151 >      }
152 >      rnemdLog_.open(rnemdFileName.c_str());
153 >
154 >      string xTempFileName;
155 >      string yTempFileName;
156 >      string zTempFileName;
157 >      if (output3DTemp_) {
158 >        xTempFileName = "temperatureX.log";
159 >        yTempFileName = "temperatureY.log";
160 >        zTempFileName = "temperatureZ.log";
161 >        xTempLog_.open(xTempFileName.c_str());
162 >        yTempLog_.open(yTempFileName.c_str());
163 >        zTempLog_.open(zTempFileName.c_str());
164 >      }
165 >
166 > #ifdef IS_MPI
167 >    }
168 > #endif
169 >
170 >    set_RNEMD_exchange_time(simParams->getRNEMD_exchangeTime());
171      set_RNEMD_nBins(simParams->getRNEMD_nBins());
172 <    exchangeSum_ = 0.0;
172 >    midBin_ = nBins_ / 2;
173 >    if (simParams->haveRNEMD_binShift()) {
174 >      if (simParams->getRNEMD_binShift()) {
175 >        zShift_ = 0.5 / (RealType)(nBins_);
176 >      } else {
177 >        zShift_ = 0.0;
178 >      }
179 >    } else {
180 >      zShift_ = 0.0;
181 >    }
182 >    //cerr << "we have zShift_ = " << zShift_ << "\n";
183 >    //shift slabs by half slab width, might be useful in heterogeneous systems
184 >    //set to 0.0 if not using it; can NOT be used in status output yet
185 >    if (simParams->haveRNEMD_logWidth()) {
186 >      set_RNEMD_logWidth(simParams->getRNEMD_logWidth());
187 >      /*arbitary rnemdLogWidth_ no checking
188 >        if (rnemdLogWidth_ != nBins_ && rnemdLogWidth_ != midBin_ + 1) {
189 >        cerr << "WARNING! RNEMD_logWidth has abnormal value!\n";
190 >        cerr << "Automaically set back to default.\n";
191 >        rnemdLogWidth_ = nBins_;
192 >        }*/
193 >    } else {
194 >      set_RNEMD_logWidth(nBins_);
195 >    }
196 >    valueHist_.resize(rnemdLogWidth_, 0.0);
197 >    valueCount_.resize(rnemdLogWidth_, 0);
198 >    xTempHist_.resize(rnemdLogWidth_, 0.0);
199 >    yTempHist_.resize(rnemdLogWidth_, 0.0);
200 >    zTempHist_.resize(rnemdLogWidth_, 0.0);
201 >    xyzTempCount_.resize(rnemdLogWidth_, 0);
202  
203 +    set_RNEMD_exchange_total(0.0);
204 +    if (simParams->haveRNEMD_targetFlux()) {
205 +      set_RNEMD_target_flux(simParams->getRNEMD_targetFlux());
206 +    } else {
207 +      set_RNEMD_target_flux(0.0);
208 +    }
209 +
210   #ifndef IS_MPI
211      if (simParams->haveSeed()) {
212        seedValue = simParams->getSeed();
# Line 123 | Line 226 | namespace oopse {
226    
227    RNEMD::~RNEMD() {
228      delete randNumGen_;
229 +    
230 + #ifdef IS_MPI
231 +    if (worldRank == 0) {
232 + #endif
233 +      
234 +      sprintf(painCave.errMsg,
235 +              "RNEMD: total failed trials: %d\n",
236 +              failTrialCount_);
237 +      painCave.isFatal = 0;
238 +      painCave.severity = OPENMD_INFO;
239 +      simError();
240 +
241 +      rnemdLog_.close();
242 +      if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPyScale) {
243 +        sprintf(painCave.errMsg,
244 +                "RNEMD: total root-checking warnings: %d\n",
245 +                failRootCount_);
246 +        painCave.isFatal = 0;
247 +        painCave.severity = OPENMD_INFO;
248 +        simError();
249 +      }
250 +      if (output3DTemp_) {
251 +        xTempLog_.close();
252 +        yTempLog_.close();
253 +        zTempLog_.close();
254 +      }
255 + #ifdef IS_MPI
256 +    }
257 + #endif
258    }
259  
260    void RNEMD::doSwap() {
129    int midBin = nBins_ / 2;
261  
262      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
263      Mat3x3d hmat = currentSnap_->getHmat();
# Line 160 | Line 291 | namespace oopse {
291        // which bin is this stuntdouble in?
292        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
293  
294 <      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
294 >      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_;
295  
296  
297        // if we're in bin 0 or the middleBin
298 <      if (binNo == 0 || binNo == midBin) {
298 >      if (binNo == 0 || binNo == midBin_) {
299          
300          RealType mass = sd->getMass();
301          Vector3d vel = sd->getVel();
302          RealType value;
303  
304          switch(rnemdType_) {
305 <        case rnemdKinetic :
305 >        case rnemdKineticSwap :
306            
307            value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
308                            vel[2]*vel[2]);
309 <          if (sd->isDirectional()) {
309 >          /*
310 >            if (sd->isDirectional()) {
311              Vector3d angMom = sd->getJ();
312              Mat3x3d I = sd->getI();
313              
314              if (sd->isLinear()) {
315 <              int i = sd->linearAxis();
316 <              int j = (i + 1) % 3;
317 <              int k = (i + 2) % 3;
318 <              value += angMom[j] * angMom[j] / I(j, j) +
319 <                angMom[k] * angMom[k] / I(k, k);
315 >            int i = sd->linearAxis();
316 >            int j = (i + 1) % 3;
317 >            int k = (i + 2) % 3;
318 >            value += angMom[j] * angMom[j] / I(j, j) +
319 >            angMom[k] * angMom[k] / I(k, k);
320              } else {                        
321 <              value += angMom[0]*angMom[0]/I(0, 0)
322 <                + angMom[1]*angMom[1]/I(1, 1)
323 <                + angMom[2]*angMom[2]/I(2, 2);
321 >            value += angMom[0]*angMom[0]/I(0, 0)
322 >            + angMom[1]*angMom[1]/I(1, 1)
323 >            + angMom[2]*angMom[2]/I(2, 2);
324              }
325 <          }
326 <          value = value * 0.5 / OOPSEConstant::energyConvert;
325 >            } no exchange of angular momenta
326 >          */
327 >          //make exchangeSum_ comparable between swap & scale
328 >          //temporarily without using energyConvert
329 >          //value = value * 0.5 / PhysicalConstants::energyConvert;
330 >          value *= 0.5;
331            break;
332          case rnemdPx :
333            value = mass * vel[0];
# Line 202 | Line 338 | namespace oopse {
338          case rnemdPz :
339            value = mass * vel[2];
340            break;
205        case rnemdUnknown :
341          default :
342            break;
343          }
# Line 218 | Line 353 | namespace oopse {
353                min_sd = sd;
354              }
355            }
356 <        } else {
356 >        } else { //midBin_
357            if (!max_found) {
358              max_val = value;
359              max_sd = sd;
# Line 243 | Line 378 | namespace oopse {
378      bool my_max_found = max_found;
379  
380      // Even if we didn't find a minimum, did someone else?
381 <    MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found,
247 <                              1, MPI::BOOL, MPI::LAND);
248 <    
381 >    MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, 1, MPI::BOOL, MPI::LOR);
382      // Even if we didn't find a maximum, did someone else?
383 <    MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found,
251 <                              1, MPI::BOOL, MPI::LAND);
252 <    
383 >    MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, 1, MPI::BOOL, MPI::LOR);
384      struct {
385        RealType val;
386        int rank;
# Line 285 | Line 416 | namespace oopse {
416   #endif
417  
418      if (max_found && min_found) {
419 <      if (min_val< max_val) {
419 >      if (min_val < max_val) {
420  
421   #ifdef IS_MPI      
422          if (max_vals.rank == worldRank && min_vals.rank == worldRank) {
423            // I have both maximum and minimum, so proceed like a single
424            // processor version:
425   #endif
426 <          // objects to be swapped: velocity & angular velocity
426 >          // objects to be swapped: velocity ONLY
427            Vector3d min_vel = min_sd->getVel();
428            Vector3d max_vel = max_sd->getVel();
429            RealType temp_vel;
430            
431            switch(rnemdType_) {
432 <          case rnemdKinetic :
432 >          case rnemdKineticSwap :
433              min_sd->setVel(max_vel);
434              max_sd->setVel(min_vel);
435 <            if (min_sd->isDirectional() && max_sd->isDirectional()) {
435 >            /*
436 >              if (min_sd->isDirectional() && max_sd->isDirectional()) {
437                Vector3d min_angMom = min_sd->getJ();
438                Vector3d max_angMom = max_sd->getJ();
439                min_sd->setJ(max_angMom);
440                max_sd->setJ(min_angMom);
441 <            }
441 >              } no angular momentum exchange
442 >            */
443              break;
444            case rnemdPx :
445              temp_vel = min_vel.x();
# Line 329 | Line 462 | namespace oopse {
462              min_sd->setVel(min_vel);
463              max_sd->setVel(max_vel);
464              break;
332          case rnemdUnknown :
465            default :
466              break;
467            }
# Line 349 | Line 481 | namespace oopse {
481                                     min_vals.rank, 0, status);
482            
483            switch(rnemdType_) {
484 <          case rnemdKinetic :
484 >          case rnemdKineticSwap :
485              max_sd->setVel(min_vel);
486 <            
486 >            //no angular momentum exchange for now
487 >            /*
488              if (max_sd->isDirectional()) {
489                Vector3d min_angMom;
490                Vector3d max_angMom = max_sd->getJ();
491 <
491 >              
492                // point-to-point swap of the angular momentum vector
493                MPI::COMM_WORLD.Sendrecv(max_angMom.getArrayPointer(), 3,
494                                         MPI::REALTYPE, min_vals.rank, 1,
495                                         min_angMom.getArrayPointer(), 3,
496                                         MPI::REALTYPE, min_vals.rank, 1,
497                                         status);
498 <
498 >              
499                max_sd->setJ(min_angMom);
500 <            }
500 >             }
501 >             */            
502              break;
503            case rnemdPx :
504              max_vel.x() = min_vel.x();
# Line 378 | Line 512 | namespace oopse {
512              max_vel.z() = min_vel.z();
513              max_sd->setVel(max_vel);
514              break;
381          case rnemdUnknown :
515            default :
516              break;
517            }
# Line 396 | Line 529 | namespace oopse {
529                                     max_vals.rank, 0, status);
530            
531            switch(rnemdType_) {
532 <          case rnemdKinetic :
532 >          case rnemdKineticSwap :
533              min_sd->setVel(max_vel);
534 <            
534 >            // no angular momentum exchange for now
535 >            /*
536              if (min_sd->isDirectional()) {
537                Vector3d min_angMom = min_sd->getJ();
538                Vector3d max_angMom;
539 <
539 >              
540                // point-to-point swap of the angular momentum vector
541                MPI::COMM_WORLD.Sendrecv(min_angMom.getArrayPointer(), 3,
542                                         MPI::REALTYPE, max_vals.rank, 1,
543                                         max_angMom.getArrayPointer(), 3,
544                                         MPI::REALTYPE, max_vals.rank, 1,
545                                         status);
546 <
546 >              
547                min_sd->setJ(max_angMom);
548              }
549 +            */
550              break;
551            case rnemdPx :
552              min_vel.x() = max_vel.x();
# Line 425 | Line 560 | namespace oopse {
560              min_vel.z() = max_vel.z();
561              min_sd->setVel(min_vel);
562              break;
428          case rnemdUnknown :
563            default :
564              break;
565            }
566          }
567   #endif
568          exchangeSum_ += max_val - min_val;
569 <      } else {
570 <        std::cerr << "exchange NOT performed.\nmin_val > max_val.\n";
569 >      } else {        
570 >        sprintf(painCave.errMsg,
571 >                "RNEMD: exchange NOT performed because min_val > max_val\n");
572 >        painCave.isFatal = 0;
573 >        painCave.severity = OPENMD_INFO;
574 >        simError();        
575 >        failTrialCount_++;
576        }
577      } else {
578 <      std::cerr << "exchange NOT performed.\none of the two slabs empty.\n";
578 >      sprintf(painCave.errMsg,
579 >              "RNEMD: exchange NOT performed because at least one\n"
580 >              "\tof the two slabs is empty\n");
581 >      painCave.isFatal = 0;
582 >      painCave.severity = OPENMD_INFO;
583 >      simError();        
584 >      failTrialCount_++;
585      }
586      
587    }
588    
589 <  void RNEMD::getStatus() {
589 >  void RNEMD::doScale() {
590  
591      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
592      Mat3x3d hmat = currentSnap_->getHmat();
448    Stats& stat = currentSnap_->statData;
449    RealType time = currentSnap_->getTime();
593  
451    stat[Stats::RNEMD_SWAP_TOTAL] = exchangeSum_;
452
594      seleMan_.setSelectionSet(evaluator_.evaluate());
595  
596      int selei;
597      StuntDouble* sd;
598      int idx;
599  
600 <    std::vector<RealType> valueHist(nBins_, 0.0); // keeps track of what's
601 <                                                  // being averaged
602 <    std::vector<int> valueCount(nBins_, 0);       // keeps track of the
603 <                                                  // number of degrees of
604 <                                                  // freedom being averaged
600 >    vector<StuntDouble*> hotBin, coldBin;
601 >
602 >    RealType Phx = 0.0;
603 >    RealType Phy = 0.0;
604 >    RealType Phz = 0.0;
605 >    RealType Khx = 0.0;
606 >    RealType Khy = 0.0;
607 >    RealType Khz = 0.0;
608 >    RealType Pcx = 0.0;
609 >    RealType Pcy = 0.0;
610 >    RealType Pcz = 0.0;
611 >    RealType Kcx = 0.0;
612 >    RealType Kcy = 0.0;
613 >    RealType Kcz = 0.0;
614 >
615 >    for (sd = seleMan_.beginSelected(selei); sd != NULL;
616 >         sd = seleMan_.nextSelected(selei)) {
617 >
618 >      idx = sd->getLocalIndex();
619 >
620 >      Vector3d pos = sd->getPos();
621 >
622 >      // wrap the stuntdouble's position back into the box:
623 >
624 >      if (usePeriodicBoundaryConditions_)
625 >        currentSnap_->wrapVector(pos);
626 >
627 >      // which bin is this stuntdouble in?
628 >      // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
629 >
630 >      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_;
631 >
632 >      // if we're in bin 0 or the middleBin
633 >      if (binNo == 0 || binNo == midBin_) {
634 >        
635 >        RealType mass = sd->getMass();
636 >        Vector3d vel = sd->getVel();
637 >      
638 >        if (binNo == 0) {
639 >          hotBin.push_back(sd);
640 >          Phx += mass * vel.x();
641 >          Phy += mass * vel.y();
642 >          Phz += mass * vel.z();
643 >          Khx += mass * vel.x() * vel.x();
644 >          Khy += mass * vel.y() * vel.y();
645 >          Khz += mass * vel.z() * vel.z();
646 >        } else { //midBin_
647 >          coldBin.push_back(sd);
648 >          Pcx += mass * vel.x();
649 >          Pcy += mass * vel.y();
650 >          Pcz += mass * vel.z();
651 >          Kcx += mass * vel.x() * vel.x();
652 >          Kcy += mass * vel.y() * vel.y();
653 >          Kcz += mass * vel.z() * vel.z();
654 >        }
655 >      }
656 >    }
657 >
658 >    Khx *= 0.5;
659 >    Khy *= 0.5;
660 >    Khz *= 0.5;
661 >    Kcx *= 0.5;
662 >    Kcy *= 0.5;
663 >    Kcz *= 0.5;
664 >
665 > #ifdef IS_MPI
666 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phx, 1, MPI::REALTYPE, MPI::SUM);
667 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phy, 1, MPI::REALTYPE, MPI::SUM);
668 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phz, 1, MPI::REALTYPE, MPI::SUM);
669 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcx, 1, MPI::REALTYPE, MPI::SUM);
670 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcy, 1, MPI::REALTYPE, MPI::SUM);
671 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcz, 1, MPI::REALTYPE, MPI::SUM);
672 >
673 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khx, 1, MPI::REALTYPE, MPI::SUM);
674 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khy, 1, MPI::REALTYPE, MPI::SUM);
675 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khz, 1, MPI::REALTYPE, MPI::SUM);
676 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcx, 1, MPI::REALTYPE, MPI::SUM);
677 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcy, 1, MPI::REALTYPE, MPI::SUM);
678 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcz, 1, MPI::REALTYPE, MPI::SUM);
679 > #endif
680 >
681 >    //use coldBin coeff's
682 >    RealType px = Pcx / Phx;
683 >    RealType py = Pcy / Phy;
684 >    RealType pz = Pcz / Phz;
685 >
686 >    RealType a000, a110, c0, a001, a111, b01, b11, c1, c;
687 >    switch(rnemdType_) {
688 >    case rnemdKineticScale :
689 >      // used hotBin coeff's & only scale x & y dimensions
690 >      /*
691 >      RealType px = Phx / Pcx;
692 >      RealType py = Phy / Pcy;
693 >      a110 = Khy;
694 >      c0 = - Khx - Khy - targetFlux_;
695 >      a000 = Khx;
696 >      a111 = Kcy * py * py;
697 >      b11 = -2.0 * Kcy * py * (1.0 + py);
698 >      c1 = Kcy * py * (2.0 + py) + Kcx * px * ( 2.0 + px) + targetFlux_;
699 >      b01 = -2.0 * Kcx * px * (1.0 + px);
700 >      a001 = Kcx * px * px;
701 >      */
702 >      //scale all three dimensions, let c_x = c_y
703 >      a000 = Kcx + Kcy;
704 >      a110 = Kcz;
705 >      c0 = targetFlux_ - Kcx - Kcy - Kcz;
706 >      a001 = Khx * px * px + Khy * py * py;
707 >      a111 = Khz * pz * pz;
708 >      b01 = -2.0 * (Khx * px * (1.0 + px) + Khy * py * (1.0 + py));
709 >      b11 = -2.0 * Khz * pz * (1.0 + pz);
710 >      c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py)
711 >        + Khz * pz * (2.0 + pz) - targetFlux_;
712 >      break;
713 >    case rnemdPxScale :
714 >      c = 1 - targetFlux_ / Pcx;
715 >      a000 = Kcy;
716 >      a110 = Kcz;
717 >      c0 = Kcx * c * c - Kcx - Kcy - Kcz;
718 >      a001 = py * py * Khy;
719 >      a111 = pz * pz * Khz;
720 >      b01 = -2.0 * Khy * py * (1.0 + py);
721 >      b11 = -2.0 * Khz * pz * (1.0 + pz);
722 >      c1 = Khy * py * (2.0 + py) + Khz * pz * (2.0 + pz)
723 >        + Khx * (fastpow(c * px - px - 1.0, 2) - 1.0);
724 >      break;
725 >    case rnemdPyScale :
726 >      c = 1 - targetFlux_ / Pcy;
727 >      a000 = Kcx;
728 >      a110 = Kcz;
729 >      c0 = Kcy * c * c - Kcx - Kcy - Kcz;
730 >      a001 = px * px * Khx;
731 >      a111 = pz * pz * Khz;
732 >      b01 = -2.0 * Khx * px * (1.0 + px);
733 >      b11 = -2.0 * Khz * pz * (1.0 + pz);
734 >      c1 = Khx * px * (2.0 + px) + Khz * pz * (2.0 + pz)
735 >        + Khy * (fastpow(c * py - py - 1.0, 2) - 1.0);
736 >      break;
737 >    case rnemdPzScale ://we don't really do this, do we?
738 >      c = 1 - targetFlux_ / Pcz;
739 >      a000 = Kcx;
740 >      a110 = Kcy;
741 >      c0 = Kcz * c * c - Kcx - Kcy - Kcz;
742 >      a001 = px * px * Khx;
743 >      a111 = py * py * Khy;
744 >      b01 = -2.0 * Khx * px * (1.0 + px);
745 >      b11 = -2.0 * Khy * py * (1.0 + py);
746 >      c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py)
747 >        + Khz * (fastpow(c * pz - pz - 1.0, 2) - 1.0);
748 >      break;      
749 >    default :
750 >      break;
751 >    }
752  
753 +    RealType v1 = a000 * a111 - a001 * a110;
754 +    RealType v2 = a000 * b01;
755 +    RealType v3 = a000 * b11;
756 +    RealType v4 = a000 * c1 - a001 * c0;
757 +    RealType v8 = a110 * b01;
758 +    RealType v10 = - b01 * c0;
759 +
760 +    RealType u0 = v2 * v10 - v4 * v4;
761 +    RealType u1 = -2.0 * v3 * v4;
762 +    RealType u2 = -v2 * v8 - v3 * v3 - 2.0 * v1 * v4;
763 +    RealType u3 = -2.0 * v1 * v3;
764 +    RealType u4 = - v1 * v1;
765 +    //rescale coefficients
766 +    RealType maxAbs = fabs(u0);
767 +    if (maxAbs < fabs(u1)) maxAbs = fabs(u1);
768 +    if (maxAbs < fabs(u2)) maxAbs = fabs(u2);
769 +    if (maxAbs < fabs(u3)) maxAbs = fabs(u3);
770 +    if (maxAbs < fabs(u4)) maxAbs = fabs(u4);
771 +    u0 /= maxAbs;
772 +    u1 /= maxAbs;
773 +    u2 /= maxAbs;
774 +    u3 /= maxAbs;
775 +    u4 /= maxAbs;
776 +    //max_element(start, end) is also available.
777 +    Polynomial<RealType> poly; //same as DoublePolynomial poly;
778 +    poly.setCoefficient(4, u4);
779 +    poly.setCoefficient(3, u3);
780 +    poly.setCoefficient(2, u2);
781 +    poly.setCoefficient(1, u1);
782 +    poly.setCoefficient(0, u0);
783 +    vector<RealType> realRoots = poly.FindRealRoots();
784 +
785 +    vector<RealType>::iterator ri;
786 +    RealType r1, r2, alpha0;
787 +    vector<pair<RealType,RealType> > rps;
788 +    for (ri = realRoots.begin(); ri !=realRoots.end(); ri++) {
789 +      r2 = *ri;
790 +      //check if FindRealRoots() give the right answer
791 +      if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) {
792 +        sprintf(painCave.errMsg,
793 +                "RNEMD Warning: polynomial solve seems to have an error!");
794 +        painCave.isFatal = 0;
795 +        simError();
796 +        failRootCount_++;
797 +      }
798 +      //might not be useful w/o rescaling coefficients
799 +      alpha0 = -c0 - a110 * r2 * r2;
800 +      if (alpha0 >= 0.0) {
801 +        r1 = sqrt(alpha0 / a000);
802 +        if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) < 1e-6)
803 +          { rps.push_back(make_pair(r1, r2)); }
804 +        if (r1 > 1e-6) { //r1 non-negative
805 +          r1 = -r1;
806 +          if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) <1e-6)
807 +            { rps.push_back(make_pair(r1, r2)); }
808 +        }
809 +      }
810 +    }
811 +    // Consider combining together the solving pair part w/ the searching
812 +    // best solution part so that we don't need the pairs vector
813 +    if (!rps.empty()) {
814 +      RealType smallestDiff = HONKING_LARGE_VALUE;
815 +      RealType diff;
816 +      pair<RealType,RealType> bestPair = make_pair(1.0, 1.0);
817 +      vector<pair<RealType,RealType> >::iterator rpi;
818 +      for (rpi = rps.begin(); rpi != rps.end(); rpi++) {
819 +        r1 = (*rpi).first;
820 +        r2 = (*rpi).second;
821 +        switch(rnemdType_) {
822 +        case rnemdKineticScale :
823 +          diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
824 +            + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2)
825 +            + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcy, 2);
826 +          break;
827 +        case rnemdPxScale :
828 +          diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
829 +            + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcy, 2);
830 +          break;
831 +        case rnemdPyScale :
832 +          diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
833 +            + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2);
834 +          break;
835 +        case rnemdPzScale :
836 +          diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
837 +            + fastpow(r1 * r1 / r2 / r2 - Kcy/Kcx, 2);
838 +        default :
839 +          break;
840 +        }
841 +        if (diff < smallestDiff) {
842 +          smallestDiff = diff;
843 +          bestPair = *rpi;
844 +        }
845 +      }
846 + #ifdef IS_MPI
847 +      if (worldRank == 0) {
848 + #endif
849 +        sprintf(painCave.errMsg,
850 +                "RNEMD: roots r1= %lf\tr2 = %lf\n",
851 +                bestPair.first, bestPair.second);
852 +        painCave.isFatal = 0;
853 +        painCave.severity = OPENMD_INFO;
854 +        simError();
855 + #ifdef IS_MPI
856 +      }
857 + #endif
858 +      
859 +      RealType x, y, z;
860 +      switch(rnemdType_) {
861 +      case rnemdKineticScale :
862 +        x = bestPair.first;
863 +        y = bestPair.first;
864 +        z = bestPair.second;
865 +        break;
866 +      case rnemdPxScale :
867 +        x = c;
868 +        y = bestPair.first;
869 +        z = bestPair.second;
870 +        break;
871 +      case rnemdPyScale :
872 +        x = bestPair.first;
873 +        y = c;
874 +        z = bestPair.second;
875 +        break;
876 +      case rnemdPzScale :
877 +        x = bestPair.first;
878 +        y = bestPair.second;
879 +        z = c;
880 +        break;          
881 +      default :
882 +        break;
883 +      }
884 +      vector<StuntDouble*>::iterator sdi;
885 +      Vector3d vel;
886 +      for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
887 +        vel = (*sdi)->getVel();
888 +        vel.x() *= x;
889 +        vel.y() *= y;
890 +        vel.z() *= z;
891 +        (*sdi)->setVel(vel);
892 +      }
893 +      //convert to hotBin coefficient
894 +      x = 1.0 + px * (1.0 - x);
895 +      y = 1.0 + py * (1.0 - y);
896 +      z = 1.0 + pz * (1.0 - z);
897 +      for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
898 +        vel = (*sdi)->getVel();
899 +        vel.x() *= x;
900 +        vel.y() *= y;
901 +        vel.z() *= z;
902 +        (*sdi)->setVel(vel);
903 +      }
904 +      exchangeSum_ += targetFlux_;
905 +      //we may want to check whether the exchange has been successful
906 +    } else {
907 +      sprintf(painCave.errMsg,
908 +              "RNEMD: exchange NOT performed!\n");
909 +      painCave.isFatal = 0;
910 +      painCave.severity = OPENMD_INFO;
911 +      simError();        
912 +      failTrialCount_++;
913 +    }
914 +
915 +  }
916 +
917 +  void RNEMD::doRNEMD() {
918 +
919 +    switch(rnemdType_) {
920 +    case rnemdKineticScale :
921 +    case rnemdPxScale :
922 +    case rnemdPyScale :
923 +    case rnemdPzScale :
924 +      doScale();
925 +      break;
926 +    case rnemdKineticSwap :
927 +    case rnemdPx :
928 +    case rnemdPy :
929 +    case rnemdPz :
930 +      doSwap();
931 +      break;
932 +    case rnemdUnknown :
933 +    default :
934 +      break;
935 +    }
936 +  }
937 +
938 +  void RNEMD::collectData() {
939 +
940 +    Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
941 +    Mat3x3d hmat = currentSnap_->getHmat();
942 +
943 +    seleMan_.setSelectionSet(evaluator_.evaluate());
944 +
945 +    int selei;
946 +    StuntDouble* sd;
947 +    int idx;
948 +
949 +    // alternative approach, track all molecules instead of only those
950 +    // selected for scaling/swapping:
951 +    /*
952 +    SimInfo::MoleculeIterator miter;
953 +    vector<StuntDouble*>::iterator iiter;
954 +    Molecule* mol;
955 +    StuntDouble* integrableObject;
956 +    for (mol = info_->beginMolecule(miter); mol != NULL;
957 +         mol = info_->nextMolecule(miter))
958 +      integrableObject is essentially sd
959 +        for (integrableObject = mol->beginIntegrableObject(iiter);
960 +             integrableObject != NULL;
961 +             integrableObject = mol->nextIntegrableObject(iiter))
962 +    */
963      for (sd = seleMan_.beginSelected(selei); sd != NULL;
964           sd = seleMan_.nextSelected(selei)) {
965        
# Line 477 | Line 975 | namespace oopse {
975        // which bin is this stuntdouble in?
976        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
977        
978 <      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;    
979 <      
978 >      int binNo = int(rnemdLogWidth_ * (pos.z() / hmat(2,2) + 0.5)) %
979 >        rnemdLogWidth_;
980 >      // no symmetrization allowed due to arbitary rnemdLogWidth_ value
981 >      /*
982 >      if (rnemdLogWidth_ == midBin_ + 1)
983 >        if (binNo > midBin_)
984 >          binNo = nBins_ - binNo;
985 >      */
986        RealType mass = sd->getMass();
987        Vector3d vel = sd->getVel();
988        RealType value;
989 +      RealType xVal, yVal, zVal;
990  
991        switch(rnemdType_) {
992 <      case rnemdKinetic :
992 >      case rnemdKineticSwap :
993 >      case rnemdKineticScale :
994          
995 <        value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
490 <                        vel[2]*vel[2]);
995 >        value = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
996          
997 <        valueCount[binNo] += 3;
997 >        valueCount_[binNo] += 3;
998          if (sd->isDirectional()) {
999            Vector3d angMom = sd->getJ();
1000            Mat3x3d I = sd->getI();
# Line 500 | Line 1005 | namespace oopse {
1005              int k = (i + 2) % 3;
1006              value += angMom[j] * angMom[j] / I(j, j) +
1007                angMom[k] * angMom[k] / I(k, k);
1008 <
1009 <            valueCount[binNo] +=2;
1010 <
1008 >            
1009 >            valueCount_[binNo] +=2;
1010 >            
1011            } else {
1012              value += angMom[0]*angMom[0]/I(0, 0)
1013                + angMom[1]*angMom[1]/I(1, 1)
1014                + angMom[2]*angMom[2]/I(2, 2);
1015 <            valueCount[binNo] +=3;
1015 >            valueCount_[binNo] +=3;
1016            }
1017          }
1018 <        value = value / OOPSEConstant::energyConvert / OOPSEConstant::kb;
1019 <
1018 >        value = value / PhysicalConstants::energyConvert / PhysicalConstants::kb;
1019 >        
1020          break;
1021        case rnemdPx :
1022 +      case rnemdPxScale :
1023          value = mass * vel[0];
1024 <        valueCount[binNo]++;
1024 >        valueCount_[binNo]++;
1025          break;
1026        case rnemdPy :
1027 +      case rnemdPyScale :
1028          value = mass * vel[1];
1029 <        valueCount[binNo]++;
1029 >        valueCount_[binNo]++;
1030          break;
1031        case rnemdPz :
1032 <        value = mass * vel[2];
1033 <        valueCount[binNo]++;
1032 >      case rnemdPzScale :
1033 >        value = pos.z(); //temporarily for homogeneous systems ONLY
1034 >        valueCount_[binNo]++;
1035          break;
1036        case rnemdUnknown :
1037        default :
1038 +        value = 1.0;
1039 +        valueCount_[binNo]++;
1040          break;
1041        }
1042 <      valueHist[binNo] += value;
1042 >      valueHist_[binNo] += value;
1043 >
1044 >      if (output3DTemp_) {
1045 >        xVal = mass * vel.x() * vel.x() / PhysicalConstants::energyConvert
1046 >          / PhysicalConstants::kb;
1047 >        yVal = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert
1048 >          / PhysicalConstants::kb;
1049 >        zVal = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert
1050 >          / PhysicalConstants::kb;
1051 >        xTempHist_[binNo] += xVal;
1052 >        yTempHist_[binNo] += yVal;
1053 >        zTempHist_[binNo] += zVal;
1054 >        xyzTempCount_[binNo]++;
1055 >      }
1056      }
1057 +  }
1058  
1059 +  void RNEMD::getStarted() {
1060 +    collectData();
1061 +    /* now should be able to output profile in step 0, but might not be useful
1062 +       Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1063 +       Stats& stat = currentSnap_->statData;
1064 +       stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
1065 +    */
1066 +    getStatus();
1067 +  }
1068 +
1069 +  void RNEMD::getStatus() {
1070 +
1071 +    Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1072 +    Stats& stat = currentSnap_->statData;
1073 +    RealType time = currentSnap_->getTime();
1074 +
1075 +    stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
1076 +    //or to be more meaningful, define another item as exchangeSum_ / time
1077 +    int j;
1078 +
1079   #ifdef IS_MPI
1080  
1081      // all processors have the same number of bins, and STL vectors pack their
1082      // arrays, so in theory, this should be safe:
1083  
1084 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueHist[0],
1085 <                              nBins_, MPI::REALTYPE, MPI::SUM);
1086 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount[0],
1087 <                              nBins_, MPI::INT, MPI::SUM);
1088 <
1084 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueHist_[0],
1085 >                              rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1086 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount_[0],
1087 >                              rnemdLogWidth_, MPI::INT, MPI::SUM);
1088 >    if (output3DTemp_) {
1089 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xTempHist_[0],
1090 >                                rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1091 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &yTempHist_[0],
1092 >                                rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1093 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &zTempHist_[0],
1094 >                                rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1095 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xyzTempCount_[0],
1096 >                                rnemdLogWidth_, MPI::INT, MPI::SUM);
1097 >    }
1098      // If we're the root node, should we print out the results
1099      int worldRank = MPI::COMM_WORLD.Get_rank();
1100      if (worldRank == 0) {
1101   #endif
1102 <      
1103 <      std::cout << time;
1104 <      for (int j = 0; j < nBins_; j++)
1105 <        std::cout << "\t" << valueHist[j] / (RealType)valueCount[j];
1106 <      std::cout << "\n";
1107 <      
1102 >      rnemdLog_ << time;
1103 >      for (j = 0; j < rnemdLogWidth_; j++) {
1104 >        rnemdLog_ << "\t" << valueHist_[j] / (RealType)valueCount_[j];
1105 >      }
1106 >      rnemdLog_ << "\n";
1107 >      if (output3DTemp_) {
1108 >        xTempLog_ << time;      
1109 >        for (j = 0; j < rnemdLogWidth_; j++) {
1110 >          xTempLog_ << "\t" << xTempHist_[j] / (RealType)xyzTempCount_[j];
1111 >        }
1112 >        xTempLog_ << "\n";
1113 >        yTempLog_ << time;
1114 >        for (j = 0; j < rnemdLogWidth_; j++) {
1115 >          yTempLog_ << "\t" << yTempHist_[j] / (RealType)xyzTempCount_[j];
1116 >        }
1117 >        yTempLog_ << "\n";
1118 >        zTempLog_ << time;
1119 >        for (j = 0; j < rnemdLogWidth_; j++) {
1120 >          zTempLog_ << "\t" << zTempHist_[j] / (RealType)xyzTempCount_[j];
1121 >        }
1122 >        zTempLog_ << "\n";
1123 >      }
1124   #ifdef IS_MPI
1125      }
1126   #endif
1127 +    for (j = 0; j < rnemdLogWidth_; j++) {
1128 +      valueCount_[j] = 0;
1129 +      valueHist_[j] = 0.0;
1130 +    }
1131 +    if (output3DTemp_)
1132 +      for (j = 0; j < rnemdLogWidth_; j++) {
1133 +        xTempHist_[j] = 0.0;
1134 +        yTempHist_[j] = 0.0;
1135 +        zTempHist_[j] = 0.0;
1136 +        xyzTempCount_[j] = 0;
1137 +      }
1138    }
1139   }

Comparing:
trunk/src/integrators/RNEMD.cpp (property svn:keywords), Revision 1350 by gezelter, Thu May 21 18:56:45 2009 UTC vs.
branches/development/src/integrators/RNEMD.cpp (property svn:keywords), Revision 1629 by gezelter, Wed Sep 14 21:15:17 2011 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines