58 |
|
#include "utils/simError.h" |
59 |
|
#include "utils/MemoryUtils.hpp" |
60 |
|
#include "utils/StringTokenizer.hpp" |
61 |
+ |
#include "brains/Thermo.hpp" |
62 |
|
|
63 |
|
#ifdef IS_MPI |
63 |
– |
|
64 |
|
#include <mpi.h> |
65 |
< |
#define TAKE_THIS_TAG_CHAR 0 |
66 |
< |
#define TAKE_THIS_TAG_INT 1 |
65 |
> |
#endif |
66 |
|
|
68 |
– |
#endif // is_mpi |
67 |
|
|
70 |
– |
|
68 |
|
namespace OpenMD { |
69 |
|
|
70 |
|
DumpReader::DumpReader(SimInfo* info, const std::string& filename) |
226 |
|
needVel_ = false; |
227 |
|
} |
228 |
|
|
229 |
< |
if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) { |
229 |
> |
if (storageLayout & DataStorage::dslAmat || |
230 |
> |
storageLayout & DataStorage::dslDipole || |
231 |
> |
storageLayout & DataStorage::dslQuadrupole) { |
232 |
|
needQuaternion_ = true; |
233 |
|
} else { |
234 |
|
needQuaternion_ = false; |
244 |
|
|
245 |
|
if (needCOMprops_) { |
246 |
|
Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot(); |
247 |
+ |
Thermo thermo(info_); |
248 |
|
Vector3d com; |
249 |
– |
Vector3d comvel; |
250 |
– |
Vector3d comw; |
251 |
– |
if (needPos_ && needVel_){ |
252 |
– |
info_->getComAll(com, comvel); |
253 |
– |
comw = info_->getAngularMomentum(); |
254 |
– |
}else{ |
255 |
– |
com = info_->getCom(); |
256 |
– |
comvel = 0.0; |
257 |
– |
comw = 0.0; |
258 |
– |
} |
259 |
– |
s->setCOMprops(com, comvel, comw); |
260 |
– |
} |
249 |
|
|
250 |
+ |
if (needPos_ && needVel_) { |
251 |
+ |
Vector3d comvel; |
252 |
+ |
Vector3d comw; |
253 |
+ |
thermo.getComAll(com, comvel); |
254 |
+ |
comw = thermo.getAngularMomentum(); |
255 |
+ |
} else { |
256 |
+ |
com = thermo.getCom(); |
257 |
+ |
} |
258 |
+ |
} |
259 |
|
} |
260 |
|
|
261 |
|
void DumpReader::readSet(int whichFrame) { |
354 |
|
|
355 |
|
int index = tokenizer.nextTokenAsInt(); |
356 |
|
|
357 |
< |
StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index); |
357 |
> |
StuntDouble* sd = info_->getIOIndexToIntegrableObject(index); |
358 |
|
|
359 |
< |
if (integrableObject == NULL) { |
359 |
> |
if (sd == NULL) { |
360 |
|
return; |
361 |
|
} |
362 |
|
std::string type = tokenizer.nextToken(); |
376 |
|
} |
377 |
|
} |
378 |
|
|
379 |
< |
if (integrableObject->isDirectional()) { |
379 |
> |
if (sd->isDirectional()) { |
380 |
|
if (needQuaternion_) { |
381 |
|
found = type.find("q"); |
382 |
|
if (found == std::string::npos) { |
399 |
|
pos[1] = tokenizer.nextTokenAsDouble(); |
400 |
|
pos[2] = tokenizer.nextTokenAsDouble(); |
401 |
|
if (needPos_) { |
402 |
< |
integrableObject->setPos(pos); |
402 |
> |
sd->setPos(pos); |
403 |
|
} |
404 |
|
break; |
405 |
|
} |
409 |
|
vel[1] = tokenizer.nextTokenAsDouble(); |
410 |
|
vel[2] = tokenizer.nextTokenAsDouble(); |
411 |
|
if (needVel_) { |
412 |
< |
integrableObject->setVel(vel); |
412 |
> |
sd->setVel(vel); |
413 |
|
} |
414 |
|
break; |
415 |
|
} |
416 |
|
|
417 |
|
case 'q' : { |
418 |
|
Quat4d q; |
419 |
< |
if (integrableObject->isDirectional()) { |
419 |
> |
if (sd->isDirectional()) { |
420 |
|
|
421 |
|
q[0] = tokenizer.nextTokenAsDouble(); |
422 |
|
q[1] = tokenizer.nextTokenAsDouble(); |
435 |
|
|
436 |
|
q.normalize(); |
437 |
|
if (needQuaternion_) { |
438 |
< |
integrableObject->setQ(q); |
438 |
> |
sd->setQ(q); |
439 |
|
} |
440 |
|
} |
441 |
|
break; |
442 |
|
} |
443 |
|
case 'j' : { |
444 |
|
Vector3d ji; |
445 |
< |
if (integrableObject->isDirectional()) { |
445 |
> |
if (sd->isDirectional()) { |
446 |
|
ji[0] = tokenizer.nextTokenAsDouble(); |
447 |
|
ji[1] = tokenizer.nextTokenAsDouble(); |
448 |
|
ji[2] = tokenizer.nextTokenAsDouble(); |
449 |
|
if (needAngMom_) { |
450 |
< |
integrableObject->setJ(ji); |
450 |
> |
sd->setJ(ji); |
451 |
|
} |
452 |
|
} |
453 |
|
break; |
458 |
|
force[0] = tokenizer.nextTokenAsDouble(); |
459 |
|
force[1] = tokenizer.nextTokenAsDouble(); |
460 |
|
force[2] = tokenizer.nextTokenAsDouble(); |
461 |
< |
integrableObject->setFrc(force); |
461 |
> |
sd->setFrc(force); |
462 |
|
break; |
463 |
|
} |
464 |
|
case 't' : { |
467 |
|
torque[0] = tokenizer.nextTokenAsDouble(); |
468 |
|
torque[1] = tokenizer.nextTokenAsDouble(); |
469 |
|
torque[2] = tokenizer.nextTokenAsDouble(); |
470 |
< |
integrableObject->setTrq(torque); |
470 |
> |
sd->setTrq(torque); |
471 |
|
break; |
472 |
|
} |
473 |
|
case 'u' : { |
474 |
|
|
475 |
|
RealType particlePot; |
476 |
|
particlePot = tokenizer.nextTokenAsDouble(); |
477 |
< |
integrableObject->setParticlePot(particlePot); |
477 |
> |
sd->setParticlePot(particlePot); |
478 |
|
break; |
479 |
|
} |
480 |
|
case 'c' : { |
481 |
|
|
482 |
|
RealType flucQPos; |
483 |
|
flucQPos = tokenizer.nextTokenAsDouble(); |
484 |
< |
integrableObject->setFlucQPos(flucQPos); |
484 |
> |
sd->setFlucQPos(flucQPos); |
485 |
|
break; |
486 |
|
} |
487 |
|
case 'w' : { |
488 |
|
|
489 |
|
RealType flucQVel; |
490 |
|
flucQVel = tokenizer.nextTokenAsDouble(); |
491 |
< |
integrableObject->setFlucQVel(flucQVel); |
491 |
> |
sd->setFlucQVel(flucQVel); |
492 |
|
break; |
493 |
|
} |
494 |
|
case 'g' : { |
495 |
|
|
496 |
|
RealType flucQFrc; |
497 |
|
flucQFrc = tokenizer.nextTokenAsDouble(); |
498 |
< |
integrableObject->setFlucQFrc(flucQFrc); |
498 |
> |
sd->setFlucQFrc(flucQFrc); |
499 |
|
break; |
500 |
|
} |
501 |
|
case 'e' : { |
504 |
|
eField[0] = tokenizer.nextTokenAsDouble(); |
505 |
|
eField[1] = tokenizer.nextTokenAsDouble(); |
506 |
|
eField[2] = tokenizer.nextTokenAsDouble(); |
507 |
< |
integrableObject->setElectricField(eField); |
507 |
> |
sd->setElectricField(eField); |
508 |
|
break; |
509 |
|
} |
510 |
|
default: { |
540 |
|
*/ |
541 |
|
|
542 |
|
int index = tokenizer.nextTokenAsInt(); |
543 |
< |
StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index); |
544 |
< |
if (integrableObject == NULL) { |
543 |
> |
StuntDouble* sd = info_->getIOIndexToIntegrableObject(index); |
544 |
> |
if (sd == NULL) { |
545 |
|
return; |
546 |
|
} |
550 |
– |
StuntDouble* sd = integrableObject; |
547 |
|
|
548 |
|
/** |
549 |
|
* Test to see if the next token is an integer or not. If not, |
557 |
|
if (i >> siteIndex) { |
558 |
|
// chew up this token and parse as an int: |
559 |
|
siteIndex = tokenizer.nextTokenAsInt(); |
560 |
< |
RigidBody* rb = static_cast<RigidBody*>(integrableObject); |
560 |
> |
RigidBody* rb = static_cast<RigidBody*>(sd); |
561 |
|
sd = rb->getAtoms()[siteIndex]; |
562 |
|
} |
563 |
|
|
710 |
|
hmat(2, 2) = tokenizer.nextTokenAsDouble(); |
711 |
|
s->setHmat(hmat); |
712 |
|
} else if (propertyName == "Thermostat") { |
713 |
< |
RealType chi = tokenizer.nextTokenAsDouble(); |
714 |
< |
RealType integralOfChiDt = tokenizer.nextTokenAsDouble(); |
715 |
< |
s->setChi(chi); |
716 |
< |
s->setIntegralOfChiDt(integralOfChiDt); |
713 |
> |
pair<RealType, RealType> thermostat; |
714 |
> |
thermostat.first = tokenizer.nextTokenAsDouble(); |
715 |
> |
thermostat.second = tokenizer.nextTokenAsDouble(); |
716 |
> |
s->setThermostat(thermostat); |
717 |
|
} else if (propertyName == "Barostat") { |
718 |
|
Mat3x3d eta; |
719 |
|
eta(0, 0) = tokenizer.nextTokenAsDouble(); |
725 |
|
eta(2, 0) = tokenizer.nextTokenAsDouble(); |
726 |
|
eta(2, 1) = tokenizer.nextTokenAsDouble(); |
727 |
|
eta(2, 2) = tokenizer.nextTokenAsDouble(); |
728 |
< |
s->setEta(eta); |
728 |
> |
s->setBarostat(eta); |
729 |
|
} else { |
730 |
|
sprintf(painCave.errMsg, |
731 |
|
"DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str()); |