mdtools/interface_implementation/TraPPE_ExFF.cpp

#include <cstdlib>
#include <cstdio>
#include <cstring>

#include <iostream>
using namespace std;

#include "ForceFields.hpp"
#include "SRI.hpp"


#ifdef IS_MPI

#include "mpiForceField.h"

int myNode;

// Declare the structures that will be passed by MPI

typedef struct{
  char name[15];
  double mass;
  double epslon;
  double sigma;
  double dipole;
  int isDipole;
  int last;      //  0  -> default
                 //  1  -> tells nodes to stop listening
                 // -1  -> an error has occured. (handled in mpiForceField)
} atomStruct;
MPI_Datatype mpiAtomStructType;

typedef struct{
  char nameA[15];
  char nameB[15];
  char type[30];
  double d0;
  int last;      //  0  -> default
                 //  1  -> tells nodes to stop listening
                 // -1  -> an error has occured. (handled in mpiForceField)
} bondStruct;
MPI_Datatype mpiBondStructType;

typedef struct{
  char nameA[15];
  char nameB[15];
  char nameC[15];
  char type[30];
  double k1, k2, k3, t0;
  int last;      //  0  -> default
                 //  1  -> tells nodes to stop listening
                 // -1  -> an error has occured. (handled in mpiForceField)
} bendStruct;
MPI_Datatype mpiBendStructType;

typedef struct{
  char nameA[15];
  char nameB[15];
  char nameC[15];
  char nameD[15];
  char type[30];
  double k1, k2, k3, k4;
  int last;      //  0  -> default
                 //  1  -> tells nodes to stop listening
                 // -1  -> an error has occured. (handled in mpiForceField)
} TorsionStruct;
MPI_Datatype mpiTorsionStructType;

#endif


