--- trunk/src/io/RestReader.cpp 2005/04/15 22:04:00 507
+++ branches/development/src/io/RestReader.cpp 2010/07/09 23:08:25 1465
@@ -1,24 +1,15 @@
-/*
- * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
- *
+/*
+ * Copyright (c) 2009 The University of Notre Dame. All Rights Reserved.
+ *
* The University of Notre Dame grants you ("Licensee") a
* non-exclusive, royalty free, license to use, modify and
* redistribute this software in source and binary code form, provided
* that the following conditions are met:
*
- * 1. Acknowledgement of the program authors must be made in any
- * publication of scientific results based in part on use of the
- * program. An acceptable form of acknowledgement is citation of
- * the article in which the program was described (Matthew
- * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
- * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
- * Parallel Simulation Engine for Molecular Dynamics,"
- * J. Comput. Chem. 26, pp. 252-271 (2005))
- *
- * 2. Redistributions of source code must retain the above copyright
+ * 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
*
- * 3. Redistributions in binary form must reproduce the above copyright
+ * 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the
* distribution.
@@ -37,799 +28,415 @@
* arising out of the use of or inability to use software, even if the
* University of Notre Dame has been advised of the possibility of
* such damages.
- */
+ *
+ * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
+ * research, please cite the appropriate papers when you publish your
+ * work. Good starting points are:
+ *
+ * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
+ * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
+ * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
+ * [4] Vardeman & Gezelter, in progress (2009).
+ */
-#define _LARGEFILE_SOURCE64
-#define _FILE_OFFSET_BITS 64
-
-#include <sys/types.h>
-#include <sys/stat.h>
-
-#include <iostream>
-#include <math.h>
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include "primitives/Molecule.hpp"
-#include "utils/MemoryUtils.hpp"
+
+#include <sys/types.h>
+#include <sys/stat.h>
+
+#include <math.h>
+#include <string>
+#include <sstream>
+#include <iostream>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "io/RestReader.hpp"
+#include "primitives/Molecule.hpp"
+#include "utils/simError.h"
+#include "utils/MemoryUtils.hpp"
#include "utils/StringTokenizer.hpp"
-#include "io/RestReader.hpp"
-#include "utils/simError.h"
+#include "restraints/ObjectRestraint.hpp"
+#include "restraints/MolecularRestraint.hpp"
+
+#ifdef IS_MPI
+
+#include <mpi.h>
+#endif
-#ifdef IS_MPI
-#include <mpi.h>
-#define TAKE_THIS_TAG_CHAR 0
-#define TAKE_THIS_TAG_INT 1
-#define TAKE_THIS_TAG_DOUBLE 2
-#endif // is_mpi
+namespace OpenMD {
-namespace oopse {
-
- RestReader::RestReader( SimInfo* info ) : info_(info){
+ void RestReader::scanFile(){
+ int lineNo = 0;
+ std::streampos prevPos;
+ std::streampos currPos;
- idealName = "idealCrystal.in";
+#ifdef IS_MPI
- isScanned = false;
-
-#ifdef IS_MPI
- if (worldRank == 0) {
-#endif
+ if (worldRank == 0) {
+#endif // is_mpi
+
+ inFile_->clear();
+ currPos = inFile_->tellg();
+ prevPos = currPos;
- inIdealFile = fopen(idealName, "r");
- if(inIdealFile == NULL){
- sprintf(painCave.errMsg,
- "RestReader: Cannot open file: %s\n", idealName);
- painCave.isFatal = 1;
- simError();
+ bool foundOpenSnapshotTag = false;
+
+ while(!foundOpenSnapshotTag && inFile_->getline(buffer, bufferSize)) {
+ ++lineNo;
+
+ std::string line = buffer;
+ currPos = inFile_->tellg();
+ if (line.find("<Snapshot>")!= std::string::npos) {
+ foundOpenSnapshotTag = true;
+ framePos_ = prevPos;
+ }
+ prevPos = currPos;
}
- inIdealFileName = idealName;
-#ifdef IS_MPI
+#ifdef IS_MPI
}
- strcpy( checkPointMsg,
- "File \"idealCrystal.in\" opened successfully for reading." );
- MPIcheckPoint();
-#endif
+ MPI_Bcast(&framePos_, 1, MPI_INT, 0, MPI_COMM_WORLD);
+#endif // is_mpi
+ }
- return;
- }
-
- RestReader :: ~RestReader( ){
-#ifdef IS_MPI
- if (worldRank == 0) {
-#endif
- int error;
- error = fclose( inIdealFile );
+
+ void RestReader::readSet(){
+ std::string line;
+
+#ifndef IS_MPI
+
+ inFile_->clear();
+ inFile_->seekg(framePos_);
+
+ std::istream& inputStream = *inFile_;
+#else
+
+ int masterNode = 0;
+ std::stringstream sstream;
+ if (worldRank == masterNode) {
+ std::string sendBuffer;
- if( error ){
- sprintf( painCave.errMsg,
- "Error closing %s\n", inIdealFileName.c_str());
- simError();
+ inFile_->clear();
+ inFile_->seekg(framePos_);
+
+ while (inFile_->getline(buffer, bufferSize)) {
+
+ line = buffer;
+ sendBuffer += line;
+ sendBuffer += '\n';
+ if (line.find("</Snapshot>") != std::string::npos) {
+ break;
+ }
}
- MemoryUtils::deletePointers(framePos);
+ int sendBufferSize = sendBuffer.size();
+ MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
+ MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
-#ifdef IS_MPI
- }
- strcpy( checkPointMsg, "Restraint file closed successfully." );
- MPIcheckPoint();
+ sstream.str(sendBuffer);
+ } else {
+ int sendBufferSize;
+ MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
+ char * recvBuffer = new char[sendBufferSize+1];
+ assert(recvBuffer);
+ recvBuffer[sendBufferSize] = '\0';
+ MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
+ sstream.str(recvBuffer);
+ delete [] recvBuffer;
+ }
+
+ std::istream& inputStream = sstream;
#endif
- return;
- }
-
-
- void RestReader :: readIdealCrystal(){
+ inputStream.getline(buffer, bufferSize);
+
+ line = buffer;
+ if (line.find("<Snapshot>") == std::string::npos) {
+ sprintf(painCave.errMsg,
+ "RestReader Error: can not find <Snapshot>\n");
+ painCave.isFatal = 1;
+ simError();
+ }
- int i;
- unsigned int j;
+ //read frameData
+ readFrameProperties(inputStream);
-#ifdef IS_MPI
- int done, which_node, which_atom; // loop counter
-#endif //is_mpi
+ //read StuntDoubles
+ readStuntDoubles(inputStream);
- const int BUFFERSIZE = 2000; // size of the read buffer
- int nTotObjs; // the number of atoms
- char read_buffer[BUFFERSIZE]; //the line buffer for reading
+ inputStream.getline(buffer, bufferSize);
+ line = buffer;
+ if (line.find("</Snapshot>") == std::string::npos) {
+ sprintf(painCave.errMsg,
+ "RestReader Error: can not find </Snapshot>\n");
+ painCave.isFatal = 1;
+ simError();
+ }
+ }
- char *eof_test; // ptr to see when we reach the end of the file
- char *parseErr;
+ void RestReader::readReferenceStructure() {
+
+ // We need temporary storage to keep track of all StuntDouble positions
+ // in case some of the restraints are molecular (i.e. if they use
+ // multiple SD positions to determine restrained orientations or positions:
+
+ all_pos_.clear();
+ all_pos_.resize(info_->getNGlobalIntegrableObjects()) ;
- std::vector<StuntDouble*> integrableObjects;
-
+ // Restraint files are just standard dump files, but with the reference
+ // structure stored in the first frame (frame 0).
+ // RestReader overloads readSet and explicitly handles all of the
+ // ObjectRestraints in that method:
+
+ scanFile();
+
+ readSet();
+
+
+ // all ObjectRestraints have been handled, now we have to worry about
+ // molecular restraints:
+
+ SimInfo::MoleculeIterator i;
+ Molecule::IntegrableObjectIterator j;
Molecule* mol;
- StuntDouble* integrableObject;
- SimInfo::MoleculeIterator mi;
- Molecule::IntegrableObjectIterator ii;
+ StuntDouble* sd;
-#ifndef IS_MPI
+ // no need to worry about parallel molecules, as molecules are not
+ // split across processor boundaries. Just loop over all molecules
+ // we know about:
- eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
- if( eof_test == NULL ){
- sprintf( painCave.errMsg,
- "RestraintReader error: error reading 1st line of \"%s\"\n",
- inIdealFileName.c_str() );
- painCave.isFatal = 1;
- simError();
- }
-
- nTotObjs = atoi( read_buffer );
-
- if( nTotObjs != info_->getNGlobalIntegrableObjects() ){
- sprintf( painCave.errMsg,
- "RestraintReader error. %s nIntegrable, %d, "
- "does not match the meta-data file's nIntegrable, %d.\n",
- inIdealFileName.c_str(), nTotObjs,
- info_->getNGlobalIntegrableObjects());
- painCave.isFatal = 1;
- simError();
- }
-
- // skip over the comment line
- eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
- if(eof_test == NULL){
- sprintf( painCave.errMsg,
- "error in reading commment in %s\n", inIdealFileName.c_str());
- painCave.isFatal = 1;
- simError();
- }
-
- // parse the ideal crystal lines
- /*
- * Note: we assume that there is a one-to-one correspondence between
- * integrable objects and lines in the idealCrystal.in file. Thermodynamic
- * integration is only supported for simple rigid bodies.
- */
- for (mol = info_->beginMolecule(mi); mol != NULL;
- mol = info_->nextMolecule(mi)) {
-
- for (integrableObject = mol->beginIntegrableObject(ii);
- integrableObject != NULL;
- integrableObject = mol->nextIntegrableObject(ii)) {
-
- eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
- if(eof_test == NULL){
- sprintf(painCave.errMsg,
- "RestReader Error: error in reading file %s\n"
- "natoms = %d; index = %d\n"
- "error reading the line from the file.\n",
- inIdealFileName.c_str(), nTotObjs,
- integrableObject->getGlobalIndex() );
- painCave.isFatal = 1;
- simError();
+ for (mol = info_->beginMolecule(i); mol != NULL;
+ mol = info_->nextMolecule(i)) {
+
+ // is this molecule restrained?
+ GenericData* data = mol->getPropertyByName("Restraint");
+
+ if (data != NULL) {
+
+ // make sure we can reinterpret the generic data as restraint data:
+
+ RestraintData* restData= dynamic_cast<RestraintData*>(data);
+
+ if (restData != NULL) {
+
+ // make sure we can reinterpet the restraint data as a
+ // pointer to a MolecularRestraint:
+
+ MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
+
+ if (mRest != NULL) {
+
+ // now we need to pack the stunt doubles for the reference
+ // structure:
+
+ std::vector<Vector3d> ref;
+ int count = 0;
+ RealType mass, mTot;
+ Vector3d COM(0.0);
+
+ mTot = 0.0;
+
+ // loop over the stunt doubles in this molecule in the order we
+ // will be looping them in the restraint code:
+
+ for (sd = mol->beginIntegrableObject(j); sd != NULL;
+ sd = mol->nextIntegrableObject(j)) {
+
+ // push back the reference positions of the stunt
+ // doubles from the *globally* sorted array of
+ // positions:
+
+ ref.push_back( all_pos_[sd->getGlobalIntegrableObjectIndex()] );
+ mass = sd->getMass();
+ COM = COM + mass * ref[count];
+ mTot = mTot + mass;
+ count = count + 1;
+ }
+ COM /= mTot;
+ mRest->setReferenceStructure(ref, COM);
+ }
}
-
- parseErr = parseIdealLine( read_buffer, integrableObject);
- if( parseErr != NULL ){
- strcpy( painCave.errMsg, parseErr );
- painCave.isFatal = 1;
- simError();
- }
}
}
+ }
+
+
+
+ void RestReader::parseDumpLine(const std::string& line) {
+
+ StringTokenizer tokenizer(line);
+ int nTokens;
+
+ nTokens = tokenizer.countTokens();
+
+ if (nTokens < 2) {
+ sprintf(painCave.errMsg,
+ "RestReader Error: Not enough Tokens.\n%s\n", line.c_str());
+ painCave.isFatal = 1;
+ simError();
+ }
+
+ int index = tokenizer.nextTokenAsInt();
+
+ StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
+
+ if (integrableObject == NULL) {
+ return;
+ }
+
+ std::string type = tokenizer.nextToken();
+ int size = type.size();
- // MPI Section of code..........
-#else //IS_MPI
-
- // first thing first, suspend fatalities.
- painCave.isEventLoop = 1;
-
- int masterNode = 0;
- int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
-
- MPI_Status istatus;
- int nCurObj;
- int nitems;
- int haveError;
+ Vector3d pos;
+ Quat4d q;
- nTotObjs = info_->getNGlobalIntegrableObjects();
- haveError = 0;
+ for(int i = 0; i < size; ++i) {
+ switch(type[i]) {
- if (worldRank == masterNode) {
- eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
- if( eof_test == NULL ){
- sprintf( painCave.errMsg,
- "Error reading 1st line of %s \n ",inIdealFileName.c_str());
- painCave.isFatal = 1;
- simError();
+ case 'p': {
+ pos[0] = tokenizer.nextTokenAsDouble();
+ pos[1] = tokenizer.nextTokenAsDouble();
+ pos[2] = tokenizer.nextTokenAsDouble();
+ break;
}
-
- nitems = atoi( read_buffer );
-
- // Check to see that the number of integrable objects in the
- // intial configuration file is the same as derived from the
- // meta-data file.
- if( nTotObjs != nitems){
- sprintf( painCave.errMsg,
- "RestraintReader Error. %s nIntegrable, %d, "
- "does not match the meta-data file's nIntegrable, %d.\n",
- inIdealFileName.c_str(), nTotObjs,
- info_->getNGlobalIntegrableObjects());
- painCave.isFatal = 1;
- simError();
+ case 'v' : {
+ Vector3d vel;
+ vel[0] = tokenizer.nextTokenAsDouble();
+ vel[1] = tokenizer.nextTokenAsDouble();
+ vel[2] = tokenizer.nextTokenAsDouble();
+ break;
}
-
- // skip over the comment line
- eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
- if(eof_test == NULL){
- sprintf( painCave.errMsg,
- "error in reading commment in %s\n", inIdealFileName.c_str());
- painCave.isFatal = 1;
- simError();
- }
- for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
- int which_node = info_->getMolToProc(i);
-
- if(which_node == masterNode){
- //molecules belong to master node
+ case 'q' : {
+ if (integrableObject->isDirectional()) {
- mol = info_->getMoleculeByGlobalIndex(i);
+ q[0] = tokenizer.nextTokenAsDouble();
+ q[1] = tokenizer.nextTokenAsDouble();
+ q[2] = tokenizer.nextTokenAsDouble();
+ q[3] = tokenizer.nextTokenAsDouble();
- if(mol == NULL) {
- sprintf(painCave.errMsg,
- "RestReader Error: Molecule not found on node %d!\n",
- worldRank);
- painCave.isFatal = 1;
- simError();
- }
-
- for (integrableObject = mol->beginIntegrableObject(ii);
- integrableObject != NULL;
- integrableObject = mol->nextIntegrableObject(ii)){
+ RealType qlen = q.length();
+ if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
- eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
-
- if(eof_test == NULL){
- sprintf(painCave.errMsg,
- "RestReader Error: error in reading file %s\n"
- "natoms = %d; index = %d\n"
- "error reading the line from the file.\n",
- inIdealFileName.c_str(), nTotObjs, i );
- painCave.isFatal = 1;
- simError();
- }
-
- parseIdealLine(read_buffer, integrableObject);
-
- }
- } else {
- //molecule belongs to slave nodes
-
- MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
- TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
-
- for(j=0; j < nCurObj; j++){
-
- eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
- if(eof_test == NULL){
- sprintf(painCave.errMsg,
- "RestReader Error: error in reading file %s\n"
- "natoms = %d; index = %d\n"
- "error reading the line from the file.\n",
- inIdealFileName.c_str(), nTotObjs, i );
- painCave.isFatal = 1;
- simError();
- }
-
- MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node,
- TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
- }
+ sprintf(painCave.errMsg,
+ "RestReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
+ painCave.isFatal = 1;
+ simError();
+ }
+
+ q.normalize();
+ }
+ break;
+ }
+ case 'j' : {
+ Vector3d ji;
+ if (integrableObject->isDirectional()) {
+ ji[0] = tokenizer.nextTokenAsDouble();
+ ji[1] = tokenizer.nextTokenAsDouble();
+ ji[2] = tokenizer.nextTokenAsDouble();
}
+ break;
+ }
+ case 'f': {
+ Vector3d force;
+ force[0] = tokenizer.nextTokenAsDouble();
+ force[1] = tokenizer.nextTokenAsDouble();
+ force[2] = tokenizer.nextTokenAsDouble();
+ break;
}
- } else {
- //actions taken at slave nodes
- for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
- int which_node = info_->getMolToProc(i);
+ case 't' : {
+ Vector3d torque;
+ torque[0] = tokenizer.nextTokenAsDouble();
+ torque[1] = tokenizer.nextTokenAsDouble();
+ torque[2] = tokenizer.nextTokenAsDouble();
+ break;
+ }
+ default: {
+ sprintf(painCave.errMsg,
+ "RestReader Error: %s is an unrecognized type\n", type.c_str());
+ painCave.isFatal = 1;
+ simError();
+ break;
+ }
+ }
+ // keep the position in case we need it for a molecular restraint:
+
+ all_pos_[index] = pos;
- if(which_node == worldRank){
- //molecule with global index i belongs to this processor
-
- mol = info_->getMoleculeByGlobalIndex(i);
-
- if(mol == NULL) {
- sprintf(painCave.errMsg,
- "RestReader Error: molecule not found on node %d\n",
- worldRank);
- painCave.isFatal = 1;
- simError();
- }
-
- nCurObj = mol->getNIntegrableObjects();
-
- MPI_Send(&nCurObj, 1, MPI_INT, masterNode,
- TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
-
- for (integrableObject = mol->beginIntegrableObject(ii);
- integrableObject != NULL;
- integrableObject = mol->nextIntegrableObject(ii)){
-
- MPI_Recv(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode,
- TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
-
- parseErr = parseIdealLine(read_buffer, integrableObject);
-
- if( parseErr != NULL ){
- strcpy( painCave.errMsg, parseErr );
- simError();
+ // is this io restrained?
+ GenericData* data = integrableObject->getPropertyByName("Restraint");
+ ObjectRestraint* oRest;
+
+ if (data != NULL) {
+ // make sure we can reinterpret the generic data as restraint data:
+ RestraintData* restData= dynamic_cast<RestraintData*>(data);
+ if (restData != NULL) {
+ // make sure we can reinterpet the restraint data as a pointer to
+ // an ObjectRestraint:
+ oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
+ if (oRest != NULL) {
+ if (integrableObject->isDirectional()) {
+ oRest->setReferenceStructure(pos, q.toRotationMatrix3());
+ } else {
+ oRest->setReferenceStructure(pos);
}
}
}
}
}
-#endif
}
- char* RestReader::parseIdealLine(char* readLine, StuntDouble* sd){
+ void RestReader::readStuntDoubles(std::istream& inputStream) {
- char *foo; // the pointer to the current string token
+ inputStream.getline(buffer, bufferSize);
+ std::string line(buffer);
- double pos[3]; // position place holders
- double q[4]; // the quaternions
- double RfromQ[3][3]; // the rotation matrix
- double normalize; // to normalize the reference unit vector
- double uX, uY, uZ; // reference unit vector place holders
- double uselessToken;
- StringTokenizer tokenizer(readLine);
- int nTokens;
-
- nTokens = tokenizer.countTokens();
-
- if (nTokens < 14) {
- sprintf(painCave.errMsg,
- "RestReader Error: Not enough Tokens.\n");
- painCave.isFatal = 1;
- simError();
+ if (line.find("<StuntDoubles>") == std::string::npos) {
+ sprintf(painCave.errMsg,
+ "RestReader Error: Missing <StuntDoubles>\n");
+ painCave.isFatal = 1;
+ simError();
}
- std::string name = tokenizer.nextToken();
-
- if (name != sd->getType()) {
+ while(inputStream.getline(buffer, bufferSize)) {
+ line = buffer;
- sprintf(painCave.errMsg,
- "RestReader Error: Atom type [%s] in %s does not "
- "match Atom Type [%s] in .md file.\n",
- name.c_str(), inIdealFileName.c_str(),
- sd->getType().c_str());
- painCave.isFatal = 1;
- simError();
+ if(line.find("</StuntDoubles>") != std::string::npos) {
+ break;
+ }
+
+ parseDumpLine(line);
}
- pos[0] = tokenizer.nextTokenAsDouble();
- pos[1] = tokenizer.nextTokenAsDouble();
- pos[2] = tokenizer.nextTokenAsDouble();
+ }
- // store the positions in the stuntdouble as generic data doubles
- DoubleGenericData* refPosX = new DoubleGenericData();
- refPosX->setID("refPosX");
- refPosX->setData(pos[0]);
- sd->addProperty(refPosX);
+
+ void RestReader::readFrameProperties(std::istream& inputStream) {
+ inputStream.getline(buffer, bufferSize);
+ std::string line(buffer);
- DoubleGenericData* refPosY = new DoubleGenericData();
- refPosY->setID("refPosY");
- refPosY->setData(pos[1]);
- sd->addProperty(refPosY);
-
- DoubleGenericData* refPosZ = new DoubleGenericData();
- refPosZ->setID("refPosZ");
- refPosZ->setData(pos[2]);
- sd->addProperty(refPosZ);
+ if (line.find("<FrameData>") == std::string::npos) {
+ sprintf(painCave.errMsg,
+ "RestReader Error: Missing <FrameData>\n");
+ painCave.isFatal = 1;
+ simError();
+ }
- // we don't need the velocities
- uselessToken = tokenizer.nextTokenAsDouble();
- uselessToken = tokenizer.nextTokenAsDouble();
- uselessToken = tokenizer.nextTokenAsDouble();
-
- if (sd->isDirectional()) {
+ // restraints don't care about frame data (unless we need to wrap
+ // coordinates, but we'll worry about that later), so
+ // we'll just scan ahead until the end of the frame data:
+
+ while(inputStream.getline(buffer, bufferSize)) {
+ line = buffer;
- q[0] = tokenizer.nextTokenAsDouble();
- q[1] = tokenizer.nextTokenAsDouble();
- q[2] = tokenizer.nextTokenAsDouble();
- q[3] = tokenizer.nextTokenAsDouble();
-
- // now build the rotation matrix and find the unit vectors
- RfromQ[0][0] = q[0]*q[0] + q[1]*q[1] - q[2]*q[2] - q[3]*q[3];
- RfromQ[0][1] = 2*(q[1]*q[2] + q[0]*q[3]);
- RfromQ[0][2] = 2*(q[1]*q[3] - q[0]*q[2]);
- RfromQ[1][0] = 2*(q[1]*q[2] - q[0]*q[3]);
- RfromQ[1][1] = q[0]*q[0] - q[1]*q[1] + q[2]*q[2] - q[3]*q[3];
- RfromQ[1][2] = 2*(q[2]*q[3] + q[0]*q[1]);
- RfromQ[2][0] = 2*(q[1]*q[3] + q[0]*q[2]);
- RfromQ[2][1] = 2*(q[2]*q[3] - q[0]*q[1]);
- RfromQ[2][2] = q[0]*q[0] - q[1]*q[1] - q[2]*q[2] + q[3]*q[3];
-
- normalize = sqrt(RfromQ[2][0]*RfromQ[2][0] + RfromQ[2][1]*RfromQ[2][1]
- + RfromQ[2][2]*RfromQ[2][2]);
- uX = RfromQ[2][0]/normalize;
- uY = RfromQ[2][1]/normalize;
- uZ = RfromQ[2][2]/normalize;
-
- // store reference unit vectors as generic data in the stuntdouble
- DoubleGenericData* refVectorX = new DoubleGenericData();
- refVectorX->setID("refVectorX");
- refVectorX->setData(uX);
- sd->addProperty(refVectorX);
-
- DoubleGenericData* refVectorY = new DoubleGenericData();
- refVectorY->setID("refVectorY");
- refVectorY->setData(uY);
- sd->addProperty(refVectorY);
-
- DoubleGenericData* refVectorZ = new DoubleGenericData();
- refVectorZ->setID("refVectorZ");
- refVectorZ->setData(uZ);
- sd->addProperty(refVectorZ);
- }
-
- // we don't need the angular velocities, so let's exit the line
- return NULL;
- }
-
- void RestReader::readZangle(){
-
- int i;
- unsigned int j;
- int isPresent;
-
- Molecule* mol;
- StuntDouble* integrableObject;
- SimInfo::MoleculeIterator mi;
- Molecule::IntegrableObjectIterator ii;
-
-#ifdef IS_MPI
- int done, which_node, which_atom; // loop counter
- int nProc;
- MPI_Status istatus;
-#endif //is_mpi
-
- const int BUFFERSIZE = 2000; // size of the read buffer
- int nTotObjs; // the number of atoms
- char read_buffer[BUFFERSIZE]; //the line buffer for reading
-
- char *eof_test; // ptr to see when we reach the end of the file
- char *parseErr;
-
- std::vector<StuntDouble*> vecParticles;
- std::vector<double> tempZangs;
-
- inAngFileName = info_->getRestFileName();
-
- inAngFileName += "0";
-
- // open the omega value file for reading
-#ifdef IS_MPI
- if (worldRank == 0) {
-#endif
- isPresent = 1;
- inAngFile = fopen(inAngFileName.c_str(), "r");
- if(!inAngFile){
- sprintf(painCave.errMsg,
- "Restraints Warning: %s file is not present\n"
- "\tAll omega values will be initialized to zero. If the\n"
- "\tsimulation is starting from the idealCrystal.in reference\n"
- "\tconfiguration, this is the desired action. If this is not\n"
- "\tthe case, the energy calculations will be incorrect.\n",
- inAngFileName.c_str());
- painCave.severity = OOPSE_WARNING;
- painCave.isFatal = 0;
- simError();
- isPresent = 0;
+ if(line.find("</FrameData>") != std::string::npos) {
+ break;
}
-#ifdef IS_MPI
- // let the other nodes know the status of the file search
- MPI_Bcast(&isPresent, 1, MPI_INT, 0, MPI_COMM_WORLD);
-#endif // is_mpi
-
- if (!isPresent) {
- zeroZangle();
- return;
- }
-
-#ifdef IS_MPI
}
- // listen to node 0 to see if we should exit this function
- if (worldRank != 0) {
- MPI_Bcast(&isPresent, 1, MPI_INT, 0, MPI_COMM_WORLD);
- if (!isPresent) {
- zeroZangle();
- return;
- }
- }
-
- strcpy( checkPointMsg, "zAngle file opened successfully for reading." );
- MPIcheckPoint();
-#endif
-
-#ifndef IS_MPI
-
- eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
- if( eof_test == NULL ){
- sprintf( painCave.errMsg,
- "RestraintReader error: error reading 1st line of \"%s\"\n",
- inAngFileName.c_str() );
- painCave.isFatal = 1;
- simError();
- }
-
- eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
- while ( eof_test != NULL ) {
- // check for and ignore blank lines
- if ( read_buffer != NULL )
- tempZangs.push_back( atof(read_buffer) );
- eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
- }
-
- nTotObjs = info_->getNGlobalIntegrableObjects();
-
- if( nTotObjs != tempZangs.size() ){
- sprintf( painCave.errMsg,
- "RestraintReader zAngle reading error. %s nIntegrable, %d, "
- "does not match the meta-data file's nIntegrable, %d.\n",
- inAngFileName.c_str(), tempZangs.size(), nTotObjs );
- painCave.isFatal = 1;
- simError();
- }
-
- // load the zAngles into the integrable objects
- i = 0;
- for (mol = info_->beginMolecule(mi); mol != NULL;
- mol = info_->nextMolecule(mi)) {
-
- for (integrableObject = mol->beginIntegrableObject(ii);
- integrableObject != NULL;
- integrableObject = mol->nextIntegrableObject(ii)) {
-
- integrableObject->setZangle(tempZangs[i]);
- i++;
- }
- }
-
- // MPI Section of code..........
-#else //IS_MPI
-
- // first thing first, suspend fatalities.
- painCave.isEventLoop = 1;
-
- int masterNode = 0;
- int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
- int haveError;
- int index;
-
- int nCurObj;
- double angleTranfer;
-
- nTotObjs = info_->getNGlobalIntegrableObjects();
- haveError = 0;
-
- if (worldRank == masterNode) {
-
- eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
- if( eof_test == NULL ){
- sprintf( painCave.errMsg,
- "Error reading 1st line of %s \n ",inAngFileName.c_str());
- haveError = 1;
- simError();
- }
-
- // let node 0 load the temporary angle vector
- eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
- while ( eof_test != NULL ) {
- // check for and ignore blank lines
- if ( read_buffer != NULL )
- tempZangs.push_back( atof(read_buffer) );
- eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
- }
-
- // Check to see that the number of integrable objects in the
- // intial configuration file is the same as derived from the
- // meta-data file.
- if( nTotObjs != tempZangs.size() ){
- sprintf( painCave.errMsg,
- "RestraintReader zAngle reading Error. %s nIntegrable, %d, "
- "does not match the meta-data file's nIntegrable, %d.\n",
- inAngFileName.c_str(), tempZangs.size(), nTotObjs);
- haveError= 1;
- simError();
- }
-
- // At this point, node 0 has a tempZangs vector completed, and
- // everyone else has nada
- index = 0;
-
- for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
- // Get the Node number which has this atom
- which_node = info_->getMolToProc(i);
-
- if (worldRank == masterNode) {
- mol = info_->getMoleculeByGlobalIndex(i);
-
- if(mol == NULL) {
- strcpy(painCave.errMsg, "Molecule not found on node 0!");
- haveError = 1;
- simError();
- }
-
- for (integrableObject = mol->beginIntegrableObject(ii);
- integrableObject != NULL;
- integrableObject = mol->nextIntegrableObject(ii)){
-
- integrableObject->setZangle(tempZangs[index]);
- index++;
- }
-
- } else {
- // I am MASTER OF THE UNIVERSE, but I don't own this molecule
-
- MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
- TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
-
- for(j=0; j < nCurObj; j++){
- angleTransfer = tempZangs[index];
- MPI_Send(&angleTransfer, 1, MPI_DOUBLE, which_node,
- TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD);
- index++;
- }
-
- }
- }
- } else {
- // I am SLAVE TO THE MASTER
- for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
- int which_node = info_->getMolToProc(i);
-
- if (which_node == worldRank) {
-
- // BUT I OWN THIS MOLECULE!!!
-
- mol = info_->getMoleculeByGlobalIndex(i);
-
- if(mol == NULL) {
- sprintf(painCave.errMsg,
- "RestReader Error: molecule not found on node %d\n",
- worldRank);
- painCave.isFatal = 1;
- simError();
- }
-
- nCurObj = mol->getNIntegrableObjects();
-
- MPI_Send(&nCurObj, 1, MPI_INT, 0,
- TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
-
- for (integrableObject = mol->beginIntegrableObject(ii);
- integrableObject != NULL;
- integrableObject = mol->nextIntegrableObject(ii)){
-
- MPI_Recv(&angleTransfer, 1, MPI_DOUBLE, 0,
- TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD, &istatus);
-
- integrableObject->setZangle(angleTransfer);
- }
- }
- }
- }
-#endif
}
-
- void RestReader :: zeroZangle(){
-
- int i;
- unsigned int j;
- int nTotObjs; // the number of atoms
-
- Molecule* mol;
- StuntDouble* integrableObject;
- SimInfo::MoleculeIterator mi;
- Molecule::IntegrableObjectIterator ii;
-
- std::vector<StuntDouble*> vecParticles;
-
-#ifndef IS_MPI
- // set all zAngles to 0.0
- for (mol = info_->beginMolecule(mi); mol != NULL;
- mol = info_->nextMolecule(mi))
-
- for (integrableObject = mol->beginIntegrableObject(ii);
- integrableObject != NULL;
- integrableObject = mol->nextIntegrableObject(ii))
- integrableObject->setZangle( 0.0 );
-
-
- // MPI Section of code..........
-#else //IS_MPI
-
- // first thing first, suspend fatalities.
- painCave.isEventLoop = 1;
-
- int masterNode = 0;
- int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
- int haveError;
- int which_node;
-
- MPI_Status istatus;
-
- int nCurObj;
- double angleTranfer;
-
- nTotObjs = info_->getNGlobalIntegrableObjects();
- haveError = 0;
- if (worldRank == masterNode) {
- for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
- // Get the Node number which has this atom
- which_node = info_->getMolToProc(i);
-
- // let's let node 0 pass out constant values to all the processors
- if (worldRank == masterNode) {
- mol = info_->getMoleculeByGlobalIndex(i);
-
- if(mol == NULL) {
- strcpy(painCave.errMsg, "Molecule not found on node 0!");
- haveError = 1;
- simError();
- }
-
- for (integrableObject = mol->beginIntegrableObject(ii);
- integrableObject != NULL;
- integrableObject = mol->nextIntegrableObject(ii)){
-
- integrableObject->setZangle( 0.0 );
-
- }
-
- } else {
- // I am MASTER OF THE UNIVERSE, but I don't own this molecule
-
- MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
- TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
-
- for(j=0; j < nCurObj; j++){
- angleTransfer = 0.0;
- MPI_Send(&angleTransfer, 1, MPI_DOUBLE, which_node,
- TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD);
-
- }
- }
- }
- } else {
- // I am SLAVE TO THE MASTER
- for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
- int which_node = info_->getMolToProc(i);
-
- if (which_node == worldRank) {
-
- // BUT I OWN THIS MOLECULE!!!
- mol = info_->getMoleculeByGlobalIndex(i);
-
- if(mol == NULL) {
- sprintf(painCave.errMsg,
- "RestReader Error: molecule not found on node %d\n",
- worldRank);
- painCave.isFatal = 1;
- simError();
- }
-
- nCurObj = mol->getNIntegrableObjects();
-
- MPI_Send(&nCurObj, 1, MPI_INT, 0,
- TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
-
- for (integrableObject = mol->beginIntegrableObject(ii);
- integrableObject != NULL;
- integrableObject = mol->nextIntegrableObject(ii)){
-
- MPI_Recv(&angleTransfer, 1, MPI_DOUBLE, 0,
- TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD, &istatus);
- vecParticles[j]->setZangle(angleTransfer);
- }
- }
- }
- }
-#endif
- }
-
-} // end namespace oopse
+
+}//end namespace OpenMD