ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DUFF.cpp
Revision: 559
Committed: Thu Jun 19 22:02:44 2003 UTC (21 years ago) by mmeineke
File size: 40852 byte(s)
Log Message: 
slowly converting to new integrator and forcefield names.

File Contents

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