ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/TraPPE_ExFF.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/TraPPE_ExFF.cpp (file contents):
Revision 378 by mmeineke, Fri Mar 21 17:42:12 2003 UTC vs.
Revision 434 by chuckv, Fri Mar 28 19:30:59 2003 UTC

# Line 85 | Line 85 | namespace TPE {  // restrict the access of the folowin
85  
86   #endif
87  
88 +  class LinkedAtomType {
89 +  public:
90 +    LinkedAtomType(){
91 +      next = NULL;
92 +      name[0] = '\0';
93 +    }
94 +    ~LinkedAtomType(){ if( next != NULL ) delete next; }
95 +
96 +    LinkedAtomType* find(char* key){
97 +      if( !strcmp(name, key) ) return this;
98 +      if( next != NULL ) return next->find(key);
99 +      return NULL;
100 +    }
101 +    
102 +    void add( atomStruct &info ){
103 +
104 +      // check for duplicates
105 +      
106 +      if( !strcmp( info.name, name ) ){
107 +        sprintf( painCave.errMsg,
108 +                 "Duplicate TraPPE_Ex atom type \"%s\" found in "
109 +                 "the TraPPE_ExFF param file./n",
110 +                 name );
111 +        painCave.isFatal = 1;
112 +        simError();
113 +      }
114 +
115 +      if( next != NULL ) next->add(info);
116 +      else{
117 +        next = new LinkedAtomType();
118 +        strcpy(next->name, info.name);
119 +        next->isDipole = info.isDipole;
120 +        next->isSSD    = info.isSSD;
121 +        next->mass     = info.mass;
122 +        next->epslon   = info.epslon;
123 +        next->sigma    = info.sigma;
124 +        next->dipole   = info.dipole;
125 +        next->w0       = info.w0;
126 +        next->v0       = info.v0;
127 +        next->ident    = info.ident;
128 +      }
129 +    }
130 +
131 + #ifdef IS_MPI
132 +    
133 +    void duplicate( atomStruct &info ){
134 +      strcpy(info.name, name);
135 +      info.isDipole = isDipole;
136 +      info.isSSD    = isSSD;
137 +      info.mass     = mass;
138 +      info.epslon   = epslon;
139 +      info.sigma    = sigma;
140 +      info.dipole   = dipole;
141 +      info.w0       = w0;
142 +      info.v0       = v0;
143 +      info.last     = 0;
144 +    }
145 +
146 +
147 + #endif
148 +
149 +    char name[15];
150 +    int isDipole;
151 +    int isSSD;
152 +    double mass;
153 +    double epslon;
154 +    double sigma;
155 +    double dipole;
156 +    double w0;
157 +    double v0;
158 +    int ident;
159 +    LinkedAtomType* next;
160 +  };
161 +
162 +  class LinkedBondType {
163 +  public:
164 +    LinkedBondType(){
165 +      next = NULL;
166 +      nameA[0] = '\0';
167 +      nameB[0] = '\0';
168 +      type[0] = '\0';
169 +    }
170 +    ~LinkedBondType(){ if( next != NULL ) delete next; }
171 +
172 +    LinkedBondType* find(char* key1, char* key2){
173 +      if( !strcmp(nameA, key1 ) && !strcmp( nameB, key2 ) ) return this;
174 +      if( !strcmp(nameA, key2 ) && !strcmp( nameB, key1 ) ) return this;
175 +      if( next != NULL ) return next->find(key1, key2);
176 +      return NULL;
177 +    }
178 +    
179 +
180 +    void add( bondStruct &info ){
181 +      
182 +      // check for duplicates
183 +      int dup = 0;
184 +
185 +      if( !strcmp(nameA, info.nameA ) && !strcmp( nameB, info.nameB ) ) dup = 1;
186 +      if( !strcmp(nameA, info.nameB ) && !strcmp( nameB, info.nameA ) ) dup = 1;
187 +      
188 +      if(dup){
189 +        sprintf( painCave.errMsg,
190 +                 "Duplicate TraPPE_Ex bond type \"%s - %s\" found in "
191 +                 "the TraPPE_ExFF param file./n",
192 +                 nameA, nameB );
193 +        painCave.isFatal = 1;
194 +        simError();
195 +      }
196 +
197 +        
198 +      if( next != NULL ) next->add(info);
199 +      else{
200 +        next = new LinkedBondType();
201 +        strcpy(next->nameA, info.nameA);
202 +        strcpy(next->nameB, info.nameB);
203 +        strcpy(next->type,  info.type);
204 +        next->d0 = info.d0;
205 +      }
206 +    }
207 +    
208 + #ifdef IS_MPI
209 +    void duplicate( bondStruct &info ){
210 +      strcpy(info.nameA, nameA);
211 +      strcpy(info.nameB, nameB);
212 +      strcpy(info.type,  type);
213 +      info.d0   = d0;
214 +      info.last = 0;
215 +    }
216 +
217 +
218 + #endif
219 +
220 +    char nameA[15];
221 +    char nameB[15];
222 +    char type[30];
223 +    double d0;
224 +
225 +    LinkedBondType* next;
226 +  };
227 +
228 +  class LinkedBendType {
229 +  public:
230 +    LinkedBendType(){
231 +      next = NULL;
232 +      nameA[0] = '\0';
233 +      nameB[0] = '\0';
234 +      nameC[0] = '\0';
235 +      type[0] = '\0';
236 +    }
237 +    ~LinkedBendType(){ if( next != NULL ) delete next; }
238 +
239 +    LinkedBendType* find( char* key1, char* key2, char* key3 ){
240 +      if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 )
241 +          && !strcmp( nameC, key3 ) ) return this;
242 +      if( !strcmp( nameA, key3 ) && !strcmp( nameB, key2 )
243 +          && !strcmp( nameC, key1 ) ) return this;
244 +      if( next != NULL ) return next->find(key1, key2, key3);
245 +      return NULL;
246 +    }
247 +    
248 +    void add( bendStruct &info ){
249 +
250 +      // check for duplicates
251 +      int dup = 0;
252 +      
253 +      if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB )
254 +          && !strcmp( nameC, info.nameC ) ) dup = 1;
255 +      if( !strcmp( nameA, info.nameC ) && !strcmp( nameB, info.nameB )
256 +          && !strcmp( nameC, info.nameA ) ) dup = 1;
257 +
258 +      if(dup){
259 +        sprintf( painCave.errMsg,
260 +                 "Duplicate TraPPE_Ex bend type \"%s - %s - %s\" found in "
261 +                 "the TraPPE_ExFF param file./n",
262 +                 nameA, nameB, nameC );
263 +        painCave.isFatal = 1;
264 +        simError();
265 +      }
266 +
267 +      if( next != NULL ) next->add(info);
268 +      else{
269 +        next = new LinkedBendType();
270 +        strcpy(next->nameA, info.nameA);
271 +        strcpy(next->nameB, info.nameB);
272 +        strcpy(next->nameC, info.nameC);
273 +        strcpy(next->type,  info.type);
274 +        next->k1 = info.k1;
275 +        next->k2 = info.k2;
276 +        next->k3 = info.k3;
277 +        next->t0 = info.t0;
278 +      }
279 +    }
280 +
281 + #ifdef IS_MPI    
282 +
283 +    void duplicate( bendStruct &info ){
284 +      strcpy(info.nameA, nameA);
285 +      strcpy(info.nameB, nameB);
286 +      strcpy(info.nameC, nameC);
287 +      strcpy(info.type,  type);
288 +      info.k1   = k1;
289 +      info.k2   = k2;
290 +      info.k3   = k3;
291 +      info.t0   = t0;
292 +      info.last = 0;
293 +    }
294 +
295 + #endif // is_mpi
296 +
297 +    char nameA[15];
298 +    char nameB[15];
299 +    char nameC[15];
300 +    char type[30];
301 +    double k1, k2, k3, t0;
302 +
303 +    LinkedBendType* next;
304 +  };
305 +
306 +  class LinkedTorsionType {
307 +  public:
308 +    LinkedTorsionType(){
309 +      next = NULL;
310 +      nameA[0] = '\0';
311 +      nameB[0] = '\0';
312 +      nameC[0] = '\0';
313 +      type[0] = '\0';
314 +    }
315 +    ~LinkedTorsionType(){ if( next != NULL ) delete next; }
316 +
317 +    LinkedTorsionType* find( char* key1, char* key2, char* key3, char* key4 ){
318 +      
319 +
320 +
321 +
322 +      if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 ) &&
323 +          !strcmp( nameC, key3 ) && !strcmp( nameD, key4 ) ) return this;
324 +
325 +      if( !strcmp( nameA, key4 ) && !strcmp( nameB, key3 ) &&
326 +          !strcmp( nameC, key2 ) && !strcmp( nameD, key1 ) ) return this;
327 +
328 +      if( next != NULL ) return next->find(key1, key2, key3, key4);
329 +      return NULL;
330 +    }
331 +
332 +    void add( torsionStruct &info ){
333 +
334 +      // check for duplicates
335 +      int dup = 0;
336 +
337 +      if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB ) &&
338 +          !strcmp( nameC, info.nameC ) && !strcmp( nameD, info.nameD ) ) dup = 1;
339 +      
340 +      if( !strcmp( nameA, info.nameD ) && !strcmp( nameB, info.nameC ) &&
341 +          !strcmp( nameC, info.nameB ) && !strcmp( nameD, info.nameA ) ) dup = 1;
342 +      
343 +      if(dup){
344 +        sprintf( painCave.errMsg,
345 +                 "Duplicate TraPPE_Ex torsion type \"%s - %s - %s - %s\" found in "
346 +                 "the TraPPE_ExFF param file./n", nameA, nameB, nameC, nameD );
347 +        painCave.isFatal = 1;
348 +        simError();
349 +      }
350 +
351 +      if( next != NULL ) next->add(info);
352 +      else{
353 +        next = new LinkedTorsionType();
354 +        strcpy(next->nameA, info.nameA);
355 +        strcpy(next->nameB, info.nameB);
356 +        strcpy(next->nameC, info.nameC);
357 +        strcpy(next->nameD, info.nameD);
358 +        strcpy(next->type,  info.type);
359 +        next->k1 = info.k1;
360 +        next->k2 = info.k2;
361 +        next->k3 = info.k3;
362 +        next->k4 = info.k4;
363 +
364 +      }
365 +    }
366 +
367 + #ifdef IS_MPI
368 +    
369 +    void duplicate( torsionStruct &info ){
370 +      strcpy(info.nameA, nameA);
371 +      strcpy(info.nameB, nameB);
372 +      strcpy(info.nameC, nameC);
373 +      strcpy(info.nameD, nameD);
374 +      strcpy(info.type,  type);
375 +      info.k1   = k1;
376 +      info.k2   = k2;
377 +      info.k3   = k3;
378 +      info.k4   = k4;
379 +      info.last = 0;
380 +    }
381 +
382 + #endif
383 +
384 +    char nameA[15];
385 +    char nameB[15];
386 +    char nameC[15];
387 +    char nameD[15];
388 +    char type[30];
389 +    double k1, k2, k3, k4;
390 +
391 +    LinkedTorsionType* next;
392 +  };
393 +
394 +
395 +  LinkedAtomType* headAtomType;
396 +  LinkedAtomType* currentAtomType;
397 +  LinkedBondType* headBondType;
398 +  LinkedBondType* currentBondType;
399 +  LinkedBendType* headBendType;
400 +  LinkedBendType* currentBendType;
401 +  LinkedTorsionType* headTorsionType;
402 +  LinkedTorsionType* currentTorsionType;
403 +
404   } // namespace
405  
406   using namespace TPE;
# Line 103 | Line 419 | TraPPE_ExFF::TraPPE_ExFF(){
419    char temp[200];
420    char errMsg[1000];
421  
422 +  headAtomType       = NULL;
423 +  currentAtomType    = NULL;
424 +  headBondType       = NULL;
425 +  currentBondType    = NULL;
426 +  headBendType       = NULL;
427 +  currentBendType    = NULL;
428 +  headTorsionType    = NULL;
429 +  currentTorsionType = NULL;
430 +
431    // do the funtion wrapping
432    wrapMeFF( this );
433  
# Line 253 | Line 578 | TraPPE_ExFF::~TraPPE_ExFF(){
578  
579   TraPPE_ExFF::~TraPPE_ExFF(){
580  
581 +  if( headAtomType != NULL ) delete headAtomType;
582 +  if( headBondType != NULL ) delete headBondType;
583 +  if( headBendType != NULL ) delete headBendType;
584 +  if( headTorsionType != NULL ) delete headTorsionType;
585 +
586   #ifdef IS_MPI
587    if( worldRank == 0 ){
588   #endif // is_mpi
# Line 264 | Line 594 | void TraPPE_ExFF::initForceField( int ljMixRule ){
594   #endif // is_mpi
595   }
596  
597 + void TraPPE_ExFF::cleanMe( void ){
598 +
599 + #ifdef IS_MPI
600 +  
601 +  // keep the linked lists in the mpi version
602 +
603 + #else // is_mpi
604 +
605 +  // delete the linked lists in the single processor version
606 +
607 +  if( headAtomType != NULL ) delete headAtomType;
608 +  if( headBondType != NULL ) delete headBondType;
609 +  if( headBendType != NULL ) delete headBendType;
610 +  if( headTorsionType != NULL ) delete headTorsionType;
611 +
612 + #endif // is_mpi
613 + }
614 +
615 +
616   void TraPPE_ExFF::initForceField( int ljMixRule ){
617    
618    initFortran( ljMixRule, entry_plug->useReactionField );
619   }
620  
621  
622 < void TraPPE_ExFF::initializeAtoms( void ){
623 <  
624 <  class LinkedType {
276 <  public:
277 <    LinkedType(){
278 <      next = NULL;
279 <      name[0] = '\0';
280 <    }
281 <    ~LinkedType(){ if( next != NULL ) delete next; }
282 <
283 <    LinkedType* find(char* key){
284 <      if( !strcmp(name, key) ) return this;
285 <      if( next != NULL ) return next->find(key);
286 <      return NULL;
287 <    }
288 <    
289 <    void add( atomStruct &info ){
290 <
291 <      // check for duplicates
292 <      
293 <      if( !strcmp( info.name, name ) ){
294 <        sprintf( painCave.errMsg,
295 <                 "Duplicate TraPPE_Ex atom type \"%s\" found in "
296 <                 "the TraPPE_ExFF param file./n",
297 <                 name );
298 <        painCave.isFatal = 1;
299 <        simError();
300 <      }
301 <
302 <      if( next != NULL ) next->add(info);
303 <      else{
304 <        next = new LinkedType();
305 <        strcpy(next->name, info.name);
306 <        next->isDipole = info.isDipole;
307 <        next->isSSD    = info.isSSD;
308 <        next->mass     = info.mass;
309 <        next->epslon   = info.epslon;
310 <        next->sigma    = info.sigma;
311 <        next->dipole   = info.dipole;
312 <        next->w0       = info.w0;
313 <        next->v0       = info.v0;
314 <        next->ident    = info.ident;
315 <      }
316 <    }
317 <
318 < #ifdef IS_MPI
319 <    
320 <    void duplicate( atomStruct &info ){
321 <      strcpy(info.name, name);
322 <      info.isDipole = isDipole;
323 <      info.isSSD    = isSSD;
324 <      info.mass     = mass;
325 <      info.epslon   = epslon;
326 <      info.sigma    = sigma;
327 <      info.dipole   = dipole;
328 <      info.w0       = w0;
329 <      info.v0       = v0;
330 <      info.last     = 0;
331 <    }
332 <
333 <
334 < #endif
335 <
336 <    char name[15];
337 <    int isDipole;
338 <    int isSSD;
339 <    double mass;
340 <    double epslon;
341 <    double sigma;
342 <    double dipole;
343 <    double w0;
344 <    double v0;
345 <    int ident;
346 <    LinkedType* next;
347 <  };
348 <  
349 <  LinkedType* headAtomType;
350 <  LinkedType* currentAtomType;
351 <  atomStruct info;
352 <  info.last = 1; // initialize last to have the last set.
353 <                 // if things go well, last will be set to 0
354 <
355 <  
356 <
357 <  int i;
622 > void TraPPE_ExFF::readParams( void ){
623 >
624 >  int i, a, b, c, d;
625    int identNum;
626 +  char* atomA;
627 +  char* atomB;
628 +  char* atomC;
629 +  char* atomD;
630    
631 <  Atom** the_atoms;
632 <  int nAtoms;
633 <  the_atoms = entry_plug->atoms;
634 <  nAtoms = entry_plug->n_atoms;
631 >  atomStruct atomInfo;
632 >  bondStruct bondInfo;
633 >  bendStruct bendInfo;
634 >  torsionStruct torsionInfo;
635    
636 <  //////////////////////////////////////////////////
366 <  // a quick water fix
636 >  bigSigma = 0.0;
637  
638 <  double waterI[3][3];
639 <  waterI[0][0] = 1.76958347772500;
640 <  waterI[0][1] = 0.0;
641 <  waterI[0][2] = 0.0;
638 >  atomInfo.last = 1;
639 >  bondInfo.last = 1;
640 >  bendInfo.last = 1;
641 >  torsionInfo.last = 1;
642  
643 <  waterI[1][0] = 0.0;
374 <  waterI[1][1] = 0.614537057924513;
375 <  waterI[1][2] = 0.0;
376 <
377 <  waterI[2][0] = 0.0;
378 <  waterI[2][1] = 0.0;
379 <  waterI[2][2] = 1.15504641980049;
380 <
381 <
382 <  double headI[3][3];
383 <  headI[0][0] = 1125;
384 <  headI[0][1] = 0.0;
385 <  headI[0][2] = 0.0;
386 <
387 <  headI[1][0] = 0.0;
388 <  headI[1][1] = 1125;
389 <  headI[1][2] = 0.0;
390 <
391 <  headI[2][0] = 0.0;
392 <  headI[2][1] = 0.0;
393 <  headI[2][2] = 250;
394 <
643 >  // read in the atom info
644    
396
397  //////////////////////////////////////////////////
398
399
645   #ifdef IS_MPI
646    if( worldRank == 0 ){
647   #endif
648      
649      // read in the atom types.
405
406    headAtomType = new LinkedType;
650      
651 +    headAtomType = new LinkedAtomType;
652 +    
653      fastForward( "AtomTypes", "initializeAtoms" );
654  
655      // we are now at the AtomTypes section.
# Line 431 | Line 676 | void TraPPE_ExFF::initializeAtoms( void ){
676        if( readLine[0] != '!' ){
677          
678          // the parser returns 0 if the line was blank
679 <        if( parseAtom( readLine, lineNum, info ) ){
680 <          info.ident = identNum;
681 <          headAtomType->add( info );;
679 >        if( parseAtom( readLine, lineNum, atomInfo ) ){
680 >          atomInfo.ident = identNum;
681 >          headAtomType->add( atomInfo );;
682            identNum++;
683          }
684        }
# Line 451 | Line 696 | void TraPPE_ExFF::initializeAtoms( void ){
696  
697      currentAtomType = headAtomType->next; //skip the first element who is a place holder.
698      while( currentAtomType != NULL ){
699 <      currentAtomType->duplicate( info );
699 >      currentAtomType->duplicate( atomInfo );
700  
701  
702  
703 <      sendFrcStruct( &info, mpiAtomStructType );
703 >      sendFrcStruct( &atomInfo, mpiAtomStructType );
704  
705        sprintf( checkPointMsg,
706                 "successfully sent TraPPE_Ex force type: \"%s\"\n",
707 <               info.name );
707 >               atomInfo.name );
708        MPIcheckPoint();
709  
710        currentAtomType = currentAtomType->next;
711      }
712 <    info.last = 1;
713 <    sendFrcStruct( &info, mpiAtomStructType );
712 >    atomInfo.last = 1;
713 >    sendFrcStruct( &atomInfo, mpiAtomStructType );
714      
715    }
716  
# Line 475 | Line 720 | void TraPPE_ExFF::initializeAtoms( void ){
720      
721      MPIcheckPoint();
722  
723 <    headAtomType = new LinkedType;
724 <    recieveFrcStruct( &info, mpiAtomStructType );
723 >    headAtomType = new LinkedAtomType;
724 >    recieveFrcStruct( &atomInfo, mpiAtomStructType );
725      
726 <    while( !info.last ){
726 >    while( !atomInfo.last ){
727  
728  
729  
730 <      headAtomType->add( info );
730 >      headAtomType->add( atomInfo );
731        
732        MPIcheckPoint();
733  
734 <      recieveFrcStruct( &info, mpiAtomStructType );
734 >      recieveFrcStruct( &atomInfo, mpiAtomStructType );
735      }
736    }
737   #endif // is_mpi
# Line 505 | Line 750 | void TraPPE_ExFF::initializeAtoms( void ){
750    currentAtomType = headAtomType;
751    while( currentAtomType != NULL ){
752      
508    if(currentAtomType->isDipole) entry_plug->useReactionField = 1;
753      if(currentAtomType->isDipole) entry_plug->useDipole = 1;
754      if(currentAtomType->isSSD)    entry_plug->useSticky = 1;
755  
# Line 546 | Line 790 | void TraPPE_ExFF::initializeAtoms( void ){
790   #endif // is_mpi
791  
792    
549  // initialize the atoms
550  
551  double bigSigma = 0.0;
552  DirectionalAtom* dAtom;
793  
794 <  for( i=0; i<nAtoms; i++ ){
555 <    
556 <    currentAtomType = headAtomType->find( the_atoms[i]->getType() );
557 <    if( currentAtomType == NULL ){
558 <      sprintf( painCave.errMsg,
559 <               "AtomType error, %s not found in force file.\n",
560 <               the_atoms[i]->getType() );
561 <      painCave.isFatal = 1;
562 <      simError();
563 <    }
564 <    
565 <    the_atoms[i]->setMass( currentAtomType->mass );
566 <    the_atoms[i]->setEpslon( currentAtomType->epslon );
567 <    the_atoms[i]->setSigma( currentAtomType->sigma );
568 <    the_atoms[i]->setIdent( currentAtomType->ident );
569 <    the_atoms[i]->setLJ();
570 <
571 <    if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
572 <
573 <    if( currentAtomType->isDipole ){
574 <      if( the_atoms[i]->isDirectional() ){
575 <        
576 <        dAtom = (DirectionalAtom *) the_atoms[i];
577 <        dAtom->setMu( currentAtomType->dipole );
578 <        dAtom->setHasDipole( 1 );
579 <        dAtom->setJx( 0.0 );
580 <        dAtom->setJy( 0.0 );
581 <        dAtom->setJz( 0.0 );
582 <        
583 <        if(!strcmp("SSD",the_atoms[i]->getType())){
584 <          dAtom->setI( waterI );
585 <          dAtom->setSSD( 1 );
586 <        }
587 <        else if(!strcmp("HEAD",the_atoms[i]->getType())){
588 <          dAtom->setI( headI );
589 <          dAtom->setSSD( 0 );
590 <        }
591 <        else{
592 <          sprintf(painCave.errMsg,
593 <                  "AtmType error, %s does not have a moment of inertia set.\n",
594 <                  the_atoms[i]->getType() );
595 <          painCave.isFatal = 1;
596 <          simError();
597 <        }
598 <        entry_plug->n_dipoles++;
599 <      }
600 <      else{
601 <        
602 <        sprintf( painCave.errMsg,
603 <                "TraPPE_ExFF error: Atom \"%s\" is a dipole, yet no standard"
604 <                 " orientation was specifed in the BASS file.\n",
605 <                 currentAtomType->name );
606 <        painCave.isFatal = 1;
607 <        simError();
608 <      }
609 <    }
610 <    else{
611 <      if( the_atoms[i]->isDirectional() ){
612 <        sprintf( painCave.errMsg,
613 <                 "TraPPE_ExFF error: Atom \"%s\" was given a standard"
614 <                 "orientation in the BASS file, yet it is not a dipole.\n",
615 <                 currentAtomType->name);
616 <        painCave.isFatal = 1;
617 <        simError();
618 <      }
619 <    }
620 <  }
621 <
622 < #ifdef IS_MPI
623 <  double tempBig = bigSigma;
624 <  MPI::COMM_WORLD.Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX );
625 < #endif  //is_mpi
626 <
627 <  //calc rCut and rList
628 <
629 <  entry_plug->rCut = 2.5 * bigSigma;
630 <  
631 <  if(entry_plug->rCut > (entry_plug->box_x / 2.0))
632 <    entry_plug->rCut = entry_plug->box_x / 2.0;
633 <  
634 <  if(entry_plug->rCut > (entry_plug->box_y / 2.0))
635 <    entry_plug->rCut = entry_plug->box_y / 2.0;
636 <  
637 <  if(entry_plug->rCut > (entry_plug->box_z / 2.0))
638 <    entry_plug->rCut = entry_plug->box_z / 2.0;
639 <
640 <  entry_plug->rList = entry_plug->rCut + 1.0;
641 <
642 <  entry_plug->useLJ = 1; // use Lennard Jones is on by default
643 <
644 <  // clean up the memory
645 <  
646 <  delete headAtomType;
647 <
648 < #ifdef IS_MPI
649 <  sprintf( checkPointMsg, "TraPPE_Ex atoms initialized succesfully" );
650 <  MPIcheckPoint();
651 < #endif // is_mpi
652 <
653 < }
654 <
655 < void TraPPE_ExFF::initializeBonds( bond_pair* the_bonds ){
656 <  
657 <  class LinkedType {
658 <  public:
659 <    LinkedType(){
660 <      next = NULL;
661 <      nameA[0] = '\0';
662 <      nameB[0] = '\0';
663 <      type[0] = '\0';
664 <    }
665 <    ~LinkedType(){ if( next != NULL ) delete next; }
666 <
667 <    LinkedType* find(char* key1, char* key2){
668 <      if( !strcmp(nameA, key1 ) && !strcmp( nameB, key2 ) ) return this;
669 <      if( !strcmp(nameA, key2 ) && !strcmp( nameB, key1 ) ) return this;
670 <      if( next != NULL ) return next->find(key1, key2);
671 <      return NULL;
672 <    }
673 <    
674 <
675 <    void add( bondStruct &info ){
676 <      
677 <      // check for duplicates
678 <      int dup = 0;
679 <
680 <      if( !strcmp(nameA, info.nameA ) && !strcmp( nameB, info.nameB ) ) dup = 1;
681 <      if( !strcmp(nameA, info.nameB ) && !strcmp( nameB, info.nameA ) ) dup = 1;
682 <      
683 <      if(dup){
684 <        sprintf( painCave.errMsg,
685 <                 "Duplicate TraPPE_Ex bond type \"%s - %s\" found in "
686 <                 "the TraPPE_ExFF param file./n",
687 <                 nameA, nameB );
688 <        painCave.isFatal = 1;
689 <        simError();
690 <      }
691 <
692 <        
693 <      if( next != NULL ) next->add(info);
694 <      else{
695 <        next = new LinkedType();
696 <        strcpy(next->nameA, info.nameA);
697 <        strcpy(next->nameB, info.nameB);
698 <        strcpy(next->type,  info.type);
699 <        next->d0 = info.d0;
700 <      }
701 <    }
702 <    
703 < #ifdef IS_MPI
704 <    void duplicate( bondStruct &info ){
705 <      strcpy(info.nameA, nameA);
706 <      strcpy(info.nameB, nameB);
707 <      strcpy(info.type,  type);
708 <      info.d0   = d0;
709 <      info.last = 0;
710 <    }
711 <
712 <
713 < #endif
714 <
715 <    char nameA[15];
716 <    char nameB[15];
717 <    char type[30];
718 <    double d0;
719 <
720 <    LinkedType* next;
721 <  };
722 <
723 <
724 <  
725 <  LinkedType* headBondType;
726 <  LinkedType* currentBondType;
727 <  bondStruct info;
728 <  info.last = 1; // initialize last to have the last set.
729 <                 // if things go well, last will be set to 0
730 <
731 <  SRI **the_sris;
732 <  Atom** the_atoms;
733 <  int nBonds;
734 <  the_sris = entry_plug->sr_interactions;
735 <  the_atoms = entry_plug->atoms;
736 <  nBonds = entry_plug->n_bonds;
737 <
738 <  int i, a, b;
739 <  char* atomA;
740 <  char* atomB;
794 >  // read in the bonds
795    
796   #ifdef IS_MPI
797    if( worldRank == 0 ){
# Line 745 | Line 799 | void TraPPE_ExFF::initializeBonds( bond_pair* the_bond
799      
800      // read in the bond types.
801      
802 <    headBondType = new LinkedType;
802 >    headBondType = new LinkedBondType;
803      
804      fastForward( "BondTypes", "initializeBonds" );
805  
# Line 772 | Line 826 | void TraPPE_ExFF::initializeBonds( bond_pair* the_bond
826        if( readLine[0] != '!' ){
827          
828          // the parser returns 0 if the line was blank
829 <        if( parseBond( readLine, lineNum, info ) ){
830 <          headBondType->add( info );
829 >        if( parseBond( readLine, lineNum, bondInfo ) ){
830 >          headBondType->add( bondInfo );
831          }
832        }
833        eof_test = fgets( readLine, sizeof(readLine), frcFile );
# Line 790 | Line 844 | void TraPPE_ExFF::initializeBonds( bond_pair* the_bond
844      
845      currentBondType = headBondType;
846      while( currentBondType != NULL ){
847 <      currentBondType->duplicate( info );
848 <      sendFrcStruct( &info, mpiBondStructType );
847 >      currentBondType->duplicate( bondInfo );
848 >      sendFrcStruct( &bondInfo, mpiBondStructType );
849        currentBondType = currentBondType->next;
850      }
851 <    info.last = 1;
852 <    sendFrcStruct( &info, mpiBondStructType );
851 >    bondInfo.last = 1;
852 >    sendFrcStruct( &bondInfo, mpiBondStructType );
853      
854    }
855  
# Line 805 | Line 859 | void TraPPE_ExFF::initializeBonds( bond_pair* the_bond
859      
860      MPIcheckPoint();
861  
862 <    headBondType = new LinkedType;
863 <    recieveFrcStruct( &info, mpiBondStructType );
864 <    while( !info.last ){
862 >    headBondType = new LinkedBondType;
863 >    recieveFrcStruct( &bondInfo, mpiBondStructType );
864 >    while( !bondInfo.last ){
865  
866 <      headBondType->add( info );
867 <      recieveFrcStruct( &info, mpiBondStructType );
866 >      headBondType->add( bondInfo );
867 >      recieveFrcStruct( &bondInfo, mpiBondStructType );
868      }
869    }
870   #endif // is_mpi
871    
818  
819  // initialize the Bonds
820  
872  
873 <  for( i=0; i<nBonds; i++ ){
823 <    
824 <    a = the_bonds[i].a;
825 <    b = the_bonds[i].b;
826 <
827 <    atomA = the_atoms[a]->getType();
828 <    atomB = the_atoms[b]->getType();
829 <    currentBondType = headBondType->find( atomA, atomB );
830 <    if( currentBondType == NULL ){
831 <      sprintf( painCave.errMsg,
832 <               "BondType error, %s - %s not found in force file.\n",
833 <               atomA, atomB );
834 <      painCave.isFatal = 1;
835 <      simError();
836 <    }
837 <    
838 <    if( !strcmp( currentBondType->type, "fixed" ) ){
839 <      
840 <      the_sris[i] = new ConstrainedBond( *the_atoms[a],
841 <                                         *the_atoms[b],
842 <                                         currentBondType->d0 );
843 <      entry_plug->n_constraints++;
844 <    }
845 <  }
846 <
847 <
848 <  // clean up the memory
849 <  
850 <  delete headBondType;
851 <
852 < #ifdef IS_MPI
853 <  sprintf( checkPointMsg, "TraPPE_Ex bonds initialized succesfully" );
854 <  MPIcheckPoint();
855 < #endif // is_mpi
856 <
857 < }
858 <
859 < void TraPPE_ExFF::initializeBends( bend_set* the_bends ){
860 <  
861 <  class LinkedType {
862 <  public:
863 <    LinkedType(){
864 <      next = NULL;
865 <      nameA[0] = '\0';
866 <      nameB[0] = '\0';
867 <      nameC[0] = '\0';
868 <      type[0] = '\0';
869 <    }
870 <    ~LinkedType(){ if( next != NULL ) delete next; }
871 <
872 <    LinkedType* find( char* key1, char* key2, char* key3 ){
873 <      if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 )
874 <          && !strcmp( nameC, key3 ) ) return this;
875 <      if( !strcmp( nameA, key3 ) && !strcmp( nameB, key2 )
876 <          && !strcmp( nameC, key1 ) ) return this;
877 <      if( next != NULL ) return next->find(key1, key2, key3);
878 <      return NULL;
879 <    }
880 <    
881 <    void add( bendStruct &info ){
882 <
883 <      // check for duplicates
884 <      int dup = 0;
885 <      
886 <      if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB )
887 <          && !strcmp( nameC, info.nameC ) ) dup = 1;
888 <      if( !strcmp( nameA, info.nameC ) && !strcmp( nameB, info.nameB )
889 <          && !strcmp( nameC, info.nameA ) ) dup = 1;
890 <
891 <      if(dup){
892 <        sprintf( painCave.errMsg,
893 <                 "Duplicate TraPPE_Ex bend type \"%s - %s - %s\" found in "
894 <                 "the TraPPE_ExFF param file./n",
895 <                 nameA, nameB, nameC );
896 <        painCave.isFatal = 1;
897 <        simError();
898 <      }
899 <
900 <      if( next != NULL ) next->add(info);
901 <      else{
902 <        next = new LinkedType();
903 <        strcpy(next->nameA, info.nameA);
904 <        strcpy(next->nameB, info.nameB);
905 <        strcpy(next->nameC, info.nameC);
906 <        strcpy(next->type,  info.type);
907 <        next->k1 = info.k1;
908 <        next->k2 = info.k2;
909 <        next->k3 = info.k3;
910 <        next->t0 = info.t0;
911 <      }
912 <    }
913 <
914 < #ifdef IS_MPI    
915 <
916 <    void duplicate( bendStruct &info ){
917 <      strcpy(info.nameA, nameA);
918 <      strcpy(info.nameB, nameB);
919 <      strcpy(info.nameC, nameC);
920 <      strcpy(info.type,  type);
921 <      info.k1   = k1;
922 <      info.k2   = k2;
923 <      info.k3   = k3;
924 <      info.t0   = t0;
925 <      info.last = 0;
926 <    }
927 <
928 < #endif // is_mpi
929 <
930 <    char nameA[15];
931 <    char nameB[15];
932 <    char nameC[15];
933 <    char type[30];
934 <    double k1, k2, k3, t0;
935 <
936 <    LinkedType* next;
937 <  };
938 <  
939 <  LinkedType* headBendType;
940 <  LinkedType* currentBendType;
941 <  bendStruct info;
942 <  info.last = 1; // initialize last to have the last set.
943 <                 // if things go well, last will be set to 0
944 <
945 <  QuadraticBend* qBend;
946 <  GhostBend* gBend;
947 <  SRI **the_sris;
948 <  Atom** the_atoms;
949 <  int nBends;
950 <  the_sris = entry_plug->sr_interactions;
951 <  the_atoms = entry_plug->atoms;
952 <  nBends = entry_plug->n_bends;
873 >  // read in the bends
874  
954  int i, a, b, c;
955  char* atomA;
956  char* atomB;
957  char* atomC;
958
959
875   #ifdef IS_MPI
876    if( worldRank == 0 ){
877   #endif
878  
879      // read in the bend types.
880  
881 <    headBendType = new LinkedType;
881 >    headBendType = new LinkedBendType;
882      
883      fastForward( "BendTypes", "initializeBends" );
884  
# Line 989 | Line 904 | void TraPPE_ExFF::initializeBends( bend_set* the_bends
904        if( readLine[0] != '!' ){
905          
906          // the parser returns 0 if the line was blank
907 <        if( parseBend( readLine, lineNum, info ) ){
908 <          headBendType->add( info );
907 >        if( parseBend( readLine, lineNum, bendInfo ) ){
908 >          headBendType->add( bendInfo );
909          }
910        }
911        eof_test = fgets( readLine, sizeof(readLine), frcFile );
# Line 1007 | Line 922 | void TraPPE_ExFF::initializeBends( bend_set* the_bends
922  
923      currentBendType = headBendType;
924      while( currentBendType != NULL ){
925 <      currentBendType->duplicate( info );
926 <      sendFrcStruct( &info, mpiBendStructType );
925 >      currentBendType->duplicate( bendInfo );
926 >      sendFrcStruct( &bendInfo, mpiBendStructType );
927        currentBendType = currentBendType->next;
928      }
929 <    info.last = 1;
930 <    sendFrcStruct( &info, mpiBendStructType );
929 >    bendInfo.last = 1;
930 >    sendFrcStruct( &bendInfo, mpiBendStructType );
931      
932    }
933  
# Line 1022 | Line 937 | void TraPPE_ExFF::initializeBends( bend_set* the_bends
937      
938      MPIcheckPoint();
939  
940 <    headBendType = new LinkedType;
941 <    recieveFrcStruct( &info, mpiBendStructType );
942 <    while( !info.last ){
940 >    headBendType = new LinkedBendType;
941 >    recieveFrcStruct( &bendInfo, mpiBendStructType );
942 >    while( !bendInfo.last ){
943  
944 <      headBendType->add( info );
945 <      recieveFrcStruct( &info, mpiBendStructType );
944 >      headBendType->add( bendInfo );
945 >      recieveFrcStruct( &bendInfo, mpiBendStructType );
946      }
947    }
948   #endif // is_mpi
949 +
950 +
951 +  // read in the torsions
952 +
953 + #ifdef IS_MPI
954 +  if( worldRank == 0 ){
955 + #endif
956 +
957 +    // read in the torsion types.
958 +
959 +    headTorsionType = new LinkedTorsionType;
960 +    
961 +    fastForward( "TorsionTypes", "initializeTorsions" );
962 +
963 +    // we are now at the torsionTypes section
964 +
965 +    eof_test =  fgets( readLine, sizeof(readLine), frcFile );
966 +    lineNum++;
967 +    
968 +    
969 +    // read a line, and start parseing out the atom types
970 +
971 +    if( eof_test == NULL ){
972 +      sprintf( painCave.errMsg,
973 +               "Error in reading torsions from force file at line %d.\n",
974 +               lineNum );
975 +      painCave.isFatal = 1;
976 +      simError();
977 +    }
978 +    
979 +    // stop reading at end of file, or at next section
980 +    while( readLine[0] != '#' && eof_test != NULL ){
981 +
982 +      // toss comment lines
983 +      if( readLine[0] != '!' ){
984 +        
985 +        // the parser returns 0 if the line was blank
986 +        if( parseTorsion( readLine, lineNum, torsionInfo ) ){
987 +          headTorsionType->add( torsionInfo );
988 +
989 +        }
990 +      }
991 +      eof_test = fgets( readLine, sizeof(readLine), frcFile );
992 +      lineNum++;
993 +    }
994 +
995 + #ifdef IS_MPI
996 +    
997 +    // send out the linked list to all the other processes
998 +    
999 +    sprintf( checkPointMsg,
1000 +             "TraPPE_Ex torsion structures read successfully." );
1001 +    MPIcheckPoint();
1002 +    
1003 +    currentTorsionType = headTorsionType;
1004 +    while( currentTorsionType != NULL ){
1005 +      currentTorsionType->duplicate( torsionInfo );
1006 +      sendFrcStruct( &torsionInfo, mpiTorsionStructType );
1007 +      currentTorsionType = currentTorsionType->next;
1008 +    }
1009 +    torsionInfo.last = 1;
1010 +    sendFrcStruct( &torsionInfo, mpiTorsionStructType );
1011 +    
1012 +  }
1013 +
1014 +  else{
1015 +    
1016 +    // listen for node 0 to send out the force params
1017 +    
1018 +    MPIcheckPoint();
1019 +
1020 +    headTorsionType = new LinkedTorsionType;
1021 +    recieveFrcStruct( &torsionInfo, mpiTorsionStructType );
1022 +    while( !torsionInfo.last ){
1023 +
1024 +      headTorsionType->add( torsionInfo );
1025 +      recieveFrcStruct( &torsionInfo, mpiTorsionStructType );
1026 +    }
1027 +  }
1028 + #endif // is_mpi
1029 +
1030 +  entry_plug->useLJ = 1;
1031 + }
1032 +
1033 +
1034 +
1035 + void TraPPE_ExFF::initializeAtoms( int nAtoms, Atom** the_atoms ){
1036    
1035  // initialize the Bends
1037    
1038 <  int index;
1038 >  //////////////////////////////////////////////////
1039 >  // a quick water fix
1040  
1041 +  double waterI[3][3];
1042 +  waterI[0][0] = 1.76958347772500;
1043 +  waterI[0][1] = 0.0;
1044 +  waterI[0][2] = 0.0;
1045 +
1046 +  waterI[1][0] = 0.0;
1047 +  waterI[1][1] = 0.614537057924513;
1048 +  waterI[1][2] = 0.0;
1049 +
1050 +  waterI[2][0] = 0.0;
1051 +  waterI[2][1] = 0.0;
1052 +  waterI[2][2] = 1.15504641980049;
1053 +
1054 +
1055 +  double headI[3][3];
1056 +  headI[0][0] = 1125;
1057 +  headI[0][1] = 0.0;
1058 +  headI[0][2] = 0.0;
1059 +
1060 +  headI[1][0] = 0.0;
1061 +  headI[1][1] = 1125;
1062 +  headI[1][2] = 0.0;
1063 +
1064 +  headI[2][0] = 0.0;
1065 +  headI[2][1] = 0.0;
1066 +  headI[2][2] = 250;
1067 +
1068 +  //////////////////////////////////////////////////
1069 +
1070 +  
1071 +  // initialize the atoms
1072 +  
1073 +  DirectionalAtom* dAtom;
1074 +
1075 +  for(int i=0; i<nAtoms; i++ ){
1076 +
1077 +    currentAtomType = headAtomType->find( the_atoms[i]->getType() );
1078 +    if( currentAtomType == NULL ){
1079 +      sprintf( painCave.errMsg,
1080 +               "AtomType error, %s not found in force file.\n",
1081 +               the_atoms[i]->getType() );
1082 +      painCave.isFatal = 1;
1083 +      simError();
1084 +    }
1085 +    
1086 +    the_atoms[i]->setMass( currentAtomType->mass );
1087 +    the_atoms[i]->setEpslon( currentAtomType->epslon );
1088 +    the_atoms[i]->setSigma( currentAtomType->sigma );
1089 +    the_atoms[i]->setIdent( currentAtomType->ident );
1090 +    the_atoms[i]->setLJ();
1091 +
1092 +    if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
1093 +
1094 +    if( currentAtomType->isDipole ){
1095 +      if( the_atoms[i]->isDirectional() ){
1096 +        
1097 +        dAtom = (DirectionalAtom *) the_atoms[i];
1098 +        dAtom->setMu( currentAtomType->dipole );
1099 +        dAtom->setHasDipole( 1 );
1100 +        dAtom->setJx( 0.0 );
1101 +        dAtom->setJy( 0.0 );
1102 +        dAtom->setJz( 0.0 );
1103 +        
1104 +        if(!strcmp("SSD",the_atoms[i]->getType())){
1105 +          dAtom->setI( waterI );
1106 +          dAtom->setSSD( 1 );
1107 +        }
1108 +        else if(!strcmp("HEAD",the_atoms[i]->getType())){
1109 +          dAtom->setI( headI );
1110 +          dAtom->setSSD( 0 );
1111 +        }
1112 +        else{
1113 +          sprintf(painCave.errMsg,
1114 +                  "AtmType error, %s does not have a moment of inertia set.\n",
1115 +                  the_atoms[i]->getType() );
1116 +          painCave.isFatal = 1;
1117 +          simError();
1118 +        }
1119 +        entry_plug->n_dipoles++;
1120 +      }
1121 +      else{
1122 +        
1123 +        sprintf( painCave.errMsg,
1124 +                "TraPPE_ExFF error: Atom \"%s\" is a dipole, yet no standard"
1125 +                 " orientation was specifed in the BASS file.\n",
1126 +                 currentAtomType->name );
1127 +        painCave.isFatal = 1;
1128 +        simError();
1129 +      }
1130 +    }
1131 +    else{
1132 +      if( the_atoms[i]->isDirectional() ){
1133 +        sprintf( painCave.errMsg,
1134 +                 "TraPPE_ExFF error: Atom \"%s\" was given a standard"
1135 +                 "orientation in the BASS file, yet it is not a dipole.\n",
1136 +                 currentAtomType->name);
1137 +        painCave.isFatal = 1;
1138 +        simError();
1139 +      }
1140 +    }
1141 +  }
1142 + }
1143 +
1144 + void TraPPE_ExFF::initializeBonds( int nBonds, Bond** bondArray,
1145 +                                   bond_pair* the_bonds ){
1146 +  int i,a,b;
1147 +  char* atomA;
1148 +  char* atomB;
1149 +  
1150 +  Atom** the_atoms;
1151 +  the_atoms = entry_plug->atoms;
1152 +  
1153 +
1154 +  // initialize the Bonds
1155 +  
1156 +  for( i=0; i<nBonds; i++ ){
1157 +    
1158 +    a = the_bonds[i].a;
1159 +    b = the_bonds[i].b;
1160 +
1161 +    atomA = the_atoms[a]->getType();
1162 +    atomB = the_atoms[b]->getType();
1163 +    currentBondType = headBondType->find( atomA, atomB );
1164 +    if( currentBondType == NULL ){
1165 +      sprintf( painCave.errMsg,
1166 +               "BondType error, %s - %s not found in force file.\n",
1167 +               atomA, atomB );
1168 +      painCave.isFatal = 1;
1169 +      simError();
1170 +    }
1171 +    
1172 +    if( !strcmp( currentBondType->type, "fixed" ) ){
1173 +      
1174 +      bondArray[i] = new ConstrainedBond( *the_atoms[a],
1175 +                                          *the_atoms[b],
1176 +                                          currentBondType->d0 );
1177 +      entry_plug->n_constraints++;
1178 +    }
1179 +  }
1180 + }
1181 +
1182 + void TraPPE_ExFF::initializeBends( int nBends, Bend** bendArray,
1183 +                                   bend_set* the_bends ){
1184 +  
1185 +  QuadraticBend* qBend;
1186 +  GhostBend* gBend;
1187 +  Atom** the_atoms;
1188 +  the_atoms = entry_plug->atoms;
1189 +  
1190 +  int i, a, b, c;
1191 +  char* atomA;
1192 +  char* atomB;
1193 +  char* atomC;
1194 +  
1195 +  // initialize the Bends
1196 +
1197    for( i=0; i<nBends; i++ ){
1198      
1199      a = the_bends[i].a;
# Line 1059 | Line 1217 | void TraPPE_ExFF::initializeBends( bend_set* the_bends
1217      
1218      if( !strcmp( currentBendType->type, "quadratic" ) ){
1219        
1062      index = i + entry_plug->n_bonds;
1063      
1220        if( the_bends[i].isGhost){
1221          
1222          if( the_bends[i].ghost == b ){
# Line 1088 | Line 1244 | void TraPPE_ExFF::initializeBends( bend_set* the_bends
1244                               currentBendType->k2,
1245                               currentBendType->k3,
1246                               currentBendType->t0 );
1247 <        the_sris[index] = gBend;
1247 >        bendArray[i] = gBend;
1248        }
1249        else{
1250          qBend = new QuadraticBend( *the_atoms[a],
# Line 1098 | Line 1254 | void TraPPE_ExFF::initializeBends( bend_set* the_bends
1254                               currentBendType->k2,
1255                               currentBendType->k3,
1256                               currentBendType->t0 );
1257 <        the_sris[index] = qBend;
1257 >        bendArray[i] = qBend;
1258        }
1259      }
1260    }
1105
1106
1107  // clean up the memory
1108  
1109  delete headBendType;
1110
1111 #ifdef IS_MPI
1112  sprintf( checkPointMsg, "TraPPE_Ex bends initialized succesfully" );
1113  MPIcheckPoint();
1114 #endif // is_mpi
1115
1261   }
1262  
1263 < void TraPPE_ExFF::initializeTorsions( torsion_set* the_torsions ){
1263 > void TraPPE_ExFF::initializeTorsions( int nTorsions, Torsion** torsionArray,
1264 >                                      torsion_set* the_torsions ){
1265  
1266 <  class LinkedType {
1121 <  public:
1122 <    LinkedType(){
1123 <      next = NULL;
1124 <      nameA[0] = '\0';
1125 <      nameB[0] = '\0';
1126 <      nameC[0] = '\0';
1127 <      type[0] = '\0';
1128 <    }
1129 <    ~LinkedType(){ if( next != NULL ) delete next; }
1130 <
1131 <    LinkedType* find( char* key1, char* key2, char* key3, char* key4 ){
1132 <      
1133 <
1134 <
1135 <
1136 <      if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 ) &&
1137 <          !strcmp( nameC, key3 ) && !strcmp( nameD, key4 ) ) return this;
1138 <
1139 <      if( !strcmp( nameA, key4 ) && !strcmp( nameB, key3 ) &&
1140 <          !strcmp( nameC, key2 ) && !strcmp( nameD, key1 ) ) return this;
1141 <
1142 <      if( next != NULL ) return next->find(key1, key2, key3, key4);
1143 <      return NULL;
1144 <    }
1145 <
1146 <    void add( torsionStruct &info ){
1147 <
1148 <      // check for duplicates
1149 <      int dup = 0;
1150 <
1151 <      if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB ) &&
1152 <          !strcmp( nameC, info.nameC ) && !strcmp( nameD, info.nameD ) ) dup = 1;
1153 <      
1154 <      if( !strcmp( nameA, info.nameD ) && !strcmp( nameB, info.nameC ) &&
1155 <          !strcmp( nameC, info.nameB ) && !strcmp( nameD, info.nameA ) ) dup = 1;
1156 <      
1157 <      if(dup){
1158 <        sprintf( painCave.errMsg,
1159 <                 "Duplicate TraPPE_Ex torsion type \"%s - %s - %s - %s\" found in "
1160 <                 "the TraPPE_ExFF param file./n", nameA, nameB, nameC, nameD );
1161 <        painCave.isFatal = 1;
1162 <        simError();
1163 <      }
1164 <
1165 <      if( next != NULL ) next->add(info);
1166 <      else{
1167 <        next = new LinkedType();
1168 <        strcpy(next->nameA, info.nameA);
1169 <        strcpy(next->nameB, info.nameB);
1170 <        strcpy(next->nameC, info.nameC);
1171 <        strcpy(next->nameD, info.nameD);
1172 <        strcpy(next->type,  info.type);
1173 <        next->k1 = info.k1;
1174 <        next->k2 = info.k2;
1175 <        next->k3 = info.k3;
1176 <        next->k4 = info.k4;
1177 <
1178 <      }
1179 <    }
1180 <
1181 < #ifdef IS_MPI
1182 <    
1183 <    void duplicate( torsionStruct &info ){
1184 <      strcpy(info.nameA, nameA);
1185 <      strcpy(info.nameB, nameB);
1186 <      strcpy(info.nameC, nameC);
1187 <      strcpy(info.nameD, nameD);
1188 <      strcpy(info.type,  type);
1189 <      info.k1   = k1;
1190 <      info.k2   = k2;
1191 <      info.k3   = k3;
1192 <      info.k4   = k4;
1193 <      info.last = 0;
1194 <    }
1195 <
1196 < #endif
1197 <
1198 <    char nameA[15];
1199 <    char nameB[15];
1200 <    char nameC[15];
1201 <    char nameD[15];
1202 <    char type[30];
1203 <    double k1, k2, k3, k4;
1204 <
1205 <    LinkedType* next;
1206 <  };
1207 <  
1208 <  LinkedType* headTorsionType;
1209 <  LinkedType* currentTorsionType;
1210 <  torsionStruct info;
1211 <  info.last = 1; // initialize last to have the last set.
1212 <                 // if things go well, last will be set to 0
1213 <
1214 <  int i, a, b, c, d, index;
1266 >  int i, a, b, c, d;
1267    char* atomA;
1268    char* atomB;
1269    char* atomC;
1270    char* atomD;
1219  CubicTorsion* cTors;
1271  
1272 <  SRI **the_sris;
1272 >  CubicTorsion* cTors;
1273    Atom** the_atoms;
1223  int nTorsions;
1224  the_sris = entry_plug->sr_interactions;
1274    the_atoms = entry_plug->atoms;
1226  nTorsions = entry_plug->n_torsions;
1275  
1228 #ifdef IS_MPI
1229  if( worldRank == 0 ){
1230 #endif
1231
1232    // read in the torsion types.
1233
1234    headTorsionType = new LinkedType;
1235    
1236    fastForward( "TorsionTypes", "initializeTorsions" );
1237
1238    // we are now at the torsionTypes section
1239
1240    eof_test =  fgets( readLine, sizeof(readLine), frcFile );
1241    lineNum++;
1242    
1243    
1244    // read a line, and start parseing out the atom types
1245
1246    if( eof_test == NULL ){
1247      sprintf( painCave.errMsg,
1248               "Error in reading torsions from force file at line %d.\n",
1249               lineNum );
1250      painCave.isFatal = 1;
1251      simError();
1252    }
1253    
1254    // stop reading at end of file, or at next section
1255    while( readLine[0] != '#' && eof_test != NULL ){
1256
1257      // toss comment lines
1258      if( readLine[0] != '!' ){
1259        
1260        // the parser returns 0 if the line was blank
1261        if( parseTorsion( readLine, lineNum, info ) ){
1262          headTorsionType->add( info );
1263
1264        }
1265      }
1266      eof_test = fgets( readLine, sizeof(readLine), frcFile );
1267      lineNum++;
1268    }
1269
1270 #ifdef IS_MPI
1271    
1272    // send out the linked list to all the other processes
1273    
1274    sprintf( checkPointMsg,
1275             "TraPPE_Ex torsion structures read successfully." );
1276    MPIcheckPoint();
1277    
1278    currentTorsionType = headTorsionType;
1279    while( currentTorsionType != NULL ){
1280      currentTorsionType->duplicate( info );
1281      sendFrcStruct( &info, mpiTorsionStructType );
1282      currentTorsionType = currentTorsionType->next;
1283    }
1284    info.last = 1;
1285    sendFrcStruct( &info, mpiTorsionStructType );
1286    
1287  }
1288
1289  else{
1290    
1291    // listen for node 0 to send out the force params
1292    
1293    MPIcheckPoint();
1294
1295    headTorsionType = new LinkedType;
1296    recieveFrcStruct( &info, mpiTorsionStructType );
1297    while( !info.last ){
1298
1299      headTorsionType->add( info );
1300      recieveFrcStruct( &info, mpiTorsionStructType );
1301    }
1302  }
1303 #endif // is_mpi
1304  
1276    // initialize the Torsions
1277  
1278    for( i=0; i<nTorsions; i++ ){
# Line 1326 | Line 1297 | void TraPPE_ExFF::initializeTorsions( torsion_set* the
1297      }
1298      
1299      if( !strcmp( currentTorsionType->type, "cubic" ) ){
1329      index = i + entry_plug->n_bonds + entry_plug->n_bends;
1300        
1301        cTors = new CubicTorsion( *the_atoms[a], *the_atoms[b],
1302                                  *the_atoms[c], *the_atoms[d] );
1303        cTors->setConstants( currentTorsionType->k1, currentTorsionType->k2,
1304                             currentTorsionType->k3, currentTorsionType->k4 );
1305 <      the_sris[index] = cTors;
1305 >      torsionArray[i] = cTors;
1306      }
1307    }
1338
1339
1340  // clean up the memory
1341  
1342  delete headTorsionType;
1343
1344 #ifdef IS_MPI
1345  sprintf( checkPointMsg, "TraPPE_Ex torsions initialized succesfully" );
1346  MPIcheckPoint();
1347 #endif // is_mpi
1348
1308   }
1309  
1310   void TraPPE_ExFF::fastForward( char* stopText, char* searchOwner ){
# Line 1560 | Line 1519 | int TPE::parseBend( char *lineBuffer, int lineNum, ben
1519  
1520      if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1521        sprintf( painCave.errMsg,
1522 <               "Error parseing BondTypes: line %d\n", lineNum );
1522 >               "Error parseing BendTypes: line %d\n", lineNum );
1523        painCave.isFatal = 1;
1524        simError();
1525      }
# Line 1569 | Line 1528 | int TPE::parseBend( char *lineBuffer, int lineNum, ben
1528  
1529      if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1530        sprintf( painCave.errMsg,
1531 <               "Error parseing BondTypes: line %d\n", lineNum );
1531 >               "Error parseing BendTypes: line %d\n", lineNum );
1532        painCave.isFatal = 1;
1533        simError();
1534      }
# Line 1578 | Line 1537 | int TPE::parseBend( char *lineBuffer, int lineNum, ben
1537  
1538      if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1539        sprintf( painCave.errMsg,
1540 <               "Error parseing BondTypes: line %d\n", lineNum );
1540 >               "Error parseing BendTypes: line %d\n", lineNum );
1541        painCave.isFatal = 1;
1542        simError();
1543      }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines