mdtools/libmdCode/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:	270
Committed:	Fri Feb 14 17:08:46 2003 UTC (21 years, 4 months ago) by mmeineke
File size:	20434 byte(s)
Log Message:	added libmdCode and a couple help scripts
#	User	Rev	Content
1	mmeineke	270	#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			}