TraPPE_ExFF::TraPPE_ExFF(){

  char fileName[200];
  char* ffPath_env = "FORCE_PARAM_PATH";
  char* ffPath;
  char temp[200];
  char errMsg[1000];

#ifdef IS_MPI
  int i;
  int mpiError;
  int cleanOpen = 0;
  
  mpiError = MPI_Comm_rank(MPI_COMM_WORLD,&myNode);
  
  // **********************************************************************
  // Init the atomStruct mpi type

  atomStruct atomProto; // mpiPrototype
  int atomBC[3] = {15,4,2};  // block counts
  MPI_Aint atomDspls[3];           // displacements
  MPI_Datatype atomMbrTypes[3];    // member mpi types

  
  MPI_Address(&atomProto.name, &atomDspls[0]);
  MPI_Address(&atomProto.mass, &atomDspls[1]);
  MPI_Address(&atomProto.isDipole, &atomDspls[2]);
  
  atomMbrTypes[0] = MPI_CHAR;
  atomMbrTypes[1] = MPI_DOUBLE;
  atomMbrTypes[2] = MPI_INT;
  
  for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
  
  MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
  MPI_Type_commit(&mpiAtomStructType);


  // **********************************************************************
  // Init the bondStruct mpi type
  
  bondStruct bondProto; // mpiPrototype
  int bondBC[3] = {60,1,1};  // block counts
  MPI_Aint bondDspls[3];           // displacements
  MPI_Datatype bondMbrTypes[3];    // member mpi types
  
  MPI_Address(&bondProto.nameA, &bondDspls[0]);
  MPI_Address(&bondProto.d0,    &bondDspls[1]);
  MPI_Address(&bondProto.last,  &bondDspls[2]);
  
  bondMbrTypes[0] = MPI_CHAR;
  bondMbrTypes[1] = MPI_DOUBLE;
  bondMbrTypes[2] = MPI_INT;
  
  for (i=2; i >= 0; i--) bondDspls[i] -= bondDspls[0];
  
  MPI_Type_struct(3, bondBC, bondDspls, bondMbrTypes, &mpiBondStructType);
  MPI_Type_commit(&mpiBondStructType);


  // **********************************************************************
  // Init the bendStruct mpi type
  
  bendStruct bendProto; // mpiPrototype
  int bendBC[3] = {75,4,1};  // block counts
  MPI_Aint bendDspls[3];           // displacements
  MPI_Datatype bendMbrTypes[3];    // member mpi types
  
  MPI_Address(&bendProto.nameA, &bendDspls[0]);
  MPI_Address(&bendProto.k1,    &bendDspls[1]);
  MPI_Address(&bendProto.last,  &bendDspls[2]);
  
  bendMbrTypes[0] = MPI_CHAR;
  bendMbrTypes[1] = MPI_DOUBLE;
  bendMbrTypes[2] = MPI_INT;
  
  for (i=2; i >= 0; i--) bendDspls[i] -= bendDspls[0];
  
  MPI_Type_struct(3, bendBC, bendDspls, bendMbrTypes, &mpiBendStructType);
  MPI_Type_commit(&mpiBendStructType);


  // **********************************************************************
  // Init the torsionStruct mpi type
  
  torsionStruct torsionProto; // mpiPrototype
  int torsionBC[3] = {90,4,1};  // block counts
  MPI_Aint torsionDspls[3];           // displacements
  MPI_Datatype torsionMbrTypes[3];    // member mpi types
  
  MPI_Address(&torsionProto.nameA, &torsionDspls[0]);
  MPI_Address(&torsionProto.k1,    &torsionDspls[1]);
  MPI_Address(&torsionProto.last,  &torsionDspls[2]);
  
  torsionMbrTypes[0] = MPI_CHAR;
  torsionMbrTypes[1] = MPI_DOUBLE;
  torsionMbrTypes[2] = MPI_INT;
  
  for (i=2; i >= 0; i--) torsionDspls[i] -= torsionDspls[0];
  
  MPI_Type_struct(3, torsionBC, torsionDspls, torsionMbrTypes, 
                  &mpiTorsionStructType);
  MPI_Type_commit(&mpiTorsionStructType);

  // ***********************************************************************
  
  if( myNode == 0 ){
#endif
    
    // generate the force file name
    
    strcpy( fileName, "TraPPE_Ex.frc" ); 
    //    fprintf( stderr,"Trying to open %s\n", fileName );
    
    // attempt to open the file in the current directory first.
    
    frcFile = fopen( fileName, "r" );
    
    if( frcFile == NULL ){
      
      // next see if the force path enviorment variable is set
      
      ffPath = getenv( ffPath_env );
      if( ffPath == NULL ) {
        sprintf( errMsg,
                 "Error opening the force field parameter file: %s\n"
                 "Have you tried setting the FORCE_PARAM_PATH environment "
                 "vairable?\n",
                 fileName );
#ifdef IS_MPI   
        cleanOpen = 0;
        mpiError = MPI_Bcast(cleanOpen,1,MPI_INT,0,MPI_COMM_WORLD);
#endif
        ffError( errMsg );
      }
      
      
      strcpy( temp, ffPath );
      strcat( temp, "/" );
      strcat( temp, fileName );
      strcpy( fileName, temp );
      
      frcFile = fopen( fileName, "r" );
      
      if( frcFile == NULL ){
        
        sprintf( errMsg,
                 "Error opening the force field parameter file: %s\n"
                 "Have you tried setting the FORCE_PARAM_PATH environment "
                 "vairable?\n",
                 fileName );
#ifdef IS_MPI   
        cleanOpen = 0;
        mpiError = MPI_Bcast(cleanOpen,1,MPI_INT,0,MPI_COMM_WORLD);
#endif
        ffError( errMsg );
      }
    }
    cleanOpen = 1;
#ifdef IS_MPI
  }

  mpiError = MPI_Bcast(cleanOpen,1,MPI_INT,0,MPI_COMM_WORLD);
  if( !cleanOpen ) fferror( NULL );
#endif
}


TraPPE_ExFF::~TraPPE_ExFF(){
  
  fclose( frcFile );
}

