ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/WATER.cpp
Revision: 1636
Committed: Fri Oct 22 22:54:01 2004 UTC (19 years, 8 months ago) by chrisfen
File size: 27594 byte(s)
Log Message:
fixey fixey the breakey breakey

File Contents

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