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:	11
Committed:	Tue Jul 9 18:40:59 2002 UTC (22 years ago) by mmeineke
File size:	20256 byte(s)
Log Message:	This commit was generated by cvs2svn to compensate for changes in r10, which included commits to RCS files with non-trunk default branches.
#	User	Rev	Content
1	mmeineke	10	#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			}