void TraPPE_ExFF::initializeAtoms( void ){
  
  class LinkedType {
  public:
    LinkedType(){ 
      next = NULL;
      name[0] = '\0';
    }
    ~LinkedType(){ if( next != NULL ) delete next; }

    LinkedType* find(char* key){ 
      if( !strcmp(name, key) ) return this;
      if( next != NULL ) return next->find(key);
      return NULL;
    }
    
#ifdef IS_MPI
    void add( atomStruct &info ){
      if( next != NULL ) next->add(info);
      else{
        next = new LinkedType();
        strcpy(next->name, info.name);
        next->isDipole = info.dipole;
        next->mass     = info.mass;
        next->epslon   = info.epslon;
        next->sigma    = info.sigma;
        next->dipole   = info.dipole;
      }
    }
    
    void duplicate( atomStruct &info ){
      strcpy(info.name, name);
      info.isDipole = dipole;
      info.mass     = mass;
      info.epslon   = epslon;
      info.sigma    = sigma;
      info.dipole   = dipole;
      info.last     = 0;
    }


#endif

    char name[15];
    int isDipole;
    double mass;
    double epslon;
    double sigma;
    double dipole;
    LinkedType* next;
  };
  
  LinkedType* headAtomType; 
  LinkedType* currentAtomType;
  LinkedType* tempAtomType;

#ifdef IS_MPI
  atomStruct info;
  info.last = -1; // initialize last to have the error set. 
                  // if things go well, last will be set to 0
#endif
  

  char readLine[500];
  char* the_token;
  char* eof_test;
  int foundAtom = 0;
  int lineNum = 0;
  int i;
  char errMsg[1000];

  
  //////////////////////////////////////////////////
  // a quick water fix 

  double waterI[3][3]; 
  waterI[0][0] = 1.76958347772500;
  waterI[0][1] = 0.0;
  waterI[0][2] = 0.0;

  waterI[1][0] = 0.0;
  waterI[1][1] = 0.614537057924513;
  waterI[1][2] = 0.0;

  waterI[2][0] = 0.0;
  waterI[2][1] = 0.0;
  waterI[2][2] = 1.15504641980049;


  double headI[3][3]; 
  headI[0][0] = 1125;
  headI[0][1] = 0.0;
  headI[0][2] = 0.0;

  headI[1][0] = 0.0;
  headI[1][1] = 1125;
  headI[1][2] = 0.0;

  headI[2][0] = 0.0;
  headI[2][1] = 0.0;
  headI[2][2] = 250;

  
  //////////////////////////////////////////////////

  Atom** the_atoms;
  int nAtoms;
  the_atoms = entry_plug->atoms;
  nAtoms = entry_plug->n_atoms;
  
  
#ifdef IS_MPI
  if( myNode == 0 ){
#endif
    
    // read in the atom types.
    
    rewind( frcFile );
    headAtomType = new LinkedType;
    currentAtomType = headAtomType;
    
    eof_test = fgets( readLine, sizeof(readLine), frcFile );
    lineNum++;
    if( eof_test == NULL ){
      sprintf( errMsg, "Error in reading Atoms from force file.\n" );
#ifdef IS_MPI
      sendFrcStruct( &info, mpiAtomStructType );
#endif
      ffError(errMsg);
    }
    
    
    while( !foundAtom ){
      while( eof_test != NULL && readLine[0] != '#' ){
        eof_test = fgets( readLine, sizeof(readLine), frcFile );
        lineNum++;
      }
      if( eof_test == NULL ){
        sprintf( errMsg, 
                 "Error in reading Atoms from force file at line %d.\n",
                 lineNum );
#ifdef IS_MPI
        sendFrcStruct( &info, mpiAtomStructType );
#endif
        ffError(errMsg);
      }
      
      the_token = strtok( readLine, " ,;\t#\n" );
      foundAtom = !strcmp( "AtomTypes", the_token );
      
      if( !foundAtom ){
        eof_test = fgets( readLine, sizeof(readLine), frcFile );
        lineNum++;
        
        if( eof_test == NULL ){
          sprintf( errMsg, 
                   "Error in reading Atoms from force file at line %d.\n",
                   lineNum );
#ifdef IS_MPI
          sendFrcStruct( &info, mpiAtomStructType );
#endif
          ffError(errMsg);
        } 
      }
    }
    
    // we are now at the AtomTypes section.
    
    eof_test = fgets( readLine, sizeof(readLine), frcFile );
    lineNum++;
    
    if( eof_test == NULL ){
      sprintf( errMsg, 
               "Error in reading Atoms from force file at line %d.\n",
               lineNum );
#ifdef IS_MPI
      sendFrcStruct( &info, mpiAtomStructType );
#endif      
      ffError(errMsg);
    }
    
    while( readLine[0] != '#' && eof_test != NULL ){
      
      if( readLine[0] != '!' ){
        
        the_token = strtok( readLine, " \n\t,;" );
        if( the_token != NULL ){
          
          strcpy( currentAtomType->name, the_token );
          
          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
            sprintf( errMsg, "Error parseing AtomTypes: line %d\n", lineNum );
#ifdef IS_MPI
            sendFrcStruct( &info, mpiAtomStructType );
#endif
            ffError(errMsg);
          }
          
          sscanf( the_token, "%d", &currentAtomType->isDipole );
          
          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
            sprintf( errMsg, "Error parseing AtomTypes: line %d\n", lineNum );
#ifdef IS_MPI
            sendFrcStruct( &info, mpiAtomStructType );
#endif
            ffError(errMsg);
          }
          
          sscanf( the_token, "%lf", &currentAtomType->mass );
          
          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
            sprintf( errMsg, "Error parseing AtomTypes: line %d\n", lineNum );
#ifdef IS_MPI
            sendFrcStruct( &info, mpiAtomStructType );
#endif
            ffError(errMsg);
          }
          
          sscanf( the_token, "%lf", &currentAtomType->epslon );
          
          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
            sprintf( errMsg, "Error parseing AtomTypes: line %d\n", lineNum );
#ifdef IS_MPI
            sendFrcStruct( &info, mpiAtomStructType );
#endif
            ffError(errMsg);
          }
          
          sscanf( the_token, "%lf", &currentAtomType->sigma );
          
          if( currentAtomType->isDipole ){
            if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
              sprintf( errMsg, "Error parseing AtomTypes: line %d\n", 
                       lineNum );
#ifdef IS_MPI
              sendFrcStruct( &info, mpiAtomStructType );
#endif
              ffError(errMsg);
            }
            
            sscanf( the_token, "%lf", &currentAtomType->dipole );
          }
        }
      }
      
      tempAtomType = new LinkedType;
      currentAtomType->next = tempAtomType;
      currentAtomType = tempAtomType;
      
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;
    }

