mdtools/interface_implementation/TraPPEFF.cpp

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

#include <iostream>
using namespace std;

#include "ForceFields.hpp"
#include "SRI.hpp"
#include "simError.h"

TraPPEFF::TraPPEFF(){

#ifdef IS_MPI
  sprintf( painCave.errMsg,
           "TraPPE ForceField does not currently support mpi" );
  painCave.isFatal = 1;
  simError();
#endif // is_mpi


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

  // generate the force file name

  strcpy( fileName, "TraPPE.frc" ); 

  // 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 );
    strcpy( temp, ffPath );
    strcat( temp, "/" );
    strcat( temp, fileName );
    strcpy( fileName, temp );

    frcFile = fopen( fileName, "r" );
    
    if( frcFile == NULL ){
      
      fprintf( stderr,
               "Error opening the force field parameter file: %s\n"
               "Have you tried setting the FORCE_PARAM_PATH environment "
               "vairable?\n",
               fileName );
      exit( 8 );
    }
  }
}

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

void TraPPEFF::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;
    }
    
    char name[15];
    double mass;
    double epslon;
    double sigma;
    LinkedType* next;
  };
  
  LinkedType* headAtomType; 
  LinkedType* currentAtomType;
  LinkedType* tempAtomType;

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

  Atom** the_atoms;
  int nAtoms;
  the_atoms = entry_plug->atoms;
  nAtoms = entry_plug->n_atoms;
  
  // read in the atom types.

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

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

    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 ){
        fprintf( stderr, 
                 "Error in reading Atoms from force file at line %d.\n",
                 lineNum );
        exit(8);
      } 
    }
  }

  // we are now at the AtomTypes section.

  eof_test = fgets( readLine, sizeof(readLine), frcFile );
  lineNum++;
  
  if( eof_test == NULL ){
    fprintf( stderr, 
             "Error in reading Atoms 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( currentAtomType->name, the_token );
        
        if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
          fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
          exit(8);
        }

        sscanf( the_token, "%lf", &currentAtomType->mass );

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

        sscanf( the_token, "%lf", &currentAtomType->sigma );
      }
    }
    
    tempAtomType = new LinkedType;
    currentAtomType->next = tempAtomType;
    currentAtomType = tempAtomType;
    
    eof_test = fgets( readLine, sizeof(readLine), frcFile );
    lineNum++;
  }

  
  // initialize the atoms

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


  // clean up the memory
  
  delete headAtomType;
}

void TraPPEFF::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;
    }
    
    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 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 TraPPEFF::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;
    }
    
    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 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 TraPPEFF::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;
    }
    
    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 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:	176
