ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/TraPPE_ExFF.cpp
Revision: 365
Committed: Tue Mar 18 22:17:31 2003 UTC (21 years, 3 months ago) by mmeineke
File size: 40957 byte(s)
Log Message: 
some initial bug fixes

File Contents

# User Rev Content
1 mmeineke 270 #include <cstdlib>
2     #include <cstdio>
3     #include <cstring>
4    
5     #include <iostream>
6     using namespace std;
7    
8     #include "ForceFields.hpp"
9     #include "SRI.hpp"
10     #include "simError.h"
11    
12 mmeineke 294 #include <fortranWrappers.hpp>
13    
14 mmeineke 291 #ifdef IS_MPI
15     #include "mpiForceField.h"
16     #endif // is_mpi
17 mmeineke 270
18 mmeineke 291 namespace TPE { // restrict the access of the folowing to this file only.
19    
20    
21     // Declare the structures that will be passed by MPI
22    
23     typedef struct{
24     char name[15];
25     double mass;
26     double epslon;
27     double sigma;
28     double dipole;
29     double w0;
30     double v0;
31     int isSSD;
32     int isDipole;
33     int ident;
34     int last; // 0 -> default
35     // 1 -> tells nodes to stop listening
36     } atomStruct;
37    
38    
39     typedef struct{
40     char nameA[15];
41     char nameB[15];
42     char type[30];
43     double d0;
44     int last; // 0 -> default
45     // 1 -> tells nodes to stop listening
46     } bondStruct;
47    
48    
49     typedef struct{
50     char nameA[15];
51     char nameB[15];
52     char nameC[15];
53     char type[30];
54     double k1, k2, k3, t0;
55     int last; // 0 -> default
56     // 1 -> tells nodes to stop listening
57     } bendStruct;
58    
59    
60     typedef struct{
61     char nameA[15];
62     char nameB[15];
63     char nameC[15];
64     char nameD[15];
65     char type[30];
66     double k1, k2, k3, k4;
67     int last; // 0 -> default
68     // 1 -> tells nodes to stop listening
69     } torsionStruct;
70    
71    
72     int parseAtom( char *lineBuffer, int lineNum, atomStruct &info );
73     int parseBond( char *lineBuffer, int lineNum, bondStruct &info );
74     int parseBend( char *lineBuffer, int lineNum, bendStruct &info );
75     int parseTorsion( char *lineBuffer, int lineNum, torsionStruct &info );
76    
77    
78 mmeineke 270 #ifdef IS_MPI
79    
80 mmeineke 291 MPI_Datatype mpiAtomStructType;
81     MPI_Datatype mpiBondStructType;
82     MPI_Datatype mpiBendStructType;
83     MPI_Datatype mpiTorsionStructType;
84 mmeineke 270
85 mmeineke 291 #endif
86 mmeineke 270
87 mmeineke 291 } // namespace
88 mmeineke 270
89 mmeineke 291 using namespace TPE;
90    
91    
92     //****************************************************************
93     // begins the actual forcefield stuff.
94     //****************************************************************
95    
96    
97 mmeineke 270 TraPPE_ExFF::TraPPE_ExFF(){
98    
99     char fileName[200];
100     char* ffPath_env = "FORCE_PARAM_PATH";
101     char* ffPath;
102     char temp[200];
103     char errMsg[1000];
104    
105 mmeineke 291 // do the funtion wrapping
106 mmeineke 294 wrapMeFF( this );
107 mmeineke 291
108    
109 mmeineke 270 #ifdef IS_MPI
110     int i;
111    
112 mmeineke 291 // **********************************************************************
113 mmeineke 270 // Init the atomStruct mpi type
114    
115     atomStruct atomProto; // mpiPrototype
116 mmeineke 291 int atomBC[3] = {15,6,4}; // block counts
117 mmeineke 270 MPI_Aint atomDspls[3]; // displacements
118     MPI_Datatype atomMbrTypes[3]; // member mpi types
119    
120     MPI_Address(&atomProto.name, &atomDspls[0]);
121     MPI_Address(&atomProto.mass, &atomDspls[1]);
122 mmeineke 291 MPI_Address(&atomProto.isSSD, &atomDspls[2]);
123 mmeineke 270
124     atomMbrTypes[0] = MPI_CHAR;
125     atomMbrTypes[1] = MPI_DOUBLE;
126     atomMbrTypes[2] = MPI_INT;
127    
128     for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
129    
130     MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
131     MPI_Type_commit(&mpiAtomStructType);
132    
133    
134     // **********************************************************************
135     // Init the bondStruct mpi type
136    
137     bondStruct bondProto; // mpiPrototype
138     int bondBC[3] = {60,1,1}; // block counts
139     MPI_Aint bondDspls[3]; // displacements
140     MPI_Datatype bondMbrTypes[3]; // member mpi types
141    
142     MPI_Address(&bondProto.nameA, &bondDspls[0]);
143     MPI_Address(&bondProto.d0, &bondDspls[1]);
144     MPI_Address(&bondProto.last, &bondDspls[2]);
145    
146     bondMbrTypes[0] = MPI_CHAR;
147     bondMbrTypes[1] = MPI_DOUBLE;
148     bondMbrTypes[2] = MPI_INT;
149    
150     for (i=2; i >= 0; i--) bondDspls[i] -= bondDspls[0];
151    
152     MPI_Type_struct(3, bondBC, bondDspls, bondMbrTypes, &mpiBondStructType);
153     MPI_Type_commit(&mpiBondStructType);
154    
155    
156     // **********************************************************************
157     // Init the bendStruct mpi type
158    
159     bendStruct bendProto; // mpiPrototype
160     int bendBC[3] = {75,4,1}; // block counts
161     MPI_Aint bendDspls[3]; // displacements
162     MPI_Datatype bendMbrTypes[3]; // member mpi types
163    
164     MPI_Address(&bendProto.nameA, &bendDspls[0]);
165     MPI_Address(&bendProto.k1, &bendDspls[1]);
166     MPI_Address(&bendProto.last, &bendDspls[2]);
167    
168     bendMbrTypes[0] = MPI_CHAR;
169     bendMbrTypes[1] = MPI_DOUBLE;
170     bendMbrTypes[2] = MPI_INT;
171    
172     for (i=2; i >= 0; i--) bendDspls[i] -= bendDspls[0];
173    
174     MPI_Type_struct(3, bendBC, bendDspls, bendMbrTypes, &mpiBendStructType);
175     MPI_Type_commit(&mpiBendStructType);
176    
177    
178     // **********************************************************************
179     // Init the torsionStruct mpi type
180    
181     torsionStruct torsionProto; // mpiPrototype
182     int torsionBC[3] = {90,4,1}; // block counts
183     MPI_Aint torsionDspls[3]; // displacements
184     MPI_Datatype torsionMbrTypes[3]; // member mpi types
185    
186     MPI_Address(&torsionProto.nameA, &torsionDspls[0]);
187     MPI_Address(&torsionProto.k1, &torsionDspls[1]);
188     MPI_Address(&torsionProto.last, &torsionDspls[2]);
189    
190     torsionMbrTypes[0] = MPI_CHAR;
191     torsionMbrTypes[1] = MPI_DOUBLE;
192     torsionMbrTypes[2] = MPI_INT;
193    
194     for (i=2; i >= 0; i--) torsionDspls[i] -= torsionDspls[0];
195    
196     MPI_Type_struct(3, torsionBC, torsionDspls, torsionMbrTypes,
197     &mpiTorsionStructType);
198     MPI_Type_commit(&mpiTorsionStructType);
199    
200     // ***********************************************************************
201    
202     if( worldRank == 0 ){
203     #endif
204    
205     // generate the force file name
206    
207     strcpy( fileName, "TraPPE_Ex.frc" );
208     // fprintf( stderr,"Trying to open %s\n", fileName );
209    
210     // attempt to open the file in the current directory first.
211    
212     frcFile = fopen( fileName, "r" );
213    
214     if( frcFile == NULL ){
215    
216     // next see if the force path enviorment variable is set
217    
218     ffPath = getenv( ffPath_env );
219     if( ffPath == NULL ) {
220 mmeineke 321 STR_DEFINE(ffPath, FRC_PATH );
221 mmeineke 270 }
222    
223    
224     strcpy( temp, ffPath );
225     strcat( temp, "/" );
226     strcat( temp, fileName );
227     strcpy( fileName, temp );
228    
229     frcFile = fopen( fileName, "r" );
230    
231     if( frcFile == NULL ){
232    
233     sprintf( painCave.errMsg,
234     "Error opening the force field parameter file: %s\n"
235     "Have you tried setting the FORCE_PARAM_PATH environment "
236     "vairable?\n",
237     fileName );
238     painCave.isFatal = 1;
239     simError();
240     }
241     }
242    
243     #ifdef IS_MPI
244     }
245    
246     sprintf( checkPointMsg, "TraPPE_ExFF file opened sucessfully." );
247     MPIcheckPoint();
248    
249     #endif // is_mpi
250     }
251    
252    
253     TraPPE_ExFF::~TraPPE_ExFF(){
254    
255     #ifdef IS_MPI
256     if( worldRank == 0 ){
257     #endif // is_mpi
258    
259     fclose( frcFile );
260    
261     #ifdef IS_MPI
262     }
263     #endif // is_mpi
264     }
265    
266 mmeineke 362 void TraPPE_ExFF::initForceField( int ljMixRule ){
267    
268     initFortran( ljMixRule, entry_plug->useReactionField );
269     }
270 mmeineke 296
271 mmeineke 362
272 mmeineke 270 void TraPPE_ExFF::initializeAtoms( void ){
273    
274     class LinkedType {
275     public:
276     LinkedType(){
277     next = NULL;
278     name[0] = '\0';
279     }
280     ~LinkedType(){ if( next != NULL ) delete next; }
281    
282     LinkedType* find(char* key){
283     if( !strcmp(name, key) ) return this;
284     if( next != NULL ) return next->find(key);
285     return NULL;
286     }
287    
288     void add( atomStruct &info ){
289 mmeineke 291
290     // check for duplicates
291    
292     if( !strcmp( info.name, name ) ){
293     sprintf( painCave.errMsg,
294     "Duplicate TraPPE_Ex atom type \"%s\" found in "
295     "the TraPPE_ExFF param file./n",
296     name );
297     painCave.isFatal = 1;
298     simError();
299     }
300    
301 mmeineke 270 if( next != NULL ) next->add(info);
302     else{
303     next = new LinkedType();
304     strcpy(next->name, info.name);
305 mmeineke 291 next->isDipole = info.isDipole;
306     next->isSSD = info.isSSD;
307 mmeineke 270 next->mass = info.mass;
308     next->epslon = info.epslon;
309     next->sigma = info.sigma;
310     next->dipole = info.dipole;
311 mmeineke 291 next->w0 = info.w0;
312     next->v0 = info.v0;
313     next->ident = info.ident;
314 mmeineke 270 }
315     }
316 mmeineke 291
317     #ifdef IS_MPI
318 mmeineke 270
319     void duplicate( atomStruct &info ){
320     strcpy(info.name, name);
321 mmeineke 291 info.isDipole = isDipole;
322     info.isSSD = isSSD;
323 mmeineke 270 info.mass = mass;
324     info.epslon = epslon;
325     info.sigma = sigma;
326     info.dipole = dipole;
327 mmeineke 291 info.w0 = w0;
328     info.v0 = v0;
329 mmeineke 270 info.last = 0;
330     }
331    
332    
333     #endif
334    
335     char name[15];
336     int isDipole;
337 mmeineke 291 int isSSD;
338 mmeineke 270 double mass;
339     double epslon;
340     double sigma;
341     double dipole;
342 mmeineke 291 double w0;
343     double v0;
344     int ident;
345 mmeineke 270 LinkedType* next;
346     };
347    
348     LinkedType* headAtomType;
349     LinkedType* currentAtomType;
350     atomStruct info;
351     info.last = 1; // initialize last to have the last set.
352     // if things go well, last will be set to 0
353 mmeineke 291
354 mmeineke 270
355    
356     int i;
357 mmeineke 291 int identNum;
358 mmeineke 270
359 mmeineke 291 Atom** the_atoms;
360     int nAtoms;
361     the_atoms = entry_plug->atoms;
362     nAtoms = entry_plug->n_atoms;
363    
364 mmeineke 270 //////////////////////////////////////////////////
365     // a quick water fix
366    
367     double waterI[3][3];
368     waterI[0][0] = 1.76958347772500;
369     waterI[0][1] = 0.0;
370     waterI[0][2] = 0.0;
371    
372     waterI[1][0] = 0.0;
373     waterI[1][1] = 0.614537057924513;
374     waterI[1][2] = 0.0;
375    
376     waterI[2][0] = 0.0;
377     waterI[2][1] = 0.0;
378     waterI[2][2] = 1.15504641980049;
379    
380    
381     double headI[3][3];
382     headI[0][0] = 1125;
383     headI[0][1] = 0.0;
384     headI[0][2] = 0.0;
385    
386     headI[1][0] = 0.0;
387     headI[1][1] = 1125;
388     headI[1][2] = 0.0;
389    
390     headI[2][0] = 0.0;
391     headI[2][1] = 0.0;
392     headI[2][2] = 250;
393    
394    
395    
396     //////////////////////////////////////////////////
397    
398 mmeineke 291
399 mmeineke 270 #ifdef IS_MPI
400     if( worldRank == 0 ){
401     #endif
402    
403     // read in the atom types.
404 mmeineke 291
405 mmeineke 270 headAtomType = new LinkedType;
406    
407 mmeineke 291 fastForward( "AtomTypes", "initializeAtoms" );
408    
409 mmeineke 270 // we are now at the AtomTypes section.
410    
411     eof_test = fgets( readLine, sizeof(readLine), frcFile );
412     lineNum++;
413    
414 mmeineke 291
415     // read a line, and start parseing out the atom types
416    
417 mmeineke 270 if( eof_test == NULL ){
418     sprintf( painCave.errMsg,
419     "Error in reading Atoms from force file at line %d.\n",
420     lineNum );
421     painCave.isFatal = 1;
422     simError();
423     }
424    
425 mmeineke 291 identNum = 1;
426     // stop reading at end of file, or at next section
427 mmeineke 270 while( readLine[0] != '#' && eof_test != NULL ){
428 mmeineke 291
429     // toss comment lines
430 mmeineke 270 if( readLine[0] != '!' ){
431    
432 mmeineke 291 // the parser returns 0 if the line was blank
433     if( parseAtom( readLine, lineNum, info ) ){
434     info.ident = identNum;
435     headAtomType->add( info );;
436     identNum++;
437 mmeineke 270 }
438     }
439     eof_test = fgets( readLine, sizeof(readLine), frcFile );
440     lineNum++;
441     }
442    
443     #ifdef IS_MPI
444    
445     // send out the linked list to all the other processes
446    
447     sprintf( checkPointMsg,
448 mmeineke 291 "TraPPE_ExFF atom structures read successfully." );
449 mmeineke 270 MPIcheckPoint();
450    
451 mmeineke 291 currentAtomType = headAtomType->next; //skip the first element who is a place holder.
452 mmeineke 270 while( currentAtomType != NULL ){
453     currentAtomType->duplicate( info );
454 mmeineke 291
455    
456    
457 mmeineke 270 sendFrcStruct( &info, mpiAtomStructType );
458 mmeineke 291
459     sprintf( checkPointMsg,
460     "successfully sent TraPPE_Ex force type: \"%s\"\n",
461     info.name );
462     MPIcheckPoint();
463    
464 mmeineke 270 currentAtomType = currentAtomType->next;
465     }
466     info.last = 1;
467     sendFrcStruct( &info, mpiAtomStructType );
468    
469     }
470    
471     else{
472    
473     // listen for node 0 to send out the force params
474    
475     MPIcheckPoint();
476    
477     headAtomType = new LinkedType;
478     recieveFrcStruct( &info, mpiAtomStructType );
479 mmeineke 291
480 mmeineke 270 while( !info.last ){
481    
482 mmeineke 291
483    
484 mmeineke 270 headAtomType->add( info );
485 mmeineke 291
486     MPIcheckPoint();
487    
488 mmeineke 270 recieveFrcStruct( &info, mpiAtomStructType );
489     }
490     }
491     #endif // is_mpi
492    
493 mmeineke 291 // call new A_types in fortran
494 mmeineke 270
495 mmeineke 291 int isError;
496 mmeineke 296
497     // dummy variables
498    
499     int isGB = 0;
500     int isLJ = 1;
501 mmeineke 362 double GB_dummy = 0.0;
502 mmeineke 296
503    
504 mmeineke 291 currentAtomType = headAtomType;
505     while( currentAtomType != NULL ){
506    
507 mmeineke 362 if(currentAtomType->isDipole) entry_plug->useReactionField = 1;
508     if(currentAtomType->isDipole) entry_plug->useDipole = 1;
509     if(currentAtomType->isSSD) entry_plug->useSticky = 1;
510    
511 mmeineke 291 if( currentAtomType->name[0] != '\0' ){
512     isError = 0;
513 mmeineke 362 makeAtype( &(currentAtomType->ident),
514     &isLJ,
515     &(currentAtomType->isSSD),
516     &(currentAtomType->isDipole),
517     &isGB,
518     &(currentAtomType->epslon),
519     &(currentAtomType->sigma),
520     &(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,
529     &isError );
530 mmeineke 291 if( isError ){
531     sprintf( painCave.errMsg,
532     "Error initializing the \"%s\" atom type in fortran\n",
533     currentAtomType->name );
534     painCave.isFatal = 1;
535     simError();
536     }
537     }
538     currentAtomType = currentAtomType->next;
539     }
540    
541     #ifdef IS_MPI
542     sprintf( checkPointMsg,
543     "TraPPE_ExFF atom structures successfully sent to fortran\n" );
544     MPIcheckPoint();
545     #endif // is_mpi
546    
547    
548 mmeineke 270 // initialize the atoms
549    
550 mmeineke 291 double bigSigma = 0.0;
551 mmeineke 270 DirectionalAtom* dAtom;
552    
553     for( i=0; i<nAtoms; i++ ){
554    
555     currentAtomType = headAtomType->find( the_atoms[i]->getType() );
556     if( currentAtomType == NULL ){
557     sprintf( painCave.errMsg,
558     "AtomType error, %s not found in force file.\n",
559     the_atoms[i]->getType() );
560     painCave.isFatal = 1;
561     simError();
562     }
563    
564     the_atoms[i]->setMass( currentAtomType->mass );
565     the_atoms[i]->setEpslon( currentAtomType->epslon );
566     the_atoms[i]->setSigma( currentAtomType->sigma );
567     the_atoms[i]->setLJ();
568    
569 mmeineke 291 if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
570    
571 mmeineke 270 if( currentAtomType->isDipole ){
572     if( the_atoms[i]->isDirectional() ){
573    
574     dAtom = (DirectionalAtom *) the_atoms[i];
575     dAtom->setMu( currentAtomType->dipole );
576     dAtom->setHasDipole( 1 );
577     dAtom->setJx( 0.0 );
578     dAtom->setJy( 0.0 );
579     dAtom->setJz( 0.0 );
580    
581     if(!strcmp("SSD",the_atoms[i]->getType())){
582     dAtom->setI( waterI );
583     dAtom->setSSD( 1 );
584     }
585     else if(!strcmp("HEAD",the_atoms[i]->getType())){
586     dAtom->setI( headI );
587     dAtom->setSSD( 0 );
588     }
589     else{
590     sprintf(painCave.errMsg,
591     "AtmType error, %s does not have a moment of inertia set.\n",
592     the_atoms[i]->getType() );
593     painCave.isFatal = 1;
594     simError();
595     }
596     entry_plug->n_dipoles++;
597     }
598     else{
599    
600     sprintf( painCave.errMsg,
601     "TraPPE_ExFF error: Atom \"%s\" is a dipole, yet no standard"
602     " orientation was specifed in the BASS file.\n",
603     currentAtomType->name );
604     painCave.isFatal = 1;
605     simError();
606     }
607     }
608     else{
609     if( the_atoms[i]->isDirectional() ){
610     sprintf( painCave.errMsg,
611     "TraPPE_ExFF error: Atom \"%s\" was given a standard"
612     "orientation in the BASS file, yet it is not a dipole.\n",
613     currentAtomType->name);
614     painCave.isFatal = 1;
615     simError();
616     }
617     }
618     }
619    
620 mmeineke 291 #ifdef IS_MPI
621     double tempBig = bigSigma;
622     MPI::COMM_WORLD.Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX );
623     #endif //is_mpi
624 mmeineke 270
625 mmeineke 291 //calc rCut and rList
626    
627     entry_plug->rCut = 2.5 * bigSigma;
628    
629     if(entry_plug->rCut > (entry_plug->box_x / 2.0))
630     entry_plug->rCut = entry_plug->box_x / 2.0;
631    
632     if(entry_plug->rCut > (entry_plug->box_y / 2.0))
633     entry_plug->rCut = entry_plug->box_y / 2.0;
634    
635     if(entry_plug->rCut > (entry_plug->box_z / 2.0))
636     entry_plug->rCut = entry_plug->box_z / 2.0;
637    
638     entry_plug->rList = entry_plug->rCut + 1.0;
639    
640 mmeineke 362 entry_plug->useLJ = 1; // use Lennard Jones is on by default
641 mmeineke 291
642 mmeineke 270 // clean up the memory
643    
644     delete headAtomType;
645    
646     #ifdef IS_MPI
647     sprintf( checkPointMsg, "TraPPE_Ex atoms initialized succesfully" );
648     MPIcheckPoint();
649     #endif // is_mpi
650    
651     }
652    
653     void TraPPE_ExFF::initializeBonds( bond_pair* the_bonds ){
654    
655     class LinkedType {
656     public:
657     LinkedType(){
658     next = NULL;
659     nameA[0] = '\0';
660     nameB[0] = '\0';
661     type[0] = '\0';
662     }
663     ~LinkedType(){ if( next != NULL ) delete next; }
664    
665     LinkedType* find(char* key1, char* key2){
666     if( !strcmp(nameA, key1 ) && !strcmp( nameB, key2 ) ) return this;
667     if( !strcmp(nameA, key2 ) && !strcmp( nameB, key1 ) ) return this;
668     if( next != NULL ) return next->find(key1, key2);
669     return NULL;
670     }
671    
672 mmeineke 291
673 mmeineke 270 void add( bondStruct &info ){
674 mmeineke 291
675     // check for duplicates
676     int dup = 0;
677    
678     if( !strcmp(nameA, info.nameA ) && !strcmp( nameB, info.nameB ) ) dup = 1;
679     if( !strcmp(nameA, info.nameB ) && !strcmp( nameB, info.nameA ) ) dup = 1;
680    
681     if(dup){
682     sprintf( painCave.errMsg,
683     "Duplicate TraPPE_Ex bond type \"%s - %s\" found in "
684     "the TraPPE_ExFF param file./n",
685     nameA, nameB );
686     painCave.isFatal = 1;
687     simError();
688     }
689    
690    
691 mmeineke 270 if( next != NULL ) next->add(info);
692     else{
693     next = new LinkedType();
694     strcpy(next->nameA, info.nameA);
695     strcpy(next->nameB, info.nameB);
696     strcpy(next->type, info.type);
697     next->d0 = info.d0;
698     }
699     }
700    
701 mmeineke 291 #ifdef IS_MPI
702 mmeineke 270 void duplicate( bondStruct &info ){
703     strcpy(info.nameA, nameA);
704     strcpy(info.nameB, nameB);
705     strcpy(info.type, type);
706     info.d0 = d0;
707     info.last = 0;
708     }
709    
710    
711     #endif
712    
713     char nameA[15];
714     char nameB[15];
715     char type[30];
716     double d0;
717    
718     LinkedType* next;
719     };
720    
721    
722    
723     LinkedType* headBondType;
724     LinkedType* currentBondType;
725     bondStruct info;
726     info.last = 1; // initialize last to have the last set.
727     // if things go well, last will be set to 0
728    
729     SRI **the_sris;
730     Atom** the_atoms;
731     int nBonds;
732     the_sris = entry_plug->sr_interactions;
733     the_atoms = entry_plug->atoms;
734     nBonds = entry_plug->n_bonds;
735    
736 mmeineke 291 int i, a, b;
737     char* atomA;
738     char* atomB;
739 mmeineke 270
740     #ifdef IS_MPI
741     if( worldRank == 0 ){
742     #endif
743    
744     // read in the bond types.
745    
746     headBondType = new LinkedType;
747    
748 mmeineke 291 fastForward( "BondTypes", "initializeBonds" );
749    
750     // we are now at the bondTypes section
751    
752     eof_test = fgets( readLine, sizeof(readLine), frcFile );
753 mmeineke 270 lineNum++;
754    
755    
756 mmeineke 291 // read a line, and start parseing out the atom types
757    
758 mmeineke 270 if( eof_test == NULL ){
759     sprintf( painCave.errMsg,
760 mmeineke 291 "Error in reading bonds from force file at line %d.\n",
761 mmeineke 270 lineNum );
762     painCave.isFatal = 1;
763     simError();
764     }
765    
766 mmeineke 291 // stop reading at end of file, or at next section
767 mmeineke 270 while( readLine[0] != '#' && eof_test != NULL ){
768 mmeineke 291
769     // toss comment lines
770 mmeineke 270 if( readLine[0] != '!' ){
771    
772 mmeineke 291 // the parser returns 0 if the line was blank
773     if( parseBond( readLine, lineNum, info ) ){
774     headBondType->add( info );
775 mmeineke 270 }
776     }
777     eof_test = fgets( readLine, sizeof(readLine), frcFile );
778     lineNum++;
779     }
780 mmeineke 291
781 mmeineke 270 #ifdef IS_MPI
782    
783     // send out the linked list to all the other processes
784    
785     sprintf( checkPointMsg,
786     "TraPPE_Ex bond structures read successfully." );
787     MPIcheckPoint();
788    
789     currentBondType = headBondType;
790     while( currentBondType != NULL ){
791     currentBondType->duplicate( info );
792     sendFrcStruct( &info, mpiBondStructType );
793     currentBondType = currentBondType->next;
794     }
795     info.last = 1;
796     sendFrcStruct( &info, mpiBondStructType );
797    
798     }
799    
800     else{
801    
802     // listen for node 0 to send out the force params
803    
804     MPIcheckPoint();
805    
806     headBondType = new LinkedType;
807     recieveFrcStruct( &info, mpiBondStructType );
808     while( !info.last ){
809    
810     headBondType->add( info );
811     recieveFrcStruct( &info, mpiBondStructType );
812     }
813     }
814     #endif // is_mpi
815    
816    
817     // initialize the Bonds
818    
819    
820     for( i=0; i<nBonds; i++ ){
821    
822     a = the_bonds[i].a;
823     b = the_bonds[i].b;
824    
825     atomA = the_atoms[a]->getType();
826     atomB = the_atoms[b]->getType();
827     currentBondType = headBondType->find( atomA, atomB );
828     if( currentBondType == NULL ){
829     sprintf( painCave.errMsg,
830     "BondType error, %s - %s not found in force file.\n",
831     atomA, atomB );
832     painCave.isFatal = 1;
833     simError();
834     }
835    
836     if( !strcmp( currentBondType->type, "fixed" ) ){
837    
838     the_sris[i] = new ConstrainedBond( *the_atoms[a],
839     *the_atoms[b],
840     currentBondType->d0 );
841     entry_plug->n_constraints++;
842     }
843     }
844    
845    
846     // clean up the memory
847    
848     delete headBondType;
849    
850     #ifdef IS_MPI
851     sprintf( checkPointMsg, "TraPPE_Ex bonds initialized succesfully" );
852     MPIcheckPoint();
853     #endif // is_mpi
854    
855     }
856    
857     void TraPPE_ExFF::initializeBends( bend_set* the_bends ){
858    
859     class LinkedType {
860     public:
861     LinkedType(){
862     next = NULL;
863     nameA[0] = '\0';
864     nameB[0] = '\0';
865     nameC[0] = '\0';
866     type[0] = '\0';
867     }
868     ~LinkedType(){ if( next != NULL ) delete next; }
869    
870     LinkedType* find( char* key1, char* key2, char* key3 ){
871     if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 )
872     && !strcmp( nameC, key3 ) ) return this;
873     if( !strcmp( nameA, key3 ) && !strcmp( nameB, key2 )
874     && !strcmp( nameC, key1 ) ) return this;
875     if( next != NULL ) return next->find(key1, key2, key3);
876     return NULL;
877     }
878    
879 mmeineke 291 void add( bendStruct &info ){
880 mmeineke 270
881 mmeineke 291 // check for duplicates
882     int dup = 0;
883    
884     if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB )
885     && !strcmp( nameC, info.nameC ) ) dup = 1;
886     if( !strcmp( nameA, info.nameC ) && !strcmp( nameB, info.nameB )
887     && !strcmp( nameC, info.nameA ) ) dup = 1;
888    
889     if(dup){
890     sprintf( painCave.errMsg,
891     "Duplicate TraPPE_Ex bend type \"%s - %s - %s\" found in "
892     "the TraPPE_ExFF param file./n",
893     nameA, nameB, nameC );
894     painCave.isFatal = 1;
895     simError();
896     }
897    
898 mmeineke 270 if( next != NULL ) next->add(info);
899     else{
900     next = new LinkedType();
901     strcpy(next->nameA, info.nameA);
902     strcpy(next->nameB, info.nameB);
903     strcpy(next->nameC, info.nameC);
904     strcpy(next->type, info.type);
905     next->k1 = info.k1;
906     next->k2 = info.k2;
907     next->k3 = info.k3;
908     next->t0 = info.t0;
909     }
910     }
911 mmeineke 291
912     #ifdef IS_MPI
913    
914 mmeineke 270 void duplicate( bendStruct &info ){
915     strcpy(info.nameA, nameA);
916     strcpy(info.nameB, nameB);
917     strcpy(info.nameC, nameC);
918     strcpy(info.type, type);
919     info.k1 = k1;
920     info.k2 = k2;
921     info.k3 = k3;
922     info.t0 = t0;
923     info.last = 0;
924     }
925    
926     #endif // is_mpi
927    
928     char nameA[15];
929     char nameB[15];
930     char nameC[15];
931     char type[30];
932     double k1, k2, k3, t0;
933    
934     LinkedType* next;
935     };
936    
937     LinkedType* headBendType;
938     LinkedType* currentBendType;
939     bendStruct info;
940     info.last = 1; // initialize last to have the last set.
941     // if things go well, last will be set to 0
942    
943     QuadraticBend* qBend;
944 mmeineke 311 GhostBend* gBend;
945 mmeineke 270 SRI **the_sris;
946     Atom** the_atoms;
947     int nBends;
948     the_sris = entry_plug->sr_interactions;
949     the_atoms = entry_plug->atoms;
950     nBends = entry_plug->n_bends;
951    
952 mmeineke 291 int i, a, b, c;
953     char* atomA;
954     char* atomB;
955     char* atomC;
956 mmeineke 270
957 mmeineke 291
958 mmeineke 270 #ifdef IS_MPI
959     if( worldRank == 0 ){
960     #endif
961    
962     // read in the bend types.
963    
964     headBendType = new LinkedType;
965    
966 mmeineke 291 fastForward( "BendTypes", "initializeBends" );
967 mmeineke 270
968 mmeineke 291 // we are now at the bendTypes section
969 mmeineke 270
970 mmeineke 291 eof_test = fgets( readLine, sizeof(readLine), frcFile );
971     lineNum++;
972    
973     // read a line, and start parseing out the bend types
974 mmeineke 270
975     if( eof_test == NULL ){
976     sprintf( painCave.errMsg,
977 mmeineke 291 "Error in reading bends from force file at line %d.\n",
978 mmeineke 270 lineNum );
979     painCave.isFatal = 1;
980     simError();
981     }
982 mmeineke 291
983     // stop reading at end of file, or at next section
984 mmeineke 270 while( readLine[0] != '#' && eof_test != NULL ){
985 mmeineke 291
986     // toss comment lines
987 mmeineke 270 if( readLine[0] != '!' ){
988    
989 mmeineke 291 // the parser returns 0 if the line was blank
990     if( parseBend( readLine, lineNum, info ) ){
991     headBendType->add( info );
992 mmeineke 270 }
993     }
994     eof_test = fgets( readLine, sizeof(readLine), frcFile );
995     lineNum++;
996     }
997 mmeineke 291
998 mmeineke 270 #ifdef IS_MPI
999    
1000     // send out the linked list to all the other processes
1001    
1002     sprintf( checkPointMsg,
1003     "TraPPE_Ex bend structures read successfully." );
1004     MPIcheckPoint();
1005    
1006     currentBendType = headBendType;
1007     while( currentBendType != NULL ){
1008     currentBendType->duplicate( info );
1009     sendFrcStruct( &info, mpiBendStructType );
1010     currentBendType = currentBendType->next;
1011     }
1012     info.last = 1;
1013     sendFrcStruct( &info, mpiBendStructType );
1014    
1015     }
1016    
1017     else{
1018    
1019     // listen for node 0 to send out the force params
1020    
1021     MPIcheckPoint();
1022    
1023     headBendType = new LinkedType;
1024     recieveFrcStruct( &info, mpiBendStructType );
1025     while( !info.last ){
1026    
1027     headBendType->add( info );
1028     recieveFrcStruct( &info, mpiBendStructType );
1029     }
1030     }
1031     #endif // is_mpi
1032    
1033     // initialize the Bends
1034 mmeineke 291
1035     int index;
1036 mmeineke 270
1037     for( i=0; i<nBends; i++ ){
1038    
1039     a = the_bends[i].a;
1040     b = the_bends[i].b;
1041     c = the_bends[i].c;
1042    
1043     atomA = the_atoms[a]->getType();
1044     atomB = the_atoms[b]->getType();
1045 mmeineke 311
1046     if( the_bends[i].isGhost ) atomC = "GHOST";
1047     else atomC = the_atoms[c]->getType();
1048    
1049 mmeineke 270 currentBendType = headBendType->find( atomA, atomB, atomC );
1050     if( currentBendType == NULL ){
1051     sprintf( painCave.errMsg, "BendType error, %s - %s - %s not found"
1052     " in force file.\n",
1053     atomA, atomB, atomC );
1054     painCave.isFatal = 1;
1055     simError();
1056     }
1057    
1058     if( !strcmp( currentBendType->type, "quadratic" ) ){
1059    
1060     index = i + entry_plug->n_bonds;
1061 mmeineke 311
1062     if( the_bends[i].isGhost){
1063    
1064     if( the_bends[i].ghost == b ){
1065     // do nothing
1066     }
1067     else if( the_bends[i].ghost == a ){
1068     c = a;
1069     a = b;
1070     b = a;
1071     }
1072     else{
1073     sprintf( painCave.errMsg,
1074     "BendType error, %s - %s - %s,\n"
1075     " --> central atom is not "
1076     "correctly identified with the "
1077     "\"ghostVectorSource = \" tag.\n",
1078     atomA, atomB, atomC );
1079     painCave.isFatal = 1;
1080     simError();
1081     }
1082    
1083     gBend = new GhostBend( *the_atoms[a],
1084     *the_atoms[b] );
1085     gBend->setConstants( currentBendType->k1,
1086     currentBendType->k2,
1087     currentBendType->k3,
1088     currentBendType->t0 );
1089     the_sris[index] = gBend;
1090     }
1091     else{
1092     qBend = new QuadraticBend( *the_atoms[a],
1093     *the_atoms[b],
1094     *the_atoms[c] );
1095     qBend->setConstants( currentBendType->k1,
1096     currentBendType->k2,
1097     currentBendType->k3,
1098     currentBendType->t0 );
1099     the_sris[index] = qBend;
1100     }
1101 mmeineke 270 }
1102     }
1103    
1104    
1105     // clean up the memory
1106    
1107     delete headBendType;
1108    
1109     #ifdef IS_MPI
1110     sprintf( checkPointMsg, "TraPPE_Ex bends initialized succesfully" );
1111     MPIcheckPoint();
1112     #endif // is_mpi
1113    
1114     }
1115    
1116     void TraPPE_ExFF::initializeTorsions( torsion_set* the_torsions ){
1117    
1118     class LinkedType {
1119     public:
1120     LinkedType(){
1121     next = NULL;
1122     nameA[0] = '\0';
1123     nameB[0] = '\0';
1124     nameC[0] = '\0';
1125     type[0] = '\0';
1126     }
1127     ~LinkedType(){ if( next != NULL ) delete next; }
1128    
1129     LinkedType* find( char* key1, char* key2, char* key3, char* key4 ){
1130     if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 ) &&
1131     !strcmp( nameC, key3 ) && !strcmp( nameD, key4 ) ) return this;
1132    
1133     if( !strcmp( nameA, key4 ) && !strcmp( nameB, key3 ) &&
1134 mmeineke 291 !strcmp( nameC, key2 ) && !strcmp( nameD, key1 ) ) return this;
1135 mmeineke 270
1136     if( next != NULL ) return next->find(key1, key2, key3, key4);
1137     return NULL;
1138     }
1139    
1140     void add( torsionStruct &info ){
1141 mmeineke 291
1142     // check for duplicates
1143     int dup = 0;
1144    
1145     if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB ) &&
1146     !strcmp( nameC, info.nameC ) && !strcmp( nameD, info.nameD ) ) dup = 1;
1147    
1148     if( !strcmp( nameA, info.nameD ) && !strcmp( nameB, info.nameC ) &&
1149     !strcmp( nameC, info.nameB ) && !strcmp( nameD, info.nameA ) ) dup = 1;
1150    
1151     if(dup){
1152     sprintf( painCave.errMsg,
1153     "Duplicate TraPPE_Ex torsion type \"%s - %s - %s - %s\" found in "
1154     "the TraPPE_ExFF param file./n", nameA, nameB, nameC, nameD );
1155     painCave.isFatal = 1;
1156     simError();
1157     }
1158    
1159 mmeineke 270 if( next != NULL ) next->add(info);
1160     else{
1161     next = new LinkedType();
1162     strcpy(next->nameA, info.nameA);
1163     strcpy(next->nameB, info.nameB);
1164     strcpy(next->nameC, info.nameC);
1165     strcpy(next->type, info.type);
1166     next->k1 = info.k1;
1167     next->k2 = info.k2;
1168     next->k3 = info.k3;
1169     next->k4 = info.k4;
1170     }
1171     }
1172 mmeineke 291
1173     #ifdef IS_MPI
1174 mmeineke 270
1175     void duplicate( torsionStruct &info ){
1176     strcpy(info.nameA, nameA);
1177     strcpy(info.nameB, nameB);
1178     strcpy(info.nameC, nameC);
1179     strcpy(info.nameD, nameD);
1180     strcpy(info.type, type);
1181     info.k1 = k1;
1182     info.k2 = k2;
1183     info.k3 = k3;
1184     info.k4 = k4;
1185     info.last = 0;
1186     }
1187    
1188     #endif
1189    
1190     char nameA[15];
1191     char nameB[15];
1192     char nameC[15];
1193     char nameD[15];
1194     char type[30];
1195     double k1, k2, k3, k4;
1196    
1197     LinkedType* next;
1198     };
1199    
1200     LinkedType* headTorsionType;
1201     LinkedType* currentTorsionType;
1202     torsionStruct info;
1203     info.last = 1; // initialize last to have the last set.
1204     // if things go well, last will be set to 0
1205    
1206     int i, a, b, c, d, index;
1207     char* atomA;
1208     char* atomB;
1209     char* atomC;
1210     char* atomD;
1211     CubicTorsion* cTors;
1212    
1213     SRI **the_sris;
1214     Atom** the_atoms;
1215     int nTorsions;
1216     the_sris = entry_plug->sr_interactions;
1217     the_atoms = entry_plug->atoms;
1218     nTorsions = entry_plug->n_torsions;
1219    
1220     #ifdef IS_MPI
1221     if( worldRank == 0 ){
1222     #endif
1223    
1224     // read in the torsion types.
1225    
1226     headTorsionType = new LinkedType;
1227 mmeineke 291
1228     fastForward( "TorsionTypes", "initializeTorsions" );
1229 mmeineke 270
1230 mmeineke 291 // we are now at the torsionTypes section
1231 mmeineke 270
1232 mmeineke 291 eof_test = fgets( readLine, sizeof(readLine), frcFile );
1233     lineNum++;
1234 mmeineke 270
1235 mmeineke 291
1236     // read a line, and start parseing out the atom types
1237 mmeineke 270
1238     if( eof_test == NULL ){
1239     sprintf( painCave.errMsg,
1240 mmeineke 291 "Error in reading torsions from force file at line %d.\n",
1241 mmeineke 270 lineNum );
1242     painCave.isFatal = 1;
1243     simError();
1244     }
1245 mmeineke 291
1246     // stop reading at end of file, or at next section
1247     while( readLine[0] != '#' && eof_test != NULL ){
1248 mmeineke 270
1249 mmeineke 291 // toss comment lines
1250 mmeineke 270 if( readLine[0] != '!' ){
1251    
1252 mmeineke 291 // the parser returns 0 if the line was blank
1253     if( parseTorsion( readLine, lineNum, info ) ){
1254     headTorsionType->add( info );
1255 mmeineke 270 }
1256     }
1257     eof_test = fgets( readLine, sizeof(readLine), frcFile );
1258     lineNum++;
1259     }
1260    
1261     #ifdef IS_MPI
1262    
1263     // send out the linked list to all the other processes
1264    
1265     sprintf( checkPointMsg,
1266     "TraPPE_Ex torsion structures read successfully." );
1267     MPIcheckPoint();
1268    
1269     currentTorsionType = headTorsionType;
1270     while( currentTorsionType != NULL ){
1271     currentTorsionType->duplicate( info );
1272     sendFrcStruct( &info, mpiTorsionStructType );
1273     currentTorsionType = currentTorsionType->next;
1274     }
1275     info.last = 1;
1276     sendFrcStruct( &info, mpiTorsionStructType );
1277    
1278     }
1279    
1280     else{
1281    
1282     // listen for node 0 to send out the force params
1283    
1284     MPIcheckPoint();
1285    
1286     headTorsionType = new LinkedType;
1287     recieveFrcStruct( &info, mpiTorsionStructType );
1288     while( !info.last ){
1289    
1290     headTorsionType->add( info );
1291     recieveFrcStruct( &info, mpiTorsionStructType );
1292     }
1293     }
1294     #endif // is_mpi
1295    
1296     // initialize the Torsions
1297    
1298     for( i=0; i<nTorsions; i++ ){
1299    
1300     a = the_torsions[i].a;
1301     b = the_torsions[i].b;
1302     c = the_torsions[i].c;
1303     d = the_torsions[i].d;
1304    
1305     atomA = the_atoms[a]->getType();
1306     atomB = the_atoms[b]->getType();
1307     atomC = the_atoms[c]->getType();
1308     atomD = the_atoms[d]->getType();
1309     currentTorsionType = headTorsionType->find( atomA, atomB, atomC, atomD );
1310     if( currentTorsionType == NULL ){
1311     sprintf( painCave.errMsg,
1312     "TorsionType error, %s - %s - %s - %s not found"
1313     " in force file.\n",
1314     atomA, atomB, atomC, atomD );
1315     painCave.isFatal = 1;
1316     simError();
1317     }
1318    
1319     if( !strcmp( currentTorsionType->type, "cubic" ) ){
1320     index = i + entry_plug->n_bonds + entry_plug->n_bends;
1321    
1322     cTors = new CubicTorsion( *the_atoms[a], *the_atoms[b],
1323     *the_atoms[c], *the_atoms[d] );
1324     cTors->setConstants( currentTorsionType->k1, currentTorsionType->k2,
1325     currentTorsionType->k3, currentTorsionType->k4 );
1326     the_sris[index] = cTors;
1327     }
1328     }
1329    
1330    
1331     // clean up the memory
1332    
1333     delete headTorsionType;
1334    
1335     #ifdef IS_MPI
1336     sprintf( checkPointMsg, "TraPPE_Ex torsions initialized succesfully" );
1337     MPIcheckPoint();
1338     #endif // is_mpi
1339    
1340     }
1341 mmeineke 291
1342     void TraPPE_ExFF::fastForward( char* stopText, char* searchOwner ){
1343    
1344     int foundText = 0;
1345     char* the_token;
1346    
1347     rewind( frcFile );
1348     lineNum = 0;
1349    
1350     eof_test = fgets( readLine, sizeof(readLine), frcFile );
1351     lineNum++;
1352     if( eof_test == NULL ){
1353     sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
1354     " file is empty.\n",
1355     searchOwner );
1356     painCave.isFatal = 1;
1357     simError();
1358     }
1359    
1360    
1361     while( !foundText ){
1362     while( eof_test != NULL && readLine[0] != '#' ){
1363     eof_test = fgets( readLine, sizeof(readLine), frcFile );
1364     lineNum++;
1365     }
1366     if( eof_test == NULL ){
1367     sprintf( painCave.errMsg,
1368     "Error fast forwarding force file for %s at "
1369     "line %d: file ended unexpectedly.\n",
1370     searchOwner,
1371     lineNum );
1372     painCave.isFatal = 1;
1373     simError();
1374     }
1375    
1376     the_token = strtok( readLine, " ,;\t#\n" );
1377     foundText = !strcmp( stopText, the_token );
1378    
1379     if( !foundText ){
1380     eof_test = fgets( readLine, sizeof(readLine), frcFile );
1381     lineNum++;
1382    
1383     if( eof_test == NULL ){
1384     sprintf( painCave.errMsg,
1385     "Error fast forwarding force file for %s at "
1386     "line %d: file ended unexpectedly.\n",
1387     searchOwner,
1388     lineNum );
1389     painCave.isFatal = 1;
1390     simError();
1391     }
1392     }
1393     }
1394     }
1395    
1396    
1397     int parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
1398    
1399     char* the_token;
1400    
1401     the_token = strtok( lineBuffer, " \n\t,;" );
1402     if( the_token != NULL ){
1403    
1404     strcpy( info.name, the_token );
1405    
1406     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1407     sprintf( painCave.errMsg,
1408     "Error parseing AtomTypes: line %d\n", lineNum );
1409     painCave.isFatal = 1;
1410     simError();
1411     }
1412    
1413     info.isDipole = atoi( the_token );
1414    
1415     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1416     sprintf( painCave.errMsg,
1417     "Error parseing AtomTypes: line %d\n", lineNum );
1418     painCave.isFatal = 1;
1419     simError();
1420     }
1421    
1422     info.isSSD = atoi( the_token );
1423    
1424     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1425     sprintf( painCave.errMsg,
1426     "Error parseing AtomTypes: line %d\n", lineNum );
1427     painCave.isFatal = 1;
1428     simError();
1429     }
1430    
1431     info.mass = atof( the_token );
1432    
1433     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1434     sprintf( painCave.errMsg,
1435     "Error parseing AtomTypes: line %d\n", lineNum );
1436     painCave.isFatal = 1;
1437     simError();
1438     }
1439    
1440     info.epslon = atof( the_token );
1441    
1442     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1443     sprintf( painCave.errMsg,
1444     "Error parseing AtomTypes: line %d\n", lineNum );
1445     painCave.isFatal = 1;
1446     simError();
1447     }
1448    
1449     info.sigma = atof( the_token );
1450    
1451     if( info.isDipole ){
1452    
1453     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1454     sprintf( painCave.errMsg,
1455     "Error parseing AtomTypes: line %d\n", lineNum );
1456     painCave.isFatal = 1;
1457     simError();
1458     }
1459    
1460     info.dipole = atof( the_token );
1461     }
1462     else info.dipole = 0.0;
1463    
1464     if( info.isSSD ){
1465    
1466     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1467     sprintf( painCave.errMsg,
1468     "Error parseing AtomTypes: line %d\n", lineNum );
1469     painCave.isFatal = 1;
1470     simError();
1471     }
1472    
1473     info.w0 = atof( the_token );
1474    
1475     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1476     sprintf( painCave.errMsg,
1477     "Error parseing AtomTypes: line %d\n", lineNum );
1478     painCave.isFatal = 1;
1479     simError();
1480     }
1481    
1482     info.v0 = atof( the_token );
1483     }
1484     else info.v0 = info.w0 = 0.0;
1485    
1486     return 1;
1487     }
1488     else return 0;
1489     }
1490    
1491     int parseBond( char *lineBuffer, int lineNum, bondStruct &info ){
1492    
1493     char* the_token;
1494    
1495     the_token = strtok( lineBuffer, " \n\t,;" );
1496     if( the_token != NULL ){
1497    
1498     strcpy( info.nameA, the_token );
1499    
1500     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1501     sprintf( painCave.errMsg,
1502     "Error parseing BondTypes: line %d\n", lineNum );
1503     painCave.isFatal = 1;
1504     simError();
1505     }
1506    
1507     strcpy( info.nameB, the_token );
1508    
1509     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1510     sprintf( painCave.errMsg,
1511     "Error parseing BondTypes: line %d\n", lineNum );
1512     painCave.isFatal = 1;
1513     simError();
1514     }
1515    
1516     strcpy( info.type, the_token );
1517    
1518     if( !strcmp( info.type, "fixed" ) ){
1519     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1520     sprintf( painCave.errMsg,
1521     "Error parseing BondTypes: line %d\n", lineNum );
1522     painCave.isFatal = 1;
1523     simError();
1524     }
1525    
1526     info.d0 = atof( the_token );
1527     }
1528     else{
1529     sprintf( painCave.errMsg,
1530     "Unknown TraPPE_Ex bond type \"%s\" at line %d\n",
1531     info.type,
1532     lineNum );
1533     painCave.isFatal = 1;
1534     simError();
1535     }
1536    
1537     return 1;
1538     }
1539     else return 0;
1540     }
1541    
1542    
1543     int parseBend( char *lineBuffer, int lineNum, bendStruct &info ){
1544    
1545     char* the_token;
1546    
1547     the_token = strtok( lineBuffer, " \n\t,;" );
1548     if( the_token != NULL ){
1549    
1550     strcpy( info.nameA, the_token );
1551    
1552     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1553     sprintf( painCave.errMsg,
1554     "Error parseing BondTypes: line %d\n", lineNum );
1555     painCave.isFatal = 1;
1556     simError();
1557     }
1558    
1559     strcpy( info.nameB, the_token );
1560    
1561     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1562     sprintf( painCave.errMsg,
1563     "Error parseing BondTypes: line %d\n", lineNum );
1564     painCave.isFatal = 1;
1565     simError();
1566     }
1567    
1568     strcpy( info.nameC, the_token );
1569    
1570     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1571     sprintf( painCave.errMsg,
1572     "Error parseing BondTypes: line %d\n", lineNum );
1573     painCave.isFatal = 1;
1574     simError();
1575     }
1576    
1577     strcpy( info.type, the_token );
1578    
1579     if( !strcmp( info.type, "quadratic" ) ){
1580     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1581     sprintf( painCave.errMsg,
1582     "Error parseing BendTypes: line %d\n", lineNum );
1583     painCave.isFatal = 1;
1584     simError();
1585     }
1586    
1587     info.k1 = atof( the_token );
1588    
1589     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1590     sprintf( painCave.errMsg,
1591     "Error parseing BendTypes: line %d\n", lineNum );
1592     painCave.isFatal = 1;
1593     simError();
1594     }
1595    
1596     info.k2 = atof( the_token );
1597    
1598     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1599     sprintf( painCave.errMsg,
1600     "Error parseing BendTypes: line %d\n", lineNum );
1601     painCave.isFatal = 1;
1602     simError();
1603     }
1604    
1605     info.k3 = atof( the_token );
1606    
1607     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1608     sprintf( painCave.errMsg,
1609     "Error parseing BendTypes: line %d\n", lineNum );
1610     painCave.isFatal = 1;
1611     simError();
1612     }
1613    
1614     info.t0 = atof( the_token );
1615     }
1616    
1617     else{
1618     sprintf( painCave.errMsg,
1619     "Unknown TraPPE_Ex bend type \"%s\" at line %d\n",
1620     info.type,
1621     lineNum );
1622     painCave.isFatal = 1;
1623     simError();
1624     }
1625    
1626     return 1;
1627     }
1628     else return 0;
1629     }
1630    
1631     int parseTorsion( char *lineBuffer, int lineNum, torsionStruct &info ){
1632    
1633     char* the_token;
1634    
1635     the_token = strtok( lineBuffer, " \n\t,;" );
1636     if( the_token != NULL ){
1637    
1638     strcpy( info.nameA, the_token );
1639    
1640     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1641     sprintf( painCave.errMsg,
1642     "Error parseing TorsionTypes: line %d\n", lineNum );
1643     painCave.isFatal = 1;
1644     simError();
1645     }
1646    
1647     strcpy( info.nameB, the_token );
1648    
1649     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1650     sprintf( painCave.errMsg,
1651     "Error parseing TorsionTypes: line %d\n", lineNum );
1652     painCave.isFatal = 1;
1653     simError();
1654     }
1655    
1656     strcpy( info.nameC, the_token );
1657    
1658     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1659     sprintf( painCave.errMsg,
1660     "Error parseing TorsionTypes: line %d\n", lineNum );
1661     painCave.isFatal = 1;
1662     simError();
1663     }
1664    
1665     strcpy( info.nameD, the_token );
1666    
1667     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1668     sprintf( painCave.errMsg,
1669     "Error parseing TorsionTypes: line %d\n", lineNum );
1670     painCave.isFatal = 1;
1671     simError();
1672     }
1673    
1674     strcpy( info.type, the_token );
1675    
1676     if( !strcmp( info.type, "cubic" ) ){
1677     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1678     sprintf( painCave.errMsg,
1679     "Error parseing TorsionTypes: line %d\n", lineNum );
1680     painCave.isFatal = 1;
1681     simError();
1682     }
1683    
1684     info.k1 = atof( the_token );
1685    
1686     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1687     sprintf( painCave.errMsg,
1688     "Error parseing TorsionTypes: line %d\n", lineNum );
1689     painCave.isFatal = 1;
1690     simError();
1691     }
1692    
1693     info.k2 = atof( the_token );
1694    
1695     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1696     sprintf( painCave.errMsg,
1697     "Error parseing TorsionTypes: line %d\n", lineNum );
1698     painCave.isFatal = 1;
1699     simError();
1700     }
1701    
1702     info.k3 = atof( the_token );
1703    
1704     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1705     sprintf( painCave.errMsg,
1706     "Error parseing TorsionTypes: line %d\n", lineNum );
1707     painCave.isFatal = 1;
1708     simError();
1709     }
1710    
1711     info.k4 = atof( the_token );
1712    
1713     }
1714    
1715     else{
1716     sprintf( painCave.errMsg,
1717     "Unknown TraPPE_Ex torsion type \"%s\" at line %d\n",
1718     info.type,
1719     lineNum );
1720     painCave.isFatal = 1;
1721     simError();
1722     }
1723    
1724     return 1;
1725     }
1726    
1727     else return 0;
1728     }