#ifdef IS_MPI
    
    // send out the linked list to all the other processes
    
    currentAtomType = headAtomType;
    while( currentAtomType != NULL ){
      currentAtomType->duplicate( info );
      sendFrcStruct( &info, mpiAtomStructType );
      currentAtomType = currentAtomType->next;
    }
    info.last = 1;
    sendFrcStruct( &info, mpiAtomStructType );
    
  }

  else{
    
    // listen for node 0 to send out the force params

    headAtomType = new LinkedType;
    recieveFrcStruct( &info, mpiAtomStructType );
    while( !info.last ){

      headAtomType->add( info );
      recieveFrcStruct( &info, mpiAtomStructType );
    }

    if( info.last < 0 ) ffError( NULL ); // an error has occured, exit quietly.

    
  }
#endif

  
  // initialize the atoms
  
  DirectionalAtom* dAtom;

  for( i=0; i<nAtoms; i++ ){
    
    currentAtomType = headAtomType->find( the_atoms[i]->getType() );
    if( currentAtomType == NULL ){
      sprintf( errMsg, "AtomType error, %s not found in force file.\n",
               the_atoms[i]->getType() );
      ffError(errMsg);
    }
    
    the_atoms[i]->setMass( currentAtomType->mass );
    the_atoms[i]->setEpslon( currentAtomType->epslon );
    the_atoms[i]->setSigma( currentAtomType->sigma );
    the_atoms[i]->setLJ();

    if( currentAtomType->isDipole ){
      if( the_atoms[i]->isDirectional() ){
        
        dAtom = (DirectionalAtom *) the_atoms[i];
        dAtom->setMu( currentAtomType->dipole );
        dAtom->setHasDipole( 1 );
        dAtom->setJx( 0.0 );
        dAtom->setJy( 0.0 );
        dAtom->setJz( 0.0 );
        
        if(!strcmp("SSD",the_atoms[i]->getType())){
          dAtom->setI( waterI );
          dAtom->setSSD( 1 );
        }
        else if(!strcmp("HEAD",the_atoms[i]->getType())){
          dAtom->setI( headI );
          dAtom->setSSD( 0 );
        }
        else{
          sprintf(errMsg,
                  "AtmType error, %s does not have a moment of inertia set.\n",
                  the_atoms[i]->getType() );
          ffError(errMsg);
        }
        entry_plug->n_dipoles++;
      }
      else{
        
        sprintf( errMsg,
                "TraPPE_ExFF error: Atom \"%s\" is a dipole, yet no standard"
                 " orientation was specifed in the BASS file.\n",
                 currentAtomType->name );
        
        ffError(errMsg);
      }
    }
    else{
      if( the_atoms[i]->isDirectional() ){
        sprintf( errMsg,
                 "TraPPE_ExFF error: Atom \"%s\" was given a standard"
                 "orientation in the BASS file, yet it is not a dipole.\n",
                 currentAtomType->name);
        ffError(errMsg);
      }
    }
  }


  // clean up the memory
  
  delete headAtomType;
}