Committed:	Thu Nov 14 22:00:44 2002 UTC (21 years, 7 months ago) by mmeineke
File size:	20434 byte(s)
Log Message:	* empty log message *
#	Content
1	#include <cstdlib>
2	#include <cstdio>
3	#include <cstring>
4
5	#include <iostream>
6	using namespace std;
7
8	#include "ForceFields.hpp"
9	#include "SRI.hpp"
10	#include "simError.h"
11
12	TraPPEFF::TraPPEFF(){
13
14	#ifdef IS_MPI
15	sprintf( painCave.errMsg,
16	"TraPPE ForceField does not currently support mpi" );
17	painCave.isFatal = 1;
18	simError();
19	#endif // is_mpi
20
21
22	char fileName[200];
23	char* ffPath_env = "FORCE_PARAM_PATH";
24	char* ffPath;
25	char temp[200];
26
27	// generate the force file name
28
29	strcpy( fileName, "TraPPE.frc" );
30
31	// attempt to open the file in the current directory first.
32
33	frcFile = fopen( fileName, "r" );
34
35	if( frcFile == NULL ){
36
37	// next see if the force path enviorment variable is set
38
39	ffPath = getenv( ffPath_env );
40	strcpy( temp, ffPath );
41	strcat( temp, "/" );
42	strcat( temp, fileName );
43	strcpy( fileName, temp );
44
45	frcFile = fopen( fileName, "r" );
46
47	if( frcFile == NULL ){
48
49	fprintf( stderr,
50	"Error opening the force field parameter file: %s\n"
51	"Have you tried setting the FORCE_PARAM_PATH environment "
52	"vairable?\n",
53	fileName );
54	exit( 8 );
55	}
56	}
57	}
58
59	TraPPEFF::~TraPPEFF(){
60
61	fclose( frcFile );
62	}
63
64	void TraPPEFF::initializeAtoms( void ){
65
66	class LinkedType {
67	public:
68	LinkedType(){
69	next = NULL;
70	name[0] = '\0';
71	}
72	~LinkedType(){ if( next != NULL ) delete next; }
73
74	LinkedType* find(char* key){
75	if( !strcmp(name, key) ) return this;
76	if( next != NULL ) return next->find(key);
77	return NULL;
78	}
79
80	char name[15];
81	double mass;
82	double epslon;
83	double sigma;
84	LinkedType* next;
85	};
86
87	LinkedType* headAtomType;
88	LinkedType* currentAtomType;
89	LinkedType* tempAtomType;
90
91	char readLine[500];
92	char* the_token;
93	char* eof_test;
94	int foundAtom = 0;
95	int lineNum = 0;
96	int i;
97
98	Atom** the_atoms;
99	int nAtoms;
100	the_atoms = entry_plug->atoms;
101	nAtoms = entry_plug->n_atoms;
102
103	// read in the atom types.
104
105	rewind( frcFile );
106	headAtomType = new LinkedType;
107	currentAtomType = headAtomType;
108
109	eof_test = fgets( readLine, sizeof(readLine), frcFile );
110	lineNum++;
111	if( eof_test == NULL ){
112	fprintf( stderr, "Error in reading Atoms from force file.\n" );
113	exit(8);
114	}
115
116
117	while( !foundAtom ){
118	while( eof_test != NULL && readLine[0] != '#' ){
119	eof_test = fgets( readLine, sizeof(readLine), frcFile );
120	lineNum++;
121	}
122	if( eof_test == NULL ){
123	fprintf( stderr,
124	"Error in reading Atoms from force file at line %d.\n",
125	lineNum );
126	exit(8);
127	}
128
129	the_token = strtok( readLine, " ,;\t#\n" );
130	foundAtom = !strcmp( "AtomTypes", the_token );
131
132	if( !foundAtom ){
133	eof_test = fgets( readLine, sizeof(readLine), frcFile );
134	lineNum++;
135
136	if( eof_test == NULL ){
137	fprintf( stderr,
138	"Error in reading Atoms from force file at line %d.\n",
139	lineNum );
140	exit(8);
141	}
142	}
143	}
144
145	// we are now at the AtomTypes section.
146
147	eof_test = fgets( readLine, sizeof(readLine), frcFile );
148	lineNum++;
149
150	if( eof_test == NULL ){
151	fprintf( stderr,
152	"Error in reading Atoms from force file at line %d.\n",
153	lineNum );
154	exit(8);
155	}
156
157	while( readLine[0] != '#' && eof_test != NULL ){
158
159	if( readLine[0] != '!' ){
160
161	the_token = strtok( readLine, " \n\t,;" );
162	if( the_token != NULL ){
163
164	strcpy( currentAtomType->name, the_token );
165
166	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
167	fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
168	exit(8);
169	}
170
171	sscanf( the_token, "%lf", &currentAtomType->mass );
172
173	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
174	fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
175	exit(8);
176	}
177
178	sscanf( the_token, "%lf", &currentAtomType->epslon );
179
180	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
181	fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
182	exit(8);
183	}
184
185	sscanf( the_token, "%lf", &currentAtomType->sigma );
186	}
187	}
188
189	tempAtomType = new LinkedType;
190	currentAtomType->next = tempAtomType;
191	currentAtomType = tempAtomType;
192
193	eof_test = fgets( readLine, sizeof(readLine), frcFile );
194	lineNum++;
195	}
196
197
198	// initialize the atoms
199
200	for( i=0; i<nAtoms; i++ ){
201
202	currentAtomType = headAtomType->find( the_atoms[i]->getType() );
203	if( currentAtomType == NULL ){
204	fprintf( stderr, "AtomType error, %s not found in force file.\n",
205	the_atoms[i]->getType() );
206	exit(8);
207	}
208
209	the_atoms[i]->setMass( currentAtomType->mass );
210	the_atoms[i]->setEpslon( currentAtomType->epslon );
211	the_atoms[i]->setSigma( currentAtomType->sigma );
212	the_atoms[i]->setLJ();
213	}
214
215
216	// clean up the memory
217
218	delete headAtomType;
219	}
220
221	void TraPPEFF::initializeBonds( bond_pair* the_bonds ){
222
223	class LinkedType {
224	public:
225	LinkedType(){
226	next = NULL;
227	nameA[0] = '\0';
228	nameB[0] = '\0';
229	type[0] = '\0';
230	}
231	~LinkedType(){ if( next != NULL ) delete next; }
232
233	LinkedType* find(char* key1, char* key2){
234	if( !strcmp(nameA, key1 ) && !strcmp( nameB, key2 ) ) return this;
235	if( !strcmp(nameA, key2 ) && !strcmp( nameB, key1 ) ) return this;
236	if( next != NULL ) return next->find(key1, key2);
237	return NULL;
238	}
239
240	char nameA[15];
241	char nameB[15];
242	char type[30];
243	double d0;
244
245	LinkedType* next;
246	};
247
248
249
250	LinkedType* headBondType;
251	LinkedType* currentBondType;
252	LinkedType* tempBondType;
253
254	char readLine[500];
255	char* the_token;
256	char* eof_test;
257	int foundBond = 0;
258	int lineNum = 0;
259	int i, a, b;
260	char* atomA;
261	char* atomB;
262
263	SRI **the_sris;
264	Atom** the_atoms;
265	int nBonds;
266	the_sris = entry_plug->sr_interactions;
267	the_atoms = entry_plug->atoms;
268	nBonds = entry_plug->n_bonds;
269
270	// read in the bond types.
271
272	rewind( frcFile );
273	headBondType = new LinkedType;
274	currentBondType = headBondType;
275
276	eof_test = fgets( readLine, sizeof(readLine), frcFile );
277	lineNum++;
278	if( eof_test == NULL ){
279	fprintf( stderr, "Error in reading Bonds from force file.\n" );
280	exit(8);
281	}
282
283
284	while( !foundBond ){
285	while( eof_test != NULL && readLine[0] != '#' ){
286	eof_test = fgets( readLine, sizeof(readLine), frcFile );
287	lineNum++;
288	}
289	if( eof_test == NULL ){
290	fprintf( stderr,
291	"Error in reading Bonds from force file at line %d.\n",
292	lineNum );
293	exit(8);
294	}
295
296	the_token = strtok( readLine, " ,;\t#\n" );
297	foundBond = !strcmp( "BondTypes", the_token );
298
299	if( !foundBond ){
300	eof_test = fgets( readLine, sizeof(readLine), frcFile );
301	lineNum++;
302
303	if( eof_test == NULL ){
304	fprintf( stderr,
305	"Error in reading Bonds from force file at line %d.\n",
306	lineNum );
307	exit(8);
308	}
309	}
310	}
311
312	// we are now at the BondTypes section.
313
314	eof_test = fgets( readLine, sizeof(readLine), frcFile );
315	lineNum++;
316
317	if( eof_test == NULL ){
318	fprintf( stderr,
319	"Error in reading Bonds from force file at line %d.\n",
320	lineNum );
321	exit(8);
322	}
323
324	while( readLine[0] != '#' && eof_test != NULL ){
325
326	if( readLine[0] != '!' ){
327
328	the_token = strtok( readLine, " \n\t,;" );
329	if( the_token != NULL ){
330
331	strcpy( currentBondType->nameA, the_token );
332
333	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
334	fprintf( stderr, "Error parseing BondTypes: line %d\n", lineNum );
335	exit(8);
336	}
337
338	strcpy( currentBondType->nameB, the_token );
339
340	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
341	fprintf( stderr, "Error parseing BondTypes: line %d\n", lineNum );
342	exit(8);
343	}
344
345	strcpy( currentBondType->type, the_token );
346
347	if( !strcmp( currentBondType->type, "fixed" ) ){
348	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
349	fprintf( stderr, "Error parseing BondTypes: line %d\n", lineNum );
350	exit(8);
351	}
352
353	sscanf( the_token, "%lf", &currentBondType->d0 );
354	}
355
356	else{
357	fprintf(stderr,
358	"Unknown TraPPE bond type \"%s\" at line %d\n",
359	currentBondType->type,
360	lineNum );
361	exit(8);
362	}
363	}
364	}
365
366	tempBondType = new LinkedType;
367	currentBondType->next = tempBondType;
368	currentBondType = tempBondType;
369
370	eof_test = fgets( readLine, sizeof(readLine), frcFile );
371	lineNum++;
372	}
373
374
375	// initialize the Bonds
376
377
378	for( i=0; i<nBonds; i++ ){
379
380	a = the_bonds[i].a;
381	b = the_bonds[i].b;
382
383	atomA = the_atoms[a]->getType();
384	atomB = the_atoms[b]->getType();
385	currentBondType = headBondType->find( atomA, atomB );
386	if( currentBondType == NULL ){
387	fprintf( stderr, "BondType error, %s - %s not found in force file.\n",
388	atomA, atomB );
389	exit(8);
390	}
391
392	if( !strcmp( currentBondType->type, "fixed" ) ){
393
394	the_sris[i] = new ConstrainedBond( *the_atoms[a],
395	*the_atoms[b],
396	currentBondType->d0 );
397	entry_plug->n_constraints++;
398	}
399	}
400
401
402	// clean up the memory
403
404	delete headBondType;
405	}
406
407	void TraPPEFF::initializeBends( bend_set* the_bends ){
408
409	class LinkedType {
410	public:
411	LinkedType(){
412	next = NULL;
413	nameA[0] = '\0';
414	nameB[0] = '\0';
415	nameC[0] = '\0';
416	type[0] = '\0';
417	}
418	~LinkedType(){ if( next != NULL ) delete next; }
419
420	LinkedType* find( char* key1, char* key2, char* key3 ){
421	if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 )
422	&& !strcmp( nameC, key3 ) ) return this;
423	if( !strcmp( nameA, key3 ) && !strcmp( nameB, key2 )
424	&& !strcmp( nameC, key1 ) ) return this;
425	if( next != NULL ) return next->find(key1, key2, key3);
426	return NULL;
427	}
428
429	char nameA[15];
430	char nameB[15];
431	char nameC[15];
432	char type[30];
433	double k1, k2, k3, t0;
434
435	LinkedType* next;
436	};
437
438	LinkedType* headBendType;
439	LinkedType* currentBendType;
440	LinkedType* tempBendType;
441
442	char readLine[500];
443	char* the_token;
444	char* eof_test;
445	int foundBend = 0;
446	int lineNum = 0;
447	int i, a, b, c, index;
448	char* atomA;
449	char* atomB;
450	char* atomC;
451	QuadraticBend* qBend;
452
453	SRI **the_sris;
454	Atom** the_atoms;
455	int nBends;
456	the_sris = entry_plug->sr_interactions;
457	the_atoms = entry_plug->atoms;
458	nBends = entry_plug->n_bends;
459
460	// read in the bend types.
461
462	rewind( frcFile );
463	headBendType = new LinkedType;
464	currentBendType = headBendType;
465
466	eof_test = fgets( readLine, sizeof(readLine), frcFile );
467	lineNum++;
468	if( eof_test == NULL ){
469	fprintf( stderr, "Error in reading Bends from force file.\n" );
470	exit(8);
471	}
472
473
474	while( !foundBend ){
475	while( eof_test != NULL && readLine[0] != '#' ){
476	eof_test = fgets( readLine, sizeof(readLine), frcFile );
477	lineNum++;
478	}
479	if( eof_test == NULL ){
480	fprintf( stderr,
481	"Error in reading Bends from force file at line %d.\n",
482	lineNum );
483	exit(8);
484	}
485
486	the_token = strtok( readLine, " ,;\t#\n" );
487	foundBend = !strcmp( "BendTypes", the_token );
488
489	if( !foundBend ){
490	eof_test = fgets( readLine, sizeof(readLine), frcFile );
491	lineNum++;
492
493	if( eof_test == NULL ){
494	fprintf( stderr,
495	"Error in reading Bends from force file at line %d.\n",
496	lineNum );
497	exit(8);
498	}
499	}
500	}
501
502	// we are now at the BendTypes section.
503
504	eof_test = fgets( readLine, sizeof(readLine), frcFile );
505	lineNum++;
506
507	if( eof_test == NULL ){
508	fprintf( stderr,
509	"Error in reading Bends from force file at line %d.\n",
510	lineNum );
511	exit(8);
512	}
513
514	while( readLine[0] != '#' && eof_test != NULL ){
515
516	if( readLine[0] != '!' ){
517
518	the_token = strtok( readLine, " \n\t,;" );
519	if( the_token != NULL ){
520
521	strcpy( currentBendType->nameA, the_token );
522
523	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
524	fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
525	exit(8);
526	}
527
528	strcpy( currentBendType->nameB, the_token );
529
530	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
531	fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
532	exit(8);
533	}
534
535	strcpy( currentBendType->nameC, the_token );
536
537	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
538	fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
539	exit(8);
540	}
541
542	strcpy( currentBendType->type, the_token );
543
544	if( !strcmp( currentBendType->type, "quadratic" ) ){
545	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
546	fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
547	exit(8);
548	}
549
550	sscanf( the_token, "%lf", &currentBendType->k1 );
551
552	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
553	fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
554	exit(8);
555	}
556
557	sscanf( the_token, "%lf", &currentBendType->k2 );
558
559	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
560	fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
561	exit(8);
562	}
563
564	sscanf( the_token, "%lf", &currentBendType->k3 );
565
566	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
567	fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
568	exit(8);
569	}
570
571	sscanf( the_token, "%lf", &currentBendType->t0 );
572	}
573
574	else{
575	fprintf(stderr,
576	"Unknown TraPPE bend type \"%s\" at line %d\n",
577	currentBendType->type,
578	lineNum );
579	exit(8);
580	}
581	}
582	}
583
584	tempBendType = new LinkedType;
585	currentBendType->next = tempBendType;
586	currentBendType = tempBendType;
587
588	eof_test = fgets( readLine, sizeof(readLine), frcFile );
589	lineNum++;
590	}
591
592
593	// initialize the Bends
594
595	for( i=0; i<nBends; i++ ){
596
597	a = the_bends[i].a;
598	b = the_bends[i].b;
599	c = the_bends[i].c;
600
601	atomA = the_atoms[a]->getType();
602	atomB = the_atoms[b]->getType();
603	atomC = the_atoms[c]->getType();
604	currentBendType = headBendType->find( atomA, atomB, atomC );
605	if( currentBendType == NULL ){
606	fprintf( stderr, "BendType error, %s - %s - %s not found"
607	" in force file.\n",
608	atomA, atomB, atomC );
609	exit(8);
610	}
611
612	if( !strcmp( currentBendType->type, "quadratic" ) ){
613
614	index = i + entry_plug->n_bonds;
615	qBend = new QuadraticBend( *the_atoms[a],
616	*the_atoms[b],
617	*the_atoms[c] );
618	qBend->setConstants( currentBendType->k1,
619	currentBendType->k2,
620	currentBendType->k3,
621	currentBendType->t0 );
622	the_sris[index] = qBend;
623	}
624	}
625
626
627	// clean up the memory
628
629	delete headBendType;
630	}
631
632	void TraPPEFF::initializeTorsions( torsion_set* the_torsions ){
633
634	class LinkedType {
635	public:
636	LinkedType(){
637	next = NULL;
638	nameA[0] = '\0';
639	nameB[0] = '\0';
640	nameC[0] = '\0';
641	type[0] = '\0';
642	}
643	~LinkedType(){ if( next != NULL ) delete next; }
644
645	LinkedType* find( char* key1, char* key2, char* key3, char* key4 ){
646	if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 ) &&
647	!strcmp( nameC, key3 ) && !strcmp( nameD, key4 ) ) return this;
648
649	if( !strcmp( nameA, key4 ) && !strcmp( nameB, key3 ) &&
650	!strcmp( nameC, key2 ) && !strcmp( nameD, key4 ) ) return this;
651
652	if( next != NULL ) return next->find(key1, key2, key3, key4);
653	return NULL;
654	}
655
656	char nameA[15];
657	char nameB[15];
658	char nameC[15];
659	char nameD[15];
660	char type[30];
661	double k1, k2, k3, k4;
662
663	LinkedType* next;
664	};
665
666	LinkedType* headTorsionType;
667	LinkedType* currentTorsionType;
668	LinkedType* tempTorsionType;
669
670	char readLine[500];
671	char* the_token;
672	char* eof_test;
673	int foundTorsion = 0;
674	int lineNum = 0;
675	int i, a, b, c, d, index;
676	char* atomA;
677	char* atomB;
678	char* atomC;
679	char* atomD;
680	CubicTorsion* cTors;
681
682	SRI **the_sris;
683	Atom** the_atoms;
684	int nTorsions;
685	the_sris = entry_plug->sr_interactions;
686	the_atoms = entry_plug->atoms;
687	nTorsions = entry_plug->n_torsions;
688
689	// read in the torsion types.
690
691	rewind( frcFile );
692	headTorsionType = new LinkedType;
693	currentTorsionType = headTorsionType;
694
695	eof_test = fgets( readLine, sizeof(readLine), frcFile );
696	lineNum++;
697	if( eof_test == NULL ){
698	fprintf( stderr, "Error in reading Torsions from force file.\n" );
699	exit(8);
700	}
701
702
703	while( !foundTorsion ){
704	while( eof_test != NULL && readLine[0] != '#' ){
705	eof_test = fgets( readLine, sizeof(readLine), frcFile );
706	lineNum++;
707	}
708	if( eof_test == NULL ){
709	fprintf( stderr,
710	"Error in reading Torsions from force file at line %d.\n",
711	lineNum );
712	exit(8);
713	}
714
715	the_token = strtok( readLine, " ,;\t#\n" );
716	foundTorsion = !strcmp( "TorsionTypes", the_token );
717
718	if( !foundTorsion ){
719	eof_test = fgets( readLine, sizeof(readLine), frcFile );
720	lineNum++;
721
722	if( eof_test == NULL ){
723	fprintf( stderr,
724	"Error in reading Torsions from force file at line %d.\n",
725	lineNum );
726	exit(8);
727	}
728	}
729	}
730
731	// we are now at the TorsionTypes section.
732
733	eof_test = fgets( readLine, sizeof(readLine), frcFile );
734	lineNum++;
735
736	if( eof_test == NULL ){
737	fprintf( stderr,
738	"Error in reading Torsions from force file at line %d.\n",
739	lineNum );
740	exit(8);
741	}
742
743	while( readLine[0] != '#' && eof_test != NULL ){
744
745	if( readLine[0] != '!' ){
746
747	the_token = strtok( readLine, " \n\t,;" );
748	if( the_token != NULL ){
749
750	strcpy( currentTorsionType->nameA, the_token );
751
752	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
753	fprintf( stderr, "Error parseing TorsionTypes: line %d\n", lineNum );
754	exit(8);
755	}
756
757	strcpy( currentTorsionType->nameB, the_token );
758
759	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
760	fprintf( stderr, "Error parseing TorsionTypes: line %d\n", lineNum );
761	exit(8);
762	}
763
764	strcpy( currentTorsionType->nameC, the_token );
765
766	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
767	fprintf( stderr, "Error parseing TorsionTypes: line %d\n", lineNum );
768	exit(8);
769	}
770
771	strcpy( currentTorsionType->nameD, the_token );
772
773	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
774	fprintf( stderr, "Error parseing TorsionTypes: line %d\n", lineNum );
775	exit(8);
776	}
777
778	strcpy( currentTorsionType->type, the_token );
779
780	if( !strcmp( currentTorsionType->type, "cubic" ) ){
781	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
782	fprintf(stderr,"Error parseing TorsionTypes: line %d\n", lineNum );
783	exit(8);
784	}
785
786	sscanf( the_token, "%lf", &currentTorsionType->k1 );
787
788	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
789	fprintf(stderr,"Error parseing TorsionTypes: line %d\n", lineNum );
790	exit(8);
791	}
792
793	sscanf( the_token, "%lf", &currentTorsionType->k2 );
794
795	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
796	fprintf(stderr,"Error parseing TorsionTypes: line %d\n", lineNum );
797	exit(8);
798	}
799
800	sscanf( the_token, "%lf", &currentTorsionType->k3 );
801
802	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
803	fprintf(stderr,"Error parseing TorsionTypes: line %d\n", lineNum );
804	exit(8);
805	}
806
807	sscanf( the_token, "%lf", &currentTorsionType->k4 );
808	}
809
810	else{
811	fprintf(stderr,
812	"Unknown TraPPE torsion type \"%s\" at line %d\n",
813	currentTorsionType->type,
814	lineNum );
815	exit(8);
816	}
817	}
818	}
819
820	tempTorsionType = new LinkedType;
821	currentTorsionType->next = tempTorsionType;
822	currentTorsionType = tempTorsionType;
823
824	eof_test = fgets( readLine, sizeof(readLine), frcFile );
825	lineNum++;
826	}
827
828
829	// initialize the Torsions
830
831	for( i=0; i<nTorsions; i++ ){
832
833	a = the_torsions[i].a;
834	b = the_torsions[i].b;
835	c = the_torsions[i].c;
836	d = the_torsions[i].d;
837
838	atomA = the_atoms[a]->getType();
839	atomB = the_atoms[b]->getType();
840	atomC = the_atoms[c]->getType();
841	atomD = the_atoms[d]->getType();
842	currentTorsionType = headTorsionType->find( atomA, atomB, atomC, atomD );
843	if( currentTorsionType == NULL ){
844	fprintf( stderr, "TorsionType error, %s - %s - %s - %s not found"
845	" in force file.\n",
846	atomA, atomB, atomC, atomD );
847	exit(8);
848	}
849
850	if( !strcmp( currentTorsionType->type, "cubic" ) ){
851	index = i + entry_plug->n_bonds + entry_plug->n_bends;
852
853	cTors = new CubicTorsion( the_atoms[a], the_atoms[b],
854	the_atoms[c], the_atoms[d] );
855	cTors->setConstants( currentTorsionType->k1, currentTorsionType->k2,
856	currentTorsionType->k3, currentTorsionType->k4 );
857	the_sris[index] = cTors;
858	}
859	}
860
861
862	// clean up the memory
863
864	delete headTorsionType;
865	}