src/io/RestraintReader.cpp

#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 "io/ReadWrite.hpp"
#include "utils/simError.h"

#ifdef IS_MPI
#include <mpi.h>
#include "brains/mpiSimulation.hpp"
#define TAKE_THIS_TAG_CHAR 0
#define TAKE_THIS_TAG_INT 1
#define TAKE_THIS_TAG_DOUBLE 2
#endif // is_mpi


RestraintReader :: RestraintReader( SimInfo* the_simnfo ){
  
  simnfo = the_simnfo;

  idealName = "idealCrystal.in";
  
  isScanned = false;

#ifdef IS_MPI
  if (worldRank == 0) {
#endif

  inIdealFile = fopen(idealName, "r");
  if(inIdealFile == NULL){
    sprintf(painCave.errMsg,
            "Cannot open file: %s\n", idealName);
    painCave.isFatal = 1;
    simError();
  }
  
  inIdealFileName = idealName;
#ifdef IS_MPI
  }
  strcpy( checkPointMsg, "Restraint file opened for reading successfully." );
  MPIcheckPoint();
#endif
  return;  
}

RestraintReader :: ~RestraintReader( ){
#ifdef IS_MPI
  if (worldRank == 0) {
#endif
  vector<fpos_t*>::iterator i;

  int error;
  error = fclose( inIdealFile );
  if( error ){
    sprintf( painCave.errMsg,
             "Error closing %s\n", inIdealFileName.c_str());
    simError();
  }

  for(i = framePos.begin(); i != framePos.end(); ++i)
    delete *i;
  framePos.clear();
  
#ifdef IS_MPI
  }
  strcpy( checkPointMsg, "Restraint file closed successfully." );
  MPIcheckPoint();
#endif

  return;
}