void TraPPE_ExFF::initializeBonds( bond_pair* the_bonds ){
  
  class LinkedType {
  public:
    LinkedType(){ 
      next = NULL;
      nameA[0] = '\0';
      nameB[0] = '\0';
      type[0] = '\0';
    }
    ~LinkedType(){ if( next != NULL ) delete next; }

    LinkedType* find(char* key1, char* key2){ 
      if( !strcmp(nameA, key1 ) && !strcmp( nameB, key2 ) ) return this;
      if( !strcmp(nameA, key2 ) && !strcmp( nameB, key1 ) ) return this;
      if( next != NULL ) return next->find(key1, key2);
      return NULL;
    }
    
#ifdef IS_MPI
    void add( bondStruct &info ){
      if( next != NULL ) next->add(info);
      else{
        next = new LinkedType();
        strcpy(next->nameA, info.nameA);
        strcpy(next->nameB, info.nameB);
        strcpy(next->type,  info.type);
        next->d0 = info.d0;
      }
    }
    
    void duplicate( bondStruct &info ){
      strcpy(info.nameA, nameA);
      strcpy(info.nameB, nameB);
      strcpy(info.type,  type);
      info.d0   = d0;
      info.last = 0;
    }


#endif

    char nameA[15];
    char nameB[15];
    char type[30];
    double d0;

    LinkedType* next;
  };


  LinkedType* headBondType; 
  LinkedType* currentBondType;
  LinkedType* tempBondType;

  char readLine[500];
  char* the_token;
  char* eof_test;
  int foundBond = 0;
  int lineNum = 0;
  int i, a, b;
  char* atomA;
  char* atomB;

  SRI **the_sris;
  Atom** the_atoms;
  int nBonds;
  the_sris = entry_plug->sr_interactions;
  the_atoms = entry_plug->atoms;
  nBonds = entry_plug->n_bonds;

  // read in the bond types.

  rewind( frcFile );
  headBondType = new LinkedType;
  currentBondType = headBondType;
  
  eof_test = fgets( readLine, sizeof(readLine), frcFile );
  lineNum++;
  if( eof_test == NULL ){
    fprintf( stderr, "Error in reading Bonds from force file.\n" );
    exit(8);
  }

  
  while( !foundBond ){
    while( eof_test != NULL && readLine[0] != '#' ){
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;
    }
    if( eof_test == NULL ){
      fprintf( stderr, 
               "Error in reading Bonds from force file at line %d.\n",
               lineNum );
      exit(8);
    }

    the_token = strtok( readLine, " ,;\t#\n" );
    foundBond = !strcmp( "BondTypes", the_token );
    
    if( !foundBond ){
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;

      if( eof_test == NULL ){
        fprintf( stderr, 
                 "Error in reading Bonds from force file at line %d.\n",
                 lineNum );
        exit(8);
      } 
    }
  }

  // we are now at the BondTypes section.

  eof_test = fgets( readLine, sizeof(readLine), frcFile );
  lineNum++;
  
  if( eof_test == NULL ){
    fprintf( stderr, 
             "Error in reading Bonds from force file at line %d.\n",
             lineNum );
    exit(8);
  }

  while( readLine[0] != '#' && eof_test != NULL ){
    
    if( readLine[0] != '!' ){
      
      the_token = strtok( readLine, " \n\t,;" );
      if( the_token != NULL ){
        
        strcpy( currentBondType->nameA, the_token );
        
        if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
          fprintf( stderr, "Error parseing BondTypes: line %d\n", lineNum );
          exit(8);
        }

        strcpy( currentBondType->nameB, the_token );

        if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
          fprintf( stderr, "Error parseing BondTypes: line %d\n", lineNum );
          exit(8);
        }

        strcpy( currentBondType->type, the_token );

        if( !strcmp( currentBondType->type, "fixed" ) ){
          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
            fprintf( stderr, "Error parseing BondTypes: line %d\n", lineNum );
            exit(8);
          }
          
          sscanf( the_token, "%lf", &currentBondType->d0 );
        }
        
        else{
          fprintf(stderr, 
                  "Unknown TraPPE_Ex bond type \"%s\" at line %d\n",
                  currentBondType->type,
                  lineNum );
          exit(8);
        }                 
      }
    }
    
    tempBondType = new LinkedType;
    currentBondType->next = tempBondType;
    currentBondType = tempBondType;
    
    eof_test = fgets( readLine, sizeof(readLine), frcFile );
    lineNum++;
  }

  
  // initialize the Bonds
  

  for( i=0; i<nBonds; i++ ){
    
    a = the_bonds[i].a;
    b = the_bonds[i].b;

    atomA = the_atoms[a]->getType();
    atomB = the_atoms[b]->getType();
    currentBondType = headBondType->find( atomA, atomB );
    if( currentBondType == NULL ){
      fprintf( stderr, "BondType error, %s - %s not found in force file.\n",
               atomA, atomB );
      exit(8);
    }
    
    if( !strcmp( currentBondType->type, "fixed" ) ){
      
      the_sris[i] = new ConstrainedBond( *the_atoms[a], 
                                         *the_atoms[b],
                                         currentBondType->d0 );
      entry_plug->n_constraints++;
    }
  }


  // clean up the memory
  
  delete headBondType;

}

