ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/WATER.cpp
Revision: 1882
Committed: Fri Dec 10 18:41:05 2004 UTC (19 years, 7 months ago) by tim
File size: 27343 byte(s)
Log Message:
fix another bug in WATER.cpp

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