void RestraintReader :: readIdealCrystal(){

  int i;
  unsigned int j;

#ifdef IS_MPI
  int done, which_node, which_atom; // loop counter
#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;

  vector<StuntDouble*> integrableObjects;


#ifndef IS_MPI

  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 != simnfo->getTotIntegrableObjects() ){
    sprintf( painCave.errMsg,
             "RestraintReader error. %s nIntegrable, %d, "
             "does not match the meta-data file's nIntegrable, %d.\n",
             inIdealFileName.c_str(), nTotObjs, simnfo->getTotIntegrableObjects());
    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( i=0; i < simnfo->n_mol; i++){

    integrableObjects = (simnfo->molecules[i]).getIntegrableObjects();

    for(j = 0; j < integrableObjects.size(); j++){

      eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
      if(eof_test == NULL){
        sprintf(painCave.errMsg,
              "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();
      }
      
      parseErr = parseIdealLine( read_buffer, integrableObjects[j]);
      if( parseErr != NULL ){
        strcpy( painCave.errMsg, parseErr );
        painCave.isFatal = 1;
        simError();
      }
    }
  }

  // MPI Section of code..........
#else //IS_MPI

  // first thing first, suspend fatalities.
  painCave.isEventLoop = 1;

  int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
  int haveError;

  MPI_Status istatus;
  int *MolToProcMap = mpiSim->getMolToProcMap();
  int localIndex;
  int nCurObj;
  int nitems;

  nTotObjs = simnfo->getTotIntegrableObjects();
  haveError = 0;
  if (worldRank == 0) {
    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());
      haveError = 1;
      simError();
    }
    
    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, simnfo->getTotIntegrableObjects());
      haveError= 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());
      haveError = 1;
      simError();
    }
    
    for (i=0 ; i < mpiSim->getNMolGlobal(); i++) {
      which_node = MolToProcMap[i];
      if(which_node == 0){
        //molecules belong to master node
        
        localIndex = mpiSim->getGlobalToLocalMol(i);
        
        if(localIndex == -1) {
          strcpy(painCave.errMsg, "Molecule not found on node 0!");
          haveError = 1;
          simError();
        }
        
        integrableObjects = (simnfo->molecules[localIndex]).getIntegrableObjects();
        for(j=0; j < integrableObjects.size(); j++){
          
          eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
          if(eof_test == NULL){
            sprintf(painCave.errMsg,
                    "error in reading file %s\n"
                    "natoms  = %d; index = %d\n"
                    "error reading the line from the file.\n",
                    inIdealFileName.c_str(), nTotObjs, i );
            haveError= 1;
            simError();
          }
          
          if(haveError) nodeZeroError();
          
          parseIdealLine(read_buffer, integrableObjects[j]);
        }
      }
      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,
                    "error in reading file %s\n"
                    "natoms  = %d; index = %d\n"
                    "error reading the line from the file.\n",
                    inIdealFileName.c_str(), nTotObjs, i );
            haveError= 1;
            simError();
          }
          
          if(haveError) nodeZeroError();
          
          MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node,
                   TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
        }
      }
    }
  }
  else{
    //actions taken at slave nodes
    for (i=0 ; i < mpiSim->getNMolGlobal(); i++) {
      which_node = MolToProcMap[i];
      
      if(which_node == worldRank){
        //molecule with global index i belongs to this processor
        
        localIndex = mpiSim->getGlobalToLocalMol(i);
        
        if(localIndex == -1) {
          sprintf(painCave.errMsg, "Molecule not found on node %d\n", worldRank);
          haveError = 1;
          simError();
        }
        
        integrableObjects = (simnfo->molecules[localIndex]).getIntegrableObjects();        
        
        nCurObj = integrableObjects.size();
        
        MPI_Send(&nCurObj, 1, MPI_INT, 0,
                 TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
        
        for(j = 0; j < integrableObjects.size(); j++){
          
          MPI_Recv(read_buffer, BUFFERSIZE, MPI_CHAR, 0,
                   TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
          
          parseErr = parseIdealLine(read_buffer, integrableObjects[j]);
          
          if( parseErr != NULL ){
            strcpy( painCave.errMsg, parseErr );
            simError();
          }
        }
      }
    }
  }
#endif
}

char* RestraintReader::parseIdealLine(char* readLine, StuntDouble* sd){

  char *foo; // the pointer to the current string token

  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

  // set the string tokenizer

  foo = strtok(readLine, " ,;\t");

  // check the atom name to the current atom

  if( strcmp( foo, sd->getType() ) ){
    sprintf( painCave.errMsg,
             "RestraintReader error.  Does not"
             " match the meta-data atom %s.\n",
             sd->getType() );
    return strdup( painCave.errMsg );
  }

  // get the positions

  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading position x from %s\n",
             inIdealFileName.c_str());
    return strdup( painCave.errMsg );
  }
  pos[0] = atof( foo );

  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading position y from %s\n",
             inIdealFileName.c_str());
    return strdup( painCave.errMsg );
  }
  pos[1] = atof( foo );

  foo = strtok(NULL, " ,;\t");
  if(foo == NULL){
    sprintf( painCave.errMsg,
             "error in reading position z from %s\n",
             inIdealFileName.c_str());
    return strdup( painCave.errMsg );
  }
  pos[2] = atof( foo );

  // store the positions in the stuntdouble as generic data doubles
  DoubleGenericData* refPosX = new DoubleGenericData();
  refPosX->setID("refPosX");
  refPosX->setData(pos[0]);
  sd->addProperty(refPosX);
  
  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);

  // we don't need the velocities
  foo = strtok(NULL, " ,;\t");
  foo = strtok(NULL, " ,;\t");
  foo = strtok(NULL, " ,;\t");

  if (!sd->isDirectional())
    return NULL;

  // get the quaternions
  if( sd->isDirectional() ){

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
                     "error in reading quaternion 0 from %s\n",
                      inIdealFileName.c_str() );
      return strdup( painCave.errMsg );
    }
    q[0] = atof( foo );

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
                     "error in reading quaternion 1 from %s\n",
                      inIdealFileName.c_str() );
      return strdup( painCave.errMsg );
    }
    q[1] = atof( foo );

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
                     "error in reading quaternion 2 from %s\n",
                      inIdealFileName.c_str() );
      return strdup( painCave.errMsg );
    }
    q[2] = atof( foo );

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      sprintf( painCave.errMsg,
                     "error in reading quaternion 3 from %s\n",
                      inIdealFileName.c_str() );
      return strdup( painCave.errMsg );
    }
    q[3] = atof( foo );

    // 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;
}