void TraPPE_ExFF::initializeBends( bend_set* the_bends ){
  
  class LinkedType {
  public:
    LinkedType(){ 
      next = NULL;
      nameA[0] = '\0';
      nameB[0] = '\0';
      nameC[0] = '\0';
      type[0] = '\0';
    }
    ~LinkedType(){ if( next != NULL ) delete next; }

    LinkedType* find( char* key1, char* key2, char* key3 ){ 
      if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 ) 
          && !strcmp( nameC, key3 ) ) return this;
      if( !strcmp( nameA, key3 ) && !strcmp( nameB, key2 ) 
          && !strcmp( nameC, key1 ) ) return this;
      if( next != NULL ) return next->find(key1, key2, key3);
      return NULL;
    }
    
#ifdef IS_MPI

    void add( bendStruct &info ){
      if( next != NULL ) next->add(info);
      else{
        next = new LinkedType();
        strcpy(next->nameA, info.nameA);
        strcpy(next->nameB, info.nameB);
        strcpy(next->nameC, info.nameC);
        strcpy(next->type,  info.type);
        next->k1 = info.k1;
        next->k2 = info.k2;
        next->k3 = info.k3;
        next->t0 = info.t0;
      }
    }
    
    void duplicate( bendStruct &info ){
      strcpy(info.nameA, nameA);
      strcpy(info.nameB, nameB);
      strcpy(info.nameC, nameC);
      strcpy(info.type,  type);
      info.k1   = k1;
      info.k2   = k2;
      info.k3   = k3;
      info.t0   = t0;
      info.last = 0;
    }

#endif

    char nameA[15];
    char nameB[15];
    char nameC[15];
    char type[30];
    double k1, k2, k3, t0;

    LinkedType* next;
  };
  
  LinkedType* headBendType; 
  LinkedType* currentBendType;
  LinkedType* tempBendType;

  char readLine[500];
  char* the_token;
  char* eof_test;
  int foundBend = 0;
  int lineNum = 0;
  int i, a, b, c, index;
  char* atomA;
  char* atomB;
  char* atomC;
  QuadraticBend* qBend;

  SRI **the_sris;
  Atom** the_atoms;
  int nBends;
  the_sris = entry_plug->sr_interactions;
  the_atoms = entry_plug->atoms;
  nBends = entry_plug->n_bends;

  // read in the bend types.

  rewind( frcFile );
  headBendType = new LinkedType;
  currentBendType = headBendType;
  
  eof_test = fgets( readLine, sizeof(readLine), frcFile );
  lineNum++;
  if( eof_test == NULL ){
    fprintf( stderr, "Error in reading Bends from force file.\n" );
    exit(8);
  }

  
  while( !foundBend ){
    while( eof_test != NULL && readLine[0] != '#' ){
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;
    }
    if( eof_test == NULL ){
      fprintf( stderr, 
               "Error in reading Bends from force file at line %d.\n",
               lineNum );
      exit(8);
    }

    the_token = strtok( readLine, " ,;\t#\n" );
    foundBend = !strcmp( "BendTypes", the_token );
    
    if( !foundBend ){
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;

      if( eof_test == NULL ){
        fprintf( stderr, 
                 "Error in reading Bends from force file at line %d.\n",
                 lineNum );
        exit(8);
      } 
    }
  }

  // we are now at the BendTypes section.

  eof_test = fgets( readLine, sizeof(readLine), frcFile );
  lineNum++;
  
  if( eof_test == NULL ){
    fprintf( stderr, 
             "Error in reading Bends from force file at line %d.\n",
             lineNum );
    exit(8);
  }

  while( readLine[0] != '#' && eof_test != NULL ){
    
    if( readLine[0] != '!' ){
      
      the_token = strtok( readLine, " \n\t,;" );
      if( the_token != NULL ){
        
        strcpy( currentBendType->nameA, the_token );
        
        if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
          fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
          exit(8);
        }

        strcpy( currentBendType->nameB, the_token );

        if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
          fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
          exit(8);
        }

        strcpy( currentBendType->nameC, the_token );

        if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
          fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
          exit(8);
        }

        strcpy( currentBendType->type, the_token );

        if( !strcmp( currentBendType->type, "quadratic" ) ){
          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
            fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
            exit(8);
          }
          
          sscanf( the_token, "%lf", &currentBendType->k1 );

          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
            fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
            exit(8);
          }
          
          sscanf( the_token, "%lf", &currentBendType->k2 );

          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
            fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
            exit(8);
          }
          
          sscanf( the_token, "%lf", &currentBendType->k3 );

          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
            fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
            exit(8);
          }
          
          sscanf( the_token, "%lf", &currentBendType->t0 );
        }
        
        else{
          fprintf(stderr, 
                  "Unknown TraPPE_Ex bend type \"%s\" at line %d\n",
                  currentBendType->type,
                  lineNum );
          exit(8);
        }                 
      }
    }
    
    tempBendType = new LinkedType;
    currentBendType->next = tempBendType;
    currentBendType = tempBendType;
    
    eof_test = fgets( readLine, sizeof(readLine), frcFile );
    lineNum++;
  }

  
  // initialize the Bends

  for( i=0; i<nBends; i++ ){
    
    a = the_bends[i].a;
    b = the_bends[i].b;
    c = the_bends[i].c;

    atomA = the_atoms[a]->getType();
    atomB = the_atoms[b]->getType();
    atomC = the_atoms[c]->getType();
    currentBendType = headBendType->find( atomA, atomB, atomC );
    if( currentBendType == NULL ){
      fprintf( stderr, "BendType error, %s - %s - %s not found"
               " in force file.\n",
               atomA, atomB, atomC );
      exit(8);
    }
    
    if( !strcmp( currentBendType->type, "quadratic" ) ){
      
      index = i + entry_plug->n_bonds;
      qBend = new QuadraticBend( *the_atoms[a], 
                                 *the_atoms[b],
                                 *the_atoms[c] );
      qBend->setConstants( currentBendType->k1,
                           currentBendType->k2,
                           currentBendType->k3,
                           currentBendType->t0 );
      the_sris[index] = qBend;
    }
  }


  // clean up the memory
  
  delete headBendType;

}

