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 394 by gezelter, Mon Mar 24 21:55:34 2003 UTC vs.
Revision 542 by mmeineke, Fri May 30 21:31:48 2003 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines