ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/WATER.cpp
Revision: 1113
Committed: Thu Apr 15 16:18:26 2004 UTC (20 years, 2 months ago) by tim
File size: 26715 byte(s)
Log Message: 
fix whole bunch of bugs :-)

File Contents

# User Rev Content
1 chrisfen 976 #include <stdlib.h>
2     #include <stdio.h>
3     #include <string.h>
4    
5     #include <iostream>
6     using namespace std;
7    
8     #ifdef IS_MPI
9     #include <mpi.h>
10     #endif //is_mpi
11    
12     #include "ForceFields.hpp"
13     #include "SRI.hpp"
14     #include "simError.h"
15    
16     #include "fortranWrappers.hpp"
17    
18     #ifdef IS_MPI
19     #include "mpiForceField.h"
20     #endif // is_mpi
21    
22    
23    
24     namespace WATER_NS{
25    
26 gezelter 989 // Declare the structures that will be passed by the parser and MPI
27 chrisfen 976
28     typedef struct{
29     char name[15];
30     double mass;
31     double epslon;
32     double sigma;
33     double charge;
34     int isDirectional;
35     int isLJ;
36     int isCharge;
37     int ident;
38     int last; // 0 -> default
39     // 1 -> in MPI: tells nodes to stop listening
40     } atomStruct;
41    
42     typedef struct{
43     char name[15];
44 gezelter 989 double Ixx;
45     double Iyy;
46     double Izz;
47 chrisfen 976 double dipole;
48     double w0;
49     double v0;
50     double v0p;
51     double rl;
52     double ru;
53     double rlp;
54     double rup;
55     int isDipole;
56     int isSticky;
57     int last; // 0 -> default
58     // 1 -> in MPI: tells nodes to stop listening
59     } directionalStruct;
60    
61     int parseAtom( char *lineBuffer, int lineNum, atomStruct &info );
62     int parseDirectional( char *lineBuffer, int lineNum, directionalStruct &info );
63    
64    
65     #ifdef IS_MPI
66    
67     MPI_Datatype mpiAtomStructType;
68     MPI_Datatype mpiDirectionalStructType;
69    
70     #endif
71    
72     class LinkedAtomType {
73     public:
74     LinkedAtomType(){
75     next = NULL;
76     name[0] = '\0';
77     }
78     ~LinkedAtomType(){ if( next != NULL ) delete next; }
79    
80     LinkedAtomType* find(char* key){
81     if( !strcmp(name, key) ) return this;
82     if( next != NULL ) return next->find(key);
83     return NULL;
84     }
85    
86    
87     void add( atomStruct &info ){
88    
89     // check for duplicates
90    
91     if( !strcmp( info.name, name ) ){
92     sprintf( painCave.errMsg,
93     "Duplicate WATER atom type \"%s\" found in "
94     "the WATER param file./n",
95     name );
96     painCave.isFatal = 1;
97     simError();
98     }
99    
100     if( next != NULL ) next->add(info);
101     else{
102     next = new LinkedAtomType();
103     strcpy(next->name, info.name);
104     next->isDirectional = info.isDirectional;
105     next->isLJ = info.isLJ;
106     next->isCharge = info.isCharge;
107     next->mass = info.mass;
108     next->epslon = info.epslon;
109     next->sigma = info.sigma;
110     next->charge = info.charge;
111     next->ident = info.ident;
112     }
113     }
114    
115     #ifdef IS_MPI
116    
117     void duplicate( atomStruct &info ){
118     strcpy(info.name, name);
119     info.isDirectional = isDirectional;
120     info.isLJ = isLJ;
121     info.isCharge = isCharge;
122     info.mass = mass;
123     info.epslon = epslon;
124     info.sigma = sigma;
125     info.charge = charge;
126     info.ident = ident;
127     info.last = 0;
128     }
129    
130     #endif
131    
132     char name[15];
133     int isDirectional;
134     int isLJ;
135     int isCharge;
136     double mass;
137     double epslon;
138     double sigma;
139     double charge;
140     int ident;
141     LinkedAtomType* next;
142     };
143    
144     class LinkedDirectionalType {
145     public:
146     LinkedDirectionalType(){
147     next = NULL;
148     name[0] = '\0';
149     }
150     ~LinkedDirectionalType(){ if( next != NULL ) delete next; }
151    
152     LinkedDirectionalType* find(char* key){
153     if( !strcmp(name, key) ) return this;
154     if( next != NULL ) return next->find(key);
155     return NULL;
156     }
157    
158    
159     void add( directionalStruct &info ){
160    
161     // check for duplicates
162    
163     if( !strcmp( info.name, name ) ){
164     sprintf( painCave.errMsg,
165     "Duplicate WATER directional type \"%s\" found in "
166     "the WATER param file./n",
167     name );
168     painCave.isFatal = 1;
169     simError();
170     }
171    
172     if( next != NULL ) next->add(info);
173     else{
174     next = new LinkedDirectionalType();
175     strcpy(next->name, info.name);
176     next->isDipole = info.isDipole;
177     next->isSticky = info.isSticky;
178 gezelter 989 next->Ixx = info.Ixx;
179     next->Iyy = info.Iyy;
180     next->Izz = info.Izz;
181     next->dipole = info.dipole;
182 chrisfen 976 next->w0 = info.w0;
183     next->v0 = info.v0;
184     next->v0p = info.v0p;
185     next->rl = info.rl;
186     next->ru = info.ru;
187     next->rlp = info.rlp;
188     next->rup = info.rup;
189     }
190     }
191    
192     #ifdef IS_MPI
193    
194     void duplicate( directionalStruct &info ){
195     strcpy(info.name, name);
196     info.isDipole = isDipole;
197     info.isSticky = isSticky;
198 gezelter 989 info.Ixx = Ixx;
199     info.Iyy = Iyy;
200     info.Izz = Izz;
201 chrisfen 976 info.dipole = dipole;
202     info.w0 = w0;
203     info.v0 = v0;
204     info.v0p = v0p;
205     info.rl = rl;
206     info.ru = ru;
207     info.rlp = rlp;
208     info.rup = rup;
209     info.last = 0;
210     }
211    
212     #endif
213    
214     char name[15];
215     int isDipole;
216     int isSticky;
217 gezelter 989 double Ixx;
218     double Iyy;
219     double Izz;
220 chrisfen 976 double dipole;
221     double w0;
222     double v0;
223     double v0p;
224     double rl;
225     double ru;
226     double rlp;
227     double rup;
228     LinkedDirectionalType* next;
229     };
230    
231     LinkedAtomType* headAtomType;
232     LinkedAtomType* currentAtomType;
233     LinkedDirectionalType* headDirectionalType;
234     LinkedDirectionalType* currentDirectionalType;
235 gezelter 989 } // namespace
236 chrisfen 976
237     using namespace WATER_NS;
238    
239     //****************************************************************
240     // begins the actual forcefield stuff.
241     //****************************************************************
242    
243    
244     WATER::WATER(){
245    
246     char fileName[200];
247     char* ffPath_env = "FORCE_PARAM_PATH";
248     char* ffPath;
249     char temp[200];
250    
251     headAtomType = NULL;
252     currentAtomType = NULL;
253     headDirectionalType = NULL;
254     currentDirectionalType = NULL;
255    
256     // do the funtion wrapping
257     wrapMeFF( this );
258    
259     #ifdef IS_MPI
260     int i;
261    
262     // **********************************************************************
263     // Init the atomStruct mpi type
264    
265     atomStruct atomProto; // mpiPrototype
266     int atomBC[3] = {15,4,5}; // block counts
267     MPI_Aint atomDspls[3]; // displacements
268     MPI_Datatype atomMbrTypes[3]; // member mpi types
269    
270     MPI_Address(&atomProto.name, &atomDspls[0]);
271     MPI_Address(&atomProto.mass, &atomDspls[1]);
272 gezelter 989 MPI_Address(&atomProto.isDirectional, &atomDspls[2]);
273 chrisfen 976
274     atomMbrTypes[0] = MPI_CHAR;
275     atomMbrTypes[1] = MPI_DOUBLE;
276     atomMbrTypes[2] = MPI_INT;
277    
278     for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
279    
280     MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
281     MPI_Type_commit(&mpiAtomStructType);
282    
283     // ***********************************************************************
284    
285     // **********************************************************************
286     // Init the directionalStruct mpi type
287    
288     directionalStruct directionalProto; // mpiPrototype
289 gezelter 989 int directionalBC[3] = {15,11,3}; // block counts
290 chrisfen 976 MPI_Aint directionalDspls[3]; // displacements
291     MPI_Datatype directionalMbrTypes[3]; // member mpi types
292    
293     MPI_Address(&directionalProto.name, &directionalDspls[0]);
294 gezelter 989 MPI_Address(&directionalProto.Ixx, &directionalDspls[1]);
295     MPI_Address(&directionalProto.isDipole, &directionalDspls[2]);
296 chrisfen 976
297     directionalMbrTypes[0] = MPI_CHAR;
298     directionalMbrTypes[1] = MPI_DOUBLE;
299     directionalMbrTypes[2] = MPI_INT;
300    
301     for (i=2; i >= 0; i--) directionalDspls[i] -= directionalDspls[0];
302    
303     MPI_Type_struct(3, directionalBC, directionalDspls, directionalMbrTypes,
304     &mpiDirectionalStructType);
305     MPI_Type_commit(&mpiDirectionalStructType);
306    
307     // ***********************************************************************
308 chrisfen 999
309 chrisfen 976 if( worldRank == 0 ){
310     #endif
311    
312     // generate the force file name
313    
314     strcpy( fileName, "WATER.frc" );
315     // fprintf( stderr,"Trying to open %s\n", fileName );
316    
317     // attempt to open the file in the current directory first.
318    
319     frcFile = fopen( fileName, "r" );
320    
321     if( frcFile == NULL ){
322    
323     // next see if the force path enviorment variable is set
324    
325     ffPath = getenv( ffPath_env );
326     if( ffPath == NULL ) {
327     STR_DEFINE(ffPath, FRC_PATH );
328     }
329    
330    
331     strcpy( temp, ffPath );
332     strcat( temp, "/" );
333     strcat( temp, fileName );
334     strcpy( fileName, temp );
335    
336     frcFile = fopen( fileName, "r" );
337    
338     if( frcFile == NULL ){
339    
340     sprintf( painCave.errMsg,
341     "Error opening the force field parameter file: %s\n"
342     "Have you tried setting the FORCE_PARAM_PATH environment "
343     "variable?\n",
344     fileName );
345     painCave.isFatal = 1;
346     simError();
347     }
348     }
349    
350     #ifdef IS_MPI
351     }
352    
353     sprintf( checkPointMsg, "WATER file opened sucessfully." );
354     MPIcheckPoint();
355    
356     #endif // is_mpi
357     }
358    
359    
360     WATER::~WATER(){
361    
362     if( headAtomType != NULL ) delete headAtomType;
363     if( headDirectionalType != NULL ) delete headDirectionalType;
364    
365     #ifdef IS_MPI
366     if( worldRank == 0 ){
367     #endif // is_mpi
368    
369     fclose( frcFile );
370    
371     #ifdef IS_MPI
372     }
373     #endif // is_mpi
374     }
375    
376     void WATER::cleanMe( void ){
377    
378     #ifdef IS_MPI
379    
380     // keep the linked list in the mpi version
381    
382     #else // is_mpi
383    
384     // delete the linked list in the single processor version
385    
386     if( headAtomType != NULL ) delete headAtomType;
387     if( headDirectionalType != NULL ) delete headDirectionalType;
388    
389     #endif // is_mpi
390     }
391    
392    
393     void WATER::initForceField( int ljMixRule ){
394    
395     initFortran( ljMixRule, entry_plug->useReactionField );
396     }
397    
398    
399     void WATER::readParams( void ){
400    
401     int identNum;
402 chrisfen 999 int tempDirect0, tempDirect1;
403 chrisfen 976
404     atomStruct atomInfo;
405     directionalStruct directionalInfo;
406 gezelter 989 fpos_t *atomPos;
407 chrisfen 976
408     atomInfo.last = 1; // initialize last to have the last set.
409     directionalInfo.last = 1; // if things go well, last will be set to 0
410    
411     atomPos = new fpos_t;
412     bigSigma = 0.0;
413 chrisfen 999
414 chrisfen 976 #ifdef IS_MPI
415     if( worldRank == 0 ){
416     #endif
417    
418     // read in the atom types.
419    
420     headAtomType = new LinkedAtomType;
421 gezelter 989 headDirectionalType = new LinkedDirectionalType;
422 chrisfen 976
423     fastForward( "AtomTypes", "initializeAtoms" );
424 chrisfen 999
425 chrisfen 976 // we are now at the AtomTypes section.
426    
427     eof_test = fgets( readLine, sizeof(readLine), frcFile );
428     lineNum++;
429    
430    
431     // read a line, and start parsing out the atom types
432 chrisfen 999
433 chrisfen 976 if( eof_test == NULL ){
434     sprintf( painCave.errMsg,
435     "Error in reading Atoms from force file at line %d.\n",
436     lineNum );
437     painCave.isFatal = 1;
438     simError();
439     }
440    
441     identNum = 1;
442     // stop reading at end of file, or at next section
443     while( readLine[0] != '#' && eof_test != NULL ){
444    
445     // toss comment lines
446     if( readLine[0] != '!' ){
447    
448     // the parser returns 0 if the line was blank
449     if( parseAtom( readLine, lineNum, atomInfo ) ){
450     atomInfo.ident = identNum;
451     headAtomType->add( atomInfo );
452     if ( atomInfo.isDirectional ) {
453    
454     // if the atom is directional, skip to the directional section
455     // and parse directional info.
456     fgetpos( frcFile, atomPos );
457     sectionSearch( "DirectionalTypes", atomInfo.name,
458     "initializeDirectional" );
459     parseDirectional( readLine, lineNum, directionalInfo );
460 gezelter 989 headDirectionalType->add( directionalInfo );
461 chrisfen 976
462     // return to the AtomTypes section
463     fsetpos( frcFile, atomPos );
464     }
465     identNum++;
466     }
467     }
468     eof_test = fgets( readLine, sizeof(readLine), frcFile );
469     lineNum++;
470     }
471    
472     #ifdef IS_MPI
473 chrisfen 999
474 chrisfen 976 // send out the linked list to all the other processes
475    
476     sprintf( checkPointMsg,
477 chrisfen 999 "WATER atom and directional structures read successfully." );
478 chrisfen 976 MPIcheckPoint();
479 chrisfen 999 currentAtomType = headAtomType->next; //skip the first element place holder
480     currentDirectionalType = headDirectionalType->next; // same w/ directional
481    
482 chrisfen 976 while( currentAtomType != NULL ){
483     currentAtomType->duplicate( atomInfo );
484 chrisfen 999
485 chrisfen 976 sendFrcStruct( &atomInfo, mpiAtomStructType );
486    
487     sprintf( checkPointMsg,
488     "successfully sent WATER force type: \"%s\"\n",
489     atomInfo.name );
490 chrisfen 999
491     if ( atomInfo.isDirectional ){
492     // send out the directional linked list to all the other processes
493    
494     currentDirectionalType->duplicate( directionalInfo );
495     sendFrcStruct( &directionalInfo, mpiDirectionalStructType );
496    
497     sprintf( checkPointMsg,
498     "successfully sent WATER directional type: \"%s\"\n",
499     directionalInfo.name );
500     }
501    
502 chrisfen 976 MPIcheckPoint();
503 chrisfen 999 tempDirect0 = atomInfo.isDirectional;
504 chrisfen 976 currentAtomType = currentAtomType->next;
505 chrisfen 999 if( tempDirect0 )
506     currentDirectionalType = currentDirectionalType->next;
507 chrisfen 976 }
508 chrisfen 999
509 chrisfen 976 atomInfo.last = 1;
510     sendFrcStruct( &atomInfo, mpiAtomStructType );
511 chrisfen 999 directionalInfo.last = 1;
512     if ( atomInfo.isDirectional )
513 chrisfen 976 sendFrcStruct( &directionalInfo, mpiDirectionalStructType );
514     }
515 chrisfen 999
516 chrisfen 976 else{
517     // listen for node 0 to send out the force params
518    
519     MPIcheckPoint();
520 gezelter 989
521 chrisfen 976 headAtomType = new LinkedAtomType;
522 chrisfen 999 headDirectionalType = new LinkedDirectionalType;
523 chrisfen 976 receiveFrcStruct( &atomInfo, mpiAtomStructType );
524 chrisfen 999
525     if ( atomInfo.isDirectional )
526     receiveFrcStruct( &directionalInfo, mpiDirectionalStructType );
527    
528 chrisfen 976 while( !atomInfo.last ){
529 chrisfen 999
530 chrisfen 976 headAtomType->add( atomInfo );
531    
532     MPIcheckPoint();
533 chrisfen 999
534 chrisfen 976 receiveFrcStruct( &atomInfo, mpiAtomStructType );
535    
536 chrisfen 999 if( atomInfo.isDirectional ){
537 chrisfen 976 headDirectionalType->add( directionalInfo );
538 chrisfen 999
539 chrisfen 976 receiveFrcStruct( &directionalInfo, mpiDirectionalStructType );
540     }
541     }
542     }
543    
544     #endif // is_mpi
545    
546     // call new A_types in fortran
547    
548     int isError;
549    
550     // dummy variables
551     int isGB = 0;
552     int isEAM = 0;
553 gezelter 989 int notDipole = 0;
554     int notSSD = 0;
555     double noDipMoment = 0.0;
556    
557    
558     currentAtomType = headAtomType->next;
559     currentDirectionalType = headDirectionalType->next;
560    
561 chrisfen 976 while( currentAtomType != NULL ){
562    
563     if( currentAtomType->isLJ ) entry_plug->useLJ = 1;
564 gezelter 989 if( currentAtomType->isCharge ) entry_plug->useCharges = 1;
565 chrisfen 976 if( currentAtomType->isDirectional ){
566 gezelter 989 // only directional atoms can have dipoles or be sticky
567     if ( currentDirectionalType->isDipole ) entry_plug->useDipoles = 1;
568     if ( currentDirectionalType->isSticky ) {
569 chrisfen 976 entry_plug->useSticky = 1;
570 gezelter 989 set_sticky_params( &(currentDirectionalType->w0),
571     &(currentDirectionalType->v0),
572     &(currentDirectionalType->v0p),
573     &(currentDirectionalType->rl),
574     &(currentDirectionalType->ru),
575     &(currentDirectionalType->rlp),
576     &(currentDirectionalType->rup));
577 chrisfen 976 }
578 gezelter 989 if( currentAtomType->name[0] != '\0' ){
579     isError = 0;
580     makeAtype( &(currentAtomType->ident),
581     &(currentAtomType->isLJ),
582     &(currentDirectionalType->isSticky),
583     &(currentDirectionalType->isDipole),
584     &isGB,
585     &isEAM,
586     &(currentAtomType->isCharge),
587     &(currentAtomType->epslon),
588     &(currentAtomType->sigma),
589     &(currentAtomType->charge),
590     &(currentDirectionalType->dipole),
591     &isError );
592     if( isError ){
593     sprintf( painCave.errMsg,
594     "Error initializing the \"%s\" atom type in fortran\n",
595     currentAtomType->name );
596     painCave.isFatal = 1;
597     simError();
598     }
599     }
600     currentDirectionalType->next;
601 chrisfen 976 }
602 gezelter 989
603     else {
604     // use all dummy variables if this is not a directional atom
605     if( currentAtomType->name[0] != '\0' ){
606     isError = 0;
607     makeAtype( &(currentAtomType->ident),
608     &(currentAtomType->isLJ),
609     &notSSD,
610     &notDipole,
611     &isGB,
612     &isEAM,
613     &(currentAtomType->isCharge),
614     &(currentAtomType->epslon),
615     &(currentAtomType->sigma),
616     &(currentAtomType->charge),
617     &noDipMoment,
618     &isError );
619     if( isError ){
620     sprintf( painCave.errMsg,
621     "Error initializing the \"%s\" atom type in fortran\n",
622     currentAtomType->name );
623     painCave.isFatal = 1;
624     simError();
625     }
626 chrisfen 976 }
627     }
628     currentAtomType = currentAtomType->next;
629     }
630    
631     #ifdef IS_MPI
632     sprintf( checkPointMsg,
633 chrisfen 999 "WATER atom and directional structures successfully"
634     "sent to fortran\n" );
635 chrisfen 976 MPIcheckPoint();
636     #endif // is_mpi
637    
638     }
639    
640    
641     void WATER::initializeAtoms( int nAtoms, Atom** the_atoms ){
642    
643 chrisfen 999 int i,j,k;
644 chrisfen 976
645     // initialize the atoms
646 gezelter 989 DirectionalAtom* dAtom;
647 gezelter 1097 double ji[3];
648 gezelter 989 double inertialMat[3][3];
649 chrisfen 976
650     for( i=0; i<nAtoms; i++ ){
651     currentAtomType = headAtomType->find( the_atoms[i]->getType() );
652     if( currentAtomType == NULL ){
653     sprintf( painCave.errMsg,
654     "AtomType error, %s not found in force file.\n",
655     the_atoms[i]->getType() );
656     painCave.isFatal = 1;
657     simError();
658     }
659     the_atoms[i]->setMass( currentAtomType->mass );
660     the_atoms[i]->setIdent( currentAtomType->ident );
661 chrisfen 999
662 gezelter 989 if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
663 chrisfen 976
664 tim 1113 the_atoms[i]->setHasCharge(currentAtomType->isCharge);
665    
666 gezelter 989 if( currentAtomType->isDirectional ){
667     currentDirectionalType =
668     headDirectionalType->find( the_atoms[i]->getType() );
669     if( currentDirectionalType == NULL ){
670     sprintf( painCave.errMsg,
671     "DirectionalType error, %s not found in force file.\n",
672     the_atoms[i]->getType() );
673     painCave.isFatal = 1;
674     simError();
675     }
676 chrisfen 976
677 gezelter 989 // zero out the moments of inertia matrix
678 chrisfen 999 for( j=0; j<3; j++ )
679     for( k=0; k<3; k++ )
680     inertialMat[j][k] = 0.0;
681 chrisfen 976
682 gezelter 989 // load the force file moments of inertia
683     inertialMat[0][0] = currentDirectionalType->Ixx;
684     inertialMat[1][1] = currentDirectionalType->Iyy;
685     inertialMat[2][2] = currentDirectionalType->Izz;
686 chrisfen 999
687 gezelter 989 dAtom = (DirectionalAtom *) the_atoms[i];
688     dAtom->setHasDipole( currentDirectionalType->isDipole );
689 chrisfen 999
690 gezelter 1097 ji[0] = 0.0;
691     ji[1] = 0.0;
692     ji[2] = 0.0;
693     dAtom->setJ( ji );
694 gezelter 989 dAtom->setI( inertialMat );
695 chrisfen 999
696 gezelter 989 entry_plug->n_dipoles++;
697     }
698 chrisfen 976 }
699     }
700    
701     void WATER::initializeBonds( int nBonds, Bond** BondArray,
702     bond_pair* the_bonds ){
703    
704     if( nBonds ){
705     sprintf( painCave.errMsg,
706     "WATER does not support bonds.\n" );
707     painCave.isFatal = 1;
708     simError();
709     }
710     }
711    
712     void WATER::initializeBends( int nBends, Bend** bendArray,
713     bend_set* the_bends ){
714    
715     if( nBends ){
716     sprintf( painCave.errMsg,
717     "WATER does not support bends.\n" );
718     painCave.isFatal = 1;
719     simError();
720     }
721     }
722    
723     void WATER::initializeTorsions( int nTorsions, Torsion** torsionArray,
724     torsion_set* the_torsions ){
725    
726     if( nTorsions ){
727     sprintf( painCave.errMsg,
728     "WATER does not support torsions.\n" );
729     painCave.isFatal = 1;
730     simError();
731     }
732     }
733    
734     void WATER::fastForward( char* stopText, char* searchOwner ){
735    
736     int foundText = 0;
737     char* the_token;
738    
739     rewind( frcFile );
740     lineNum = 0;
741    
742     eof_test = fgets( readLine, sizeof(readLine), frcFile );
743     lineNum++;
744     if( eof_test == NULL ){
745     sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
746     " file is empty.\n",
747     searchOwner );
748     painCave.isFatal = 1;
749     simError();
750     }
751    
752    
753     while( !foundText ){
754     while( eof_test != NULL && readLine[0] != '#' ){
755     eof_test = fgets( readLine, sizeof(readLine), frcFile );
756     lineNum++;
757     }
758     if( eof_test == NULL ){
759     sprintf( painCave.errMsg,
760     "Error fast forwarding force file for %s at "
761     "line %d: file ended unexpectedly.\n",
762     searchOwner,
763     lineNum );
764     painCave.isFatal = 1;
765     simError();
766     }
767    
768     the_token = strtok( readLine, " ,;\t#\n" );
769     foundText = !strcmp( stopText, the_token );
770    
771     if( !foundText ){
772     eof_test = fgets( readLine, sizeof(readLine), frcFile );
773     lineNum++;
774    
775     if( eof_test == NULL ){
776     sprintf( painCave.errMsg,
777     "Error fast forwarding force file for %s at "
778     "line %d: file ended unexpectedly.\n",
779     searchOwner,
780     lineNum );
781     painCave.isFatal = 1;
782     simError();
783     }
784     }
785     }
786     }
787    
788     void WATER::sectionSearch( char* secHead, char* stopText, char* searchOwner ){
789    
790     int foundSection = 0;
791     int foundText = 0;
792     char* the_token;
793 gezelter 989 fpos_t *tempPos;
794 chrisfen 976
795     rewind( frcFile );
796     lineNum = 0;
797 gezelter 989 tempPos = new fpos_t;
798 chrisfen 976
799     eof_test = fgets( readLine, sizeof(readLine), frcFile );
800     lineNum++;
801     if( eof_test == NULL ){
802     sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
803     " file is empty.\n",
804     searchOwner );
805     painCave.isFatal = 1;
806     simError();
807     }
808    
809     while( !foundSection ){
810     while( eof_test != NULL && readLine[0] != '#' ){
811     eof_test = fgets( readLine, sizeof(readLine), frcFile );
812     lineNum++;
813     }
814     if( eof_test == NULL ){
815     sprintf( painCave.errMsg,
816     "Error fast forwarding force file for %s at "
817     "line %d: file ended unexpectedly.\n",
818     searchOwner,
819     lineNum );
820     painCave.isFatal = 1;
821     simError();
822     }
823     the_token = strtok( readLine, " ,;\t#\n" );
824     foundSection = !strcmp( secHead, the_token );
825    
826     if( !foundSection ){
827     eof_test = fgets( readLine, sizeof(readLine), frcFile );
828     lineNum++;
829    
830     if( eof_test == NULL ){
831     sprintf( painCave.errMsg,
832     "Error section searching force file for %s at "
833     "line %d: file ended unexpectedly.\n",
834     searchOwner,
835     lineNum );
836     painCave.isFatal = 1;
837     simError();
838     }
839     }
840 gezelter 989 }
841 chrisfen 976
842     while( !foundText ){
843 gezelter 989
844     fgetpos( frcFile, tempPos );
845     eof_test = fgets( readLine, sizeof(readLine), frcFile );
846     lineNum++;
847    
848     if( eof_test == NULL ){
849     sprintf( painCave.errMsg,
850     "Error fast forwarding force file for %s at "
851     "line %d: file ended unexpectedly.\n",
852     searchOwner,
853     lineNum );
854     painCave.isFatal = 1;
855 chrisfen 976 simError();
856 gezelter 989 }
857 chrisfen 976
858     the_token = strtok( readLine, " ,;\t#\n" );
859 gezelter 989 if( the_token != NULL ){
860     foundText = !strcmp( stopText, the_token );
861     }
862 chrisfen 976 }
863 gezelter 989 fsetpos( frcFile, tempPos );
864     eof_test = fgets( readLine, sizeof(readLine), frcFile );
865     lineNum++;
866 chrisfen 976 }
867    
868    
869     int WATER_NS::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
870    
871     char* the_token;
872    
873     the_token = strtok( lineBuffer, " \n\t,;" );
874     if( the_token != NULL ){
875    
876     strcpy( info.name, the_token );
877    
878     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
879     sprintf( painCave.errMsg,
880     "Error parseing AtomTypes: line %d\n", lineNum );
881     painCave.isFatal = 1;
882     simError();
883     }
884    
885     info.isDirectional = atoi( the_token );
886    
887     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
888     sprintf( painCave.errMsg,
889     "Error parseing AtomTypes: line %d\n", lineNum );
890     painCave.isFatal = 1;
891     simError();
892     }
893    
894     info.isLJ = atoi( the_token );
895    
896     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
897     sprintf( painCave.errMsg,
898     "Error parseing AtomTypes: line %d\n", lineNum );
899     painCave.isFatal = 1;
900     simError();
901     }
902    
903     info.isCharge = atoi( the_token );
904    
905     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
906     sprintf( painCave.errMsg,
907     "Error parseing AtomTypes: line %d\n", lineNum );
908     painCave.isFatal = 1;
909     simError();
910     }
911    
912     info.mass = atof( the_token );
913    
914     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
915     sprintf( painCave.errMsg,
916     "Error parseing AtomTypes: line %d\n", lineNum );
917     painCave.isFatal = 1;
918     simError();
919     }
920    
921     info.epslon = atof( the_token );
922    
923     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
924     sprintf( painCave.errMsg,
925     "Error parseing AtomTypes: line %d\n", lineNum );
926     painCave.isFatal = 1;
927     simError();
928     }
929    
930     info.sigma = atof( the_token );
931    
932     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
933     sprintf( painCave.errMsg,
934     "Error parseing AtomTypes: line %d\n", lineNum );
935     painCave.isFatal = 1;
936     simError();
937     }
938    
939     info.charge = atof( the_token );
940    
941     return 1;
942     }
943     else return 0;
944     }
945    
946     int WATER_NS::parseDirectional( char *lineBuffer, int lineNum, directionalStruct &info ){
947    
948     char* the_token;
949    
950     the_token = strtok( lineBuffer, " \n\t,;" );
951     if( the_token != NULL ){
952    
953     strcpy( info.name, the_token );
954    
955     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
956     sprintf( painCave.errMsg,
957     "Error parseing DirectionalTypes: line %d\n", lineNum );
958     painCave.isFatal = 1;
959     simError();
960     }
961    
962     info.isDipole = atoi( the_token );
963    
964     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
965     sprintf( painCave.errMsg,
966     "Error parseing DirectionalTypes: line %d\n", lineNum );
967     painCave.isFatal = 1;
968     simError();
969     }
970    
971     info.isSticky = atoi( the_token );
972    
973     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
974     sprintf( painCave.errMsg,
975     "Error parseing DirectionalTypes: line %d\n", lineNum );
976     painCave.isFatal = 1;
977     simError();
978     }
979    
980 gezelter 989 info.Ixx = atof( the_token );
981 chrisfen 976
982     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
983     sprintf( painCave.errMsg,
984     "Error parseing DirectionalTypes: line %d\n", lineNum );
985     painCave.isFatal = 1;
986     simError();
987     }
988    
989 gezelter 989 info.Iyy = atof( the_token );
990 chrisfen 976
991     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
992     sprintf( painCave.errMsg,
993     "Error parseing DirectionalTypes: line %d\n", lineNum );
994     painCave.isFatal = 1;
995     simError();
996     }
997    
998 gezelter 989 info.Izz = atof( the_token );
999 chrisfen 976
1000     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1001     sprintf( painCave.errMsg,
1002     "Error parseing DirectionalTypes: line %d\n", lineNum );
1003     painCave.isFatal = 1;
1004     simError();
1005     }
1006    
1007    
1008     info.dipole = atof( the_token );
1009    
1010     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1011     sprintf( painCave.errMsg,
1012     "Error parseing DirectionalTypes: line %d\n", lineNum );
1013     painCave.isFatal = 1;
1014     simError();
1015     }
1016    
1017     info.w0 = atof( the_token );
1018    
1019     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1020     sprintf( painCave.errMsg,
1021     "Error parseing DirectionalTypes: line %d\n", lineNum );
1022     painCave.isFatal = 1;
1023     simError();
1024     }
1025    
1026     info.v0 = atof( the_token );
1027    
1028     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1029     sprintf( painCave.errMsg,
1030     "Error parseing DirectionalTypes: line %d\n", lineNum );
1031     painCave.isFatal = 1;
1032     simError();
1033     }
1034    
1035     info.v0p = atof( the_token );
1036    
1037     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1038     sprintf( painCave.errMsg,
1039     "Error parseing DirectionalTypes: line %d\n", lineNum );
1040     painCave.isFatal = 1;
1041     simError();
1042     }
1043    
1044     info.rl = atof( the_token );
1045    
1046     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1047     sprintf( painCave.errMsg,
1048     "Error parseing DirectionalTypes: line %d\n", lineNum );
1049     painCave.isFatal = 1;
1050     simError();
1051     }
1052    
1053     info.ru = atof( the_token );
1054    
1055     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1056     sprintf( painCave.errMsg,
1057     "Error parseing DirectionalTypes: line %d\n", lineNum );
1058     painCave.isFatal = 1;
1059     simError();
1060     }
1061    
1062     info.rlp = atof( the_token );
1063    
1064     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1065     sprintf( painCave.errMsg,
1066     "Error parseing DirectionalTypes: line %d\n", lineNum );
1067     painCave.isFatal = 1;
1068     simError();
1069     }
1070    
1071     info.rup = atof( the_token );
1072    
1073     return 1;
1074     }
1075     else return 0;
1076     }