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