# | Line 35 | Line 35 | |
---|---|---|
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). |
38 | > | * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). |
39 | * [4] Vardeman & Gezelter, in progress (2009). | |
40 | */ | |
41 | ||
42 | #include <cmath> | |
43 | + | #include <sstream> |
44 | + | #include <string> |
45 | + | |
46 | #include "rnemd/RNEMD.hpp" | |
47 | #include "math/Vector3.hpp" | |
48 | #include "math/Vector.hpp" | |
# | Line 49 | Line 52 | |
52 | #include "primitives/StuntDouble.hpp" | |
53 | #include "utils/PhysicalConstants.hpp" | |
54 | #include "utils/Tuple.hpp" | |
55 | + | #include "brains/Thermo.hpp" |
56 | + | #include "math/ConvexHull.hpp" |
57 | #ifdef IS_MPI | |
58 | #include <mpi.h> | |
59 | #endif | |
60 | ||
61 | + | #ifdef _MSC_VER |
62 | + | #define isnan(x) _isnan((x)) |
63 | + | #define isinf(x) (!_finite(x) && !_isnan(x)) |
64 | + | #endif |
65 | + | |
66 | #define HONKING_LARGE_VALUE 1.0e10 | |
67 | ||
68 | using namespace std; | |
69 | namespace OpenMD { | |
70 | ||
71 | RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), | |
72 | + | evaluatorA_(info), seleManA_(info), |
73 | + | commonA_(info), evaluatorB_(info), |
74 | + | seleManB_(info), commonB_(info), |
75 | usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) { | |
76 | ||
77 | trialCount_ = 0; | |
78 | failTrialCount_ = 0; | |
79 | failRootCount_ = 0; | |
80 | ||
81 | < | int seedValue; |
69 | < | Globals * simParams = info->getSimParams(); |
81 | > | Globals* simParams = info->getSimParams(); |
82 | RNEMDParameters* rnemdParams = simParams->getRNEMDParameters(); | |
83 | ||
84 | + | doRNEMD_ = rnemdParams->getUseRNEMD(); |
85 | + | if (!doRNEMD_) return; |
86 | + | |
87 | stringToMethod_["Swap"] = rnemdSwap; | |
88 | stringToMethod_["NIVS"] = rnemdNIVS; | |
89 | stringToMethod_["VSS"] = rnemdVSS; | |
# | Line 77 | Line 92 | namespace OpenMD { | |
92 | stringToFluxType_["Px"] = rnemdPx; | |
93 | stringToFluxType_["Py"] = rnemdPy; | |
94 | stringToFluxType_["Pz"] = rnemdPz; | |
95 | + | stringToFluxType_["Pvector"] = rnemdPvector; |
96 | + | stringToFluxType_["Lx"] = rnemdLx; |
97 | + | stringToFluxType_["Ly"] = rnemdLy; |
98 | + | stringToFluxType_["Lz"] = rnemdLz; |
99 | + | stringToFluxType_["Lvector"] = rnemdLvector; |
100 | stringToFluxType_["KE+Px"] = rnemdKePx; | |
101 | stringToFluxType_["KE+Py"] = rnemdKePy; | |
102 | stringToFluxType_["KE+Pvector"] = rnemdKePvector; | |
103 | + | stringToFluxType_["KE+Lx"] = rnemdKeLx; |
104 | + | stringToFluxType_["KE+Ly"] = rnemdKeLy; |
105 | + | stringToFluxType_["KE+Lz"] = rnemdKeLz; |
106 | + | stringToFluxType_["KE+Lvector"] = rnemdKeLvector; |
107 | ||
108 | runTime_ = simParams->getRunTime(); | |
109 | statusTime_ = simParams->getStatusTime(); | |
110 | ||
87 | – | rnemdObjectSelection_ = rnemdParams->getObjectSelection(); |
88 | – | evaluator_.loadScriptString(rnemdObjectSelection_); |
89 | – | seleMan_.setSelectionSet(evaluator_.evaluate()); |
90 | – | |
111 | const string methStr = rnemdParams->getMethod(); | |
112 | bool hasFluxType = rnemdParams->haveFluxType(); | |
113 | ||
114 | + | rnemdObjectSelection_ = rnemdParams->getObjectSelection(); |
115 | + | |
116 | string fluxStr; | |
117 | if (hasFluxType) { | |
118 | fluxStr = rnemdParams->getFluxType(); | |
# | Line 98 | Line 120 | namespace OpenMD { | |
120 | sprintf(painCave.errMsg, | |
121 | "RNEMD: No fluxType was set in the md file. This parameter,\n" | |
122 | "\twhich must be one of the following values:\n" | |
123 | < | "\tKE, Px, Py, Pz, KE+Px, KE+Py, KE+Pvector, must be set to\n" |
124 | < | "\tuse RNEMD\n"); |
123 | > | "\tKE, Px, Py, Pz, Pvector, Lx, Ly, Lz, Lvector,\n" |
124 | > | "\tKE+Px, KE+Py, KE+Pvector, KE+Lx, KE+Ly, KE+Lz, KE+Lvector\n" |
125 | > | "\tmust be set to use RNEMD\n"); |
126 | painCave.isFatal = 1; | |
127 | painCave.severity = OPENMD_ERROR; | |
128 | simError(); | |
# | Line 108 | Line 131 | namespace OpenMD { | |
131 | bool hasKineticFlux = rnemdParams->haveKineticFlux(); | |
132 | bool hasMomentumFlux = rnemdParams->haveMomentumFlux(); | |
133 | bool hasMomentumFluxVector = rnemdParams->haveMomentumFluxVector(); | |
134 | + | bool hasAngularMomentumFlux = rnemdParams->haveAngularMomentumFlux(); |
135 | + | bool hasAngularMomentumFluxVector = rnemdParams->haveAngularMomentumFluxVector(); |
136 | + | hasSelectionA_ = rnemdParams->haveSelectionA(); |
137 | + | hasSelectionB_ = rnemdParams->haveSelectionB(); |
138 | bool hasSlabWidth = rnemdParams->haveSlabWidth(); | |
139 | bool hasSlabACenter = rnemdParams->haveSlabACenter(); | |
140 | bool hasSlabBCenter = rnemdParams->haveSlabBCenter(); | |
141 | + | bool hasSphereARadius = rnemdParams->haveSphereARadius(); |
142 | + | hasSphereBRadius_ = rnemdParams->haveSphereBRadius(); |
143 | + | bool hasCoordinateOrigin = rnemdParams->haveCoordinateOrigin(); |
144 | bool hasOutputFileName = rnemdParams->haveOutputFileName(); | |
145 | bool hasOutputFields = rnemdParams->haveOutputFields(); | |
146 | ||
# | Line 195 | Line 225 | namespace OpenMD { | |
225 | case rnemdPz: | |
226 | hasCorrectFlux = hasMomentumFlux; | |
227 | break; | |
228 | + | case rnemdLx: |
229 | + | case rnemdLy: |
230 | + | case rnemdLz: |
231 | + | hasCorrectFlux = hasAngularMomentumFlux; |
232 | + | break; |
233 | case rnemdPvector: | |
234 | hasCorrectFlux = hasMomentumFluxVector; | |
235 | + | break; |
236 | + | case rnemdLvector: |
237 | + | hasCorrectFlux = hasAngularMomentumFluxVector; |
238 | + | break; |
239 | case rnemdKePx: | |
240 | case rnemdKePy: | |
241 | hasCorrectFlux = hasMomentumFlux && hasKineticFlux; | |
242 | break; | |
243 | + | case rnemdKeLx: |
244 | + | case rnemdKeLy: |
245 | + | case rnemdKeLz: |
246 | + | hasCorrectFlux = hasAngularMomentumFlux && hasKineticFlux; |
247 | + | break; |
248 | case rnemdKePvector: | |
249 | hasCorrectFlux = hasMomentumFluxVector && hasKineticFlux; | |
250 | break; | |
251 | + | case rnemdKeLvector: |
252 | + | hasCorrectFlux = hasAngularMomentumFluxVector && hasKineticFlux; |
253 | + | break; |
254 | default: | |
255 | methodFluxMismatch = true; | |
256 | break; | |
# | Line 224 | Line 271 | namespace OpenMD { | |
271 | } | |
272 | if (!hasCorrectFlux) { | |
273 | sprintf(painCave.errMsg, | |
274 | < | "RNEMD: The current method,\n" |
228 | < | "\t%s, and flux type %s\n" |
274 | > | "RNEMD: The current method, %s, and flux type, %s,\n" |
275 | "\tdid not have the correct flux value specified. Options\n" | |
276 | < | "\tinclude: kineticFlux, momentumFlux, and momentumFluxVector\n", |
276 | > | "\tinclude: kineticFlux, momentumFlux, angularMomentumFlux,\n" |
277 | > | "\tmomentumFluxVector, and angularMomentumFluxVector.\n", |
278 | methStr.c_str(), fluxStr.c_str()); | |
279 | painCave.isFatal = 1; | |
280 | painCave.severity = OPENMD_ERROR; | |
# | Line 235 | Line 282 | namespace OpenMD { | |
282 | } | |
283 | ||
284 | if (hasKineticFlux) { | |
285 | < | kineticFlux_ = rnemdParams->getKineticFlux(); |
285 | > | // convert the kcal / mol / Angstroms^2 / fs values in the md file |
286 | > | // into amu / fs^3: |
287 | > | kineticFlux_ = rnemdParams->getKineticFlux() |
288 | > | * PhysicalConstants::energyConvert; |
289 | } else { | |
290 | kineticFlux_ = 0.0; | |
291 | } | |
# | Line 264 | Line 314 | namespace OpenMD { | |
314 | default: | |
315 | break; | |
316 | } | |
317 | < | } |
318 | < | } |
319 | < | |
320 | < | // do some sanity checking |
321 | < | |
322 | < | int selectionCount = seleMan_.getSelectionCount(); |
323 | < | int nIntegrable = info->getNGlobalIntegrableObjects(); |
324 | < | |
325 | < | if (selectionCount > nIntegrable) { |
326 | < | sprintf(painCave.errMsg, |
327 | < | "RNEMD: The current objectSelection,\n" |
328 | < | "\t\t%s\n" |
329 | < | "\thas resulted in %d selected objects. However,\n" |
330 | < | "\tthe total number of integrable objects in the system\n" |
331 | < | "\tis only %d. This is almost certainly not what you want\n" |
332 | < | "\tto do. A likely cause of this is forgetting the _RB_0\n" |
333 | < | "\tselector in the selection script!\n", |
334 | < | rnemdObjectSelection_.c_str(), |
335 | < | selectionCount, nIntegrable); |
336 | < | painCave.isFatal = 0; |
337 | < | painCave.severity = OPENMD_WARNING; |
338 | < | simError(); |
339 | < | } |
340 | < | |
341 | < | nBins_ = rnemdParams->getOutputBins(); |
317 | > | } |
318 | > | if (hasAngularMomentumFluxVector) { |
319 | > | angularMomentumFluxVector_ = rnemdParams->getAngularMomentumFluxVector(); |
320 | > | } else { |
321 | > | angularMomentumFluxVector_ = V3Zero; |
322 | > | if (hasAngularMomentumFlux) { |
323 | > | RealType angularMomentumFlux = rnemdParams->getAngularMomentumFlux(); |
324 | > | switch (rnemdFluxType_) { |
325 | > | case rnemdLx: |
326 | > | angularMomentumFluxVector_.x() = angularMomentumFlux; |
327 | > | break; |
328 | > | case rnemdLy: |
329 | > | angularMomentumFluxVector_.y() = angularMomentumFlux; |
330 | > | break; |
331 | > | case rnemdLz: |
332 | > | angularMomentumFluxVector_.z() = angularMomentumFlux; |
333 | > | break; |
334 | > | case rnemdKeLx: |
335 | > | angularMomentumFluxVector_.x() = angularMomentumFlux; |
336 | > | break; |
337 | > | case rnemdKeLy: |
338 | > | angularMomentumFluxVector_.y() = angularMomentumFlux; |
339 | > | break; |
340 | > | case rnemdKeLz: |
341 | > | angularMomentumFluxVector_.z() = angularMomentumFlux; |
342 | > | break; |
343 | > | default: |
344 | > | break; |
345 | > | } |
346 | > | } |
347 | > | } |
348 | ||
349 | < | data_.resize(RNEMD::ENDINDEX); |
350 | < | OutputData z; |
351 | < | z.units = "Angstroms"; |
352 | < | z.title = "Z"; |
353 | < | z.dataType = "RealType"; |
298 | < | z.accumulator.reserve(nBins_); |
299 | < | for (unsigned int i = 0; i < nBins_; i++) |
300 | < | z.accumulator.push_back( new Accumulator() ); |
301 | < | data_[Z] = z; |
302 | < | outputMap_["Z"] = Z; |
349 | > | if (hasCoordinateOrigin) { |
350 | > | coordinateOrigin_ = rnemdParams->getCoordinateOrigin(); |
351 | > | } else { |
352 | > | coordinateOrigin_ = V3Zero; |
353 | > | } |
354 | ||
355 | < | OutputData temperature; |
305 | < | temperature.units = "K"; |
306 | < | temperature.title = "Temperature"; |
307 | < | temperature.dataType = "RealType"; |
308 | < | temperature.accumulator.reserve(nBins_); |
309 | < | for (unsigned int i = 0; i < nBins_; i++) |
310 | < | temperature.accumulator.push_back( new Accumulator() ); |
311 | < | data_[TEMPERATURE] = temperature; |
312 | < | outputMap_["TEMPERATURE"] = TEMPERATURE; |
355 | > | // do some sanity checking |
356 | ||
357 | < | OutputData velocity; |
315 | < | velocity.units = "amu/fs"; |
316 | < | velocity.title = "Velocity"; |
317 | < | velocity.dataType = "Vector3d"; |
318 | < | velocity.accumulator.reserve(nBins_); |
319 | < | for (unsigned int i = 0; i < nBins_; i++) |
320 | < | velocity.accumulator.push_back( new VectorAccumulator() ); |
321 | < | data_[VELOCITY] = velocity; |
322 | < | outputMap_["VELOCITY"] = VELOCITY; |
357 | > | int selectionCount = seleMan_.getSelectionCount(); |
358 | ||
359 | < | OutputData density; |
360 | < | density.units = "g cm^-3"; |
361 | < | density.title = "Density"; |
362 | < | density.dataType = "RealType"; |
363 | < | density.accumulator.reserve(nBins_); |
364 | < | for (unsigned int i = 0; i < nBins_; i++) |
365 | < | density.accumulator.push_back( new Accumulator() ); |
366 | < | data_[DENSITY] = density; |
367 | < | outputMap_["DENSITY"] = DENSITY; |
359 | > | int nIntegrable = info->getNGlobalIntegrableObjects(); |
360 | > | |
361 | > | if (selectionCount > nIntegrable) { |
362 | > | sprintf(painCave.errMsg, |
363 | > | "RNEMD: The current objectSelection,\n" |
364 | > | "\t\t%s\n" |
365 | > | "\thas resulted in %d selected objects. However,\n" |
366 | > | "\tthe total number of integrable objects in the system\n" |
367 | > | "\tis only %d. This is almost certainly not what you want\n" |
368 | > | "\tto do. A likely cause of this is forgetting the _RB_0\n" |
369 | > | "\tselector in the selection script!\n", |
370 | > | rnemdObjectSelection_.c_str(), |
371 | > | selectionCount, nIntegrable); |
372 | > | painCave.isFatal = 0; |
373 | > | painCave.severity = OPENMD_WARNING; |
374 | > | simError(); |
375 | > | } |
376 | > | |
377 | > | areaAccumulator_ = new Accumulator(); |
378 | > | |
379 | > | nBins_ = rnemdParams->getOutputBins(); |
380 | > | binWidth_ = rnemdParams->getOutputBinWidth(); |
381 | > | |
382 | > | data_.resize(RNEMD::ENDINDEX); |
383 | > | OutputData z; |
384 | > | z.units = "Angstroms"; |
385 | > | z.title = "Z"; |
386 | > | z.dataType = "RealType"; |
387 | > | z.accumulator.reserve(nBins_); |
388 | > | for (int i = 0; i < nBins_; i++) |
389 | > | z.accumulator.push_back( new Accumulator() ); |
390 | > | data_[Z] = z; |
391 | > | outputMap_["Z"] = Z; |
392 | ||
393 | < | if (hasOutputFields) { |
394 | < | parseOutputFileFormat(rnemdParams->getOutputFields()); |
395 | < | } else { |
396 | < | outputMask_.set(Z); |
397 | < | switch (rnemdFluxType_) { |
398 | < | case rnemdKE: |
399 | < | case rnemdRotKE: |
400 | < | case rnemdFullKE: |
401 | < | outputMask_.set(TEMPERATURE); |
343 | < | break; |
344 | < | case rnemdPx: |
345 | < | case rnemdPy: |
346 | < | outputMask_.set(VELOCITY); |
347 | < | break; |
348 | < | case rnemdPz: |
349 | < | case rnemdPvector: |
350 | < | outputMask_.set(VELOCITY); |
351 | < | outputMask_.set(DENSITY); |
352 | < | break; |
353 | < | case rnemdKePx: |
354 | < | case rnemdKePy: |
355 | < | outputMask_.set(TEMPERATURE); |
356 | < | outputMask_.set(VELOCITY); |
357 | < | break; |
358 | < | case rnemdKePvector: |
359 | < | outputMask_.set(TEMPERATURE); |
360 | < | outputMask_.set(VELOCITY); |
361 | < | outputMask_.set(DENSITY); |
362 | < | break; |
363 | < | default: |
364 | < | break; |
365 | < | } |
366 | < | } |
367 | < | |
368 | < | if (hasOutputFileName) { |
369 | < | rnemdFileName_ = rnemdParams->getOutputFileName(); |
370 | < | } else { |
371 | < | rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd"; |
372 | < | } |
393 | > | OutputData r; |
394 | > | r.units = "Angstroms"; |
395 | > | r.title = "R"; |
396 | > | r.dataType = "RealType"; |
397 | > | r.accumulator.reserve(nBins_); |
398 | > | for (int i = 0; i < nBins_; i++) |
399 | > | r.accumulator.push_back( new Accumulator() ); |
400 | > | data_[R] = r; |
401 | > | outputMap_["R"] = R; |
402 | ||
403 | < | exchangeTime_ = rnemdParams->getExchangeTime(); |
403 | > | OutputData temperature; |
404 | > | temperature.units = "K"; |
405 | > | temperature.title = "Temperature"; |
406 | > | temperature.dataType = "RealType"; |
407 | > | temperature.accumulator.reserve(nBins_); |
408 | > | for (int i = 0; i < nBins_; i++) |
409 | > | temperature.accumulator.push_back( new Accumulator() ); |
410 | > | data_[TEMPERATURE] = temperature; |
411 | > | outputMap_["TEMPERATURE"] = TEMPERATURE; |
412 | ||
413 | < | Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot(); |
414 | < | Mat3x3d hmat = currentSnap_->getHmat(); |
415 | < | |
416 | < | // Target exchange quantities (in each exchange) = 2 Lx Ly dt flux |
417 | < | // Lx, Ly = box dimensions in x & y |
418 | < | // dt = exchange time interval |
419 | < | // flux = target flux |
413 | > | OutputData velocity; |
414 | > | velocity.units = "angstroms/fs"; |
415 | > | velocity.title = "Velocity"; |
416 | > | velocity.dataType = "Vector3d"; |
417 | > | velocity.accumulator.reserve(nBins_); |
418 | > | for (int i = 0; i < nBins_; i++) |
419 | > | velocity.accumulator.push_back( new VectorAccumulator() ); |
420 | > | data_[VELOCITY] = velocity; |
421 | > | outputMap_["VELOCITY"] = VELOCITY; |
422 | ||
423 | < | kineticTarget_ = 2.0*kineticFlux_*exchangeTime_*hmat(0,0)*hmat(1,1); |
424 | < | momentumTarget_ = 2.0*momentumFluxVector_*exchangeTime_*hmat(0,0)*hmat(1,1); |
423 | > | OutputData angularVelocity; |
424 | > | angularVelocity.units = "angstroms^2/fs"; |
425 | > | angularVelocity.title = "AngularVelocity"; |
426 | > | angularVelocity.dataType = "Vector3d"; |
427 | > | angularVelocity.accumulator.reserve(nBins_); |
428 | > | for (int i = 0; i < nBins_; i++) |
429 | > | angularVelocity.accumulator.push_back( new VectorAccumulator() ); |
430 | > | data_[ANGULARVELOCITY] = angularVelocity; |
431 | > | outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY; |
432 | ||
433 | < | // total exchange sums are zeroed out at the beginning: |
433 | > | OutputData density; |
434 | > | density.units = "g cm^-3"; |
435 | > | density.title = "Density"; |
436 | > | density.dataType = "RealType"; |
437 | > | density.accumulator.reserve(nBins_); |
438 | > | for (int i = 0; i < nBins_; i++) |
439 | > | density.accumulator.push_back( new Accumulator() ); |
440 | > | data_[DENSITY] = density; |
441 | > | outputMap_["DENSITY"] = DENSITY; |
442 | ||
443 | < | kineticExchange_ = 0.0; |
444 | < | momentumExchange_ = V3Zero; |
443 | > | if (hasOutputFields) { |
444 | > | parseOutputFileFormat(rnemdParams->getOutputFields()); |
445 | > | } else { |
446 | > | if (usePeriodicBoundaryConditions_) |
447 | > | outputMask_.set(Z); |
448 | > | else |
449 | > | outputMask_.set(R); |
450 | > | switch (rnemdFluxType_) { |
451 | > | case rnemdKE: |
452 | > | case rnemdRotKE: |
453 | > | case rnemdFullKE: |
454 | > | outputMask_.set(TEMPERATURE); |
455 | > | break; |
456 | > | case rnemdPx: |
457 | > | case rnemdPy: |
458 | > | outputMask_.set(VELOCITY); |
459 | > | break; |
460 | > | case rnemdPz: |
461 | > | case rnemdPvector: |
462 | > | outputMask_.set(VELOCITY); |
463 | > | outputMask_.set(DENSITY); |
464 | > | break; |
465 | > | case rnemdLx: |
466 | > | case rnemdLy: |
467 | > | case rnemdLz: |
468 | > | case rnemdLvector: |
469 | > | outputMask_.set(ANGULARVELOCITY); |
470 | > | break; |
471 | > | case rnemdKeLx: |
472 | > | case rnemdKeLy: |
473 | > | case rnemdKeLz: |
474 | > | case rnemdKeLvector: |
475 | > | outputMask_.set(TEMPERATURE); |
476 | > | outputMask_.set(ANGULARVELOCITY); |
477 | > | break; |
478 | > | case rnemdKePx: |
479 | > | case rnemdKePy: |
480 | > | outputMask_.set(TEMPERATURE); |
481 | > | outputMask_.set(VELOCITY); |
482 | > | break; |
483 | > | case rnemdKePvector: |
484 | > | outputMask_.set(TEMPERATURE); |
485 | > | outputMask_.set(VELOCITY); |
486 | > | outputMask_.set(DENSITY); |
487 | > | break; |
488 | > | default: |
489 | > | break; |
490 | > | } |
491 | > | } |
492 | > | |
493 | > | if (hasOutputFileName) { |
494 | > | rnemdFileName_ = rnemdParams->getOutputFileName(); |
495 | > | } else { |
496 | > | rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd"; |
497 | > | } |
498 | ||
499 | < | if (hasSlabWidth) |
500 | < | slabWidth_ = rnemdParams->getSlabWidth(); |
501 | < | else |
502 | < | slabWidth_ = hmat(2,2) / 10.0; |
503 | < | |
504 | < | if (hasSlabACenter) |
505 | < | slabACenter_ = rnemdParams->getSlabACenter(); |
506 | < | else |
507 | < | slabACenter_ = 0.0; |
499 | > | exchangeTime_ = rnemdParams->getExchangeTime(); |
500 | > | |
501 | > | Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot(); |
502 | > | // total exchange sums are zeroed out at the beginning: |
503 | > | |
504 | > | kineticExchange_ = 0.0; |
505 | > | momentumExchange_ = V3Zero; |
506 | > | angularMomentumExchange_ = V3Zero; |
507 | > | |
508 | > | std::ostringstream selectionAstream; |
509 | > | std::ostringstream selectionBstream; |
510 | ||
511 | < | if (hasSlabBCenter) |
512 | < | slabBCenter_ = rnemdParams->getSlabBCenter(); |
513 | < | else |
514 | < | slabBCenter_ = hmat(2,2) / 2.0; |
511 | > | if (hasSelectionA_) { |
512 | > | selectionA_ = rnemdParams->getSelectionA(); |
513 | > | } else { |
514 | > | if (usePeriodicBoundaryConditions_) { |
515 | > | Mat3x3d hmat = currentSnap_->getHmat(); |
516 | > | |
517 | > | if (hasSlabWidth) |
518 | > | slabWidth_ = rnemdParams->getSlabWidth(); |
519 | > | else |
520 | > | slabWidth_ = hmat(2,2) / 10.0; |
521 | > | |
522 | > | if (hasSlabACenter) |
523 | > | slabACenter_ = rnemdParams->getSlabACenter(); |
524 | > | else |
525 | > | slabACenter_ = 0.0; |
526 | > | |
527 | > | selectionAstream << "select wrappedz > " |
528 | > | << slabACenter_ - 0.5*slabWidth_ |
529 | > | << " && wrappedz < " |
530 | > | << slabACenter_ + 0.5*slabWidth_; |
531 | > | selectionA_ = selectionAstream.str(); |
532 | > | } else { |
533 | > | if (hasSphereARadius) |
534 | > | sphereARadius_ = rnemdParams->getSphereARadius(); |
535 | > | else { |
536 | > | // use an initial guess to the size of the inner slab to be 1/10 the |
537 | > | // radius of an approximately spherical hull: |
538 | > | Thermo thermo(info); |
539 | > | RealType hVol = thermo.getHullVolume(); |
540 | > | sphereARadius_ = 0.1 * pow((3.0 * hVol / (4.0 * M_PI)), 1.0/3.0); |
541 | > | } |
542 | > | selectionAstream << "select r < " << sphereARadius_; |
543 | > | selectionA_ = selectionAstream.str(); |
544 | > | } |
545 | > | } |
546 | ||
547 | + | if (hasSelectionB_) { |
548 | + | selectionB_ = rnemdParams->getSelectionB(); |
549 | + | } else { |
550 | + | if (usePeriodicBoundaryConditions_) { |
551 | + | Mat3x3d hmat = currentSnap_->getHmat(); |
552 | + | |
553 | + | if (hasSlabWidth) |
554 | + | slabWidth_ = rnemdParams->getSlabWidth(); |
555 | + | else |
556 | + | slabWidth_ = hmat(2,2) / 10.0; |
557 | + | |
558 | + | if (hasSlabBCenter) |
559 | + | slabBCenter_ = rnemdParams->getSlabACenter(); |
560 | + | else |
561 | + | slabBCenter_ = hmat(2,2) / 2.0; |
562 | + | |
563 | + | selectionBstream << "select wrappedz > " |
564 | + | << slabBCenter_ - 0.5*slabWidth_ |
565 | + | << " && wrappedz < " |
566 | + | << slabBCenter_ + 0.5*slabWidth_; |
567 | + | selectionB_ = selectionBstream.str(); |
568 | + | } else { |
569 | + | if (hasSphereBRadius_) { |
570 | + | sphereBRadius_ = rnemdParams->getSphereBRadius(); |
571 | + | selectionBstream << "select r > " << sphereBRadius_; |
572 | + | selectionB_ = selectionBstream.str(); |
573 | + | } else { |
574 | + | selectionB_ = "select hull"; |
575 | + | hasSelectionB_ = true; |
576 | + | } |
577 | + | } |
578 | + | } |
579 | + | } |
580 | + | // object evaluator: |
581 | + | evaluator_.loadScriptString(rnemdObjectSelection_); |
582 | + | seleMan_.setSelectionSet(evaluator_.evaluate()); |
583 | + | |
584 | + | evaluatorA_.loadScriptString(selectionA_); |
585 | + | evaluatorB_.loadScriptString(selectionB_); |
586 | + | |
587 | + | seleManA_.setSelectionSet(evaluatorA_.evaluate()); |
588 | + | seleManB_.setSelectionSet(evaluatorB_.evaluate()); |
589 | + | |
590 | + | commonA_ = seleManA_ & seleMan_; |
591 | + | commonB_ = seleManB_ & seleMan_; |
592 | } | |
593 | < | |
409 | < | RNEMD::~RNEMD() { |
593 | > | |
594 | ||
595 | + | RNEMD::~RNEMD() { |
596 | + | if (!doRNEMD_) return; |
597 | #ifdef IS_MPI | |
598 | if (worldRank == 0) { | |
599 | #endif | |
# | Line 421 | Line 607 | namespace OpenMD { | |
607 | #endif | |
608 | } | |
609 | ||
610 | < | bool RNEMD::inSlabA(Vector3d pos) { |
611 | < | return (abs(pos.z() - slabACenter_) < 0.5*slabWidth_); |
612 | < | } |
613 | < | bool RNEMD::inSlabB(Vector3d pos) { |
428 | < | return (abs(pos.z() - slabBCenter_) < 0.5*slabWidth_); |
429 | < | } |
610 | > | void RNEMD::doSwap(SelectionManager& smanA, SelectionManager& smanB) { |
611 | > | if (!doRNEMD_) return; |
612 | > | int selei; |
613 | > | int selej; |
614 | ||
431 | – | void RNEMD::doSwap() { |
432 | – | |
615 | Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); | |
616 | Mat3x3d hmat = currentSnap_->getHmat(); | |
617 | ||
436 | – | seleMan_.setSelectionSet(evaluator_.evaluate()); |
437 | – | |
438 | – | int selei; |
618 | StuntDouble* sd; | |
440 | – | int idx; |
619 | ||
620 | RealType min_val; | |
621 | bool min_found = false; | |
# | Line 447 | Line 625 | namespace OpenMD { | |
625 | bool max_found = false; | |
626 | StuntDouble* max_sd; | |
627 | ||
628 | < | for (sd = seleMan_.beginSelected(selei); sd != NULL; |
629 | < | sd = seleMan_.nextSelected(selei)) { |
628 | > | for (sd = seleManA_.beginSelected(selei); sd != NULL; |
629 | > | sd = seleManA_.nextSelected(selei)) { |
630 | ||
453 | – | idx = sd->getLocalIndex(); |
454 | – | |
631 | Vector3d pos = sd->getPos(); | |
632 | < | |
632 | > | |
633 | // wrap the stuntdouble's position back into the box: | |
634 | < | |
634 | > | |
635 | if (usePeriodicBoundaryConditions_) | |
636 | currentSnap_->wrapVector(pos); | |
637 | < | bool inA = inSlabA(pos); |
638 | < | bool inB = inSlabB(pos); |
639 | < | |
640 | < | if (inA || inB) { |
637 | > | |
638 | > | RealType mass = sd->getMass(); |
639 | > | Vector3d vel = sd->getVel(); |
640 | > | RealType value; |
641 | > | |
642 | > | switch(rnemdFluxType_) { |
643 | > | case rnemdKE : |
644 | ||
645 | < | RealType mass = sd->getMass(); |
646 | < | Vector3d vel = sd->getVel(); |
647 | < | RealType value; |
648 | < | |
649 | < | switch(rnemdFluxType_) { |
471 | < | case rnemdKE : |
645 | > | value = mass * vel.lengthSquare(); |
646 | > | |
647 | > | if (sd->isDirectional()) { |
648 | > | Vector3d angMom = sd->getJ(); |
649 | > | Mat3x3d I = sd->getI(); |
650 | ||
651 | < | value = mass * vel.lengthSquare(); |
652 | < | |
653 | < | if (sd->isDirectional()) { |
654 | < | Vector3d angMom = sd->getJ(); |
655 | < | Mat3x3d I = sd->getI(); |
656 | < | |
657 | < | if (sd->isLinear()) { |
658 | < | int i = sd->linearAxis(); |
659 | < | int j = (i + 1) % 3; |
660 | < | int k = (i + 2) % 3; |
661 | < | value += angMom[j] * angMom[j] / I(j, j) + |
662 | < | angMom[k] * angMom[k] / I(k, k); |
663 | < | } else { |
664 | < | value += angMom[0]*angMom[0]/I(0, 0) |
665 | < | + angMom[1]*angMom[1]/I(1, 1) |
666 | < | + angMom[2]*angMom[2]/I(2, 2); |
667 | < | } |
668 | < | } //angular momenta exchange enabled |
669 | < | //energyConvert temporarily disabled |
670 | < | //make kineticExchange_ comparable between swap & scale |
671 | < | //value = value * 0.5 / PhysicalConstants::energyConvert; |
672 | < | value *= 0.5; |
673 | < | break; |
674 | < | case rnemdPx : |
675 | < | value = mass * vel[0]; |
676 | < | break; |
677 | < | case rnemdPy : |
678 | < | value = mass * vel[1]; |
679 | < | break; |
680 | < | case rnemdPz : |
681 | < | value = mass * vel[2]; |
682 | < | break; |
683 | < | default : |
684 | < | break; |
651 | > | if (sd->isLinear()) { |
652 | > | int i = sd->linearAxis(); |
653 | > | int j = (i + 1) % 3; |
654 | > | int k = (i + 2) % 3; |
655 | > | value += angMom[j] * angMom[j] / I(j, j) + |
656 | > | angMom[k] * angMom[k] / I(k, k); |
657 | > | } else { |
658 | > | value += angMom[0]*angMom[0]/I(0, 0) |
659 | > | + angMom[1]*angMom[1]/I(1, 1) |
660 | > | + angMom[2]*angMom[2]/I(2, 2); |
661 | > | } |
662 | > | } //angular momenta exchange enabled |
663 | > | value *= 0.5; |
664 | > | break; |
665 | > | case rnemdPx : |
666 | > | value = mass * vel[0]; |
667 | > | break; |
668 | > | case rnemdPy : |
669 | > | value = mass * vel[1]; |
670 | > | break; |
671 | > | case rnemdPz : |
672 | > | value = mass * vel[2]; |
673 | > | break; |
674 | > | default : |
675 | > | break; |
676 | > | } |
677 | > | if (!max_found) { |
678 | > | max_val = value; |
679 | > | max_sd = sd; |
680 | > | max_found = true; |
681 | > | } else { |
682 | > | if (max_val < value) { |
683 | > | max_val = value; |
684 | > | max_sd = sd; |
685 | } | |
686 | + | } |
687 | + | } |
688 | ||
689 | < | if (inA == 0) { |
690 | < | if (!min_found) { |
691 | < | min_val = value; |
692 | < | min_sd = sd; |
693 | < | min_found = true; |
694 | < | } else { |
695 | < | if (min_val > value) { |
696 | < | min_val = value; |
697 | < | min_sd = sd; |
698 | < | } |
699 | < | } |
700 | < | } else { |
701 | < | if (!max_found) { |
702 | < | max_val = value; |
703 | < | max_sd = sd; |
704 | < | max_found = true; |
705 | < | } else { |
706 | < | if (max_val < value) { |
707 | < | max_val = value; |
708 | < | max_sd = sd; |
709 | < | } |
710 | < | } |
711 | < | } |
689 | > | for (sd = seleManB_.beginSelected(selej); sd != NULL; |
690 | > | sd = seleManB_.nextSelected(selej)) { |
691 | > | |
692 | > | Vector3d pos = sd->getPos(); |
693 | > | |
694 | > | // wrap the stuntdouble's position back into the box: |
695 | > | |
696 | > | if (usePeriodicBoundaryConditions_) |
697 | > | currentSnap_->wrapVector(pos); |
698 | > | |
699 | > | RealType mass = sd->getMass(); |
700 | > | Vector3d vel = sd->getVel(); |
701 | > | RealType value; |
702 | > | |
703 | > | switch(rnemdFluxType_) { |
704 | > | case rnemdKE : |
705 | > | |
706 | > | value = mass * vel.lengthSquare(); |
707 | > | |
708 | > | if (sd->isDirectional()) { |
709 | > | Vector3d angMom = sd->getJ(); |
710 | > | Mat3x3d I = sd->getI(); |
711 | > | |
712 | > | if (sd->isLinear()) { |
713 | > | int i = sd->linearAxis(); |
714 | > | int j = (i + 1) % 3; |
715 | > | int k = (i + 2) % 3; |
716 | > | value += angMom[j] * angMom[j] / I(j, j) + |
717 | > | angMom[k] * angMom[k] / I(k, k); |
718 | > | } else { |
719 | > | value += angMom[0]*angMom[0]/I(0, 0) |
720 | > | + angMom[1]*angMom[1]/I(1, 1) |
721 | > | + angMom[2]*angMom[2]/I(2, 2); |
722 | > | } |
723 | > | } //angular momenta exchange enabled |
724 | > | value *= 0.5; |
725 | > | break; |
726 | > | case rnemdPx : |
727 | > | value = mass * vel[0]; |
728 | > | break; |
729 | > | case rnemdPy : |
730 | > | value = mass * vel[1]; |
731 | > | break; |
732 | > | case rnemdPz : |
733 | > | value = mass * vel[2]; |
734 | > | break; |
735 | > | default : |
736 | > | break; |
737 | > | } |
738 | > | |
739 | > | if (!min_found) { |
740 | > | min_val = value; |
741 | > | min_sd = sd; |
742 | > | min_found = true; |
743 | > | } else { |
744 | > | if (min_val > value) { |
745 | > | min_val = value; |
746 | > | min_sd = sd; |
747 | > | } |
748 | } | |
749 | } | |
750 | ||
751 | < | #ifdef IS_MPI |
752 | < | int nProc, worldRank; |
751 | > | #ifdef IS_MPI |
752 | > | int worldRank = MPI::COMM_WORLD.Get_rank(); |
753 | ||
538 | – | nProc = MPI::COMM_WORLD.Get_size(); |
539 | – | worldRank = MPI::COMM_WORLD.Get_rank(); |
540 | – | |
754 | bool my_min_found = min_found; | |
755 | bool my_max_found = max_found; | |
756 | ||
# | Line 728 | Line 941 | namespace OpenMD { | |
941 | ||
942 | switch(rnemdFluxType_) { | |
943 | case rnemdKE: | |
731 | – | cerr << "KE\n"; |
944 | kineticExchange_ += max_val - min_val; | |
945 | break; | |
946 | case rnemdPx: | |
# | Line 741 | Line 953 | namespace OpenMD { | |
953 | momentumExchange_.z() += max_val - min_val; | |
954 | break; | |
955 | default: | |
744 | – | cerr << "default\n"; |
956 | break; | |
957 | } | |
958 | } else { | |
# | Line 763 | Line 974 | namespace OpenMD { | |
974 | } | |
975 | } | |
976 | ||
977 | < | void RNEMD::doNIVS() { |
977 | > | void RNEMD::doNIVS(SelectionManager& smanA, SelectionManager& smanB) { |
978 | > | if (!doRNEMD_) return; |
979 | > | int selei; |
980 | > | int selej; |
981 | ||
982 | Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); | |
983 | + | RealType time = currentSnap_->getTime(); |
984 | Mat3x3d hmat = currentSnap_->getHmat(); | |
985 | ||
771 | – | seleMan_.setSelectionSet(evaluator_.evaluate()); |
772 | – | |
773 | – | int selei; |
986 | StuntDouble* sd; | |
775 | – | int idx; |
987 | ||
988 | vector<StuntDouble*> hotBin, coldBin; | |
989 | ||
# | Line 791 | Line 1002 | namespace OpenMD { | |
1002 | RealType Kcz = 0.0; | |
1003 | RealType Kcw = 0.0; | |
1004 | ||
1005 | < | for (sd = seleMan_.beginSelected(selei); sd != NULL; |
1006 | < | sd = seleMan_.nextSelected(selei)) { |
796 | < | |
797 | < | idx = sd->getLocalIndex(); |
798 | < | |
799 | < | Vector3d pos = sd->getPos(); |
1005 | > | for (sd = smanA.beginSelected(selei); sd != NULL; |
1006 | > | sd = smanA.nextSelected(selei)) { |
1007 | ||
1008 | + | Vector3d pos = sd->getPos(); |
1009 | + | |
1010 | // wrap the stuntdouble's position back into the box: | |
1011 | < | |
1011 | > | |
1012 | if (usePeriodicBoundaryConditions_) | |
1013 | currentSnap_->wrapVector(pos); | |
1014 | < | |
1015 | < | // which bin is this stuntdouble in? |
1016 | < | bool inA = inSlabA(pos); |
1017 | < | bool inB = inSlabB(pos); |
1014 | > | |
1015 | > | |
1016 | > | RealType mass = sd->getMass(); |
1017 | > | Vector3d vel = sd->getVel(); |
1018 | > | |
1019 | > | hotBin.push_back(sd); |
1020 | > | Phx += mass * vel.x(); |
1021 | > | Phy += mass * vel.y(); |
1022 | > | Phz += mass * vel.z(); |
1023 | > | Khx += mass * vel.x() * vel.x(); |
1024 | > | Khy += mass * vel.y() * vel.y(); |
1025 | > | Khz += mass * vel.z() * vel.z(); |
1026 | > | if (sd->isDirectional()) { |
1027 | > | Vector3d angMom = sd->getJ(); |
1028 | > | Mat3x3d I = sd->getI(); |
1029 | > | if (sd->isLinear()) { |
1030 | > | int i = sd->linearAxis(); |
1031 | > | int j = (i + 1) % 3; |
1032 | > | int k = (i + 2) % 3; |
1033 | > | Khw += angMom[j] * angMom[j] / I(j, j) + |
1034 | > | angMom[k] * angMom[k] / I(k, k); |
1035 | > | } else { |
1036 | > | Khw += angMom[0]*angMom[0]/I(0, 0) |
1037 | > | + angMom[1]*angMom[1]/I(1, 1) |
1038 | > | + angMom[2]*angMom[2]/I(2, 2); |
1039 | > | } |
1040 | > | } |
1041 | > | } |
1042 | > | for (sd = smanB.beginSelected(selej); sd != NULL; |
1043 | > | sd = smanB.nextSelected(selej)) { |
1044 | > | Vector3d pos = sd->getPos(); |
1045 | > | |
1046 | > | // wrap the stuntdouble's position back into the box: |
1047 | > | |
1048 | > | if (usePeriodicBoundaryConditions_) |
1049 | > | currentSnap_->wrapVector(pos); |
1050 | > | |
1051 | > | RealType mass = sd->getMass(); |
1052 | > | Vector3d vel = sd->getVel(); |
1053 | ||
1054 | < | if (inA || inB) { |
1055 | < | |
1056 | < | RealType mass = sd->getMass(); |
1057 | < | Vector3d vel = sd->getVel(); |
1058 | < | |
1059 | < | if (inA) { |
1060 | < | hotBin.push_back(sd); |
1061 | < | Phx += mass * vel.x(); |
1062 | < | Phy += mass * vel.y(); |
1063 | < | Phz += mass * vel.z(); |
1064 | < | Khx += mass * vel.x() * vel.x(); |
1065 | < | Khy += mass * vel.y() * vel.y(); |
1066 | < | Khz += mass * vel.z() * vel.z(); |
1067 | < | if (sd->isDirectional()) { |
1068 | < | Vector3d angMom = sd->getJ(); |
1069 | < | Mat3x3d I = sd->getI(); |
1070 | < | if (sd->isLinear()) { |
1071 | < | int i = sd->linearAxis(); |
1072 | < | int j = (i + 1) % 3; |
1073 | < | int k = (i + 2) % 3; |
1074 | < | Khw += angMom[j] * angMom[j] / I(j, j) + |
831 | < | angMom[k] * angMom[k] / I(k, k); |
832 | < | } else { |
833 | < | Khw += angMom[0]*angMom[0]/I(0, 0) |
834 | < | + angMom[1]*angMom[1]/I(1, 1) |
835 | < | + angMom[2]*angMom[2]/I(2, 2); |
836 | < | } |
837 | < | } |
838 | < | } else { |
839 | < | coldBin.push_back(sd); |
840 | < | Pcx += mass * vel.x(); |
841 | < | Pcy += mass * vel.y(); |
842 | < | Pcz += mass * vel.z(); |
843 | < | Kcx += mass * vel.x() * vel.x(); |
844 | < | Kcy += mass * vel.y() * vel.y(); |
845 | < | Kcz += mass * vel.z() * vel.z(); |
846 | < | if (sd->isDirectional()) { |
847 | < | Vector3d angMom = sd->getJ(); |
848 | < | Mat3x3d I = sd->getI(); |
849 | < | if (sd->isLinear()) { |
850 | < | int i = sd->linearAxis(); |
851 | < | int j = (i + 1) % 3; |
852 | < | int k = (i + 2) % 3; |
853 | < | Kcw += angMom[j] * angMom[j] / I(j, j) + |
854 | < | angMom[k] * angMom[k] / I(k, k); |
855 | < | } else { |
856 | < | Kcw += angMom[0]*angMom[0]/I(0, 0) |
857 | < | + angMom[1]*angMom[1]/I(1, 1) |
858 | < | + angMom[2]*angMom[2]/I(2, 2); |
859 | < | } |
860 | < | } |
861 | < | } |
1054 | > | coldBin.push_back(sd); |
1055 | > | Pcx += mass * vel.x(); |
1056 | > | Pcy += mass * vel.y(); |
1057 | > | Pcz += mass * vel.z(); |
1058 | > | Kcx += mass * vel.x() * vel.x(); |
1059 | > | Kcy += mass * vel.y() * vel.y(); |
1060 | > | Kcz += mass * vel.z() * vel.z(); |
1061 | > | if (sd->isDirectional()) { |
1062 | > | Vector3d angMom = sd->getJ(); |
1063 | > | Mat3x3d I = sd->getI(); |
1064 | > | if (sd->isLinear()) { |
1065 | > | int i = sd->linearAxis(); |
1066 | > | int j = (i + 1) % 3; |
1067 | > | int k = (i + 2) % 3; |
1068 | > | Kcw += angMom[j] * angMom[j] / I(j, j) + |
1069 | > | angMom[k] * angMom[k] / I(k, k); |
1070 | > | } else { |
1071 | > | Kcw += angMom[0]*angMom[0]/I(0, 0) |
1072 | > | + angMom[1]*angMom[1]/I(1, 1) |
1073 | > | + angMom[2]*angMom[2]/I(2, 2); |
1074 | > | } |
1075 | } | |
1076 | } | |
1077 | ||
# | Line 908 | Line 1121 | namespace OpenMD { | |
1121 | ||
1122 | if ((c > 0.81) && (c < 1.21)) {//restrict scaling coefficients | |
1123 | c = sqrt(c); | |
1124 | < | //std::cerr << "cold slab scaling coefficient: " << c << endl; |
912 | < | //now convert to hotBin coefficient |
1124 | > | |
1125 | RealType w = 0.0; | |
1126 | if (rnemdFluxType_ == rnemdFullKE) { | |
1127 | x = 1.0 + px * (1.0 - c); | |
# | Line 947 | Line 1159 | namespace OpenMD { | |
1159 | } | |
1160 | } | |
1161 | w = sqrt(w); | |
950 | – | // std::cerr << "xh= " << x << "\tyh= " << y << "\tzh= " << z |
951 | – | // << "\twh= " << w << endl; |
1162 | for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) { | |
1163 | if (rnemdFluxType_ == rnemdFullKE) { | |
1164 | vel = (*sdi)->getVel(); | |
# | Line 1211 | Line 1421 | namespace OpenMD { | |
1421 | failTrialCount_++; | |
1422 | } | |
1423 | } | |
1424 | < | |
1425 | < | void RNEMD::doVSS() { |
1424 | > | |
1425 | > | void RNEMD::doVSS(SelectionManager& smanA, SelectionManager& smanB) { |
1426 | > | if (!doRNEMD_) return; |
1427 | > | int selei; |
1428 | > | int selej; |
1429 | ||
1430 | Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); | |
1431 | RealType time = currentSnap_->getTime(); | |
1432 | Mat3x3d hmat = currentSnap_->getHmat(); | |
1433 | ||
1221 | – | seleMan_.setSelectionSet(evaluator_.evaluate()); |
1222 | – | |
1223 | – | int selei; |
1434 | StuntDouble* sd; | |
1225 | – | int idx; |
1435 | ||
1436 | vector<StuntDouble*> hotBin, coldBin; | |
1437 | ||
1438 | Vector3d Ph(V3Zero); | |
1439 | + | Vector3d Lh(V3Zero); |
1440 | RealType Mh = 0.0; | |
1441 | + | Mat3x3d Ih(0.0); |
1442 | RealType Kh = 0.0; | |
1443 | Vector3d Pc(V3Zero); | |
1444 | + | Vector3d Lc(V3Zero); |
1445 | RealType Mc = 0.0; | |
1446 | + | Mat3x3d Ic(0.0); |
1447 | RealType Kc = 0.0; | |
1448 | ||
1449 | + | for (sd = smanA.beginSelected(selei); sd != NULL; |
1450 | + | sd = smanA.nextSelected(selei)) { |
1451 | ||
1237 | – | for (sd = seleMan_.beginSelected(selei); sd != NULL; |
1238 | – | sd = seleMan_.nextSelected(selei)) { |
1239 | – | |
1240 | – | idx = sd->getLocalIndex(); |
1241 | – | |
1452 | Vector3d pos = sd->getPos(); | |
1453 | ||
1454 | // wrap the stuntdouble's position back into the box: | |
1455 | + | |
1456 | + | if (usePeriodicBoundaryConditions_) |
1457 | + | currentSnap_->wrapVector(pos); |
1458 | + | |
1459 | + | RealType mass = sd->getMass(); |
1460 | + | Vector3d vel = sd->getVel(); |
1461 | + | Vector3d rPos = sd->getPos() - coordinateOrigin_; |
1462 | + | RealType r2; |
1463 | + | |
1464 | + | hotBin.push_back(sd); |
1465 | + | Ph += mass * vel; |
1466 | + | Mh += mass; |
1467 | + | Kh += mass * vel.lengthSquare(); |
1468 | + | Lh += mass * cross(rPos, vel); |
1469 | + | Ih -= outProduct(rPos, rPos) * mass; |
1470 | + | r2 = rPos.lengthSquare(); |
1471 | + | Ih(0, 0) += mass * r2; |
1472 | + | Ih(1, 1) += mass * r2; |
1473 | + | Ih(2, 2) += mass * r2; |
1474 | + | |
1475 | + | if (rnemdFluxType_ == rnemdFullKE) { |
1476 | + | if (sd->isDirectional()) { |
1477 | + | Vector3d angMom = sd->getJ(); |
1478 | + | Mat3x3d I = sd->getI(); |
1479 | + | if (sd->isLinear()) { |
1480 | + | int i = sd->linearAxis(); |
1481 | + | int j = (i + 1) % 3; |
1482 | + | int k = (i + 2) % 3; |
1483 | + | Kh += angMom[j] * angMom[j] / I(j, j) + |
1484 | + | angMom[k] * angMom[k] / I(k, k); |
1485 | + | } else { |
1486 | + | Kh += angMom[0] * angMom[0] / I(0, 0) + |
1487 | + | angMom[1] * angMom[1] / I(1, 1) + |
1488 | + | angMom[2] * angMom[2] / I(2, 2); |
1489 | + | } |
1490 | + | } |
1491 | + | } |
1492 | + | } |
1493 | + | for (sd = smanB.beginSelected(selej); sd != NULL; |
1494 | + | sd = smanB.nextSelected(selej)) { |
1495 | ||
1496 | + | Vector3d pos = sd->getPos(); |
1497 | + | |
1498 | + | // wrap the stuntdouble's position back into the box: |
1499 | + | |
1500 | if (usePeriodicBoundaryConditions_) | |
1501 | currentSnap_->wrapVector(pos); | |
1502 | + | |
1503 | + | RealType mass = sd->getMass(); |
1504 | + | Vector3d vel = sd->getVel(); |
1505 | + | Vector3d rPos = sd->getPos() - coordinateOrigin_; |
1506 | + | RealType r2; |
1507 | ||
1508 | < | // which bin is this stuntdouble in? |
1509 | < | bool inA = inSlabA(pos); |
1510 | < | bool inB = inSlabB(pos); |
1508 | > | coldBin.push_back(sd); |
1509 | > | Pc += mass * vel; |
1510 | > | Mc += mass; |
1511 | > | Kc += mass * vel.lengthSquare(); |
1512 | > | Lc += mass * cross(rPos, vel); |
1513 | > | Ic -= outProduct(rPos, rPos) * mass; |
1514 | > | r2 = rPos.lengthSquare(); |
1515 | > | Ic(0, 0) += mass * r2; |
1516 | > | Ic(1, 1) += mass * r2; |
1517 | > | Ic(2, 2) += mass * r2; |
1518 | ||
1519 | < | if (inA || inB) { |
1520 | < | |
1521 | < | RealType mass = sd->getMass(); |
1522 | < | Vector3d vel = sd->getVel(); |
1523 | < | |
1524 | < | if (inA) { |
1525 | < | hotBin.push_back(sd); |
1526 | < | //std::cerr << "before, velocity = " << vel << endl; |
1527 | < | Ph += mass * vel; |
1528 | < | //std::cerr << "after, velocity = " << vel << endl; |
1529 | < | Mh += mass; |
1530 | < | Kh += mass * vel.lengthSquare(); |
1531 | < | if (rnemdFluxType_ == rnemdFullKE) { |
1532 | < | if (sd->isDirectional()) { |
1533 | < | Vector3d angMom = sd->getJ(); |
1534 | < | Mat3x3d I = sd->getI(); |
1269 | < | if (sd->isLinear()) { |
1270 | < | int i = sd->linearAxis(); |
1271 | < | int j = (i + 1) % 3; |
1272 | < | int k = (i + 2) % 3; |
1273 | < | Kh += angMom[j] * angMom[j] / I(j, j) + |
1274 | < | angMom[k] * angMom[k] / I(k, k); |
1275 | < | } else { |
1276 | < | Kh += angMom[0] * angMom[0] / I(0, 0) + |
1277 | < | angMom[1] * angMom[1] / I(1, 1) + |
1278 | < | angMom[2] * angMom[2] / I(2, 2); |
1279 | < | } |
1280 | < | } |
1281 | < | } |
1282 | < | } else { //midBin_ |
1283 | < | coldBin.push_back(sd); |
1284 | < | Pc += mass * vel; |
1285 | < | Mc += mass; |
1286 | < | Kc += mass * vel.lengthSquare(); |
1287 | < | if (rnemdFluxType_ == rnemdFullKE) { |
1288 | < | if (sd->isDirectional()) { |
1289 | < | Vector3d angMom = sd->getJ(); |
1290 | < | Mat3x3d I = sd->getI(); |
1291 | < | if (sd->isLinear()) { |
1292 | < | int i = sd->linearAxis(); |
1293 | < | int j = (i + 1) % 3; |
1294 | < | int k = (i + 2) % 3; |
1295 | < | Kc += angMom[j] * angMom[j] / I(j, j) + |
1296 | < | angMom[k] * angMom[k] / I(k, k); |
1297 | < | } else { |
1298 | < | Kc += angMom[0] * angMom[0] / I(0, 0) + |
1299 | < | angMom[1] * angMom[1] / I(1, 1) + |
1300 | < | angMom[2] * angMom[2] / I(2, 2); |
1301 | < | } |
1302 | < | } |
1303 | < | } |
1304 | < | } |
1519 | > | if (rnemdFluxType_ == rnemdFullKE) { |
1520 | > | if (sd->isDirectional()) { |
1521 | > | Vector3d angMom = sd->getJ(); |
1522 | > | Mat3x3d I = sd->getI(); |
1523 | > | if (sd->isLinear()) { |
1524 | > | int i = sd->linearAxis(); |
1525 | > | int j = (i + 1) % 3; |
1526 | > | int k = (i + 2) % 3; |
1527 | > | Kc += angMom[j] * angMom[j] / I(j, j) + |
1528 | > | angMom[k] * angMom[k] / I(k, k); |
1529 | > | } else { |
1530 | > | Kc += angMom[0] * angMom[0] / I(0, 0) + |
1531 | > | angMom[1] * angMom[1] / I(1, 1) + |
1532 | > | angMom[2] * angMom[2] / I(2, 2); |
1533 | > | } |
1534 | > | } |
1535 | } | |
1536 | } | |
1537 | ||
1538 | Kh *= 0.5; | |
1539 | Kc *= 0.5; | |
1310 | – | |
1311 | – | // std::cerr << "Mh= " << Mh << "\tKh= " << Kh << "\tMc= " << Mc |
1312 | – | // << "\tKc= " << Kc << endl; |
1313 | – | // std::cerr << "Ph= " << Ph << "\tPc= " << Pc << endl; |
1540 | ||
1541 | #ifdef IS_MPI | |
1542 | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Ph[0], 3, MPI::REALTYPE, MPI::SUM); | |
1543 | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pc[0], 3, MPI::REALTYPE, MPI::SUM); | |
1544 | + | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Lh[0], 3, MPI::REALTYPE, MPI::SUM); |
1545 | + | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Lc[0], 3, MPI::REALTYPE, MPI::SUM); |
1546 | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mh, 1, MPI::REALTYPE, MPI::SUM); | |
1547 | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kh, 1, MPI::REALTYPE, MPI::SUM); | |
1548 | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mc, 1, MPI::REALTYPE, MPI::SUM); | |
1549 | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kc, 1, MPI::REALTYPE, MPI::SUM); | |
1550 | + | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, Ih.getArrayPointer(), 9, |
1551 | + | MPI::REALTYPE, MPI::SUM); |
1552 | + | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, Ic.getArrayPointer(), 9, |
1553 | + | MPI::REALTYPE, MPI::SUM); |
1554 | #endif | |
1555 | < | |
1555 | > | |
1556 | bool successfulExchange = false; | |
1557 | if ((Mh > 0.0) && (Mc > 0.0)) {//both slabs are not empty | |
1558 | Vector3d vc = Pc / Mc; | |
1559 | Vector3d ac = -momentumTarget_ / Mc + vc; | |
1560 | Vector3d acrec = -momentumTarget_ / Mc; | |
1561 | < | RealType cNumerator = Kc - kineticTarget_ - 0.5 * Mc * ac.lengthSquare(); |
1561 | > | |
1562 | > | // We now need the inverse of the inertia tensor to calculate the |
1563 | > | // angular velocity of the cold slab; |
1564 | > | Mat3x3d Ici = Ic.inverse(); |
1565 | > | Vector3d omegac = Ici * Lc; |
1566 | > | Vector3d bc = -(Ici * angularMomentumTarget_) + omegac; |
1567 | > | Vector3d bcrec = bc - omegac; |
1568 | > | |
1569 | > | RealType cNumerator = Kc - kineticTarget_ |
1570 | > | - 0.5 * Mc * ac.lengthSquare() - 0.5 * ( dot(bc, Ic * bc)); |
1571 | if (cNumerator > 0.0) { | |
1572 | < | RealType cDenominator = Kc - 0.5 * Mc * vc.lengthSquare(); |
1572 | > | |
1573 | > | RealType cDenominator = Kc - 0.5 * Mc * vc.lengthSquare() |
1574 | > | - 0.5*(dot(omegac, Ic * omegac)); |
1575 | > | |
1576 | if (cDenominator > 0.0) { | |
1577 | RealType c = sqrt(cNumerator / cDenominator); | |
1578 | if ((c > 0.9) && (c < 1.1)) {//restrict scaling coefficients | |
1579 | + | |
1580 | Vector3d vh = Ph / Mh; | |
1581 | Vector3d ah = momentumTarget_ / Mh + vh; | |
1582 | Vector3d ahrec = momentumTarget_ / Mh; | |
1583 | < | RealType hNumerator = Kh + kineticTarget_ |
1584 | < | - 0.5 * Mh * ah.lengthSquare(); |
1585 | < | if (hNumerator > 0.0) { |
1586 | < | RealType hDenominator = Kh - 0.5 * Mh * vh.lengthSquare(); |
1583 | > | |
1584 | > | // We now need the inverse of the inertia tensor to |
1585 | > | // calculate the angular velocity of the hot slab; |
1586 | > | Mat3x3d Ihi = Ih.inverse(); |
1587 | > | Vector3d omegah = Ihi * Lh; |
1588 | > | Vector3d bh = (Ihi * angularMomentumTarget_) + omegah; |
1589 | > | Vector3d bhrec = bh - omegah; |
1590 | > | |
1591 | > | RealType hNumerator = Kh + kineticTarget_ |
1592 | > | - 0.5 * Mh * ah.lengthSquare() - 0.5 * ( dot(bh, Ih * bh));; |
1593 | > | if (hNumerator > 0.0) { |
1594 | > | |
1595 | > | RealType hDenominator = Kh - 0.5 * Mh * vh.lengthSquare() |
1596 | > | - 0.5*(dot(omegah, Ih * omegah)); |
1597 | > | |
1598 | if (hDenominator > 0.0) { | |
1599 | RealType h = sqrt(hNumerator / hDenominator); | |
1600 | if ((h > 0.9) && (h < 1.1)) { | |
1601 | < | // std::cerr << "cold slab scaling coefficient: " << c << "\n"; |
1346 | < | // std::cerr << "hot slab scaling coefficient: " << h << "\n"; |
1601 | > | |
1602 | vector<StuntDouble*>::iterator sdi; | |
1603 | Vector3d vel; | |
1604 | + | Vector3d rPos; |
1605 | + | |
1606 | for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) { | |
1607 | //vel = (*sdi)->getVel(); | |
1608 | < | vel = ((*sdi)->getVel() - vc) * c + ac; |
1608 | > | rPos = (*sdi)->getPos() - coordinateOrigin_; |
1609 | > | vel = ((*sdi)->getVel() - vc - cross(omegac, rPos)) * c |
1610 | > | + ac + cross(bc, rPos); |
1611 | (*sdi)->setVel(vel); | |
1612 | if (rnemdFluxType_ == rnemdFullKE) { | |
1613 | if ((*sdi)->isDirectional()) { | |
# | Line 1359 | Line 1618 | namespace OpenMD { | |
1618 | } | |
1619 | for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) { | |
1620 | //vel = (*sdi)->getVel(); | |
1621 | < | vel = ((*sdi)->getVel() - vh) * h + ah; |
1621 | > | rPos = (*sdi)->getPos() - coordinateOrigin_; |
1622 | > | vel = ((*sdi)->getVel() - vh - cross(omegah, rPos)) * h |
1623 | > | + ah + cross(bh, rPos); |
1624 | > | cerr << "setting vel to " << vel << "\n"; |
1625 | (*sdi)->setVel(vel); | |
1626 | if (rnemdFluxType_ == rnemdFullKE) { | |
1627 | if ((*sdi)->isDirectional()) { | |
# | Line 1371 | Line 1633 | namespace OpenMD { | |
1633 | successfulExchange = true; | |
1634 | kineticExchange_ += kineticTarget_; | |
1635 | momentumExchange_ += momentumTarget_; | |
1636 | + | angularMomentumExchange_ += angularMomentumTarget_; |
1637 | } | |
1638 | } | |
1639 | } | |
# | Line 1390 | Line 1653 | namespace OpenMD { | |
1653 | } | |
1654 | } | |
1655 | ||
1656 | < | void RNEMD::doRNEMD() { |
1656 | > | RealType RNEMD::getDividingArea() { |
1657 | ||
1658 | + | if (hasDividingArea_) return dividingArea_; |
1659 | + | |
1660 | + | RealType areaA, areaB; |
1661 | + | Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); |
1662 | + | |
1663 | + | if (hasSelectionA_) { |
1664 | + | int isd; |
1665 | + | StuntDouble* sd; |
1666 | + | vector<StuntDouble*> aSites; |
1667 | + | ConvexHull* surfaceMeshA = new ConvexHull(); |
1668 | + | seleManA_.setSelectionSet(evaluatorA_.evaluate()); |
1669 | + | for (sd = seleManA_.beginSelected(isd); sd != NULL; |
1670 | + | sd = seleManA_.nextSelected(isd)) { |
1671 | + | aSites.push_back(sd); |
1672 | + | } |
1673 | + | surfaceMeshA->computeHull(aSites); |
1674 | + | areaA = surfaceMeshA->getArea(); |
1675 | + | } else { |
1676 | + | if (usePeriodicBoundaryConditions_) { |
1677 | + | // in periodic boundaries, the surface area is twice the x-y |
1678 | + | // area of the current box: |
1679 | + | areaA = 2.0 * snap->getXYarea(); |
1680 | + | } else { |
1681 | + | // in non-periodic simulations, without explicitly setting |
1682 | + | // selections, the sphere radius sets the surface area of the |
1683 | + | // dividing surface: |
1684 | + | areaA = 4.0 * M_PI * pow(sphereARadius_, 2); |
1685 | + | } |
1686 | + | } |
1687 | + | |
1688 | + | if (hasSelectionB_) { |
1689 | + | int isd; |
1690 | + | StuntDouble* sd; |
1691 | + | vector<StuntDouble*> bSites; |
1692 | + | ConvexHull* surfaceMeshB = new ConvexHull(); |
1693 | + | seleManB_.setSelectionSet(evaluatorB_.evaluate()); |
1694 | + | for (sd = seleManB_.beginSelected(isd); sd != NULL; |
1695 | + | sd = seleManB_.nextSelected(isd)) { |
1696 | + | bSites.push_back(sd); |
1697 | + | } |
1698 | + | surfaceMeshB->computeHull(bSites); |
1699 | + | areaB = surfaceMeshB->getArea(); |
1700 | + | } else { |
1701 | + | if (usePeriodicBoundaryConditions_) { |
1702 | + | // in periodic boundaries, the surface area is twice the x-y |
1703 | + | // area of the current box: |
1704 | + | areaB = 2.0 * snap->getXYarea(); |
1705 | + | } else { |
1706 | + | // in non-periodic simulations, without explicitly setting |
1707 | + | // selections, but if a sphereBradius has been set, just use that: |
1708 | + | areaB = 4.0 * M_PI * pow(sphereBRadius_, 2); |
1709 | + | } |
1710 | + | } |
1711 | + | |
1712 | + | dividingArea_ = min(areaA, areaB); |
1713 | + | hasDividingArea_ = true; |
1714 | + | return dividingArea_; |
1715 | + | } |
1716 | + | |
1717 | + | void RNEMD::doRNEMD() { |
1718 | + | if (!doRNEMD_) return; |
1719 | trialCount_++; | |
1720 | + | |
1721 | + | cerr << "trialCount = " << trialCount_ << "\n"; |
1722 | + | // object evaluator: |
1723 | + | evaluator_.loadScriptString(rnemdObjectSelection_); |
1724 | + | seleMan_.setSelectionSet(evaluator_.evaluate()); |
1725 | + | |
1726 | + | evaluatorA_.loadScriptString(selectionA_); |
1727 | + | evaluatorB_.loadScriptString(selectionB_); |
1728 | + | |
1729 | + | seleManA_.setSelectionSet(evaluatorA_.evaluate()); |
1730 | + | seleManB_.setSelectionSet(evaluatorB_.evaluate()); |
1731 | + | |
1732 | + | commonA_ = seleManA_ & seleMan_; |
1733 | + | commonB_ = seleManB_ & seleMan_; |
1734 | + | |
1735 | + | // Target exchange quantities (in each exchange) = dividingArea * dt * flux |
1736 | + | // dt = exchange time interval |
1737 | + | // flux = target flux |
1738 | + | // dividingArea = smallest dividing surface between the two regions |
1739 | + | |
1740 | + | hasDividingArea_ = false; |
1741 | + | RealType area = getDividingArea(); |
1742 | + | |
1743 | + | kineticTarget_ = kineticFlux_ * exchangeTime_ * area; |
1744 | + | momentumTarget_ = momentumFluxVector_ * exchangeTime_ * area; |
1745 | + | angularMomentumTarget_ = angularMomentumFluxVector_ * exchangeTime_ * area; |
1746 | + | |
1747 | switch(rnemdMethod_) { | |
1748 | case rnemdSwap: | |
1749 | < | doSwap(); |
1749 | > | doSwap(commonA_, commonB_); |
1750 | break; | |
1751 | case rnemdNIVS: | |
1752 | < | doNIVS(); |
1752 | > | doNIVS(commonA_, commonB_); |
1753 | break; | |
1754 | case rnemdVSS: | |
1755 | < | doVSS(); |
1755 | > | doVSS(commonA_, commonB_); |
1756 | break; | |
1757 | case rnemdUnkownMethod: | |
1758 | default : | |
# | Line 1410 | Line 1761 | namespace OpenMD { | |
1761 | } | |
1762 | ||
1763 | void RNEMD::collectData() { | |
1764 | < | |
1764 | > | if (!doRNEMD_) return; |
1765 | Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); | |
1766 | + | |
1767 | + | cerr << "collecting data\n"; |
1768 | + | // collectData can be called more frequently than the doRNEMD, so use the |
1769 | + | // computed area from the last exchange time: |
1770 | + | RealType area = getDividingArea(); |
1771 | + | areaAccumulator_->add(area); |
1772 | Mat3x3d hmat = currentSnap_->getHmat(); | |
1416 | – | |
1773 | seleMan_.setSelectionSet(evaluator_.evaluate()); | |
1774 | ||
1775 | < | int selei; |
1775 | > | int selei(0); |
1776 | StuntDouble* sd; | |
1777 | < | int idx; |
1777 | > | int binNo; |
1778 | ||
1779 | vector<RealType> binMass(nBins_, 0.0); | |
1780 | vector<RealType> binPx(nBins_, 0.0); | |
1781 | vector<RealType> binPy(nBins_, 0.0); | |
1782 | vector<RealType> binPz(nBins_, 0.0); | |
1783 | + | vector<RealType> binOmegax(nBins_, 0.0); |
1784 | + | vector<RealType> binOmegay(nBins_, 0.0); |
1785 | + | vector<RealType> binOmegaz(nBins_, 0.0); |
1786 | vector<RealType> binKE(nBins_, 0.0); | |
1787 | vector<int> binDOF(nBins_, 0); | |
1788 | vector<int> binCount(nBins_, 0); | |
# | Line 1431 | Line 1790 | namespace OpenMD { | |
1790 | // alternative approach, track all molecules instead of only those | |
1791 | // selected for scaling/swapping: | |
1792 | /* | |
1793 | < | SimInfo::MoleculeIterator miter; |
1794 | < | vector<StuntDouble*>::iterator iiter; |
1795 | < | Molecule* mol; |
1796 | < | StuntDouble* sd; |
1797 | < | for (mol = info_->beginMolecule(miter); mol != NULL; |
1793 | > | SimInfo::MoleculeIterator miter; |
1794 | > | vector<StuntDouble*>::iterator iiter; |
1795 | > | Molecule* mol; |
1796 | > | StuntDouble* sd; |
1797 | > | for (mol = info_->beginMolecule(miter); mol != NULL; |
1798 | mol = info_->nextMolecule(miter)) | |
1799 | sd is essentially sd | |
1800 | < | for (sd = mol->beginIntegrableObject(iiter); |
1801 | < | sd != NULL; |
1802 | < | sd = mol->nextIntegrableObject(iiter)) |
1800 | > | for (sd = mol->beginIntegrableObject(iiter); |
1801 | > | sd != NULL; |
1802 | > | sd = mol->nextIntegrableObject(iiter)) |
1803 | */ | |
1804 | + | |
1805 | for (sd = seleMan_.beginSelected(selei); sd != NULL; | |
1806 | < | sd = seleMan_.nextSelected(selei)) { |
1807 | < | |
1448 | < | idx = sd->getLocalIndex(); |
1449 | < | |
1806 | > | sd = seleMan_.nextSelected(selei)) { |
1807 | > | |
1808 | Vector3d pos = sd->getPos(); | |
1809 | ||
1810 | // wrap the stuntdouble's position back into the box: | |
1811 | ||
1812 | < | if (usePeriodicBoundaryConditions_) |
1812 | > | if (usePeriodicBoundaryConditions_) { |
1813 | currentSnap_->wrapVector(pos); | |
1814 | + | // which bin is this stuntdouble in? |
1815 | + | // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)] |
1816 | + | // Shift molecules by half a box to have bins start at 0 |
1817 | + | // The modulo operator is used to wrap the case when we are |
1818 | + | // beyond the end of the bins back to the beginning. |
1819 | + | binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_; |
1820 | + | } else { |
1821 | + | Vector3d rPos = pos - coordinateOrigin_; |
1822 | + | binNo = int(rPos.length() / binWidth_); |
1823 | + | } |
1824 | ||
1457 | – | // which bin is this stuntdouble in? |
1458 | – | // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)] |
1459 | – | // Shift molecules by half a box to have bins start at 0 |
1460 | – | // The modulo operator is used to wrap the case when we are |
1461 | – | // beyond the end of the bins back to the beginning. |
1462 | – | int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_; |
1463 | – | |
1825 | RealType mass = sd->getMass(); | |
1826 | Vector3d vel = sd->getVel(); | |
1827 | < | |
1828 | < | binCount[binNo]++; |
1829 | < | binMass[binNo] += mass; |
1830 | < | binPx[binNo] += mass*vel.x(); |
1831 | < | binPy[binNo] += mass*vel.y(); |
1832 | < | binPz[binNo] += mass*vel.z(); |
1833 | < | binKE[binNo] += 0.5 * (mass * vel.lengthSquare()); |
1834 | < | binDOF[binNo] += 3; |
1835 | < | |
1836 | < | if (sd->isDirectional()) { |
1837 | < | Vector3d angMom = sd->getJ(); |
1838 | < | Mat3x3d I = sd->getI(); |
1839 | < | if (sd->isLinear()) { |
1840 | < | int i = sd->linearAxis(); |
1841 | < | int j = (i + 1) % 3; |
1842 | < | int k = (i + 2) % 3; |
1843 | < | binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / I(j, j) + |
1844 | < | angMom[k] * angMom[k] / I(k, k)); |
1845 | < | binDOF[binNo] += 2; |
1846 | < | } else { |
1847 | < | binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) + |
1848 | < | angMom[1] * angMom[1] / I(1, 1) + |
1849 | < | angMom[2] * angMom[2] / I(2, 2)); |
1850 | < | binDOF[binNo] += 3; |
1827 | > | Vector3d rPos = sd->getPos() - coordinateOrigin_; |
1828 | > | Vector3d aVel = cross(rPos, vel); |
1829 | > | |
1830 | > | if (binNo >= 0 && binNo < nBins_) { |
1831 | > | binCount[binNo]++; |
1832 | > | binMass[binNo] += mass; |
1833 | > | binPx[binNo] += mass*vel.x(); |
1834 | > | binPy[binNo] += mass*vel.y(); |
1835 | > | binPz[binNo] += mass*vel.z(); |
1836 | > | binOmegax[binNo] += aVel.x(); |
1837 | > | binOmegay[binNo] += aVel.y(); |
1838 | > | binOmegaz[binNo] += aVel.z(); |
1839 | > | binKE[binNo] += 0.5 * (mass * vel.lengthSquare()); |
1840 | > | binDOF[binNo] += 3; |
1841 | > | |
1842 | > | if (sd->isDirectional()) { |
1843 | > | Vector3d angMom = sd->getJ(); |
1844 | > | Mat3x3d I = sd->getI(); |
1845 | > | if (sd->isLinear()) { |
1846 | > | int i = sd->linearAxis(); |
1847 | > | int j = (i + 1) % 3; |
1848 | > | int k = (i + 2) % 3; |
1849 | > | binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / I(j, j) + |
1850 | > | angMom[k] * angMom[k] / I(k, k)); |
1851 | > | binDOF[binNo] += 2; |
1852 | > | } else { |
1853 | > | binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) + |
1854 | > | angMom[1] * angMom[1] / I(1, 1) + |
1855 | > | angMom[2] * angMom[2] / I(2, 2)); |
1856 | > | binDOF[binNo] += 3; |
1857 | > | } |
1858 | } | |
1859 | } | |
1860 | } | |
1861 | ||
1494 | – | |
1862 | #ifdef IS_MPI | |
1863 | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0], | |
1864 | nBins_, MPI::INT, MPI::SUM); | |
# | Line 1503 | Line 1870 | namespace OpenMD { | |
1870 | nBins_, MPI::REALTYPE, MPI::SUM); | |
1871 | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPz[0], | |
1872 | nBins_, MPI::REALTYPE, MPI::SUM); | |
1873 | + | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegax[0], |
1874 | + | nBins_, MPI::REALTYPE, MPI::SUM); |
1875 | + | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegay[0], |
1876 | + | nBins_, MPI::REALTYPE, MPI::SUM); |
1877 | + | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegaz[0], |
1878 | + | nBins_, MPI::REALTYPE, MPI::SUM); |
1879 | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binKE[0], | |
1880 | nBins_, MPI::REALTYPE, MPI::SUM); | |
1881 | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binDOF[0], | |
# | Line 1510 | Line 1883 | namespace OpenMD { | |
1883 | #endif | |
1884 | ||
1885 | Vector3d vel; | |
1886 | + | Vector3d aVel; |
1887 | RealType den; | |
1888 | RealType temp; | |
1889 | RealType z; | |
1890 | + | RealType r; |
1891 | for (int i = 0; i < nBins_; i++) { | |
1892 | < | z = (((RealType)i + 0.5) / (RealType)nBins_) * hmat(2,2); |
1892 | > | if (usePeriodicBoundaryConditions_) { |
1893 | > | z = (((RealType)i + 0.5) / (RealType)nBins_) * hmat(2,2); |
1894 | > | den = binMass[i] * nBins_ * PhysicalConstants::densityConvert |
1895 | > | / currentSnap_->getVolume() ; |
1896 | > | } else { |
1897 | > | r = (((RealType)i + 0.5) * binWidth_); |
1898 | > | RealType rinner = (RealType)i * binWidth_; |
1899 | > | RealType router = (RealType)(i+1) * binWidth_; |
1900 | > | den = binMass[i] * 3.0 * PhysicalConstants::densityConvert |
1901 | > | / (4.0 * M_PI * (pow(router,3) - pow(rinner,3))); |
1902 | > | } |
1903 | vel.x() = binPx[i] / binMass[i]; | |
1904 | vel.y() = binPy[i] / binMass[i]; | |
1905 | vel.z() = binPz[i] / binMass[i]; | |
1906 | < | den = binCount[i] * nBins_ / (hmat(0,0) * hmat(1,1) * hmat(2,2)); |
1907 | < | temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb * |
1908 | < | PhysicalConstants::energyConvert); |
1906 | > | aVel.x() = binOmegax[i]; |
1907 | > | aVel.y() = binOmegay[i]; |
1908 | > | aVel.z() = binOmegaz[i]; |
1909 | ||
1910 | < | for (unsigned int j = 0; j < outputMask_.size(); ++j) { |
1911 | < | if(outputMask_[j]) { |
1912 | < | switch(j) { |
1913 | < | case Z: |
1914 | < | (data_[j].accumulator[i])->add(z); |
1915 | < | break; |
1916 | < | case TEMPERATURE: |
1917 | < | data_[j].accumulator[i]->add(temp); |
1918 | < | break; |
1919 | < | case VELOCITY: |
1920 | < | dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel); |
1921 | < | break; |
1922 | < | case DENSITY: |
1923 | < | data_[j].accumulator[i]->add(den); |
1924 | < | break; |
1910 | > | if (binCount[i] > 0) { |
1911 | > | // only add values if there are things to add |
1912 | > | temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb * |
1913 | > | PhysicalConstants::energyConvert); |
1914 | > | |
1915 | > | for (unsigned int j = 0; j < outputMask_.size(); ++j) { |
1916 | > | if(outputMask_[j]) { |
1917 | > | switch(j) { |
1918 | > | case Z: |
1919 | > | dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(z); |
1920 | > | break; |
1921 | > | case R: |
1922 | > | dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(r); |
1923 | > | break; |
1924 | > | case TEMPERATURE: |
1925 | > | dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(temp); |
1926 | > | break; |
1927 | > | case VELOCITY: |
1928 | > | dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel); |
1929 | > | break; |
1930 | > | case ANGULARVELOCITY: |
1931 | > | dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(aVel); |
1932 | > | break; |
1933 | > | case DENSITY: |
1934 | > | dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den); |
1935 | > | break; |
1936 | > | } |
1937 | } | |
1938 | } | |
1939 | } | |
# | Line 1544 | Line 1941 | namespace OpenMD { | |
1941 | } | |
1942 | ||
1943 | void RNEMD::getStarted() { | |
1944 | + | if (!doRNEMD_) return; |
1945 | + | hasDividingArea_ = false; |
1946 | collectData(); | |
1947 | writeOutputFile(); | |
1948 | } | |
1949 | ||
1950 | void RNEMD::parseOutputFileFormat(const std::string& format) { | |
1951 | + | if (!doRNEMD_) return; |
1952 | StringTokenizer tokenizer(format, " ,;|\t\n\r"); | |
1953 | ||
1954 | while(tokenizer.hasMoreTokens()) { | |
# | Line 1569 | Line 1969 | namespace OpenMD { | |
1969 | } | |
1970 | ||
1971 | void RNEMD::writeOutputFile() { | |
1972 | + | if (!doRNEMD_) return; |
1973 | ||
1974 | #ifdef IS_MPI | |
1975 | // If we're the root node, should we print out the results | |
# | Line 1588 | Line 1989 | namespace OpenMD { | |
1989 | Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); | |
1990 | ||
1991 | RealType time = currentSnap_->getTime(); | |
1992 | < | |
1993 | < | |
1992 | > | RealType avgArea; |
1993 | > | areaAccumulator_->getAverage(avgArea); |
1994 | > | RealType Jz = kineticExchange_ / (time * avgArea) |
1995 | > | / PhysicalConstants::energyConvert; |
1996 | > | Vector3d JzP = momentumExchange_ / (time * avgArea); |
1997 | > | Vector3d JzL = angularMomentumExchange_ / (time * avgArea); |
1998 | > | |
1999 | rnemdFile_ << "#######################################################\n"; | |
2000 | rnemdFile_ << "# RNEMD {\n"; | |
2001 | ||
2002 | map<string, RNEMDMethod>::iterator mi; | |
2003 | for(mi = stringToMethod_.begin(); mi != stringToMethod_.end(); ++mi) { | |
2004 | if ( (*mi).second == rnemdMethod_) | |
2005 | < | rnemdFile_ << "# exchangeMethod = " << (*mi).first << "\n"; |
2005 | > | rnemdFile_ << "# exchangeMethod = \"" << (*mi).first << "\";\n"; |
2006 | } | |
2007 | map<string, RNEMDFluxType>::iterator fi; | |
2008 | for(fi = stringToFluxType_.begin(); fi != stringToFluxType_.end(); ++fi) { | |
2009 | if ( (*fi).second == rnemdFluxType_) | |
2010 | < | rnemdFile_ << "# fluxType = " << (*fi).first << "\n"; |
2010 | > | rnemdFile_ << "# fluxType = \"" << (*fi).first << "\";\n"; |
2011 | } | |
2012 | ||
2013 | < | rnemdFile_ << "# exchangeTime = " << exchangeTime_ << " fs\n"; |
2013 | > | rnemdFile_ << "# exchangeTime = " << exchangeTime_ << ";\n"; |
2014 | ||
2015 | rnemdFile_ << "# objectSelection = \"" | |
2016 | < | << rnemdObjectSelection_ << "\"\n"; |
2017 | < | rnemdFile_ << "# slabWidth = " << slabWidth_ << " angstroms\n"; |
2018 | < | rnemdFile_ << "# slabAcenter = " << slabACenter_ << " angstroms\n"; |
1613 | < | rnemdFile_ << "# slabBcenter = " << slabBCenter_ << " angstroms\n"; |
2016 | > | << rnemdObjectSelection_ << "\";\n"; |
2017 | > | rnemdFile_ << "# selectionA = \"" << selectionA_ << "\";\n"; |
2018 | > | rnemdFile_ << "# selectionB = \"" << selectionB_ << "\";\n"; |
2019 | rnemdFile_ << "# }\n"; | |
2020 | rnemdFile_ << "#######################################################\n"; | |
2021 | < | |
2022 | < | rnemdFile_ << "# running time = " << time << " fs\n"; |
2023 | < | rnemdFile_ << "# target kinetic flux = " << kineticFlux_ << "\n"; |
2024 | < | rnemdFile_ << "# target momentum flux = " << momentumFluxVector_ << "\n"; |
2025 | < | |
2026 | < | rnemdFile_ << "# target one-time kinetic exchange = " << kineticTarget_ |
2027 | < | << "\n"; |
2028 | < | rnemdFile_ << "# target one-time momentum exchange = " << momentumTarget_ |
2029 | < | << "\n"; |
2030 | < | |
2031 | < | rnemdFile_ << "# actual kinetic exchange = " << kineticExchange_ << "\n"; |
2032 | < | rnemdFile_ << "# actual momentum exchange = " << momentumExchange_ |
2033 | < | << "\n"; |
2034 | < | |
2035 | < | rnemdFile_ << "# attempted exchanges: " << trialCount_ << "\n"; |
2036 | < | rnemdFile_ << "# failed exchanges: " << failTrialCount_ << "\n"; |
2037 | < | |
2038 | < | |
2021 | > | rnemdFile_ << "# RNEMD report:\n"; |
2022 | > | rnemdFile_ << "# running time = " << time << " fs\n"; |
2023 | > | rnemdFile_ << "# Target flux:\n"; |
2024 | > | rnemdFile_ << "# kinetic = " |
2025 | > | << kineticFlux_ / PhysicalConstants::energyConvert |
2026 | > | << " (kcal/mol/A^2/fs)\n"; |
2027 | > | rnemdFile_ << "# momentum = " << momentumFluxVector_ |
2028 | > | << " (amu/A/fs^2)\n"; |
2029 | > | rnemdFile_ << "# angular momentum = " << angularMomentumFluxVector_ |
2030 | > | << " (amu/A^2/fs^2)\n"; |
2031 | > | rnemdFile_ << "# Target one-time exchanges:\n"; |
2032 | > | rnemdFile_ << "# kinetic = " |
2033 | > | << kineticTarget_ / PhysicalConstants::energyConvert |
2034 | > | << " (kcal/mol)\n"; |
2035 | > | rnemdFile_ << "# momentum = " << momentumTarget_ |
2036 | > | << " (amu*A/fs)\n"; |
2037 | > | rnemdFile_ << "# angular momentum = " << angularMomentumTarget_ |
2038 | > | << " (amu*A^2/fs)\n"; |
2039 | > | rnemdFile_ << "# Actual exchange totals:\n"; |
2040 | > | rnemdFile_ << "# kinetic = " |
2041 | > | << kineticExchange_ / PhysicalConstants::energyConvert |
2042 | > | << " (kcal/mol)\n"; |
2043 | > | rnemdFile_ << "# momentum = " << momentumExchange_ |
2044 | > | << " (amu*A/fs)\n"; |
2045 | > | rnemdFile_ << "# angular momentum = " << angularMomentumExchange_ |
2046 | > | << " (amu*A^2/fs)\n"; |
2047 | > | rnemdFile_ << "# Actual flux:\n"; |
2048 | > | rnemdFile_ << "# kinetic = " << Jz |
2049 | > | << " (kcal/mol/A^2/fs)\n"; |
2050 | > | rnemdFile_ << "# momentum = " << JzP |
2051 | > | << " (amu/A/fs^2)\n"; |
2052 | > | rnemdFile_ << "# angular momentum = " << JzL |
2053 | > | << " (amu/A^2/fs^2)\n"; |
2054 | > | rnemdFile_ << "# Exchange statistics:\n"; |
2055 | > | rnemdFile_ << "# attempted = " << trialCount_ << "\n"; |
2056 | > | rnemdFile_ << "# failed = " << failTrialCount_ << "\n"; |
2057 | if (rnemdMethod_ == rnemdNIVS) { | |
2058 | < | rnemdFile_ << "# NIVS root-check warnings: " << failRootCount_ << "\n"; |
2058 | > | rnemdFile_ << "# NIVS root-check errors = " |
2059 | > | << failRootCount_ << "\n"; |
2060 | } | |
1637 | – | |
2061 | rnemdFile_ << "#######################################################\n"; | |
2062 | ||
2063 | ||
# | Line 1645 | Line 2068 | namespace OpenMD { | |
2068 | if (outputMask_[i]) { | |
2069 | rnemdFile_ << "\t" << data_[i].title << | |
2070 | "(" << data_[i].units << ")"; | |
2071 | + | // add some extra tabs for column alignment |
2072 | + | if (data_[i].dataType == "Vector3d") rnemdFile_ << "\t\t"; |
2073 | } | |
2074 | } | |
2075 | rnemdFile_ << std::endl; | |
2076 | ||
2077 | rnemdFile_.precision(8); | |
2078 | ||
2079 | < | for (unsigned int j = 0; j < nBins_; j++) { |
2079 | > | for (int j = 0; j < nBins_; j++) { |
2080 | ||
2081 | for (unsigned int i = 0; i < outputMask_.size(); ++i) { | |
2082 | if (outputMask_[i]) { | |
2083 | if (data_[i].dataType == "RealType") | |
2084 | writeReal(i,j); | |
2085 | < | else if (data_[i].dataType == "Vector3d") |
2085 | > | else if (data_[i].dataType == "Vector3d") |
2086 | writeVector(i,j); | |
2087 | else { | |
2088 | sprintf( painCave.errMsg, | |
# | Line 1671 | Line 2096 | namespace OpenMD { | |
2096 | rnemdFile_ << std::endl; | |
2097 | ||
2098 | } | |
2099 | + | |
2100 | + | rnemdFile_ << "#######################################################\n"; |
2101 | + | rnemdFile_ << "# Standard Deviations in those quantities follow:\n"; |
2102 | + | rnemdFile_ << "#######################################################\n"; |
2103 | + | |
2104 | + | |
2105 | + | for (int j = 0; j < nBins_; j++) { |
2106 | + | rnemdFile_ << "#"; |
2107 | + | for (unsigned int i = 0; i < outputMask_.size(); ++i) { |
2108 | + | if (outputMask_[i]) { |
2109 | + | if (data_[i].dataType == "RealType") |
2110 | + | writeRealStdDev(i,j); |
2111 | + | else if (data_[i].dataType == "Vector3d") |
2112 | + | writeVectorStdDev(i,j); |
2113 | + | else { |
2114 | + | sprintf( painCave.errMsg, |
2115 | + | "RNEMD found an unknown data type for: %s ", |
2116 | + | data_[i].title.c_str()); |
2117 | + | painCave.isFatal = 1; |
2118 | + | simError(); |
2119 | + | } |
2120 | + | } |
2121 | + | } |
2122 | + | rnemdFile_ << std::endl; |
2123 | + | |
2124 | + | } |
2125 | ||
2126 | rnemdFile_.flush(); | |
2127 | rnemdFile_.close(); | |
# | Line 1682 | Line 2133 | namespace OpenMD { | |
2133 | } | |
2134 | ||
2135 | void RNEMD::writeReal(int index, unsigned int bin) { | |
2136 | + | if (!doRNEMD_) return; |
2137 | assert(index >=0 && index < ENDINDEX); | |
2138 | < | assert(bin >=0 && bin < nBins_); |
2138 | > | assert(int(bin) < nBins_); |
2139 | RealType s; | |
2140 | + | int count; |
2141 | ||
2142 | < | data_[index].accumulator[bin]->getAverage(s); |
2142 | > | count = data_[index].accumulator[bin]->count(); |
2143 | > | if (count == 0) return; |
2144 | ||
2145 | + | dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getAverage(s); |
2146 | + | |
2147 | if (! isinf(s) && ! isnan(s)) { | |
2148 | rnemdFile_ << "\t" << s; | |
2149 | } else{ | |
# | Line 1700 | Line 2156 | namespace OpenMD { | |
2156 | } | |
2157 | ||
2158 | void RNEMD::writeVector(int index, unsigned int bin) { | |
2159 | + | if (!doRNEMD_) return; |
2160 | assert(index >=0 && index < ENDINDEX); | |
2161 | < | assert(bin >=0 && bin < nBins_); |
2161 | > | assert(int(bin) < nBins_); |
2162 | Vector3d s; | |
2163 | + | int count; |
2164 | + | |
2165 | + | count = data_[index].accumulator[bin]->count(); |
2166 | + | |
2167 | + | if (count == 0) return; |
2168 | + | |
2169 | dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getAverage(s); | |
2170 | if (isinf(s[0]) || isnan(s[0]) || | |
2171 | isinf(s[1]) || isnan(s[1]) || | |
# | Line 1716 | Line 2179 | namespace OpenMD { | |
2179 | rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2]; | |
2180 | } | |
2181 | } | |
2182 | + | |
2183 | + | void RNEMD::writeRealStdDev(int index, unsigned int bin) { |
2184 | + | if (!doRNEMD_) return; |
2185 | + | assert(index >=0 && index < ENDINDEX); |
2186 | + | assert(int(bin) < nBins_); |
2187 | + | RealType s; |
2188 | + | int count; |
2189 | + | |
2190 | + | count = data_[index].accumulator[bin]->count(); |
2191 | + | if (count == 0) return; |
2192 | + | |
2193 | + | dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getStdDev(s); |
2194 | + | |
2195 | + | if (! isinf(s) && ! isnan(s)) { |
2196 | + | rnemdFile_ << "\t" << s; |
2197 | + | } else{ |
2198 | + | sprintf( painCave.errMsg, |
2199 | + | "RNEMD detected a numerical error writing: %s std. dev. for bin %d", |
2200 | + | data_[index].title.c_str(), bin); |
2201 | + | painCave.isFatal = 1; |
2202 | + | simError(); |
2203 | + | } |
2204 | + | } |
2205 | + | |
2206 | + | void RNEMD::writeVectorStdDev(int index, unsigned int bin) { |
2207 | + | if (!doRNEMD_) return; |
2208 | + | assert(index >=0 && index < ENDINDEX); |
2209 | + | assert(int(bin) < nBins_); |
2210 | + | Vector3d s; |
2211 | + | int count; |
2212 | + | |
2213 | + | count = data_[index].accumulator[bin]->count(); |
2214 | + | if (count == 0) return; |
2215 | + | |
2216 | + | dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s); |
2217 | + | if (isinf(s[0]) || isnan(s[0]) || |
2218 | + | isinf(s[1]) || isnan(s[1]) || |
2219 | + | isinf(s[2]) || isnan(s[2]) ) { |
2220 | + | sprintf( painCave.errMsg, |
2221 | + | "RNEMD detected a numerical error writing: %s std. dev. for bin %d", |
2222 | + | data_[index].title.c_str(), bin); |
2223 | + | painCave.isFatal = 1; |
2224 | + | simError(); |
2225 | + | } else { |
2226 | + | rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2]; |
2227 | + | } |
2228 | + | } |
2229 | } | |
2230 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |