ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/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 1338 by skuang, Thu Apr 23 18:22:30 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 59 | Line 62 | namespace oopse {
62  
63   namespace oopse {
64    
65 <  RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info) {
65 >  RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
66      
67      int seedValue;
68      Globals * simParams = info->getSimParams();
# Line 74 | Line 77 | namespace oopse {
77  
78      std::cerr << "calling  evaluator with " << rnemdObjectSelection_ << "\n";
79      evaluator_.loadScriptString(rnemdObjectSelection_);
80 <    std::cerr << "done";
80 >    std::cerr << "done\n";
81      
82      const std::string st = simParams->getRNEMD_swapType();
83  
# Line 85 | Line 88 | namespace oopse {
88      set_RNEMD_swapTime(simParams->getRNEMD_swapTime());
89      set_RNEMD_nBins(simParams->getRNEMD_nBins());
90      exchangeSum_ = 0.0;
91 <    
91 >    counter_ = 0; //added by shenyu
92 >    //profile_.open("profile", std::ios::out);
93 >
94   #ifndef IS_MPI
95      if (simParams->haveSeed()) {
96        seedValue = simParams->getSeed();
# Line 105 | Line 110 | namespace oopse {
110    
111    RNEMD::~RNEMD() {
112      delete randNumGen_;
113 +    //profile_.close();
114    }
115  
116    void RNEMD::doSwap() {
117 <    std::cerr << "in RNEMD!\n";  
118 <    std::cerr << "nBins = " << nBins_ << "\n";
119 <    std::cerr << "swapTime = " << swapTime_ << "\n";
120 <    std::cerr << "exchangeSum = " << exchangeSum_ << "\n";
121 <    std::cerr << "swapType = " << rnemdType_ << "\n";
122 <    std::cerr << "selection = " << rnemdObjectSelection_ << "\n";
117 >    //std::cerr << "in RNEMD!\n";  
118 >    //std::cerr << "nBins = " << nBins_ << "\n";
119 >    int midBin = nBins_ / 2;
120 >    //std::cerr << "midBin = " << midBin << "\n";
121 >    //std::cerr << "swapTime = " << swapTime_ << "\n";
122 >    //std::cerr << "swapType = " << rnemdType_ << "\n";
123 >    //std::cerr << "selection = " << rnemdObjectSelection_ << "\n";
124  
125 +    Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
126 +    Mat3x3d hmat = currentSnap_->getHmat();
127 +
128 +    //std::cerr << "hmat = " << hmat << "\n";
129 +
130      seleMan_.setSelectionSet(evaluator_.evaluate());
131  
132 <    std::cerr << "selectionCount = " << seleMan_.getSelectionCount() << "\n\n";
132 >    //std::cerr << "selectionCount = " << seleMan_.getSelectionCount() << "\n\n";
133  
134 <    int i;
134 >    int selei;
135      StuntDouble* sd;
136 +    int idx;
137  
138 <    for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
138 >    RealType min_val;
139 >    bool min_found = false;  
140 >    StuntDouble* min_sd;
141 >
142 >    RealType max_val;
143 >    bool max_found = false;
144 >    StuntDouble* max_sd;
145 >
146 >    for (sd = seleMan_.beginSelected(selei); sd != NULL;
147 >         sd = seleMan_.nextSelected(selei)) {
148 >
149 >      idx = sd->getLocalIndex();
150 >
151        Vector3d pos = sd->getPos();
152 <      //wrap the stuntdoubles into a cell      
152 >
153 >      // wrap the stuntdouble's position back into the box:
154 >
155        if (usePeriodicBoundaryConditions_)
156 <        info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
157 <      int binNo = int(nBins_ * (pos.z()) / hmat(2,2));
158 <      sliceSDLists_[binNo].push_back(sd);
156 >        currentSnap_->wrapVector(pos);
157 >
158 >      // which bin is this stuntdouble in?
159 >      // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
160 >
161 >      int binNo = midBin + int(nBins_ * (pos.z()) / hmat(2,2));
162 >
163 >      //std::cerr << "pos.z() = " << pos.z() << " bin = " << binNo << "\n";
164 >
165 >      // if we're in bin 0 or the middleBin
166 >      if (binNo == 0 || binNo == midBin) {
167 >        
168 >        RealType mass = sd->getMass();
169 >        Vector3d vel = sd->getVel();
170 >        RealType value;
171 >
172 >        switch(rnemdType_) {
173 >        case rnemdKinetic :
174 >          
175 >          value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
176 >                          vel[2]*vel[2]);
177 >          
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 +    //std::cerr << "smallest value = " << min_val  << "\n";
236 +    //std::cerr << "largest value = " << max_val << "\n";
237 +    
238 +    // missing:  swap information in parallel
239  
240 <  }  
240 >    if (max_found && min_found) {
241 >      if (min_val< max_val) {
242 >        Vector3d min_vel = min_sd->getVel();
243 >        Vector3d max_vel = max_sd->getVel();
244 >        RealType temp_vel;
245 >        switch(rnemdType_) {
246 >        case rnemdKinetic :
247 >          min_sd->setVel(max_vel);
248 >          max_sd->setVel(min_vel);                
249 >          if (min_sd->isDirectional() && max_sd->isDirectional()) {
250 >            Vector3d min_angMom = min_sd->getJ();
251 >            Vector3d max_angMom = max_sd->getJ();
252 >            min_sd->setJ(max_angMom);
253 >            max_sd->setJ(min_angMom);
254 >          }
255 >          break;
256 >        case rnemdPx :
257 >          temp_vel = min_vel.x();
258 >          min_vel.x() = max_vel.x();
259 >          max_vel.x() = temp_vel;
260 >          min_sd->setVel(min_vel);
261 >          max_sd->setVel(max_vel);
262 >          break;
263 >        case rnemdPy :
264 >          temp_vel = min_vel.y();
265 >          min_vel.y() = max_vel.y();
266 >          max_vel.y() = temp_vel;
267 >          min_sd->setVel(min_vel);
268 >          max_sd->setVel(max_vel);
269 >          break;
270 >        case rnemdPz :
271 >          temp_vel = min_vel.z();
272 >          min_vel.z() = max_vel.z();
273 >          max_vel.z() = temp_vel;
274 >          min_sd->setVel(min_vel);
275 >          max_sd->setVel(max_vel);
276 >          break;
277 >        case rnemdUnknown :
278 >        default :
279 >          break;
280 >        }
281 >      exchangeSum_ += max_val - min_val;
282 >      } else {
283 >        std::cerr << "exchange NOT performed.\nmin_val > max_val.\n";
284 >      }
285 >    } else {
286 >      std::cerr << "exchange NOT performed.\none of the two slabs empty.\n";
287 >    }
288 >    std::cerr << "exchangeSum = " << exchangeSum_ << "\n";
289 >  }
290 >
291 >  void RNEMD::getStatus() {
292 >    //std::cerr << "in RNEMD!\n";  
293 >    //std::cerr << "nBins = " << nBins_ << "\n";
294 >    int midBin = nBins_ / 2;
295 >    //std::cerr << "midBin = " << midBin << "\n";
296 >    //std::cerr << "swapTime = " << swapTime_ << "\n";
297 >    //std::cerr << "exchangeSum = " << exchangeSum_ << "\n";
298 >    //std::cerr << "swapType = " << rnemdType_ << "\n";
299 >    //std::cerr << "selection = " << rnemdObjectSelection_ << "\n";
300 >
301 >    Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
302 >    Mat3x3d hmat = currentSnap_->getHmat();
303 >
304 >    //std::cerr << "hmat = " << hmat << "\n";
305 >
306 >    seleMan_.setSelectionSet(evaluator_.evaluate());
307 >
308 >    //std::cerr << "selectionCount = " << seleMan_.getSelectionCount() << "\n\n";
309 >
310 >    int selei;
311 >    StuntDouble* sd;
312 >    int idx;
313 >    /*
314 >    RealType min_val;
315 >    bool min_found = false;  
316 >    StuntDouble* min_sd;
317 >
318 >    RealType max_val;
319 >    bool max_found = false;
320 >    StuntDouble* max_sd;
321 >    */
322 >    std::vector<RealType> valueHist;
323 >    std::vector<int> valueCount;
324 >    valueHist.resize(nBins_);
325 >    valueCount.resize(nBins_);
326 >    //do they initialize themselves to zero automatically?
327 >    for (sd = seleMan_.beginSelected(selei); sd != NULL;
328 >         sd = seleMan_.nextSelected(selei)) {
329 >      
330 >      idx = sd->getLocalIndex();
331 >      
332 >      Vector3d pos = sd->getPos();
333 >
334 >      // wrap the stuntdouble's position back into the box:
335 >      
336 >      if (usePeriodicBoundaryConditions_)
337 >        currentSnap_->wrapVector(pos);
338 >      
339 >      // which bin is this stuntdouble in?
340 >      // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
341 >      
342 >      int binNo = midBin + int(nBins_ * (pos.z()) / hmat(2,2));
343 >      
344 >      //std::cerr << "pos.z() = " << pos.z() << " bin = " << binNo << "\n";
345 >      
346 >      RealType mass = sd->getMass();
347 >      Vector3d vel = sd->getVel();
348 >      //std::cerr << "mass = " << mass << " vel = " << vel << "\n";
349 >      RealType value;
350 >
351 >      switch(rnemdType_) {
352 >      case rnemdKinetic :
353 >        
354 >        value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
355 >                        vel[2]*vel[2]);
356 >        
357 >        if (sd->isDirectional()) {
358 >          Vector3d angMom = sd->getJ();
359 >          Mat3x3d I = sd->getI();
360 >          
361 >          if (sd->isLinear()) {
362 >            int i = sd->linearAxis();
363 >            int j = (i + 1) % 3;
364 >            int k = (i + 2) % 3;
365 >            value += angMom[j] * angMom[j] / I(j, j) +
366 >              angMom[k] * angMom[k] / I(k, k);
367 >          } else {                        
368 >            value += angMom[0]*angMom[0]/I(0, 0)
369 >              + angMom[1]*angMom[1]/I(1, 1)
370 >              + angMom[2]*angMom[2]/I(2, 2);
371 >          }
372 >        }
373 >        //std::cerr <<"this value = " << value << "\n";
374 >        value = value * 0.5 / OOPSEConstant::energyConvert;
375 >        //std::cerr <<"this value = " << value << "\n";
376 >        break;
377 >      case rnemdPx :
378 >        value = mass * vel[0];
379 >        break;
380 >      case rnemdPy :
381 >        value = mass * vel[1];
382 >        break;
383 >      case rnemdPz :
384 >        value = mass * vel[2];
385 >        break;
386 >      case rnemdUnknown :
387 >      default :
388 >        break;
389 >      }
390 >      //std::cerr << "bin = " << binNo << " value = " << value ;
391 >      valueHist[binNo] += value;
392 >      valueCount[binNo]++;
393 >      //std::cerr << " hist = " << valueHist[binNo] << " count = " << valueCount[binNo] << "\n";
394 >    }
395 >    
396 >    std::cout << counter_++;
397 >    for (int j = 0; j < nBins_; j++)
398 >      std::cout << "\t" << valueHist[j] / (RealType)valueCount[j];
399 >    std::cout << "\n";
400 >  }
401   }

Diff Legend

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