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) |
75 |
|
if (worldRank == 0) { |
76 |
|
#endif |
77 |
|
|
78 |
< |
inFile_ = new std::ifstream(filename_.c_str()); |
78 |
> |
inFile_ = new std::ifstream(filename_.c_str(), |
79 |
> |
ifstream::in | ifstream::binary); |
80 |
|
|
81 |
|
if (inFile_->fail()) { |
82 |
|
sprintf(painCave.errMsg, |
141 |
|
prevPos = currPos; |
142 |
|
bool foundOpenSnapshotTag = false; |
143 |
|
bool foundClosedSnapshotTag = false; |
144 |
+ |
|
145 |
|
while(inFile_->getline(buffer, bufferSize)) { |
146 |
|
++lineNo; |
147 |
|
|
202 |
|
#ifdef IS_MPI |
203 |
|
} |
204 |
|
|
205 |
< |
MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD); |
205 |
> |
MPI::COMM_WORLD.Bcast(&nframes_, 1, MPI::INT, 0); |
206 |
|
|
207 |
|
#endif // is_mpi |
208 |
|
|
227 |
|
needVel_ = false; |
228 |
|
} |
229 |
|
|
230 |
< |
if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) { |
230 |
> |
if (storageLayout & DataStorage::dslAmat || |
231 |
> |
storageLayout & DataStorage::dslDipole || |
232 |
> |
storageLayout & DataStorage::dslQuadrupole) { |
233 |
|
needQuaternion_ = true; |
234 |
|
} else { |
235 |
|
needQuaternion_ = false; |
245 |
|
|
246 |
|
if (needCOMprops_) { |
247 |
|
Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot(); |
248 |
+ |
Thermo thermo(info_); |
249 |
|
Vector3d com; |
248 |
– |
Vector3d comvel; |
249 |
– |
Vector3d comw; |
250 |
– |
if (needPos_ && needVel_){ |
251 |
– |
info_->getComAll(com, comvel); |
252 |
– |
comw = info_->getAngularMomentum(); |
253 |
– |
}else{ |
254 |
– |
com = info_->getCom(); |
255 |
– |
comvel = 0.0; |
256 |
– |
comw = 0.0; |
257 |
– |
} |
258 |
– |
s->setCOMprops(com, comvel, comw); |
259 |
– |
} |
250 |
|
|
251 |
+ |
if (needPos_ && needVel_) { |
252 |
+ |
Vector3d comvel; |
253 |
+ |
Vector3d comw; |
254 |
+ |
thermo.getComAll(com, comvel); |
255 |
+ |
comw = thermo.getAngularMomentum(); |
256 |
+ |
} else { |
257 |
+ |
com = thermo.getCom(); |
258 |
+ |
} |
259 |
+ |
} |
260 |
|
} |
261 |
|
|
262 |
|
void DumpReader::readSet(int whichFrame) { |
288 |
|
} |
289 |
|
|
290 |
|
int sendBufferSize = sendBuffer.size(); |
291 |
< |
MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
292 |
< |
MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
291 |
> |
MPI::COMM_WORLD.Bcast(&sendBufferSize, 1, MPI::INT, masterNode); |
292 |
> |
MPI::COMM_WORLD.Bcast((void *)sendBuffer.c_str(), sendBufferSize, |
293 |
> |
MPI::CHAR, masterNode); |
294 |
|
|
295 |
|
sstream.str(sendBuffer); |
296 |
|
} else { |
297 |
|
int sendBufferSize; |
298 |
< |
MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
298 |
> |
MPI::COMM_WORLD.Bcast(&sendBufferSize, 1, MPI::INT, masterNode); |
299 |
|
char * recvBuffer = new char[sendBufferSize+1]; |
300 |
|
assert(recvBuffer); |
301 |
|
recvBuffer[sendBufferSize] = '\0'; |
302 |
< |
MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
302 |
> |
MPI::COMM_WORLD.Bcast(recvBuffer, sendBufferSize, MPI::CHAR, masterNode); |
303 |
|
sstream.str(recvBuffer); |
304 |
|
delete [] recvBuffer; |
305 |
|
} |
325 |
|
|
326 |
|
inputStream.getline(buffer, bufferSize); |
327 |
|
line = buffer; |
328 |
< |
if (line.find("</Snapshot>") == std::string::npos) { |
329 |
< |
sprintf(painCave.errMsg, |
330 |
< |
"DumpReader Error: can not find </Snapshot>\n"); |
331 |
< |
painCave.isFatal = 1; |
332 |
< |
simError(); |
333 |
< |
} |
334 |
< |
|
328 |
> |
|
329 |
> |
if (line.find("<SiteData>") != std::string::npos) { |
330 |
> |
//read SiteData |
331 |
> |
readSiteData(inputStream); |
332 |
> |
} else { |
333 |
> |
if (line.find("</Snapshot>") == std::string::npos) { |
334 |
> |
sprintf(painCave.errMsg, |
335 |
> |
"DumpReader Error: can not find </Snapshot>\n"); |
336 |
> |
painCave.isFatal = 1; |
337 |
> |
simError(); |
338 |
> |
} |
339 |
> |
} |
340 |
|
} |
341 |
|
|
342 |
|
void DumpReader::parseDumpLine(const std::string& line) { |
356 |
|
|
357 |
|
int index = tokenizer.nextTokenAsInt(); |
358 |
|
|
359 |
< |
StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index); |
359 |
> |
StuntDouble* sd = info_->getIOIndexToIntegrableObject(index); |
360 |
|
|
361 |
< |
if (integrableObject == NULL) { |
361 |
> |
if (sd == NULL) { |
362 |
|
return; |
363 |
|
} |
364 |
|
std::string type = tokenizer.nextToken(); |
378 |
|
} |
379 |
|
} |
380 |
|
|
381 |
< |
if (integrableObject->isDirectional()) { |
381 |
> |
if (sd->isDirectional()) { |
382 |
|
if (needQuaternion_) { |
383 |
|
found = type.find("q"); |
384 |
|
if (found == std::string::npos) { |
401 |
|
pos[1] = tokenizer.nextTokenAsDouble(); |
402 |
|
pos[2] = tokenizer.nextTokenAsDouble(); |
403 |
|
if (needPos_) { |
404 |
< |
integrableObject->setPos(pos); |
404 |
> |
sd->setPos(pos); |
405 |
|
} |
406 |
|
break; |
407 |
|
} |
411 |
|
vel[1] = tokenizer.nextTokenAsDouble(); |
412 |
|
vel[2] = tokenizer.nextTokenAsDouble(); |
413 |
|
if (needVel_) { |
414 |
< |
integrableObject->setVel(vel); |
414 |
> |
sd->setVel(vel); |
415 |
|
} |
416 |
|
break; |
417 |
|
} |
418 |
|
|
419 |
|
case 'q' : { |
420 |
|
Quat4d q; |
421 |
< |
if (integrableObject->isDirectional()) { |
421 |
> |
if (sd->isDirectional()) { |
422 |
|
|
423 |
|
q[0] = tokenizer.nextTokenAsDouble(); |
424 |
|
q[1] = tokenizer.nextTokenAsDouble(); |
437 |
|
|
438 |
|
q.normalize(); |
439 |
|
if (needQuaternion_) { |
440 |
< |
integrableObject->setQ(q); |
440 |
> |
sd->setQ(q); |
441 |
|
} |
442 |
|
} |
443 |
|
break; |
444 |
|
} |
445 |
|
case 'j' : { |
446 |
|
Vector3d ji; |
447 |
< |
if (integrableObject->isDirectional()) { |
447 |
> |
if (sd->isDirectional()) { |
448 |
|
ji[0] = tokenizer.nextTokenAsDouble(); |
449 |
|
ji[1] = tokenizer.nextTokenAsDouble(); |
450 |
|
ji[2] = tokenizer.nextTokenAsDouble(); |
451 |
|
if (needAngMom_) { |
452 |
< |
integrableObject->setJ(ji); |
452 |
> |
sd->setJ(ji); |
453 |
|
} |
454 |
|
} |
455 |
|
break; |
460 |
|
force[0] = tokenizer.nextTokenAsDouble(); |
461 |
|
force[1] = tokenizer.nextTokenAsDouble(); |
462 |
|
force[2] = tokenizer.nextTokenAsDouble(); |
463 |
< |
integrableObject->setFrc(force); |
463 |
> |
sd->setFrc(force); |
464 |
|
break; |
465 |
|
} |
466 |
|
case 't' : { |
469 |
|
torque[0] = tokenizer.nextTokenAsDouble(); |
470 |
|
torque[1] = tokenizer.nextTokenAsDouble(); |
471 |
|
torque[2] = tokenizer.nextTokenAsDouble(); |
472 |
< |
integrableObject->setTrq(torque); |
472 |
> |
sd->setTrq(torque); |
473 |
|
break; |
474 |
|
} |
475 |
|
case 'u' : { |
476 |
|
|
477 |
|
RealType particlePot; |
478 |
|
particlePot = tokenizer.nextTokenAsDouble(); |
479 |
< |
integrableObject->setParticlePot(particlePot); |
479 |
> |
sd->setParticlePot(particlePot); |
480 |
|
break; |
481 |
|
} |
482 |
|
case 'c' : { |
483 |
|
|
484 |
|
RealType flucQPos; |
485 |
|
flucQPos = tokenizer.nextTokenAsDouble(); |
486 |
< |
integrableObject->setFlucQPos(flucQPos); |
486 |
> |
sd->setFlucQPos(flucQPos); |
487 |
|
break; |
488 |
|
} |
489 |
|
case 'w' : { |
490 |
|
|
491 |
|
RealType flucQVel; |
492 |
|
flucQVel = tokenizer.nextTokenAsDouble(); |
493 |
< |
integrableObject->setFlucQVel(flucQVel); |
493 |
> |
sd->setFlucQVel(flucQVel); |
494 |
|
break; |
495 |
|
} |
496 |
|
case 'g' : { |
497 |
|
|
498 |
|
RealType flucQFrc; |
499 |
|
flucQFrc = tokenizer.nextTokenAsDouble(); |
500 |
< |
integrableObject->setFlucQFrc(flucQFrc); |
500 |
> |
sd->setFlucQFrc(flucQFrc); |
501 |
|
break; |
502 |
|
} |
503 |
|
case 'e' : { |
506 |
|
eField[0] = tokenizer.nextTokenAsDouble(); |
507 |
|
eField[1] = tokenizer.nextTokenAsDouble(); |
508 |
|
eField[2] = tokenizer.nextTokenAsDouble(); |
509 |
< |
integrableObject->setElectricField(eField); |
509 |
> |
sd->setElectricField(eField); |
510 |
|
break; |
511 |
|
} |
512 |
|
default: { |
523 |
|
} |
524 |
|
|
525 |
|
|
526 |
< |
void DumpReader::readStuntDoubles(std::istream& inputStream) { |
526 |
> |
void DumpReader::parseSiteLine(const std::string& line) { |
527 |
|
|
528 |
+ |
StringTokenizer tokenizer(line); |
529 |
+ |
int nTokens; |
530 |
+ |
|
531 |
+ |
nTokens = tokenizer.countTokens(); |
532 |
+ |
|
533 |
+ |
if (nTokens < 2) { |
534 |
+ |
sprintf(painCave.errMsg, |
535 |
+ |
"DumpReader Error: Not enough Tokens.\n%s\n", line.c_str()); |
536 |
+ |
painCave.isFatal = 1; |
537 |
+ |
simError(); |
538 |
+ |
} |
539 |
+ |
|
540 |
+ |
/** |
541 |
+ |
* The first token is the global integrable object index. |
542 |
+ |
*/ |
543 |
+ |
|
544 |
+ |
int index = tokenizer.nextTokenAsInt(); |
545 |
+ |
StuntDouble* sd = info_->getIOIndexToIntegrableObject(index); |
546 |
+ |
if (sd == NULL) { |
547 |
+ |
return; |
548 |
+ |
} |
549 |
+ |
|
550 |
+ |
/** |
551 |
+ |
* Test to see if the next token is an integer or not. If not, |
552 |
+ |
* we've got data on the integrable object itself. If there is an |
553 |
+ |
* integer, we're parsing data for a site on a rigid body. |
554 |
+ |
*/ |
555 |
+ |
|
556 |
+ |
std::string indexTest = tokenizer.peekNextToken(); |
557 |
+ |
std::istringstream i(indexTest); |
558 |
+ |
int siteIndex; |
559 |
+ |
if (i >> siteIndex) { |
560 |
+ |
// chew up this token and parse as an int: |
561 |
+ |
siteIndex = tokenizer.nextTokenAsInt(); |
562 |
+ |
RigidBody* rb = static_cast<RigidBody*>(sd); |
563 |
+ |
sd = rb->getAtoms()[siteIndex]; |
564 |
+ |
} |
565 |
+ |
|
566 |
+ |
/** |
567 |
+ |
* The next token contains information on what follows. |
568 |
+ |
*/ |
569 |
+ |
std::string type = tokenizer.nextToken(); |
570 |
+ |
int size = type.size(); |
571 |
+ |
|
572 |
+ |
for(int i = 0; i < size; ++i) { |
573 |
+ |
switch(type[i]) { |
574 |
+ |
|
575 |
+ |
case 'u' : { |
576 |
+ |
|
577 |
+ |
RealType particlePot; |
578 |
+ |
particlePot = tokenizer.nextTokenAsDouble(); |
579 |
+ |
sd->setParticlePot(particlePot); |
580 |
+ |
break; |
581 |
+ |
} |
582 |
+ |
case 'c' : { |
583 |
+ |
|
584 |
+ |
RealType flucQPos; |
585 |
+ |
flucQPos = tokenizer.nextTokenAsDouble(); |
586 |
+ |
sd->setFlucQPos(flucQPos); |
587 |
+ |
break; |
588 |
+ |
} |
589 |
+ |
case 'w' : { |
590 |
+ |
|
591 |
+ |
RealType flucQVel; |
592 |
+ |
flucQVel = tokenizer.nextTokenAsDouble(); |
593 |
+ |
sd->setFlucQVel(flucQVel); |
594 |
+ |
break; |
595 |
+ |
} |
596 |
+ |
case 'g' : { |
597 |
+ |
|
598 |
+ |
RealType flucQFrc; |
599 |
+ |
flucQFrc = tokenizer.nextTokenAsDouble(); |
600 |
+ |
sd->setFlucQFrc(flucQFrc); |
601 |
+ |
break; |
602 |
+ |
} |
603 |
+ |
case 'e' : { |
604 |
+ |
|
605 |
+ |
Vector3d eField; |
606 |
+ |
eField[0] = tokenizer.nextTokenAsDouble(); |
607 |
+ |
eField[1] = tokenizer.nextTokenAsDouble(); |
608 |
+ |
eField[2] = tokenizer.nextTokenAsDouble(); |
609 |
+ |
sd->setElectricField(eField); |
610 |
+ |
break; |
611 |
+ |
} |
612 |
+ |
default: { |
613 |
+ |
sprintf(painCave.errMsg, |
614 |
+ |
"DumpReader Error: %s is an unrecognized type\n", type.c_str()); |
615 |
+ |
painCave.isFatal = 1; |
616 |
+ |
simError(); |
617 |
+ |
break; |
618 |
+ |
} |
619 |
+ |
} |
620 |
+ |
} |
621 |
+ |
} |
622 |
+ |
|
623 |
+ |
|
624 |
+ |
void DumpReader::readStuntDoubles(std::istream& inputStream) { |
625 |
+ |
|
626 |
|
inputStream.getline(buffer, bufferSize); |
627 |
|
std::string line(buffer); |
628 |
|
|
645 |
|
|
646 |
|
} |
647 |
|
|
648 |
+ |
void DumpReader::readSiteData(std::istream& inputStream) { |
649 |
+ |
|
650 |
+ |
inputStream.getline(buffer, bufferSize); |
651 |
+ |
std::string line(buffer); |
652 |
+ |
|
653 |
+ |
if (line.find("<SiteData>") == std::string::npos) { |
654 |
+ |
// site data isn't required for a simulation, so skip |
655 |
+ |
return; |
656 |
+ |
} |
657 |
+ |
|
658 |
+ |
while(inputStream.getline(buffer, bufferSize)) { |
659 |
+ |
line = buffer; |
660 |
+ |
|
661 |
+ |
if(line.find("</SiteData>") != std::string::npos) { |
662 |
+ |
break; |
663 |
+ |
} |
664 |
+ |
|
665 |
+ |
parseSiteLine(line); |
666 |
+ |
} |
667 |
+ |
|
668 |
+ |
} |
669 |
+ |
|
670 |
|
void DumpReader::readFrameProperties(std::istream& inputStream) { |
671 |
|
|
672 |
|
Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot(); |
712 |
|
hmat(2, 2) = tokenizer.nextTokenAsDouble(); |
713 |
|
s->setHmat(hmat); |
714 |
|
} else if (propertyName == "Thermostat") { |
715 |
< |
RealType chi = tokenizer.nextTokenAsDouble(); |
716 |
< |
RealType integralOfChiDt = tokenizer.nextTokenAsDouble(); |
717 |
< |
s->setChi(chi); |
718 |
< |
s->setIntegralOfChiDt(integralOfChiDt); |
715 |
> |
pair<RealType, RealType> thermostat; |
716 |
> |
thermostat.first = tokenizer.nextTokenAsDouble(); |
717 |
> |
thermostat.second = tokenizer.nextTokenAsDouble(); |
718 |
> |
s->setThermostat(thermostat); |
719 |
|
} else if (propertyName == "Barostat") { |
720 |
|
Mat3x3d eta; |
721 |
|
eta(0, 0) = tokenizer.nextTokenAsDouble(); |
727 |
|
eta(2, 0) = tokenizer.nextTokenAsDouble(); |
728 |
|
eta(2, 1) = tokenizer.nextTokenAsDouble(); |
729 |
|
eta(2, 2) = tokenizer.nextTokenAsDouble(); |
730 |
< |
s->setEta(eta); |
730 |
> |
s->setBarostat(eta); |
731 |
|
} else { |
732 |
|
sprintf(painCave.errMsg, |
733 |
|
"DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str()); |