void TraPPE_ExFF::initializeTorsions( torsion_set* the_torsions ){

  class LinkedType {
  public:
    LinkedType(){ 
      next = NULL;
      nameA[0] = '\0';
      nameB[0] = '\0';
      nameC[0] = '\0';
      type[0] = '\0';
    }
    ~LinkedType(){ if( next != NULL ) delete next; }

    LinkedType* find( char* key1, char* key2, char* key3, char* key4 ){ 
      if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 ) &&
          !strcmp( nameC, key3 ) && !strcmp( nameD, key4 ) ) return this;

      if( !strcmp( nameA, key4 ) && !strcmp( nameB, key3 ) &&
          !strcmp( nameC, key2 ) && !strcmp( nameD, key4 ) ) return this;

      if( next != NULL ) return next->find(key1, key2, key3, key4);
      return NULL;
    }
    
#ifdef IS_MPI

    void add( torsionStruct &info ){
      if( next != NULL ) next->add(info);
      else{
        next = new LinkedType();
        strcpy(next->nameA, info.nameA);
        strcpy(next->nameB, info.nameB);
        strcpy(next->nameC, info.nameC);
        strcpy(next->type,  info.type);
        next->k1 = info.k1;
        next->k2 = info.k2;
        next->k3 = info.k3;
        next->k4 = info.k4;
      }
    }
    
    void duplicate( torsionStruct &info ){
      strcpy(info.nameA, nameA);
      strcpy(info.nameB, nameB);
      strcpy(info.nameC, nameC);
      strcpy(info.nameD, nameD);
      strcpy(info.type,  type);
      info.k1   = k1;
      info.k2   = k2;
      info.k3   = k3;
      info.k4   = k4;
      info.last = 0;
    }