#ifdef IS_MPI
void RestraintReader::nodeZeroError( void ){
  int j, myStatus;

  myStatus = 0;
  for (j = 0; j < mpiSim->getNProcessors(); j++) {
    MPI_Send( &myStatus, 1, MPI_INT, j,
              TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
  }


  MPI_Finalize();
  exit (0);

}

void RestraintReader::anonymousNodeDie( void ){

  MPI_Finalize();
  exit (0);
}
#endif

void RestraintReader::readZangle( const char *in_name ){

  int i;
  unsigned int j;

#ifdef IS_MPI
  int done, which_node, which_atom; // loop counter

  int *potatoes;
  int myPotato;

  int nProc;
#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;

  vector<StuntDouble*> vecParticles;
  vector<double> tempZangs;

  // open the omega value file for reading
#ifdef IS_MPI
  if (worldRank == 0) {
#endif

  inAngFile = fopen(in_name, "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",
            in_name);
    painCave.severity = OOPSE_WARNING;
    painCave.isFatal = 0;
    simError();   
    return;
  }

  inAngFileName = in_name;
#ifdef IS_MPI
  }
  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 = simnfo->getTotIntegrableObjects();

  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
  vecParticles = simnfo->integrableObjects;

  for ( i=0; i<vecParticles.size(); i++ ) {
    vecParticles[i]->setZangle(tempZangs[i]);
  }

  // MPI Section of code..........
#else //IS_MPI

  // first thing first, suspend fatalities.
  painCave.isEventLoop = 1;

  int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
  int haveError, index;

  MPI_Status istatus;
  int *MolToProcMap = mpiSim->getMolToProcMap();
  int localIndex;
  int nCurObj;
  double angleTranfer;

  nTotObjs = simnfo->getTotIntegrableObjects();
  haveError = 0;
  if (worldRank == 0) {

    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 < mpiSim->getNMolGlobal(); i++) {
    // Get the Node number which has this atom
    which_node = MolToProcMap[i];
    
    if (worldRank == 0) {
      if (which_node == 0) {
        localIndex = mpiSim->getGlobalToLocalMol(i);
        
        if(localIndex == -1) {
          strcpy(painCave.errMsg, "Molecule not found on node 0!");
          haveError = 1;
          simError();
        }

        vecParticles = (simnfo->molecules[localIndex]).getIntegrableObjects();  
        for(j = 0; j < vecParticles.size(); j++){         
          vecParticles[j]->setZangle(tempZangs[index]);
          index++;
        }       
        
        // restraints is limited to a single zAngle per molecule
        vecParticles = (simnfo->molecules[localIndex]).getIntegrableObjects();
        for(j=0; j < vecParticles.size(); j++)
          vecParticles[j]->setZangle(tempZangs[i]);
        
      } 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

      if (which_node == worldRank) {

        // BUT I OWN THIS MOLECULE!!!

        localIndex = mpiSim->getGlobalToLocalMol(i);
        vecParticles = (simnfo->molecules[localIndex]).getIntegrableObjects();  
        nCurObj = vecParticles.size();
        
        MPI_Send(&nCurObj, 1, MPI_INT, 0,
                        TAKE_THIS_TAG_INT, MPI_COMM_WORLD);

        for(j = 0; j < vecParticles.size(); j++){
          
          MPI_Recv(&angleTransfer, 1, MPI_DOUBLE, 0,
                   TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD, &istatus);
          vecParticles[j]->setZangle(angleTransfer);
        }       
      }
    }
  }
#endif
}

