[ViewVC] Diff of: group/trunk/OOPSE/libmdtools/TraPPE

Comparing trunk/OOPSE/libmdtools/TraPPE_ExFF.cpp (file contents):
Revision 388 by chuckv, Fri Mar 21 22:11:50 2003 UTC vs.
Revision 542 by mmeineke, Fri May 30 21:31:48 2003 UTC

+#include "fortranWrappers.hpp"
+#ifdef IS_MPI
-–
+#include <mpi++.h>
+#include "mpiForceField.h"
+#endif // is_mpi
+  MPI_Datatype mpiBondStructType;
+  MPI_Datatype mpiBendStructType;
+  MPI_Datatype mpiTorsionStructType;
-+
-+
+#endif
-+
-+
+  class LinkedAtomType {
-+
+  public:
-+
+    LinkedAtomType(){
-+
+      next = NULL;
-+
+      name[0] = '\0';
-+
+    }
-+
+    ~LinkedAtomType(){ if( next != NULL ) delete next; }
-+
-+
+    LinkedAtomType* find(char* key){
-+
+      if( !strcmp(name, key) ) return this;
-+
+      if( next != NULL ) return next->find(key);
-+
+      return NULL;
-+
+    }
-+
-+
+    void printMe( void ){
-+
-+
+      std::cerr << "LinkedAtype " << name << ": ident = " << ident << "\n";
-+
+      if( next != NULL ) next->printMe();
-+
-+
+    }
-+
-+
+    void add( atomStruct &info ){
-+
-+
+      // check for duplicates
-+
-+
+      if( !strcmp( info.name, name ) ){
-+
+        sprintf( painCave.errMsg,
-+
+                 "Duplicate TraPPE_Ex atom type \"%s\" found in "
-+
+                 "the TraPPE_ExFF param file./n",
-+
+                 name );
-+
+        painCave.isFatal = 1;
-+
+        simError();
-+
+      }
-+
-+
+      if( next != NULL ) next->add(info);
-+
+      else{
-+
+        next = new LinkedAtomType();
-+
+        strcpy(next->name, info.name);
-+
+        next->isDipole = info.isDipole;
-+
+        next->isSSD    = info.isSSD;
-+
+        next->mass     = info.mass;
-+
+        next->epslon   = info.epslon;
-+
+        next->sigma    = info.sigma;
-+
+        next->dipole   = info.dipole;
-+
+        next->w0       = info.w0;
-+
+        next->v0       = info.v0;
-+
+        next->ident    = info.ident;
-+
+      }
-+
+    }
-+
-+
+#ifdef IS_MPI
-+
-+
+    void duplicate( atomStruct &info ){
-+
+      strcpy(info.name, name);
-+
+      info.isDipole = isDipole;
-+
+      info.isSSD    = isSSD;
-+
+      info.mass     = mass;
-+
+      info.epslon   = epslon;
-+
+      info.sigma    = sigma;
-+
+      info.dipole   = dipole;
-+
+      info.w0       = w0;
-+
+      info.v0       = v0;
-+
+      info.ident    = ident;
-+
+      info.last     = 0;
-+
+    }
-+
-+
-+
+#endif
-+
-+
+    char name[15];
-+
+    int isDipole;
-+
+    int isSSD;
-+
+    double mass;
-+
+    double epslon;
-+
+    double sigma;
-+
+    double dipole;
-+
+    double w0;
-+
+    double v0;
-+
+    int ident;
-+
+    LinkedAtomType* next;
-+
+  };
-+
-+
+  class LinkedBondType {
-+
+  public:
-+
+    LinkedBondType(){
-+
+      next = NULL;
-+
+      nameA[0] = '\0';
-+
+      nameB[0] = '\0';
-+
+      type[0] = '\0';
-+
+    }
-+
+    ~LinkedBondType(){ if( next != NULL ) delete next; }
-+
-+
+    LinkedBondType* 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;
-+
+    }
-+
-+
-+
+    void add( bondStruct &info ){
-+
-+
+      // check for duplicates
-+
+      int dup = 0;
-+
-+
+      if( !strcmp(nameA, info.nameA ) && !strcmp( nameB, info.nameB ) ) dup = 1;
-+
+      if( !strcmp(nameA, info.nameB ) && !strcmp( nameB, info.nameA ) ) dup = 1;
-+
-+
+      if(dup){
-+
+        sprintf( painCave.errMsg,
-+
+                 "Duplicate TraPPE_Ex bond type \"%s - %s\" found in "
-+
+                 "the TraPPE_ExFF param file./n",
-+
+                 nameA, nameB );
-+
+        painCave.isFatal = 1;
-+
+        simError();
-+
+      }
-+
-+
-+
+      if( next != NULL ) next->add(info);
-+
+      else{
-+
+        next = new LinkedBondType();
-+
+        strcpy(next->nameA, info.nameA);
-+
+        strcpy(next->nameB, info.nameB);
-+
+        strcpy(next->type,  info.type);
-+
+        next->d0 = info.d0;
-+
+      }
-+
+    }
-+
-+
+#ifdef IS_MPI
-+
+    void duplicate( bondStruct &info ){
-+
+      strcpy(info.nameA, nameA);
-+
+      strcpy(info.nameB, nameB);
-+
+      strcpy(info.type,  type);
-+
+      info.d0   = d0;
-+
+      info.last = 0;
-+
+    }
-+
-+
-+
+#endif
-+
-+
+    char nameA[15];
-+
+    char nameB[15];
-+
+    char type[30];
-+
+    double d0;
-+
-+
+    LinkedBondType* next;
-+
+  };
-+
-+
+  class LinkedBendType {
-+
+  public:
-+
+    LinkedBendType(){
-+
+      next = NULL;
-+
+      nameA[0] = '\0';
-+
+      nameB[0] = '\0';
-+
+      nameC[0] = '\0';
-+
+      type[0] = '\0';
-+
+    }
-+
+    ~LinkedBendType(){ if( next != NULL ) delete next; }
-+
-+
+    LinkedBendType* 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;
-+
+    }
-+
-+
+    void add( bendStruct &info ){
-+
-+
+      // check for duplicates
-+
+      int dup = 0;
-+
-+
+      if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB )
-+
+          && !strcmp( nameC, info.nameC ) ) dup = 1;
-+
+      if( !strcmp( nameA, info.nameC ) && !strcmp( nameB, info.nameB )
-+
+          && !strcmp( nameC, info.nameA ) ) dup = 1;
-+
-+
+      if(dup){
-+
+        sprintf( painCave.errMsg,
-+
+                 "Duplicate TraPPE_Ex bend type \"%s - %s - %s\" found in "
-+
+                 "the TraPPE_ExFF param file./n",
-+
+                 nameA, nameB, nameC );
-+
+        painCave.isFatal = 1;
-+
+        simError();
-+
+      }
-+
+      if( next != NULL ) next->add(info);
-+
+      else{
-+
+        next = new LinkedBendType();
-+
+        strcpy(next->nameA, info.nameA);
-+
+        strcpy(next->nameB, info.nameB);
-+
+        strcpy(next->nameC, info.nameC);
-+
+        strcpy(next->type,  info.type);
-+
+        next->k1 = info.k1;
-+
+        next->k2 = info.k2;
-+
+        next->k3 = info.k3;
-+
+        next->t0 = info.t0;
-+
+      }
-+
+    }
-+
-+
+#ifdef IS_MPI
-+
-+
+    void duplicate( bendStruct &info ){
-+
+      strcpy(info.nameA, nameA);
-+
+      strcpy(info.nameB, nameB);
-+
+      strcpy(info.nameC, nameC);
-+
+      strcpy(info.type,  type);
-+
+      info.k1   = k1;
-+
+      info.k2   = k2;
-+
+      info.k3   = k3;
-+
+      info.t0   = t0;
-+
+      info.last = 0;
-+
+    }
-+
-+
+#endif // is_mpi
-+
-+
+    char nameA[15];
-+
+    char nameB[15];
-+
+    char nameC[15];
-+
+    char type[30];
-+
+    double k1, k2, k3, t0;
-+
-+
+    LinkedBendType* next;
-+
+  };
-+
-+
+  class LinkedTorsionType {
-+
+  public:
-+
+    LinkedTorsionType(){
-+
+      next = NULL;
-+
+      nameA[0] = '\0';
-+
+      nameB[0] = '\0';
-+
+      nameC[0] = '\0';
-+
+      type[0] = '\0';
-+
+    }
-+
+    ~LinkedTorsionType(){ if( next != NULL ) delete next; }
-+
-+
+    LinkedTorsionType* 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, key1 ) ) return this;
-+
-+
+      if( next != NULL ) return next->find(key1, key2, key3, key4);
-+
+      return NULL;
-+
+    }
-+
-+
+    void add( torsionStruct &info ){
-+
-+
+      // check for duplicates
-+
+      int dup = 0;
-+
-+
+      if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB ) &&
-+
+          !strcmp( nameC, info.nameC ) && !strcmp( nameD, info.nameD ) ) dup = 1;
-+
-+
+      if( !strcmp( nameA, info.nameD ) && !strcmp( nameB, info.nameC ) &&
-+
+          !strcmp( nameC, info.nameB ) && !strcmp( nameD, info.nameA ) ) dup = 1;
-+
-+
+      if(dup){
-+
+        sprintf( painCave.errMsg,
-+
+                 "Duplicate TraPPE_Ex torsion type \"%s - %s - %s - %s\" found in "
-+
+                 "the TraPPE_ExFF param file./n", nameA, nameB, nameC, nameD );
-+
+        painCave.isFatal = 1;
-+
+        simError();
-+
+      }
-+
-+
+      if( next != NULL ) next->add(info);
-+
+      else{
-+
+        next = new LinkedTorsionType();
-+
+        strcpy(next->nameA, info.nameA);
-+
+        strcpy(next->nameB, info.nameB);
-+
+        strcpy(next->nameC, info.nameC);
-+
+        strcpy(next->nameD, info.nameD);
-+
+        strcpy(next->type,  info.type);
-+
+        next->k1 = info.k1;
-+
+        next->k2 = info.k2;
-+
+        next->k3 = info.k3;
-+
+        next->k4 = info.k4;
-+
-+
+      }
-+
+    }
-+
-+
+#ifdef IS_MPI
-+
-+
+    void duplicate( torsionStruct &info ){
-+
+      strcpy(info.nameA, nameA);
-+
+      strcpy(info.nameB, nameB);
-+
+      strcpy(info.nameC, nameC);
-+
+      strcpy(info.nameD, nameD);
-+
+      strcpy(info.type,  type);
-+
+      info.k1   = k1;
-+
+      info.k2   = k2;
-+
+      info.k3   = k3;
-+
+      info.k4   = k4;
-+
+      info.last = 0;
-+
+    }
-+
+#endif
-+
+    char nameA[15];
-+
+    char nameB[15];
-+
+    char nameC[15];
-+
+    char nameD[15];
-+
+    char type[30];
-+
+    double k1, k2, k3, k4;
-+
-+
+    LinkedTorsionType* next;
-+
+  };
-+
-+
-+
+  LinkedAtomType* headAtomType;
-+
+  LinkedAtomType* currentAtomType;
-+
+  LinkedBondType* headBondType;
-+
+  LinkedBondType* currentBondType;
-+
+  LinkedBendType* headBendType;
-+
+  LinkedBendType* currentBendType;
-+
+  LinkedTorsionType* headTorsionType;
-+
+  LinkedTorsionType* currentTorsionType;
-+
+} // namespace
+using namespace TPE;
+  char temp[200];
+  char errMsg[1000];
-+
+  headAtomType       = NULL;
-+
+  currentAtomType    = NULL;
-+
+  headBondType       = NULL;
-+
+  currentBondType    = NULL;
-+
+  headBendType       = NULL;
-+
+  currentBendType    = NULL;
-+
+  headTorsionType    = NULL;
-+
+  currentTorsionType = NULL;
-+
+  // do the funtion wrapping
+  wrapMeFF( this );
+TraPPE_ExFF::~TraPPE_ExFF(){
-+
+  if( headAtomType != NULL ) delete headAtomType;
-+
+  if( headBondType != NULL ) delete headBondType;
-+
+  if( headBendType != NULL ) delete headBendType;
-+
+  if( headTorsionType != NULL ) delete headTorsionType;
-+
+#ifdef IS_MPI
+  if( worldRank == 0 ){
+#endif // is_mpi
+#endif // is_mpi
-+
+void TraPPE_ExFF::cleanMe( void ){
-+
-+
+#ifdef IS_MPI
-+
-+
+  // keep the linked lists in the mpi version
-+
-+
+#else // is_mpi
-+
-+
+  // delete the linked lists in the single processor version
-+
-+
+  if( headAtomType != NULL ) delete headAtomType;
-+
+  if( headBondType != NULL ) delete headBondType;
-+
+  if( headBendType != NULL ) delete headBendType;
-+
+  if( headTorsionType != NULL ) delete headTorsionType;
-+
-+
+#endif // is_mpi
-+
+}
-+
-+
+void TraPPE_ExFF::initForceField( int ljMixRule ){
+  initFortran( ljMixRule, entry_plug->useReactionField );
-<
+void TraPPE_ExFF::initializeAtoms( void ){
->
+void TraPPE_ExFF::readParams( void ){
->
->
+  int i, a, b, c, d;
->
+  int identNum;
->
+  char* atomA;
->
+  char* atomB;
->
+  char* atomC;
->
+  char* atomD;
-<
+  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;
-<
+    }
-<
-<
+    void add( atomStruct &info ){
-<
-<
+      // check for duplicates
-<
-<
+      if( !strcmp( info.name, name ) ){
-<
+        sprintf( painCave.errMsg,
-<
+                 "Duplicate TraPPE_Ex atom type \"%s\" found in "
-<
+                 "the TraPPE_ExFF param file./n",
-<
+                 name );
-<
+        painCave.isFatal = 1;
-<
+        simError();
-<
+      }
-<
-<
+      if( next != NULL ) next->add(info);
-<
+      else{
-<
+        next = new LinkedType();
-<
+        strcpy(next->name, info.name);
-<
+        next->isDipole = info.isDipole;
-<
+        next->isSSD    = info.isSSD;
-<
+        next->mass     = info.mass;
-<
+        next->epslon   = info.epslon;
-<
+        next->sigma    = info.sigma;
-<
+        next->dipole   = info.dipole;
-<
+        next->w0       = info.w0;
-<
+        next->v0       = info.v0;
-<
+        next->ident    = info.ident;
-<
+      }
-<
+    }
-<
-<
+#ifdef IS_MPI
-<
-<
+    void duplicate( atomStruct &info ){
-<
+      strcpy(info.name, name);
-<
+      info.isDipole = isDipole;
-<
+      info.isSSD    = isSSD;
-<
+      info.mass     = mass;
-<
+      info.epslon   = epslon;
-<
+      info.sigma    = sigma;
-<
+      info.dipole   = dipole;
-<
+      info.w0       = w0;
-<
+      info.v0       = v0;
-<
+      info.last     = 0;
-<
+    }
-<
-<
-<
+#endif
-<
-<
+    char name[15];
-<
+    int isDipole;
-<
+    int isSSD;
-<
+    double mass;
-<
+    double epslon;
-<
+    double sigma;
-<
+    double dipole;
-<
+    double w0;
-<
+    double v0;
-<
+    int ident;
-<
+    LinkedType* next;
-<
+  };
->
+  atomStruct atomInfo;
->
+  bondStruct bondInfo;
->
+  bendStruct bendInfo;
->
+  torsionStruct torsionInfo;
-<
+  LinkedType* headAtomType;
-<
+  LinkedType* currentAtomType;
-<
+  atomStruct info;
-<
+  info.last = 1; // initialize last to have the last set.
-<
+                 // if things go well, last will be set to 0
->
+  bigSigma = 0.0;
-<
->
+  atomInfo.last = 1;
->
+  bondInfo.last = 1;
->
+  bendInfo.last = 1;
->
+  torsionInfo.last = 1;
-<
+  int i;
-<
+  int identNum;
->
+  // read in the atom info
-–
+  Atom** the_atoms;
-–
+  int nAtoms;
-–
+  the_atoms = entry_plug->atoms;
-–
+  nAtoms = entry_plug->n_atoms;
-–
-–
+  //////////////////////////////////////////////////
-–
+  // a quick water fix
-–
-–
+  double waterI[3][3];
-–
+  waterI[0][0] = 1.76958347772500;
-–
+  waterI[0][1] = 0.0;
-–
+  waterI[0][2] = 0.0;
-–
-–
+  waterI[1][0] = 0.0;
-–
+  waterI[1][1] = 0.614537057924513;
-–
+  waterI[1][2] = 0.0;
-–
-–
+  waterI[2][0] = 0.0;
-–
+  waterI[2][1] = 0.0;
-–
+  waterI[2][2] = 1.15504641980049;
-–
-–
-–
+  double headI[3][3];
-–
+  headI[0][0] = 1125;
-–
+  headI[0][1] = 0.0;
-–
+  headI[0][2] = 0.0;
-–
-–
+  headI[1][0] = 0.0;
-–
+  headI[1][1] = 1125;
-–
+  headI[1][2] = 0.0;
-–
-–
+  headI[2][0] = 0.0;
-–
+  headI[2][1] = 0.0;
-–
+  headI[2][2] = 250;
-–
-–
-–
-–
+  //////////////////////////////////////////////////
-–
-–
+#ifdef IS_MPI
+  if( worldRank == 0 ){
+#endif
+    // read in the atom types.
-–
-–
+    headAtomType = new LinkedType;
-+
+    headAtomType = new LinkedAtomType;
-+
+    fastForward( "AtomTypes", "initializeAtoms" );
+    // we are now at the AtomTypes section.
+      if( readLine[0] != '!' ){
+        // the parser returns 0 if the line was blank
-<
+        if( parseAtom( readLine, lineNum, info ) ){
-<
+          info.ident = identNum;
-<
+          headAtomType->add( info );;
->
+        if( parseAtom( readLine, lineNum, atomInfo ) ){
->
+          atomInfo.ident = identNum;
->
+          headAtomType->add( atomInfo );;
+          identNum++;
+    currentAtomType = headAtomType->next; //skip the first element who is a place holder.
+    while( currentAtomType != NULL ){
-<
+      currentAtomType->duplicate( info );
->
+      currentAtomType->duplicate( atomInfo );
-<
+      sendFrcStruct( &info, mpiAtomStructType );
->
+      sendFrcStruct( &atomInfo, mpiAtomStructType );
+      sprintf( checkPointMsg,
+               "successfully sent TraPPE_Ex force type: \"%s\"\n",
-<
+               info.name );
->
+               atomInfo.name );
+      MPIcheckPoint();
+      currentAtomType = currentAtomType->next;
-<
+    info.last = 1;
-<
+    sendFrcStruct( &info, mpiAtomStructType );
->
+    atomInfo.last = 1;
->
+    sendFrcStruct( &atomInfo, mpiAtomStructType );
+    MPIcheckPoint();
-<
+    headAtomType = new LinkedType;
-<
+    recieveFrcStruct( &info, mpiAtomStructType );
->
+    headAtomType = new LinkedAtomType;
->
+    recieveFrcStruct( &atomInfo, mpiAtomStructType );
-<
+    while( !info.last ){
->
+    while( !atomInfo.last ){
-<
+      headAtomType->add( info );
->
+      headAtomType->add( atomInfo );
+      MPIcheckPoint();
-<
+      recieveFrcStruct( &info, mpiAtomStructType );
->
+      recieveFrcStruct( &atomInfo, mpiAtomStructType );
-+
+#endif // is_mpi
-+
-+
+  // call new A_types in fortran
+  int isError;
+  double GB_dummy = 0.0;
-<
+  currentAtomType = headAtomType;
->
+  currentAtomType = headAtomType->next;;
+  while( currentAtomType != NULL ){
-–
+    if(currentAtomType->isDipole) entry_plug->useReactionField = 1;
+    if(currentAtomType->isDipole) entry_plug->useDipole = 1;
-<
+    if(currentAtomType->isSSD)    entry_plug->useSticky = 1;
->
+    if(currentAtomType->isSSD) {
->
+      entry_plug->useSticky = 1;
->
+      set_sticky_params( &(currentAtomType->w0), &(currentAtomType->v0));
->
+    }
+    if( currentAtomType->name[0] != '\0' ){
+      isError = 0;
+                 &(currentAtomType->epslon),
+                 &(currentAtomType->sigma),
+                 &(currentAtomType->dipole),
-–
+                 &(currentAtomType->w0),
-–
+                 &(currentAtomType->v0),
-–
+                 &GB_dummy,
-–
+                 &GB_dummy,
-–
+                 &GB_dummy,
-–
+                 &GB_dummy,
-–
+                 &GB_dummy,
-–
+                 &GB_dummy,
+                 &isError );
+      if( isError ){
+        sprintf( painCave.errMsg,
+#endif // is_mpi
-–
+  // initialize the atoms
-–
-–
+  double bigSigma = 0.0;
-–
+  DirectionalAtom* dAtom;
-<
+  for( i=0; i<nAtoms; i++ ){
-<
-<
+    currentAtomType = headAtomType->find( the_atoms[i]->getType() );
-<
+    if( currentAtomType == NULL ){
-<
+      sprintf( painCave.errMsg,
-<
+               "AtomType error, %s not found in force file.\n",
-<
+               the_atoms[i]->getType() );
-<
+      painCave.isFatal = 1;
-<
+      simError();
-<
+    }
-<
-<
+    the_atoms[i]->setMass( currentAtomType->mass );
-<
+    the_atoms[i]->setEpslon( currentAtomType->epslon );
-<
+    the_atoms[i]->setSigma( currentAtomType->sigma );
-<
+    the_atoms[i]->setIdent( currentAtomType->ident );
-<
+    the_atoms[i]->setLJ();
-<
-<
+    if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
-<
-<
+    if( currentAtomType->isDipole ){
-<
+      if( the_atoms[i]->isDirectional() ){
-<
-<
+        dAtom = (DirectionalAtom *) the_atoms[i];
-<
+        dAtom->setMu( currentAtomType->dipole );
-<
+        dAtom->setHasDipole( 1 );
-<
+        dAtom->setJx( 0.0 );
-<
+        dAtom->setJy( 0.0 );
-<
+        dAtom->setJz( 0.0 );
-<
-<
+        if(!strcmp("SSD",the_atoms[i]->getType())){
-<
+          dAtom->setI( waterI );
-<
+          dAtom->setSSD( 1 );
-<
+        }
-<
+        else if(!strcmp("HEAD",the_atoms[i]->getType())){
-<
+          dAtom->setI( headI );
-<
+          dAtom->setSSD( 0 );
-<
+        }
-<
+        else{
-<
+          sprintf(painCave.errMsg,
-<
+                  "AtmType error, %s does not have a moment of inertia set.\n",
-<
+                  the_atoms[i]->getType() );
-<
+          painCave.isFatal = 1;
-<
+          simError();
-<
+        }
-<
+        entry_plug->n_dipoles++;
-<
+      }
-<
+      else{
-<
-<
+        sprintf( painCave.errMsg,
-<
+                "TraPPE_ExFF error: Atom \"%s\" is a dipole, yet no standard"
-<
+                 " orientation was specifed in the BASS file.\n",
-<
+                 currentAtomType->name );
-<
+        painCave.isFatal = 1;
-<
+        simError();
-<
+      }
-<
+    }
-<
+    else{
-<
+      if( the_atoms[i]->isDirectional() ){
-<
+        sprintf( painCave.errMsg,
-<
+                 "TraPPE_ExFF error: Atom \"%s\" was given a standard"
-<
+                 "orientation in the BASS file, yet it is not a dipole.\n",
-<
+                 currentAtomType->name);
-<
+        painCave.isFatal = 1;
-<
+        simError();
-<
+      }
-<
+    }
-<
+  }
-<
-<
+#ifdef IS_MPI
-<
+  double tempBig = bigSigma;
-<
+  MPI::COMM_WORLD.Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX );
-<
+#endif  //is_mpi
-<
-<
+  //calc rCut and rList
-<
-<
+  entry_plug->rCut = 2.5 * bigSigma;
-<
-<
+  if(entry_plug->rCut > (entry_plug->box_x / 2.0))
-<
+    entry_plug->rCut = entry_plug->box_x / 2.0;
-<
-<
+  if(entry_plug->rCut > (entry_plug->box_y / 2.0))
-<
+    entry_plug->rCut = entry_plug->box_y / 2.0;
-<
-<
+  if(entry_plug->rCut > (entry_plug->box_z / 2.0))
-<
+    entry_plug->rCut = entry_plug->box_z / 2.0;
-<
-<
+  entry_plug->rList = entry_plug->rCut + 1.0;
-<
-<
+  entry_plug->useLJ = 1; // use Lennard Jones is on by default
-<
-<
+  // clean up the memory
-<
-<
+  delete headAtomType;
-<
-<
+#ifdef IS_MPI
-<
+  sprintf( checkPointMsg, "TraPPE_Ex atoms initialized succesfully" );
-<
+  MPIcheckPoint();
-<
+#endif // is_mpi
-<
-<
+}
-<
-<
+void TraPPE_ExFF::initializeBonds( bond_pair* the_bonds ){
-<
-<
+  class LinkedType {
-<
+  public:
-<
+    LinkedType(){
-<
+      next = NULL;
-<
+      nameA[0] = '\0';
-<
+      nameB[0] = '\0';
-<
+      type[0] = '\0';
-<
+    }
-<
+    ~LinkedType(){ if( next != NULL ) delete next; }
-<
-<
+    LinkedType* find(char* key1, char* key2){
-<
+      if( !strcmp(nameA, key1 ) && !strcmp( nameB, key2 ) ) return this;
-<
+      if( !strcmp(nameA, key2 ) && !strcmp( nameB, key1 ) ) return this;
-<
+      if( next != NULL ) return next->find(key1, key2);
-<
+      return NULL;
-<
+    }
-<
-<
-<
+    void add( bondStruct &info ){
-<
-<
+      // check for duplicates
-<
+      int dup = 0;
-<
-<
+      if( !strcmp(nameA, info.nameA ) && !strcmp( nameB, info.nameB ) ) dup = 1;
-<
+      if( !strcmp(nameA, info.nameB ) && !strcmp( nameB, info.nameA ) ) dup = 1;
-<
-<
+      if(dup){
-<
+        sprintf( painCave.errMsg,
-<
+                 "Duplicate TraPPE_Ex bond type \"%s - %s\" found in "
-<
+                 "the TraPPE_ExFF param file./n",
-<
+                 nameA, nameB );
-<
+        painCave.isFatal = 1;
-<
+        simError();
-<
+      }
-<
-<
-<
+      if( next != NULL ) next->add(info);
-<
+      else{
-<
+        next = new LinkedType();
-<
+        strcpy(next->nameA, info.nameA);
-<
+        strcpy(next->nameB, info.nameB);
-<
+        strcpy(next->type,  info.type);
-<
+        next->d0 = info.d0;
-<
+      }
-<
+    }
-<
-<
+#ifdef IS_MPI
-<
+    void duplicate( bondStruct &info ){
-<
+      strcpy(info.nameA, nameA);
-<
+      strcpy(info.nameB, nameB);
-<
+      strcpy(info.type,  type);
-<
+      info.d0   = d0;
-<
+      info.last = 0;
-<
+    }
-<
-<
-<
+#endif
-<
-<
+    char nameA[15];
-<
+    char nameB[15];
-<
+    char type[30];
-<
+    double d0;
-<
-<
+    LinkedType* next;
-<
+  };
-<
-<
->
+  // read in the bonds
-–
+  LinkedType* headBondType;
-–
+  LinkedType* currentBondType;
-–
+  bondStruct info;
-–
+  info.last = 1; // initialize last to have the last set.
-–
+                 // if things go well, last will be set to 0
-–
-–
+  SRI **the_sris;
-–
+  Atom** the_atoms;
-–
+  int nBonds;
-–
+  the_sris = entry_plug->sr_interactions;
-–
+  the_atoms = entry_plug->atoms;
-–
+  nBonds = entry_plug->n_bonds;
-–
-–
+  int i, a, b;
-–
+  char* atomA;
-–
+  char* atomB;
-–
+#ifdef IS_MPI
+  if( worldRank == 0 ){
+#endif
+    // read in the bond types.
-<
+    headBondType = new LinkedType;
->
+    headBondType = new LinkedBondType;
+    fastForward( "BondTypes", "initializeBonds" );
+      if( readLine[0] != '!' ){
+        // the parser returns 0 if the line was blank
-<
+        if( parseBond( readLine, lineNum, info ) ){
-<
+          headBondType->add( info );
->
+        if( parseBond( readLine, lineNum, bondInfo ) ){
->
+          headBondType->add( bondInfo );
+      eof_test = fgets( readLine, sizeof(readLine), frcFile );
+             "TraPPE_Ex bond structures read successfully." );
+    MPIcheckPoint();
-<
+    currentBondType = headBondType;
->
+    currentBondType = headBondType->next;
+    while( currentBondType != NULL ){
-<
+      currentBondType->duplicate( info );
-<
+      sendFrcStruct( &info, mpiBondStructType );
->
+      currentBondType->duplicate( bondInfo );
->
+      sendFrcStruct( &bondInfo, mpiBondStructType );
+      currentBondType = currentBondType->next;
-<
+    info.last = 1;
-<
+    sendFrcStruct( &info, mpiBondStructType );
->
+    bondInfo.last = 1;
->
+    sendFrcStruct( &bondInfo, mpiBondStructType );
+    // listen for node 0 to send out the force params
-<
+    MPIcheckPoint();
->
+    MPIcheckPoint();
-<
+    headBondType = new LinkedType;
-<
+    recieveFrcStruct( &info, mpiBondStructType );
-<
+    while( !info.last ){
->
+    headBondType = new LinkedBondType;
->
+    recieveFrcStruct( &bondInfo, mpiBondStructType );
->
+    while( !bondInfo.last ){
-<
+      headBondType->add( info );
-<
+      recieveFrcStruct( &info, mpiBondStructType );
->
+      headBondType->add( bondInfo );
->
+      recieveFrcStruct( &bondInfo, mpiBondStructType );
-–
+#endif // is_mpi
-–
-–
-–
+  // 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 ){
-<
+      sprintf( painCave.errMsg,
-<
+               "BondType error, %s - %s not found in force file.\n",
-<
+               atomA, atomB );
-<
+      painCave.isFatal = 1;
-<
+      simError();
-<
+    }
-<
-<
+    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;
-<
-<
+#ifdef IS_MPI
-<
+  sprintf( checkPointMsg, "TraPPE_Ex bonds initialized succesfully" );
->
+  sprintf( checkPointMsg,
->
+           "TraPPE_ExFF bond structures broadcast successfully." );
+  MPIcheckPoint();
-–
+#endif // is_mpi
-–
+}
-–
-–
+void TraPPE_ExFF::initializeBends( bend_set* the_bends ){
-–
-–
+  class LinkedType {
-–
+  public:
-–
+    LinkedType(){
-–
+      next = NULL;
-–
+      nameA[0] = '\0';
-–
+      nameB[0] = '\0';
-–
+      nameC[0] = '\0';
-–
+      type[0] = '\0';
-–
+    }
-–
+    ~LinkedType(){ if( next != NULL ) delete next; }
-–
-–
+    LinkedType* find( char* key1, char* key2, char* key3 ){
-–
+      if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 )
-–
+          && !strcmp( nameC, key3 ) ) return this;
-–
+      if( !strcmp( nameA, key3 ) && !strcmp( nameB, key2 )
-–
+          && !strcmp( nameC, key1 ) ) return this;
-–
+      if( next != NULL ) return next->find(key1, key2, key3);
-–
+      return NULL;
-–
+    }
-–
-–
+    void add( bendStruct &info ){
-–
-–
+      // check for duplicates
-–
+      int dup = 0;
-–
-–
+      if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB )
-–
+          && !strcmp( nameC, info.nameC ) ) dup = 1;
-–
+      if( !strcmp( nameA, info.nameC ) && !strcmp( nameB, info.nameB )
-–
+          && !strcmp( nameC, info.nameA ) ) dup = 1;
-–
-–
+      if(dup){
-–
+        sprintf( painCave.errMsg,
-–
+                 "Duplicate TraPPE_Ex bend type \"%s - %s - %s\" found in "
-–
+                 "the TraPPE_ExFF param file./n",
-–
+                 nameA, nameB, nameC );
-–
+        painCave.isFatal = 1;
-–
+        simError();
-–
+      }
-–
-–
+      if( next != NULL ) next->add(info);
-–
+      else{
-–
+        next = new LinkedType();
-–
+        strcpy(next->nameA, info.nameA);
-–
+        strcpy(next->nameB, info.nameB);
-–
+        strcpy(next->nameC, info.nameC);
-–
+        strcpy(next->type,  info.type);
-–
+        next->k1 = info.k1;
-–
+        next->k2 = info.k2;
-–
+        next->k3 = info.k3;
-–
+        next->t0 = info.t0;
-–
+      }
-–
+    }
-–
-–
+#ifdef IS_MPI
-–
-–
+    void duplicate( bendStruct &info ){
-–
+      strcpy(info.nameA, nameA);
-–
+      strcpy(info.nameB, nameB);
-–
+      strcpy(info.nameC, nameC);
-–
+      strcpy(info.type,  type);
-–
+      info.k1   = k1;
-–
+      info.k2   = k2;
-–
+      info.k3   = k3;
-–
+      info.t0   = t0;
-–
+      info.last = 0;
-–
+    }
-–
+#endif // is_mpi
-–
-–
+    char nameA[15];
-–
+    char nameB[15];
-–
+    char nameC[15];
-–
+    char type[30];
-–
+    double k1, k2, k3, t0;
-–
-–
+    LinkedType* next;
-–
+  };
-–
+  LinkedType* headBendType;
-–
+  LinkedType* currentBendType;
-–
+  bendStruct info;
-–
+  info.last = 1; // initialize last to have the last set.
-–
+                 // if things go well, last will be set to 0
-<
+  QuadraticBend* qBend;
-<
+  GhostBend* gBend;
-<
+  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 bends
-–
+  int i, a, b, c;
-–
+  char* atomA;
-–
+  char* atomB;
-–
+  char* atomC;
-–
-–
+#ifdef IS_MPI
+  if( worldRank == 0 ){
+#endif
+    // read in the bend types.
-<
+    headBendType = new LinkedType;
->
+    headBendType = new LinkedBendType;
+    fastForward( "BendTypes", "initializeBends" );
+      if( readLine[0] != '!' ){
+        // the parser returns 0 if the line was blank
-<
+        if( parseBend( readLine, lineNum, info ) ){
-<
+          headBendType->add( info );
->
+        if( parseBend( readLine, lineNum, bendInfo ) ){
->
+          headBendType->add( bendInfo );
+      eof_test = fgets( readLine, sizeof(readLine), frcFile );
+             "TraPPE_Ex bend structures read successfully." );
+    MPIcheckPoint();
-<
+    currentBendType = headBendType;
->
+    currentBendType = headBendType->next;
+    while( currentBendType != NULL ){
-<
+      currentBendType->duplicate( info );
-<
+      sendFrcStruct( &info, mpiBendStructType );
->
+      currentBendType->duplicate( bendInfo );
->
+      sendFrcStruct( &bendInfo, mpiBendStructType );
+      currentBendType = currentBendType->next;
-<
+    info.last = 1;
-<
+    sendFrcStruct( &info, mpiBendStructType );
->
+    bendInfo.last = 1;
->
+    sendFrcStruct( &bendInfo, mpiBendStructType );
+    MPIcheckPoint();
-<
+    headBendType = new LinkedType;
-<
+    recieveFrcStruct( &info, mpiBendStructType );
-<
+    while( !info.last ){
->
+    headBendType = new LinkedBendType;
->
+    recieveFrcStruct( &bendInfo, mpiBendStructType );
->
+    while( !bendInfo.last ){
-<
+      headBendType->add( info );
-<
+      recieveFrcStruct( &info, mpiBendStructType );
->
+      headBendType->add( bendInfo );
->
+      recieveFrcStruct( &bendInfo, mpiBendStructType );
-+
-+
+  sprintf( checkPointMsg,
-+
+           "TraPPE_ExFF bend structures broadcast successfully." );
-+
+  MPIcheckPoint();
-+
+#endif // is_mpi
-+
-+
-+
+  // read in the torsions
-+
-+
+#ifdef IS_MPI
-+
+  if( worldRank == 0 ){
-+
+#endif
-+
-+
+    // read in the torsion types.
-+
-+
+    headTorsionType = new LinkedTorsionType;
-+
-+
+    fastForward( "TorsionTypes", "initializeTorsions" );
-+
-+
+    // we are now at the torsionTypes section
-+
-+
+    eof_test =  fgets( readLine, sizeof(readLine), frcFile );
-+
+    lineNum++;
-+
-+
-+
+    // read a line, and start parseing out the atom types
-+
-+
+    if( eof_test == NULL ){
-+
+      sprintf( painCave.errMsg,
-+
+               "Error in reading torsions from force file at line %d.\n",
-+
+               lineNum );
-+
+      painCave.isFatal = 1;
-+
+      simError();
-+
+    }
-+
-+
+    // stop reading at end of file, or at next section
-+
+    while( readLine[0] != '#' && eof_test != NULL ){
-+
-+
+      // toss comment lines
-+
+      if( readLine[0] != '!' ){
-+
-+
+        // the parser returns 0 if the line was blank
-+
+        if( parseTorsion( readLine, lineNum, torsionInfo ) ){
-+
+          headTorsionType->add( torsionInfo );
-+
-+
+        }
-+
+      }
-+
+      eof_test = fgets( readLine, sizeof(readLine), frcFile );
-+
+      lineNum++;
-+
+    }
-+
-+
+#ifdef IS_MPI
-+
-+
+    // send out the linked list to all the other processes
-+
-+
+    sprintf( checkPointMsg,
-+
+             "TraPPE_Ex torsion structures read successfully." );
-+
+    MPIcheckPoint();
-+
-+
+    currentTorsionType = headTorsionType->next;
-+
+    while( currentTorsionType != NULL ){
-+
+      currentTorsionType->duplicate( torsionInfo );
-+
+      sendFrcStruct( &torsionInfo, mpiTorsionStructType );
-+
+      currentTorsionType = currentTorsionType->next;
-+
+    }
-+
+    torsionInfo.last = 1;
-+
+    sendFrcStruct( &torsionInfo, mpiTorsionStructType );
-+
-+
+  }
-+
-+
+  else{
-+
-+
+    // listen for node 0 to send out the force params
-+
-+
+    MPIcheckPoint();
-+
-+
+    headTorsionType = new LinkedTorsionType;
-+
+    recieveFrcStruct( &torsionInfo, mpiTorsionStructType );
-+
+    while( !torsionInfo.last ){
-+
-+
+      headTorsionType->add( torsionInfo );
-+
+      recieveFrcStruct( &torsionInfo, mpiTorsionStructType );
-+
+    }
-+
+  }
-+
-+
+  sprintf( checkPointMsg,
-+
+           "TraPPE_ExFF torsion structures broadcast successfully." );
-+
+  MPIcheckPoint();
-+
-+
+#endif // is_mpi
-+
-+
+  entry_plug->useLJ = 1;
-+
+}
-+
-+
-+
-+
+void TraPPE_ExFF::initializeAtoms( int nAtoms, Atom** the_atoms ){
-–
+  // initialize the Bends
-<
+  int index;
->
+  //////////////////////////////////////////////////
->
+  // a quick water fix
->
->
+  double waterI[3][3];
->
+  waterI[0][0] = 1.76958347772500;
->
+  waterI[0][1] = 0.0;
->
+  waterI[0][2] = 0.0;
->
->
+  waterI[1][0] = 0.0;
->
+  waterI[1][1] = 0.614537057924513;
->
+  waterI[1][2] = 0.0;
->
->
+  waterI[2][0] = 0.0;
->
+  waterI[2][1] = 0.0;
->
+  waterI[2][2] = 1.15504641980049;
->
->
->
+  double headI[3][3];
->
+  headI[0][0] = 1125;
->
+  headI[0][1] = 0.0;
->
+  headI[0][2] = 0.0;
->
->
+  headI[1][0] = 0.0;
->
+  headI[1][1] = 1125;
->
+  headI[1][2] = 0.0;
->
->
+  headI[2][0] = 0.0;
->
+  headI[2][1] = 0.0;
->
+  headI[2][2] = 250;
->
->
+  //////////////////////////////////////////////////
->
->
->
+  // initialize the atoms
->
->
+  DirectionalAtom* dAtom;
->
->
+  for(int i=0; i<nAtoms; i++ ){
-+
+    currentAtomType = headAtomType->find( the_atoms[i]->getType() );
-+
+    if( currentAtomType == NULL ){
-+
+      sprintf( painCave.errMsg,
-+
+               "AtomType error, %s not found in force file.\n",
-+
+               the_atoms[i]->getType() );
-+
+      painCave.isFatal = 1;
-+
+      simError();
-+
+    }
-+
-+
+    the_atoms[i]->setMass( currentAtomType->mass );
-+
+    the_atoms[i]->setEpslon( currentAtomType->epslon );
-+
+    the_atoms[i]->setSigma( currentAtomType->sigma );
-+
+    the_atoms[i]->setIdent( currentAtomType->ident );
-+
+    the_atoms[i]->setLJ();
-+
-+
+    if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
-+
-+
+    if( currentAtomType->isDipole ){
-+
+      if( the_atoms[i]->isDirectional() ){
-+
-+
+        dAtom = (DirectionalAtom *) the_atoms[i];
-+
+        dAtom->setMu( currentAtomType->dipole );
-+
+        dAtom->setHasDipole( 1 );
-+
+        dAtom->setJx( 0.0 );
-+
+        dAtom->setJy( 0.0 );
-+
+        dAtom->setJz( 0.0 );
-+
-+
+        if(!strcmp("SSD",the_atoms[i]->getType())){
-+
+          dAtom->setI( waterI );
-+
+          dAtom->setSSD( 1 );
-+
+        }
-+
+        else if(!strcmp("HEAD",the_atoms[i]->getType())){
-+
+          dAtom->setI( headI );
-+
+          dAtom->setSSD( 0 );
-+
+        }
-+
+        else{
-+
+          sprintf(painCave.errMsg,
-+
+                  "AtmType error, %s does not have a moment of inertia set.\n",
-+
+                  the_atoms[i]->getType() );
-+
+          painCave.isFatal = 1;
-+
+          simError();
-+
+        }
-+
+        entry_plug->n_dipoles++;
-+
+      }
-+
+      else{
-+
-+
+        sprintf( painCave.errMsg,
-+
+                "TraPPE_ExFF error: Atom \"%s\" is a dipole, yet no standard"
-+
+                 " orientation was specifed in the BASS file.\n",
-+
+                 currentAtomType->name );
-+
+        painCave.isFatal = 1;
-+
+        simError();
-+
+      }
-+
+    }
-+
+    else{
-+
+      if( the_atoms[i]->isDirectional() ){
-+
+        sprintf( painCave.errMsg,
-+
+                 "TraPPE_ExFF error: Atom \"%s\" was given a standard"
-+
+                 "orientation in the BASS file, yet it is not a dipole.\n",
-+
+                 currentAtomType->name);
-+
+        painCave.isFatal = 1;
-+
+        simError();
-+
+      }
-+
+    }
-+
+  }
-+
+}
-+
-+
+void TraPPE_ExFF::initializeBonds( int nBonds, Bond** bondArray,
-+
+                                   bond_pair* the_bonds ){
-+
+  int i,a,b;
-+
+  char* atomA;
-+
+  char* atomB;
-+
-+
+  Atom** the_atoms;
-+
+  the_atoms = entry_plug->atoms;
-+
-+
-+
+  // 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 ){
-+
+      sprintf( painCave.errMsg,
-+
+               "BondType error, %s - %s not found in force file.\n",
-+
+               atomA, atomB );
-+
+      painCave.isFatal = 1;
-+
+      simError();
-+
+    }
-+
-+
+    if( !strcmp( currentBondType->type, "fixed" ) ){
-+
-+
+      bondArray[i] = new ConstrainedBond( *the_atoms[a],
-+
+                                          *the_atoms[b],
-+
+                                          currentBondType->d0 );
-+
+      entry_plug->n_constraints++;
-+
+    }
-+
+  }
-+
+}
-+
-+
+void TraPPE_ExFF::initializeBends( int nBends, Bend** bendArray,
-+
+                                   bend_set* the_bends ){
-+
-+
+  QuadraticBend* qBend;
-+
+  GhostBend* gBend;
-+
+  Atom** the_atoms;
-+
+  the_atoms = entry_plug->atoms;
-+
-+
+  int i, a, b, c;
-+
+  char* atomA;
-+
+  char* atomB;
-+
+  char* atomC;
-+
-+
+  // initialize the Bends
-+
+  for( i=0; i<nBends; i++ ){
+    a = the_bends[i].a;
+    if( !strcmp( currentBendType->type, "quadratic" ) ){
-–
+      index = i + entry_plug->n_bonds;
-–
+      if( the_bends[i].isGhost){
+        if( the_bends[i].ghost == b ){
+        else if( the_bends[i].ghost == a ){
+          c = a;
+          a = b;
-<
+          b = a;
->
+          b = c;
+        else{
+          sprintf( painCave.errMsg,
+                             currentBendType->k2,
+                             currentBendType->k3,
+                             currentBendType->t0 );
-<
+        the_sris[index] = gBend;
->
+        bendArray[i] = gBend;
+      else{
+        qBend = new QuadraticBend( *the_atoms[a],
+                             currentBendType->k2,
+                             currentBendType->k3,
+                             currentBendType->t0 );
-<
+        the_sris[index] = qBend;
->
+        bendArray[i] = qBend;
-–
-–
-–
+  // clean up the memory
-–
-–
+  delete headBendType;
-–
-–
+#ifdef IS_MPI
-–
+  sprintf( checkPointMsg, "TraPPE_Ex bends initialized succesfully" );
-–
+  MPIcheckPoint();
-–
+#endif // is_mpi
-–
-<
+void TraPPE_ExFF::initializeTorsions( torsion_set* the_torsions ){
-<
-<
+  class LinkedType {
-<
+  public:
-<
+    LinkedType(){
-<
+      next = NULL;
-<
+      nameA[0] = '\0';
-<
+      nameB[0] = '\0';
-<
+      nameC[0] = '\0';
-<
+      type[0] = '\0';
-<
+    }
-<
+    ~LinkedType(){ if( next != NULL ) delete next; }
-<
-<
+    LinkedType* find( char* key1, char* key2, char* key3, char* key4 ){
-<
-<
-<
-<
-<
+      if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 ) &&
-<
+          !strcmp( nameC, key3 ) && !strcmp( nameD, key4 ) ) return this;
-<
-<
+      if( !strcmp( nameA, key4 ) && !strcmp( nameB, key3 ) &&
-<
+          !strcmp( nameC, key2 ) && !strcmp( nameD, key1 ) ) return this;
-<
-<
+      if( next != NULL ) return next->find(key1, key2, key3, key4);
-<
+      return NULL;
-<
+    }
-<
-<
+    void add( torsionStruct &info ){
->
+void TraPPE_ExFF::initializeTorsions( int nTorsions, Torsion** torsionArray,
->
+                                      torsion_set* the_torsions ){
-<
+      // check for duplicates
-<
+      int dup = 0;
-<
-<
+      if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB ) &&
-<
+          !strcmp( nameC, info.nameC ) && !strcmp( nameD, info.nameD ) ) dup = 1;
-<
-<
+      if( !strcmp( nameA, info.nameD ) && !strcmp( nameB, info.nameC ) &&
-<
+          !strcmp( nameC, info.nameB ) && !strcmp( nameD, info.nameA ) ) dup = 1;
-<
-<
+      if(dup){
-<
+        sprintf( painCave.errMsg,
-<
+                 "Duplicate TraPPE_Ex torsion type \"%s - %s - %s - %s\" found in "
-<
+                 "the TraPPE_ExFF param file./n", nameA, nameB, nameC, nameD );
-<
+        painCave.isFatal = 1;
-<
+        simError();
-<
+      }
-<
-<
+      if( next != NULL ) next->add(info);
-<
+      else{
-<
+        next = new LinkedType();
-<
+        strcpy(next->nameA, info.nameA);
-<
+        strcpy(next->nameB, info.nameB);
-<
+        strcpy(next->nameC, info.nameC);
-<
+        strcpy(next->nameD, info.nameD);
-<
+        strcpy(next->type,  info.type);
-<
+        next->k1 = info.k1;
-<
+        next->k2 = info.k2;
-<
+        next->k3 = info.k3;
-<
+        next->k4 = info.k4;
-<
-<
+      }
-<
+    }
-<
-<
+#ifdef IS_MPI
-<
-<
+    void duplicate( torsionStruct &info ){
-<
+      strcpy(info.nameA, nameA);
-<
+      strcpy(info.nameB, nameB);
-<
+      strcpy(info.nameC, nameC);
-<
+      strcpy(info.nameD, nameD);
-<
+      strcpy(info.type,  type);
-<
+      info.k1   = k1;
-<
+      info.k2   = k2;
-<
+      info.k3   = k3;
-<
+      info.k4   = k4;
-<
+      info.last = 0;
-<
+    }
-<
-<
+#endif
-<
-<
+    char nameA[15];
-<
+    char nameB[15];
-<
+    char nameC[15];
-<
+    char nameD[15];
-<
+    char type[30];
-<
+    double k1, k2, k3, k4;
-<
-<
+    LinkedType* next;
-<
+  };
-<
-<
+  LinkedType* headTorsionType;
-<
+  LinkedType* currentTorsionType;
-<
+  torsionStruct info;
-<
+  info.last = 1; // initialize last to have the last set.
-<
+                 // if things go well, last will be set to 0
-<
-<
+  int i, a, b, c, d, index;
->
+  int i, a, b, c, d;
+  char* atomA;
+  char* atomB;
+  char* atomC;
+  char* atomD;
-–
+  CubicTorsion* cTors;
-<
+  SRI **the_sris;
->
+  CubicTorsion* cTors;
+  Atom** the_atoms;
-–
+  int nTorsions;
-–
+  the_sris = entry_plug->sr_interactions;
+  the_atoms = entry_plug->atoms;
-–
+  nTorsions = entry_plug->n_torsions;
-–
+#ifdef IS_MPI
-–
+  if( worldRank == 0 ){
-–
+#endif
-–
-–
+    // read in the torsion types.
-–
-–
+    headTorsionType = new LinkedType;
-–
-–
+    fastForward( "TorsionTypes", "initializeTorsions" );
-–
-–
+    // we are now at the torsionTypes section
-–
-–
+    eof_test =  fgets( readLine, sizeof(readLine), frcFile );
-–
+    lineNum++;
-–
-–
-–
+    // read a line, and start parseing out the atom types
-–
-–
+    if( eof_test == NULL ){
-–
+      sprintf( painCave.errMsg,
-–
+               "Error in reading torsions from force file at line %d.\n",
-–
+               lineNum );
-–
+      painCave.isFatal = 1;
-–
+      simError();
-–
+    }
-–
-–
+    // stop reading at end of file, or at next section
-–
+    while( readLine[0] != '#' && eof_test != NULL ){
-–
-–
+      // toss comment lines
-–
+      if( readLine[0] != '!' ){
-–
-–
+        // the parser returns 0 if the line was blank
-–
+        if( parseTorsion( readLine, lineNum, info ) ){
-–
+          headTorsionType->add( info );
-–
-–
+        }
-–
+      }
-–
+      eof_test = fgets( readLine, sizeof(readLine), frcFile );
-–
+      lineNum++;
-–
+    }
-–
-–
+#ifdef IS_MPI
-–
-–
+    // send out the linked list to all the other processes
-–
-–
+    sprintf( checkPointMsg,
-–
+             "TraPPE_Ex torsion structures read successfully." );
-–
+    MPIcheckPoint();
-–
-–
+    currentTorsionType = headTorsionType;
-–
+    while( currentTorsionType != NULL ){
-–
+      currentTorsionType->duplicate( info );
-–
+      sendFrcStruct( &info, mpiTorsionStructType );
-–
+      currentTorsionType = currentTorsionType->next;
-–
+    }
-–
+    info.last = 1;
-–
+    sendFrcStruct( &info, mpiTorsionStructType );
-–
-–
+  }
-–
-–
+  else{
-–
-–
+    // listen for node 0 to send out the force params
-–
-–
+    MPIcheckPoint();
-–
-–
+    headTorsionType = new LinkedType;
-–
+    recieveFrcStruct( &info, mpiTorsionStructType );
-–
+    while( !info.last ){
-–
-–
+      headTorsionType->add( info );
-–
+      recieveFrcStruct( &info, mpiTorsionStructType );
-–
+    }
-–
+  }
-–
+#endif // is_mpi
-–
+  // initialize the Torsions
+  for( i=0; i<nTorsions; i++ ){
+    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;
->
+      torsionArray[i] = cTors;
-–
-–
-–
+  // clean up the memory
-–
-–
+  delete headTorsionType;
-–
-–
+#ifdef IS_MPI
-–
+  sprintf( checkPointMsg, "TraPPE_Ex torsions initialized succesfully" );
-–
+  MPIcheckPoint();
-–
+#endif // is_mpi
-–
+void TraPPE_ExFF::fastForward( char* stopText, char* searchOwner ){

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines

Comparing trunk/OOPSE/libmdtools/TraPPE_ExFF.cpp (file contents): Revision 388 by chuckv, Fri Mar 21 22:11:50 2003 UTC vs. Revision 542 by mmeineke, Fri May 30 21:31:48 2003 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/TraPPE_ExFF.cpp (file contents):
Revision 388 by chuckv, Fri Mar 21 22:11:50 2003 UTC vs.
Revision 542 by mmeineke, Fri May 30 21:31:48 2003 UTC