#endif

    char nameA[15];
    char nameB[15];
    char nameC[15];
    char nameD[15];
    char type[30];
    double k1, k2, k3, k4;

    LinkedType* next;
  };
  
  LinkedType* headTorsionType; 
  LinkedType* currentTorsionType;
  LinkedType* tempTorsionType;

  char readLine[500];
  char* the_token;
  char* eof_test;
  int foundTorsion = 0;
  int lineNum = 0;
  int i, a, b, c, d, index;
  char* atomA;
  char* atomB;
  char* atomC;
  char* atomD;
  CubicTorsion* cTors;

  SRI **the_sris;
  Atom** the_atoms;
  int nTorsions;
  the_sris = entry_plug->sr_interactions;
  the_atoms = entry_plug->atoms;
  nTorsions = entry_plug->n_torsions;

  // read in the torsion types.

  rewind( frcFile );
  headTorsionType = new LinkedType;
  currentTorsionType = headTorsionType;
  
  eof_test = fgets( readLine, sizeof(readLine), frcFile );
  lineNum++;
  if( eof_test == NULL ){
    fprintf( stderr, "Error in reading Torsions from force file.\n" );
    exit(8);
  }

  
  while( !foundTorsion ){
    while( eof_test != NULL && readLine[0] != '#' ){
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;
    }
    if( eof_test == NULL ){
      fprintf( stderr, 
               "Error in reading Torsions from force file at line %d.\n",
               lineNum );
      exit(8);
    }

    the_token = strtok( readLine, " ,;\t#\n" );
    foundTorsion = !strcmp( "TorsionTypes", the_token );
    
    if( !foundTorsion ){
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;

      if( eof_test == NULL ){
        fprintf( stderr, 
                 "Error in reading Torsions from force file at line %d.\n",
                 lineNum );
        exit(8);
      } 
    }
  }

  // we are now at the TorsionTypes section.

  eof_test = fgets( readLine, sizeof(readLine), frcFile );
  lineNum++;
  
  if( eof_test == NULL ){
    fprintf( stderr, 
             "Error in reading Torsions from force file at line %d.\n",
             lineNum );
    exit(8);
  }

  while( readLine[0] != '#' && eof_test != NULL ){
    
    if( readLine[0] != '!' ){
      
      the_token = strtok( readLine, " \n\t,;" );
      if( the_token != NULL ){
        
        strcpy( currentTorsionType->nameA, the_token );
        
        if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
          fprintf( stderr, "Error parseing TorsionTypes: line %d\n", lineNum );
          exit(8);
        }

        strcpy( currentTorsionType->nameB, the_token );

        if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
          fprintf( stderr, "Error parseing TorsionTypes: line %d\n", lineNum );
          exit(8);
        }

        strcpy( currentTorsionType->nameC, the_token );

        if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
          fprintf( stderr, "Error parseing TorsionTypes: line %d\n", lineNum );
          exit(8);
        }

        strcpy( currentTorsionType->nameD, the_token );

        if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
          fprintf( stderr, "Error parseing TorsionTypes: line %d\n", lineNum );
          exit(8);
        }

        strcpy( currentTorsionType->type, the_token );

        if( !strcmp( currentTorsionType->type, "cubic" ) ){
          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
            fprintf(stderr,"Error parseing TorsionTypes: line %d\n", lineNum );
            exit(8);
          }
          
          sscanf( the_token, "%lf", &currentTorsionType->k1 );

          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
            fprintf(stderr,"Error parseing TorsionTypes: line %d\n", lineNum );
            exit(8);
          }
          
          sscanf( the_token, "%lf", &currentTorsionType->k2 );

          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
            fprintf(stderr,"Error parseing TorsionTypes: line %d\n", lineNum );
            exit(8);
          }
          
          sscanf( the_token, "%lf", &currentTorsionType->k3 );

          if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
            fprintf(stderr,"Error parseing TorsionTypes: line %d\n", lineNum );
            exit(8);
          }
          
          sscanf( the_token, "%lf", &currentTorsionType->k4 );
        }
        
        else{
          fprintf(stderr, 
                  "Unknown TraPPE_Ex torsion type \"%s\" at line %d\n",
                  currentTorsionType->type,
                  lineNum );
          exit(8);
        }                 
      }
    }
    
    tempTorsionType = new LinkedType;
    currentTorsionType->next = tempTorsionType;
    currentTorsionType = tempTorsionType;
    
    eof_test = fgets( readLine, sizeof(readLine), frcFile );
    lineNum++;
  }

  
  // initialize the Torsions

  for( i=0; i<nTorsions; i++ ){
    
    a = the_torsions[i].a;
    b = the_torsions[i].b;
    c = the_torsions[i].c;
    d = the_torsions[i].d;

    atomA = the_atoms[a]->getType();
    atomB = the_atoms[b]->getType();
    atomC = the_atoms[c]->getType();
    atomD = the_atoms[d]->getType();
    currentTorsionType = headTorsionType->find( atomA, atomB, atomC, atomD );
    if( currentTorsionType == NULL ){
      fprintf( stderr, "TorsionType error, %s - %s - %s - %s not found"
               " in force file.\n",
               atomA, atomB, atomC, atomD );
      exit(8);
    }
    
    if( !strcmp( currentTorsionType->type, "cubic" ) ){
      index = i + entry_plug->n_bonds + entry_plug->n_bends;
      
      cTors = new CubicTorsion( *the_atoms[a], *the_atoms[b],
                                *the_atoms[c], *the_atoms[d] );
      cTors->setConstants( currentTorsionType->k1, currentTorsionType->k2,
                           currentTorsionType->k3, currentTorsionType->k4 );
      the_sris[index] = cTors;
    }
  }


  // clean up the memory
  
  delete headTorsionType;

}
Revision:	152
Committed:	Mon Oct 21 22:03:00 2002 UTC (21 years, 8 months ago) by mmeineke
File size:	32105 byte(s)
Log Message:	* empty log message *