| 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). |
| 39 |
< |
* [4] Vardeman & Gezelter, in progress (2009). |
| 38 |
> |
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). |
| 39 |
> |
* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
| 40 |
> |
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
| 41 |
|
*/ |
| 41 |
– |
|
| 42 |
|
|
| 43 |
+ |
#ifdef IS_MPI |
| 44 |
+ |
#include <mpi.h> |
| 45 |
+ |
#endif |
| 46 |
+ |
|
| 47 |
|
#include <sys/types.h> |
| 48 |
|
#include <sys/stat.h> |
| 49 |
|
|
| 62 |
|
#include "utils/StringTokenizer.hpp" |
| 63 |
|
#include "restraints/ObjectRestraint.hpp" |
| 64 |
|
#include "restraints/MolecularRestraint.hpp" |
| 61 |
– |
|
| 62 |
– |
#ifdef IS_MPI |
| 63 |
– |
|
| 64 |
– |
#include <mpi.h> |
| 65 |
– |
#endif |
| 65 |
|
|
| 66 |
|
namespace OpenMD { |
| 67 |
|
|
| 68 |
|
void RestReader::scanFile(){ |
| 69 |
< |
int lineNo = 0; |
| 69 |
> |
|
| 70 |
|
std::streampos prevPos; |
| 71 |
< |
std::streampos currPos; |
| 71 |
> |
std::streampos currPos; |
| 72 |
|
|
| 73 |
|
#ifdef IS_MPI |
| 74 |
|
|
| 80 |
|
prevPos = currPos; |
| 81 |
|
|
| 82 |
|
bool foundOpenSnapshotTag = false; |
| 83 |
< |
|
| 83 |
> |
int lineNo = 0; |
| 84 |
|
while(!foundOpenSnapshotTag && inFile_->getline(buffer, bufferSize)) { |
| 85 |
|
++lineNo; |
| 86 |
|
|
| 88 |
|
currPos = inFile_->tellg(); |
| 89 |
|
if (line.find("<Snapshot>")!= std::string::npos) { |
| 90 |
|
foundOpenSnapshotTag = true; |
| 91 |
< |
framePos_ = prevPos; |
| 91 |
> |
framePos_ = (long long) prevPos; |
| 92 |
|
} |
| 93 |
|
prevPos = currPos; |
| 94 |
|
} |
| 95 |
|
|
| 96 |
|
#ifdef IS_MPI |
| 97 |
|
} |
| 98 |
< |
MPI_Bcast(&framePos_, 1, MPI_INT, 0, MPI_COMM_WORLD); |
| 98 |
> |
MPI_Bcast(&framePos_, 1, MPI_LONG_LONG, 0, MPI_COMM_WORLD); |
| 99 |
|
#endif // is_mpi |
| 100 |
|
} |
| 101 |
|
|
| 130 |
|
} |
| 131 |
|
|
| 132 |
|
int sendBufferSize = sendBuffer.size(); |
| 133 |
< |
MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
| 134 |
< |
MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
| 133 |
> |
MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
| 134 |
> |
MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, |
| 135 |
> |
MPI_CHAR, masterNode, MPI_COMM_WORLD); |
| 136 |
|
|
| 137 |
|
sstream.str(sendBuffer); |
| 138 |
|
} else { |
| 139 |
|
int sendBufferSize; |
| 140 |
< |
MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
| 140 |
> |
MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
| 141 |
|
char * recvBuffer = new char[sendBufferSize+1]; |
| 142 |
|
assert(recvBuffer); |
| 143 |
|
recvBuffer[sendBufferSize] = '\0'; |
| 144 |
< |
MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
| 144 |
> |
MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, |
| 145 |
> |
MPI_COMM_WORLD); |
| 146 |
|
sstream.str(recvBuffer); |
| 147 |
|
delete [] recvBuffer; |
| 148 |
|
} |
| 192 |
|
|
| 193 |
|
scanFile(); |
| 194 |
|
|
| 194 |
– |
|
| 195 |
|
readSet(); |
| 196 |
|
|
| 197 |
|
|
| 280 |
|
|
| 281 |
|
int index = tokenizer.nextTokenAsInt(); |
| 282 |
|
|
| 283 |
+ |
StuntDouble* sd = info_->getIOIndexToIntegrableObject(index); |
| 284 |
+ |
|
| 285 |
+ |
if (sd == NULL) { |
| 286 |
+ |
return; |
| 287 |
+ |
} |
| 288 |
|
|
| 289 |
< |
if (index < 1116){ |
| 290 |
< |
std::string type = tokenizer.nextToken(); |
| 291 |
< |
|
| 292 |
< |
Vector3d pos; |
| 293 |
< |
|
| 294 |
< |
pos[0] = tokenizer.nextTokenAsDouble(); |
| 295 |
< |
pos[1] = tokenizer.nextTokenAsDouble(); |
| 296 |
< |
pos[2] = tokenizer.nextTokenAsDouble(); |
| 297 |
< |
|
| 298 |
< |
all_pos_[index] = pos; |
| 294 |
< |
}else{ |
| 295 |
< |
|
| 296 |
< |
bool exist = false; |
| 297 |
< |
int indexCount = 0; |
| 298 |
< |
|
| 299 |
< |
while ( (!exist) && (indexCount < stuntDoubleIndex_.size())){ |
| 300 |
< |
if (index == stuntDoubleIndex_[indexCount]) |
| 301 |
< |
exist = true; |
| 302 |
< |
indexCount++; |
| 303 |
< |
} |
| 304 |
< |
|
| 305 |
< |
StuntDouble* integrableObject; |
| 306 |
< |
|
| 307 |
< |
if (exist){ |
| 308 |
< |
|
| 309 |
< |
integrableObject = info_->getIOIndexToIntegrableObject(index); |
| 310 |
< |
|
| 311 |
< |
int compareindex = integrableObject->getGlobalIntegrableObjectIndex(); |
| 312 |
< |
|
| 313 |
< |
if (integrableObject == NULL) { |
| 314 |
< |
return; |
| 315 |
< |
} |
| 316 |
< |
|
| 317 |
< |
std::string type = tokenizer.nextToken(); |
| 318 |
< |
|
| 319 |
< |
Vector3d pos; |
| 320 |
< |
|
| 289 |
> |
std::string type = tokenizer.nextToken(); |
| 290 |
> |
int size = type.size(); |
| 291 |
> |
|
| 292 |
> |
Vector3d pos; |
| 293 |
> |
Quat4d q; |
| 294 |
> |
|
| 295 |
> |
for(int i = 0; i < size; ++i) { |
| 296 |
> |
switch(type[i]) { |
| 297 |
> |
|
| 298 |
> |
case 'p': { |
| 299 |
|
pos[0] = tokenizer.nextTokenAsDouble(); |
| 300 |
|
pos[1] = tokenizer.nextTokenAsDouble(); |
| 301 |
|
pos[2] = tokenizer.nextTokenAsDouble(); |
| 302 |
< |
|
| 302 |
> |
break; |
| 303 |
> |
} |
| 304 |
> |
case 'v' : { |
| 305 |
|
Vector3d vel; |
| 306 |
|
vel[0] = tokenizer.nextTokenAsDouble(); |
| 307 |
|
vel[1] = tokenizer.nextTokenAsDouble(); |
| 308 |
|
vel[2] = tokenizer.nextTokenAsDouble(); |
| 309 |
< |
|
| 310 |
< |
|
| 311 |
< |
Quat4d q; |
| 312 |
< |
if (integrableObject->isDirectional()) { |
| 309 |
> |
break; |
| 310 |
> |
} |
| 311 |
> |
|
| 312 |
> |
case 'q' : { |
| 313 |
> |
if (sd->isDirectional()) { |
| 314 |
|
|
| 315 |
|
q[0] = tokenizer.nextTokenAsDouble(); |
| 316 |
|
q[1] = tokenizer.nextTokenAsDouble(); |
| 317 |
|
q[2] = tokenizer.nextTokenAsDouble(); |
| 318 |
|
q[3] = tokenizer.nextTokenAsDouble(); |
| 319 |
< |
} |
| 320 |
< |
// keep the position in case we need it for a molecular restraint: |
| 319 |
> |
|
| 320 |
> |
RealType qlen = q.length(); |
| 321 |
> |
if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0 |
| 322 |
> |
|
| 323 |
> |
sprintf(painCave.errMsg, |
| 324 |
> |
"RestReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n"); |
| 325 |
> |
painCave.isFatal = 1; |
| 326 |
> |
simError(); |
| 327 |
> |
} |
| 328 |
> |
|
| 329 |
> |
q.normalize(); |
| 330 |
> |
} |
| 331 |
> |
break; |
| 332 |
> |
} |
| 333 |
> |
case 'j' : { |
| 334 |
> |
Vector3d ji; |
| 335 |
> |
if (sd->isDirectional()) { |
| 336 |
> |
ji[0] = tokenizer.nextTokenAsDouble(); |
| 337 |
> |
ji[1] = tokenizer.nextTokenAsDouble(); |
| 338 |
> |
ji[2] = tokenizer.nextTokenAsDouble(); |
| 339 |
> |
} |
| 340 |
> |
break; |
| 341 |
> |
} |
| 342 |
> |
case 'f': { |
| 343 |
> |
Vector3d force; |
| 344 |
> |
force[0] = tokenizer.nextTokenAsDouble(); |
| 345 |
> |
force[1] = tokenizer.nextTokenAsDouble(); |
| 346 |
> |
force[2] = tokenizer.nextTokenAsDouble(); |
| 347 |
> |
break; |
| 348 |
> |
} |
| 349 |
> |
case 't' : { |
| 350 |
> |
Vector3d torque; |
| 351 |
> |
torque[0] = tokenizer.nextTokenAsDouble(); |
| 352 |
> |
torque[1] = tokenizer.nextTokenAsDouble(); |
| 353 |
> |
torque[2] = tokenizer.nextTokenAsDouble(); |
| 354 |
> |
break; |
| 355 |
> |
} |
| 356 |
> |
default: { |
| 357 |
> |
sprintf(painCave.errMsg, |
| 358 |
> |
"RestReader Error: %s is an unrecognized type\n", type.c_str()); |
| 359 |
> |
painCave.isFatal = 1; |
| 360 |
> |
simError(); |
| 361 |
> |
break; |
| 362 |
> |
} |
| 363 |
> |
} |
| 364 |
> |
// keep the position in case we need it for a molecular restraint: |
| 365 |
> |
|
| 366 |
> |
all_pos_[index] = pos; |
| 367 |
|
|
| 368 |
< |
all_pos_[index] = pos; |
| 369 |
< |
|
| 370 |
< |
// is this io restrained? |
| 371 |
< |
GenericData* data = integrableObject->getPropertyByName("Restraint"); |
| 372 |
< |
ObjectRestraint* oRest; |
| 373 |
< |
|
| 374 |
< |
if (data != NULL) { |
| 375 |
< |
// make sure we can reinterpret the generic data as restraint data: |
| 349 |
< |
RestraintData* restData= dynamic_cast<RestraintData*>(data); |
| 350 |
< |
if (restData != NULL) { |
| 351 |
< |
// make sure we can reinterpet the restraint data as a pointer to |
| 368 |
> |
// is this io restrained? |
| 369 |
> |
GenericData* data = sd->getPropertyByName("Restraint"); |
| 370 |
> |
|
| 371 |
> |
if (data != NULL) { |
| 372 |
> |
// make sure we can reinterpret the generic data as restraint data: |
| 373 |
> |
RestraintData* restData= dynamic_cast<RestraintData*>(data); |
| 374 |
> |
if (restData != NULL) { |
| 375 |
> |
// make sure we can reinterpet the restraint data as a pointer to |
| 376 |
|
// an ObjectRestraint: |
| 377 |
< |
oRest = dynamic_cast<ObjectRestraint*>(restData->getData()); |
| 378 |
< |
if (oRest != NULL) { |
| 379 |
< |
if (integrableObject->isDirectional()) { |
| 380 |
< |
oRest->setReferenceStructure(pos, q.toRotationMatrix3()); |
| 381 |
< |
} else { |
| 382 |
< |
oRest->setReferenceStructure(pos); |
| 359 |
< |
} |
| 377 |
> |
ObjectRestraint* oRest = dynamic_cast<ObjectRestraint*>(restData->getData()); |
| 378 |
> |
if (oRest != NULL) { |
| 379 |
> |
if (sd->isDirectional()) { |
| 380 |
> |
oRest->setReferenceStructure(pos, q.toRotationMatrix3()); |
| 381 |
> |
} else { |
| 382 |
> |
oRest->setReferenceStructure(pos); |
| 383 |
|
} |
| 384 |
|
} |
| 385 |
|
} |