void RestraintReader :: zeroZangle(){

  int i;
  unsigned int j;

  vector<StuntDouble*> vecParticles;

#ifndef IS_MPI
  // set all zAngles to 0.0
  vecParticles = simnfo->integrableObjects;

  for ( i=0; i<vecParticles.size(); i++ ) {
    vecParticles[i]->setZangle( 0.0 );
  }

  // MPI Section of code..........
#else //IS_MPI

  // first thing first, suspend fatalities.
  painCave.isEventLoop = 1;

  int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
  int haveError;

  MPI_Status istatus;
  int *MolToProcMap = mpiSim->getMolToProcMap();
  int localIndex;
  int which_node;

  haveError = 0;

  for (i=0 ; i < mpiSim->getNMolGlobal(); i++) {
    which_node = MolToProcMap[i];
    
    if(which_node == worldRank){
      //molecule with global index i belongs to this processor
      
      localIndex = mpiSim->getGlobalToLocalMol(i);
      
      if(localIndex == -1) {
        sprintf(painCave.errMsg, "Molecule not found on node %d\n", worldRank);
        haveError = 1;
        simError();
      }
      
      vecParticles = (simnfo->molecules[localIndex]).getIntegrableObjects();
      
      // set zAngle to 0.0 for all integrable objects
      for(j = 0; j < vecParticles.size(); j++){
        vecParticles[j]->setZangle( 0.0 );
      }
    }
  }
#endif
}
Revision:	1772
Committed:	Tue Nov 23 22:48:31 2004 UTC (19 years, 7 months ago) by chrisfen
File size:	19981 byte(s)
Log Message:	Improvements to restraints
#	User	Rev	Content
1	chrisfen	1772	#define _LARGEFILE_SOURCE64
2			#define _FILE_OFFSET_BITS 64
3
4			#include <sys/types.h>
5			#include <sys/stat.h>
6
7			#include <iostream>
8			#include <math.h>
9
10			#include <stdio.h>
11			#include <stdlib.h>
12			#include <string.h>
13
14
15			#include "io/ReadWrite.hpp"
16			#include "utils/simError.h"
17
18			#ifdef IS_MPI
19			#include <mpi.h>
20			#include "brains/mpiSimulation.hpp"
21			#define TAKE_THIS_TAG_CHAR 0
22			#define TAKE_THIS_TAG_INT 1
23			#define TAKE_THIS_TAG_DOUBLE 2
24			#endif // is_mpi
25
26
27			RestraintReader :: RestraintReader( SimInfo* the_simnfo ){
28
29			simnfo = the_simnfo;
30
31			idealName = "idealCrystal.in";
32
33			isScanned = false;
34
35			#ifdef IS_MPI
36			if (worldRank == 0) {
37			#endif
38
39			inIdealFile = fopen(idealName, "r");
40			if(inIdealFile == NULL){
41			sprintf(painCave.errMsg,
42			"Cannot open file: %s\n", idealName);
43			painCave.isFatal = 1;
44			simError();
45			}
46
47			inIdealFileName = idealName;
48			#ifdef IS_MPI
49			}
50			strcpy( checkPointMsg, "Restraint file opened for reading successfully." );
51			MPIcheckPoint();
52			#endif
53			return;
54			}
55
56			RestraintReader :: ~RestraintReader( ){
57			#ifdef IS_MPI
58			if (worldRank == 0) {
59			#endif
60			vector<fpos_t*>::iterator i;
61
62			int error;
63			error = fclose( inIdealFile );
64			if( error ){
65			sprintf( painCave.errMsg,
66			"Error closing %s\n", inIdealFileName.c_str());
67			simError();
68			}
69
70			for(i = framePos.begin(); i != framePos.end(); ++i)
71			delete *i;
72			framePos.clear();
73
74			#ifdef IS_MPI
75			}
76			strcpy( checkPointMsg, "Restraint file closed successfully." );
77			MPIcheckPoint();
78			#endif
79
80			return;
81			}
82
83
84			void RestraintReader :: readIdealCrystal(){
85
86			int i;
87			unsigned int j;
88
89			#ifdef IS_MPI
90			int done, which_node, which_atom; // loop counter
91			#endif //is_mpi
92
93			const int BUFFERSIZE = 2000; // size of the read buffer
94			int nTotObjs; // the number of atoms
95			char read_buffer[BUFFERSIZE]; //the line buffer for reading
96
97			char *eof_test; // ptr to see when we reach the end of the file
98			char *parseErr;
99
100			vector<StuntDouble*> integrableObjects;
101
102
103			#ifndef IS_MPI
104
105			eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
106			if( eof_test == NULL ){
107			sprintf( painCave.errMsg,
108			"RestraintReader error: error reading 1st line of \"%s\"\n",
109			inIdealFileName.c_str() );
110			painCave.isFatal = 1;
111			simError();
112			}
113
114			nTotObjs = atoi( read_buffer );
115
116			if( nTotObjs != simnfo->getTotIntegrableObjects() ){
117			sprintf( painCave.errMsg,
118			"RestraintReader error. %s nIntegrable, %d, "
119			"does not match the meta-data file's nIntegrable, %d.\n",
120			inIdealFileName.c_str(), nTotObjs, simnfo->getTotIntegrableObjects());
121			painCave.isFatal = 1;
122			simError();
123			}
124
125			// skip over the comment line
126			eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
127			if(eof_test == NULL){
128			sprintf( painCave.errMsg,
129			"error in reading commment in %s\n", inIdealFileName.c_str());
130			painCave.isFatal = 1;
131			simError();
132			}
133
134			// parse the ideal crystal lines
135			/*
136			Note: we assume that there is a one-to-one correspondence between
137			integrable objects and lines in the idealCrystal.in file. Thermodynamic
138			integration is only supported for simple rigid bodies.
139			*/
140			for( i=0; i < simnfo->n_mol; i++){
141
142			integrableObjects = (simnfo->molecules[i]).getIntegrableObjects();
143
144			for(j = 0; j < integrableObjects.size(); j++){
145
146			eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
147			if(eof_test == NULL){
148			sprintf(painCave.errMsg,
149			"error in reading file %s\n"
150			"natoms = %d; index = %d\n"
151			"error reading the line from the file.\n",
152			inIdealFileName.c_str(), nTotObjs, i );
153			painCave.isFatal = 1;
154			simError();
155			}
156
157			parseErr = parseIdealLine( read_buffer, integrableObjects[j]);
158			if( parseErr != NULL ){
159			strcpy( painCave.errMsg, parseErr );
160			painCave.isFatal = 1;
161			simError();
162			}
163			}
164			}
165
166			// MPI Section of code..........
167			#else //IS_MPI
168
169			// first thing first, suspend fatalities.
170			painCave.isEventLoop = 1;
171
172			int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
173			int haveError;
174
175			MPI_Status istatus;
176			int *MolToProcMap = mpiSim->getMolToProcMap();
177			int localIndex;
178			int nCurObj;
179			int nitems;
180
181			nTotObjs = simnfo->getTotIntegrableObjects();
182			haveError = 0;
183			if (worldRank == 0) {
184			eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
185			if( eof_test == NULL ){
186			sprintf( painCave.errMsg,
187			"Error reading 1st line of %s \n ",inIdealFileName.c_str());
188			haveError = 1;
189			simError();
190			}
191
192			nitems = atoi( read_buffer );
193
194			// Check to see that the number of integrable objects in the
195			// intial configuration file is the same as derived from the
196			// meta-data file.
197			if( nTotObjs != nitems){
198			sprintf( painCave.errMsg,
199			"RestraintReader Error. %s nIntegrable, %d, "
200			"does not match the meta-data file's nIntegrable, %d.\n",
201			inIdealFileName.c_str(), nTotObjs, simnfo->getTotIntegrableObjects());
202			haveError= 1;
203			simError();
204			}
205
206			// skip over the comment line
207			eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
208			if(eof_test == NULL){
209			sprintf( painCave.errMsg,
210			"error in reading commment in %s\n", inIdealFileName.c_str());
211			haveError = 1;
212			simError();
213			}
214
215			for (i=0 ; i < mpiSim->getNMolGlobal(); i++) {
216			which_node = MolToProcMap[i];
217			if(which_node == 0){
218			//molecules belong to master node
219
220			localIndex = mpiSim->getGlobalToLocalMol(i);
221
222			if(localIndex == -1) {
223			strcpy(painCave.errMsg, "Molecule not found on node 0!");
224			haveError = 1;
225			simError();
226			}
227
228			integrableObjects = (simnfo->molecules[localIndex]).getIntegrableObjects();
229			for(j=0; j < integrableObjects.size(); j++){
230
231			eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
232			if(eof_test == NULL){
233			sprintf(painCave.errMsg,
234			"error in reading file %s\n"
235			"natoms = %d; index = %d\n"
236			"error reading the line from the file.\n",
237			inIdealFileName.c_str(), nTotObjs, i );
238			haveError= 1;
239			simError();
240			}
241
242			if(haveError) nodeZeroError();
243
244			parseIdealLine(read_buffer, integrableObjects[j]);
245			}
246			}
247			else{
248			//molecule belongs to slave nodes
249
250			MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
251			TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
252
253			for(j=0; j < nCurObj; j++){
254
255			eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
256			if(eof_test == NULL){
257			sprintf(painCave.errMsg,
258			"error in reading file %s\n"
259			"natoms = %d; index = %d\n"
260			"error reading the line from the file.\n",
261			inIdealFileName.c_str(), nTotObjs, i );
262			haveError= 1;
263			simError();
264			}
265
266			if(haveError) nodeZeroError();
267
268			MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node,
269			TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
270			}
271			}
272			}
273			}
274			else{
275			//actions taken at slave nodes
276			for (i=0 ; i < mpiSim->getNMolGlobal(); i++) {
277			which_node = MolToProcMap[i];
278
279			if(which_node == worldRank){
280			//molecule with global index i belongs to this processor
281
282			localIndex = mpiSim->getGlobalToLocalMol(i);
283
284			if(localIndex == -1) {
285			sprintf(painCave.errMsg, "Molecule not found on node %d\n", worldRank);
286			haveError = 1;
287			simError();
288			}
289
290			integrableObjects = (simnfo->molecules[localIndex]).getIntegrableObjects();
291
292			nCurObj = integrableObjects.size();
293
294			MPI_Send(&nCurObj, 1, MPI_INT, 0,
295			TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
296
297			for(j = 0; j < integrableObjects.size(); j++){
298
299			MPI_Recv(read_buffer, BUFFERSIZE, MPI_CHAR, 0,
300			TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
301
302			parseErr = parseIdealLine(read_buffer, integrableObjects[j]);
303
304			if( parseErr != NULL ){
305			strcpy( painCave.errMsg, parseErr );
306			simError();
307			}
308			}
309			}
310			}
311			}
312			#endif
313			}
314
315			char* RestraintReader::parseIdealLine(char* readLine, StuntDouble* sd){
316
317			char *foo; // the pointer to the current string token
318
319			double pos[3]; // position place holders
320			double q[4]; // the quaternions
321			double RfromQ[3][3]; // the rotation matrix
322			double normalize; // to normalize the reference unit vector
323			double uX, uY, uZ; // reference unit vector place holders
324
325			// set the string tokenizer
326
327			foo = strtok(readLine, " ,;\t");
328
329			// check the atom name to the current atom
330
331			if( strcmp( foo, sd->getType() ) ){
332			sprintf( painCave.errMsg,
333			"RestraintReader error. Does not"
334			" match the meta-data atom %s.\n",
335			sd->getType() );
336			return strdup( painCave.errMsg );
337			}
338
339			// get the positions
340
341			foo = strtok(NULL, " ,;\t");
342			if(foo == NULL){
343			sprintf( painCave.errMsg,
344			"error in reading position x from %s\n",
345			inIdealFileName.c_str());
346			return strdup( painCave.errMsg );
347			}
348			pos[0] = atof( foo );
349
350			foo = strtok(NULL, " ,;\t");
351			if(foo == NULL){
352			sprintf( painCave.errMsg,
353			"error in reading position y from %s\n",
354			inIdealFileName.c_str());
355			return strdup( painCave.errMsg );
356			}
357			pos[1] = atof( foo );
358
359			foo = strtok(NULL, " ,;\t");
360			if(foo == NULL){
361			sprintf( painCave.errMsg,
362			"error in reading position z from %s\n",
363			inIdealFileName.c_str());
364			return strdup( painCave.errMsg );
365			}
366			pos[2] = atof( foo );
367
368			// store the positions in the stuntdouble as generic data doubles
369			DoubleGenericData* refPosX = new DoubleGenericData();
370			refPosX->setID("refPosX");
371			refPosX->setData(pos[0]);
372			sd->addProperty(refPosX);
373
374			DoubleGenericData* refPosY = new DoubleGenericData();
375			refPosY->setID("refPosY");
376			refPosY->setData(pos[1]);
377			sd->addProperty(refPosY);
378
379			DoubleGenericData* refPosZ = new DoubleGenericData();
380			refPosZ->setID("refPosZ");
381			refPosZ->setData(pos[2]);
382			sd->addProperty(refPosZ);
383
384			// we don't need the velocities
385			foo = strtok(NULL, " ,;\t");
386			foo = strtok(NULL, " ,;\t");
387			foo = strtok(NULL, " ,;\t");
388
389			if (!sd->isDirectional())
390			return NULL;
391
392			// get the quaternions
393			if( sd->isDirectional() ){
394
395			foo = strtok(NULL, " ,;\t");
396			if(foo == NULL){
397			sprintf( painCave.errMsg,
398			"error in reading quaternion 0 from %s\n",
399			inIdealFileName.c_str() );
400			return strdup( painCave.errMsg );
401			}
402			q[0] = atof( foo );
403
404			foo = strtok(NULL, " ,;\t");
405			if(foo == NULL){
406			sprintf( painCave.errMsg,
407			"error in reading quaternion 1 from %s\n",
408			inIdealFileName.c_str() );
409			return strdup( painCave.errMsg );
410			}
411			q[1] = atof( foo );
412
413			foo = strtok(NULL, " ,;\t");
414			if(foo == NULL){
415			sprintf( painCave.errMsg,
416			"error in reading quaternion 2 from %s\n",
417			inIdealFileName.c_str() );
418			return strdup( painCave.errMsg );
419			}
420			q[2] = atof( foo );
421
422			foo = strtok(NULL, " ,;\t");
423			if(foo == NULL){
424			sprintf( painCave.errMsg,
425			"error in reading quaternion 3 from %s\n",
426			inIdealFileName.c_str() );
427			return strdup( painCave.errMsg );
428			}
429			q[3] = atof( foo );
430
431			// now build the rotation matrix and find the unit vectors
432			RfromQ[0][0] = q[0]q[0] + q[1]q[1] - q[2]q[2] - q[3]q[3];
433			RfromQ[0][1] = 2(q[1]q[2] + q[0]*q[3]);
434			RfromQ[0][2] = 2(q[1]q[3] - q[0]*q[2]);
435			RfromQ[1][0] = 2(q[1]q[2] - q[0]*q[3]);
436			RfromQ[1][1] = q[0]q[0] - q[1]q[1] + q[2]q[2] - q[3]q[3];
437			RfromQ[1][2] = 2(q[2]q[3] + q[0]*q[1]);
438			RfromQ[2][0] = 2(q[1]q[3] + q[0]*q[2]);
439			RfromQ[2][1] = 2(q[2]q[3] - q[0]*q[1]);
440			RfromQ[2][2] = q[0]q[0] - q[1]q[1] - q[2]q[2] + q[3]q[3];
441
442			normalize = sqrt(RfromQ[2][0]RfromQ[2][0] + RfromQ[2][1]RfromQ[2][1]
443			+ RfromQ[2][2]*RfromQ[2][2]);
444			uX = RfromQ[2][0]/normalize;
445			uY = RfromQ[2][1]/normalize;
446			uZ = RfromQ[2][2]/normalize;
447
448			// store reference unit vectors as generic data in the stuntdouble
449			DoubleGenericData* refVectorX = new DoubleGenericData();
450			refVectorX->setID("refVectorX");
451			refVectorX->setData(uX);
452			sd->addProperty(refVectorX);
453
454			DoubleGenericData* refVectorY = new DoubleGenericData();
455			refVectorY->setID("refVectorY");
456			refVectorY->setData(uY);
457			sd->addProperty(refVectorY);
458
459			DoubleGenericData* refVectorZ = new DoubleGenericData();
460			refVectorZ->setID("refVectorZ");
461			refVectorZ->setData(uZ);
462			sd->addProperty(refVectorZ);
463			}
464
465			// we don't need the angular velocities, so let's exit the line
466			return NULL;
467			}
468
469			#ifdef IS_MPI
470			void RestraintReader::nodeZeroError( void ){
471			int j, myStatus;
472
473			myStatus = 0;
474			for (j = 0; j < mpiSim->getNProcessors(); j++) {
475			MPI_Send( &myStatus, 1, MPI_INT, j,
476			TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
477			}
478
479
480			MPI_Finalize();
481			exit (0);
482
483			}
484
485			void RestraintReader::anonymousNodeDie( void ){
486
487			MPI_Finalize();
488			exit (0);
489			}
490			#endif
491
492			void RestraintReader::readZangle( const char *in_name ){
493
494			int i;
495			unsigned int j;
496
497			#ifdef IS_MPI
498			int done, which_node, which_atom; // loop counter
499
500			int *potatoes;
501			int myPotato;
502
503			int nProc;
504			#endif //is_mpi
505
506			const int BUFFERSIZE = 2000; // size of the read buffer
507			int nTotObjs; // the number of atoms
508			char read_buffer[BUFFERSIZE]; //the line buffer for reading
509
510			char *eof_test; // ptr to see when we reach the end of the file
511			char *parseErr;
512
513			vector<StuntDouble*> vecParticles;
514			vector<double> tempZangs;
515
516			// open the omega value file for reading
517			#ifdef IS_MPI
518			if (worldRank == 0) {
519			#endif
520
521			inAngFile = fopen(in_name, "r");
522			if(!inAngFile){
523			sprintf(painCave.errMsg,
524			"Restraints Warning: %s file is not present\n"
525			"\tAll omega values will be initialized to zero. If the\n"
526			"\tsimulation is starting from the idealCrystal.in reference\n"
527			"\tconfiguration, this is the desired action. If this is not\n"
528			"\tthe case, the energy calculations will be incorrect.\n",
529			in_name);
530			painCave.severity = OOPSE_WARNING;
531			painCave.isFatal = 0;
532			simError();
533			return;
534			}
535
536			inAngFileName = in_name;
537			#ifdef IS_MPI
538			}
539			strcpy( checkPointMsg, "zAngle file opened successfully for reading." );
540			MPIcheckPoint();
541			#endif
542
543			#ifndef IS_MPI
544
545			eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
546			if( eof_test == NULL ){
547			sprintf( painCave.errMsg,
548			"RestraintReader error: error reading 1st line of \"%s\"\n",
549			inAngFileName.c_str() );
550			painCave.isFatal = 1;
551			simError();
552			}
553
554			eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
555			while ( eof_test != NULL ) {
556			// check for and ignore blank lines
557			if ( read_buffer != NULL )
558			tempZangs.push_back( atof(read_buffer) );
559			eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
560			}
561
562			nTotObjs = simnfo->getTotIntegrableObjects();
563
564			if( nTotObjs != tempZangs.size() ){
565			sprintf( painCave.errMsg,
566			"RestraintReader zAngle reading error. %s nIntegrable, %d, "
567			"does not match the meta-data file's nIntegrable, %d.\n",
568			inAngFileName.c_str(), tempZangs.size(), nTotObjs );
569			painCave.isFatal = 1;
570			simError();
571			}
572
573			// load the zAngles into the integrable objects
574			vecParticles = simnfo->integrableObjects;
575
576			for ( i=0; i<vecParticles.size(); i++ ) {
577			vecParticles[i]->setZangle(tempZangs[i]);
578			}
579
580			// MPI Section of code..........
581			#else //IS_MPI
582
583			// first thing first, suspend fatalities.
584			painCave.isEventLoop = 1;
585
586			int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
587			int haveError, index;
588
589			MPI_Status istatus;
590			int *MolToProcMap = mpiSim->getMolToProcMap();
591			int localIndex;
592			int nCurObj;
593			double angleTranfer;
594
595			nTotObjs = simnfo->getTotIntegrableObjects();
596			haveError = 0;
597			if (worldRank == 0) {
598
599			eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
600			if( eof_test == NULL ){
601			sprintf( painCave.errMsg,
602			"Error reading 1st line of %s \n ",inAngFileName.c_str());
603			haveError = 1;
604			simError();
605			}
606
607			// let node 0 load the temporary angle vector
608			eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
609			while ( eof_test != NULL ) {
610			// check for and ignore blank lines
611			if ( read_buffer != NULL )
612			tempZangs.push_back( atof(read_buffer) );
613			eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
614			}
615
616			// Check to see that the number of integrable objects in the
617			// intial configuration file is the same as derived from the
618			// meta-data file.
619			if( nTotObjs != tempZangs.size() ){
620			sprintf( painCave.errMsg,
621			"RestraintReader zAngle reading Error. %s nIntegrable, %d, "
622			"does not match the meta-data file's nIntegrable, %d.\n",
623			inAngFileName.c_str(), tempZangs.size(), nTotObjs);
624			haveError= 1;
625			simError();
626			}
627
628			}
629			// At this point, node 0 has a tempZangs vector completed, and
630			// everyone else has nada
631			index = 0;
632
633			for (i=0 ; i < mpiSim->getNMolGlobal(); i++) {
634			// Get the Node number which has this atom
635			which_node = MolToProcMap[i];
636
637			if (worldRank == 0) {
638			if (which_node == 0) {
639			localIndex = mpiSim->getGlobalToLocalMol(i);
640
641			if(localIndex == -1) {
642			strcpy(painCave.errMsg, "Molecule not found on node 0!");
643			haveError = 1;
644			simError();
645			}
646
647			vecParticles = (simnfo->molecules[localIndex]).getIntegrableObjects();
648			for(j = 0; j < vecParticles.size(); j++){
649			vecParticles[j]->setZangle(tempZangs[index]);
650			index++;
651			}
652
653			// restraints is limited to a single zAngle per molecule
654			vecParticles = (simnfo->molecules[localIndex]).getIntegrableObjects();
655			for(j=0; j < vecParticles.size(); j++)
656			vecParticles[j]->setZangle(tempZangs[i]);
657
658			} else {
659			// I am MASTER OF THE UNIVERSE, but I don't own this molecule
660
661			MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
662			TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
663
664			for(j=0; j < nCurObj; j++){
665			angleTransfer = tempZangs[index];
666			MPI_Send(&angleTransfer, 1, MPI_DOUBLE, which_node,
667			TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD);
668			index++;
669			}
670
671			}
672
673			} else {
674			// I am SLAVE TO THE MASTER
675
676			if (which_node == worldRank) {
677
678			// BUT I OWN THIS MOLECULE!!!
679
680			localIndex = mpiSim->getGlobalToLocalMol(i);
681			vecParticles = (simnfo->molecules[localIndex]).getIntegrableObjects();
682			nCurObj = vecParticles.size();
683
684			MPI_Send(&nCurObj, 1, MPI_INT, 0,
685			TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
686
687			for(j = 0; j < vecParticles.size(); j++){
688
689			MPI_Recv(&angleTransfer, 1, MPI_DOUBLE, 0,
690			TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD, &istatus);
691			vecParticles[j]->setZangle(angleTransfer);
692			}
693			}
694			}
695			}
696			#endif
697			}
698
699			void RestraintReader :: zeroZangle(){
700
701			int i;
702			unsigned int j;
703
704			vector<StuntDouble*> vecParticles;
705
706			#ifndef IS_MPI
707			// set all zAngles to 0.0
708			vecParticles = simnfo->integrableObjects;
709
710			for ( i=0; i<vecParticles.size(); i++ ) {
711			vecParticles[i]->setZangle( 0.0 );
712			}
713
714			// MPI Section of code..........
715			#else //IS_MPI
716
717			// first thing first, suspend fatalities.
718			painCave.isEventLoop = 1;
719
720			int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
721			int haveError;
722
723			MPI_Status istatus;
724			int *MolToProcMap = mpiSim->getMolToProcMap();
725			int localIndex;
726			int which_node;
727
728			haveError = 0;
729
730			for (i=0 ; i < mpiSim->getNMolGlobal(); i++) {
731			which_node = MolToProcMap[i];
732
733			if(which_node == worldRank){
734			//molecule with global index i belongs to this processor
735
736			localIndex = mpiSim->getGlobalToLocalMol(i);
737
738			if(localIndex == -1) {
739			sprintf(painCave.errMsg, "Molecule not found on node %d\n", worldRank);
740			haveError = 1;
741			simError();
742			}
743
744			vecParticles = (simnfo->molecules[localIndex]).getIntegrableObjects();
745
746			// set zAngle to 0.0 for all integrable objects
747			for(j = 0; j < vecParticles.size(); j++){
748			vecParticles[j]->setZangle( 0.0 );
749			}
750			}
751			}
752			#endif
753			}