--- trunk/src/io/RestReader.cpp 2006/06/07 18:05:19 985 +++ branches/development/src/io/RestReader.cpp 2013/02/20 15:39:39 1850 @@ -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,417 @@ * 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, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). + */ -#define _LARGEFILE_SOURCE64 -#define _FILE_OFFSET_BITS 64 - -#include -#include - -#include -#include - -#include -#include -#include - -#include "primitives/Molecule.hpp" -#include "utils/MemoryUtils.hpp" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#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 +#endif -#ifdef IS_MPI -#include -#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("")!= 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::COMM_WORLD.Bcast(&framePos_, 1, MPI::INT, 0); +#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("") != std::string::npos) { + break; + } } - MemoryUtils::deletePointers(framePos); + int sendBufferSize = sendBuffer.size(); + MPI::COMM_WORLD.Bcast(&sendBufferSize, 1, MPI::INT, masterNode); + MPI::COMM_WORLD.Bcast((void *)sendBuffer.c_str(), sendBufferSize, + MPI::CHAR, masterNode); -#ifdef IS_MPI - } - strcpy( checkPointMsg, "Restraint file closed successfully." ); - MPIcheckPoint(); + sstream.str(sendBuffer); + } else { + int sendBufferSize; + MPI::COMM_WORLD.Bcast(&sendBufferSize, 1, MPI::INT, masterNode); + char * recvBuffer = new char[sendBufferSize+1]; + assert(recvBuffer); + recvBuffer[sendBufferSize] = '\0'; + MPI::COMM_WORLD.Bcast(recvBuffer, sendBufferSize, MPI::CHAR, masterNode); + sstream.str(recvBuffer); + delete [] recvBuffer; + } + + std::istream& inputStream = sstream; #endif - return; - } - - - void RestReader :: readIdealCrystal(){ + inputStream.getline(buffer, bufferSize); + + line = buffer; + if (line.find("") == std::string::npos) { + sprintf(painCave.errMsg, + "RestReader Error: can not find \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("") == std::string::npos) { + sprintf(painCave.errMsg, + "RestReader Error: can not find \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 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(data); + + if (restData != NULL) { + + // make sure we can reinterpet the restraint data as a + // pointer to a MolecularRestraint: + + MolecularRestraint* mRest = dynamic_cast(restData->getData()); + + if (mRest != NULL) { + + // now we need to pack the stunt doubles for the reference + // structure: + + std::vector 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* sd = info_->getIOIndexToIntegrableObject(index); + + if (sd == 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 (sd->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 (sd->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 = sd->getPropertyByName("Restraint"); + ObjectRestraint* oRest; + + if (data != NULL) { + // make sure we can reinterpret the generic data as restraint data: + RestraintData* restData= dynamic_cast(data); + if (restData != NULL) { + // make sure we can reinterpet the restraint data as a pointer to + // an ObjectRestraint: + oRest = dynamic_cast(restData->getData()); + if (oRest != NULL) { + if (sd->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); - RealType pos[3]; // position place holders - RealType q[4]; // the quaternions - RealType RfromQ[3][3]; // the rotation matrix - RealType normalize; // to normalize the reference unit vector - RealType uX, uY, uZ; // reference unit vector place holders - RealType 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("") == std::string::npos) { + sprintf(painCave.errMsg, + "RestReader Error: Missing \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("") != 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("") == std::string::npos) { + sprintf(painCave.errMsg, + "RestReader Error: Missing \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 vecParticles; - std::vector 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("") != 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; - RealType 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 (which_node == 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_REALTYPE, 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_REALTYPE, 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 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; - RealType 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_REALTYPE, 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_REALTYPE, 0, - TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD, &istatus); - vecParticles[j]->setZangle(angleTransfer); - } - } - } - } -#endif - } - -} // end namespace oopse + +}//end namespace OpenMD