mdtools/interface_implementation/TraPPEFF.cpp

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

#include <iostream>
using namespace std;

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


TraPPEFF::TraPPEFF(){

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