ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/interface_implementation/TraPPE_ExFF.cpp
Revision: 137
Committed: Wed Oct 16 20:02:05 2002 UTC (21 years, 11 months ago) by chuckv
File size: 22821 byte(s)
Log Message: 
*** empty log message ***

File Contents

# User Rev Content
1 mmeineke 10 #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    
11    
12     TraPPE_ExFF::TraPPE_ExFF(){
13    
14     char fileName[200];
15     char* ffPath_env = "FORCE_PARAM_PATH";
16     char* ffPath;
17     char temp[200];
18    
19     // generate the force file name
20    
21 chuckv 137
22 mmeineke 10 strcpy( fileName, "TraPPE_Ex.frc" );
23 chuckv 137 fprintf( stderr,"Trying to open %s\n", fileName );
24    
25 mmeineke 10 // attempt to open the file in the current directory first.
26    
27     frcFile = fopen( fileName, "r" );
28    
29     if( frcFile == NULL ){
30    
31     // next see if the force path enviorment variable is set
32    
33     ffPath = getenv( ffPath_env );
34 chuckv 137 if( ffPath == NULL ) {
35     fprintf( stderr,
36     "Error opening the force field parameter file: %s\n"
37     "Have you tried setting the FORCE_PARAM_PATH environment "
38     "vairable?\n",
39     fileName );
40     exit( 8 );
41     }
42    
43    
44 mmeineke 10 strcpy( temp, ffPath );
45     strcat( temp, "/" );
46     strcat( temp, fileName );
47     strcpy( fileName, temp );
48    
49     frcFile = fopen( fileName, "r" );
50    
51     if( frcFile == NULL ){
52    
53     fprintf( stderr,
54     "Error opening the force field parameter file: %s\n"
55     "Have you tried setting the FORCE_PARAM_PATH environment "
56     "vairable?\n",
57     fileName );
58     exit( 8 );
59     }
60     }
61     }
62    
63     TraPPE_ExFF::~TraPPE_ExFF(){
64    
65     fclose( frcFile );
66     }
67    
68     void TraPPE_ExFF::initializeAtoms( void ){
69    
70     class LinkedType {
71     public:
72     LinkedType(){
73     next = NULL;
74     name[0] = '\0';
75     }
76     ~LinkedType(){ if( next != NULL ) delete next; }
77    
78     LinkedType* find(char* key){
79     if( !strcmp(name, key) ) return this;
80     if( next != NULL ) return next->find(key);
81     return NULL;
82     }
83    
84     char name[15];
85     int isDipole;
86     double mass;
87     double epslon;
88     double sigma;
89     double dipole;
90     LinkedType* next;
91     };
92    
93     LinkedType* headAtomType;
94     LinkedType* currentAtomType;
95     LinkedType* tempAtomType;
96    
97     char readLine[500];
98     char* the_token;
99     char* eof_test;
100     int foundAtom = 0;
101     int lineNum = 0;
102     int i;
103    
104     //////////////////////////////////////////////////
105     // a quick water fix
106    
107     double waterI[3][3];
108     waterI[0][0] = 1.76958347772500;
109     waterI[0][1] = 0.0;
110     waterI[0][2] = 0.0;
111    
112     waterI[1][0] = 0.0;
113     waterI[1][1] = 0.614537057924513;
114     waterI[1][2] = 0.0;
115    
116     waterI[2][0] = 0.0;
117     waterI[2][1] = 0.0;
118     waterI[2][2] = 1.15504641980049;
119    
120    
121     double headI[3][3];
122     headI[0][0] = 1125;
123     headI[0][1] = 0.0;
124     headI[0][2] = 0.0;
125    
126     headI[1][0] = 0.0;
127     headI[1][1] = 1125;
128     headI[1][2] = 0.0;
129    
130     headI[2][0] = 0.0;
131     headI[2][1] = 0.0;
132     headI[2][2] = 250;
133    
134    
135    
136     //////////////////////////////////////////////////
137    
138     Atom** the_atoms;
139     int nAtoms;
140     the_atoms = entry_plug->atoms;
141     nAtoms = entry_plug->n_atoms;
142    
143     // read in the atom types.
144    
145     rewind( frcFile );
146     headAtomType = new LinkedType;
147     currentAtomType = headAtomType;
148    
149     eof_test = fgets( readLine, sizeof(readLine), frcFile );
150     lineNum++;
151     if( eof_test == NULL ){
152     fprintf( stderr, "Error in reading Atoms from force file.\n" );
153     exit(8);
154     }
155    
156    
157     while( !foundAtom ){
158     while( eof_test != NULL && readLine[0] != '#' ){
159     eof_test = fgets( readLine, sizeof(readLine), frcFile );
160     lineNum++;
161     }
162     if( eof_test == NULL ){
163     fprintf( stderr,
164     "Error in reading Atoms from force file at line %d.\n",
165     lineNum );
166     exit(8);
167     }
168    
169     the_token = strtok( readLine, " ,;\t#\n" );
170     foundAtom = !strcmp( "AtomTypes", the_token );
171    
172     if( !foundAtom ){
173     eof_test = fgets( readLine, sizeof(readLine), frcFile );
174     lineNum++;
175    
176     if( eof_test == NULL ){
177     fprintf( stderr,
178     "Error in reading Atoms from force file at line %d.\n",
179     lineNum );
180     exit(8);
181     }
182     }
183     }
184    
185     // we are now at the AtomTypes section.
186    
187     eof_test = fgets( readLine, sizeof(readLine), frcFile );
188     lineNum++;
189    
190     if( eof_test == NULL ){
191     fprintf( stderr,
192     "Error in reading Atoms from force file at line %d.\n",
193     lineNum );
194     exit(8);
195     }
196    
197     while( readLine[0] != '#' && eof_test != NULL ){
198    
199     if( readLine[0] != '!' ){
200    
201     the_token = strtok( readLine, " \n\t,;" );
202     if( the_token != NULL ){
203    
204     strcpy( currentAtomType->name, the_token );
205    
206     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
207     fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
208     exit(8);
209     }
210    
211     sscanf( the_token, "%d", &currentAtomType->isDipole );
212    
213     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
214     fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
215     exit(8);
216     }
217    
218     sscanf( the_token, "%lf", &currentAtomType->mass );
219    
220     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
221     fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
222     exit(8);
223     }
224    
225     sscanf( the_token, "%lf", &currentAtomType->epslon );
226    
227     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
228     fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
229     exit(8);
230     }
231    
232     sscanf( the_token, "%lf", &currentAtomType->sigma );
233    
234     if( currentAtomType->isDipole ){
235     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
236     fprintf( stderr, "Error parseing AtomTypes: line %d\n", lineNum );
237     exit(8);
238     }
239    
240     sscanf( the_token, "%lf", &currentAtomType->dipole );
241     }
242     }
243     }
244    
245     tempAtomType = new LinkedType;
246     currentAtomType->next = tempAtomType;
247     currentAtomType = tempAtomType;
248    
249     eof_test = fgets( readLine, sizeof(readLine), frcFile );
250     lineNum++;
251     }
252    
253    
254     // initialize the atoms
255    
256     DirectionalAtom* dAtom;
257    
258     for( i=0; i<nAtoms; i++ ){
259    
260     currentAtomType = headAtomType->find( the_atoms[i]->getType() );
261     if( currentAtomType == NULL ){
262     fprintf( stderr, "AtomType error, %s not found in force file.\n",
263     the_atoms[i]->getType() );
264     exit(8);
265     }
266    
267     the_atoms[i]->setMass( currentAtomType->mass );
268     the_atoms[i]->setEpslon( currentAtomType->epslon );
269     the_atoms[i]->setSigma( currentAtomType->sigma );
270     the_atoms[i]->setLJ();
271    
272     if( currentAtomType->isDipole ){
273     if( the_atoms[i]->isDirectional() ){
274    
275     dAtom = (DirectionalAtom *) the_atoms[i];
276     dAtom->setMu( currentAtomType->dipole );
277     dAtom->setHasDipole( 1 );
278     dAtom->setJx( 0.0 );
279     dAtom->setJy( 0.0 );
280     dAtom->setJz( 0.0 );
281    
282     if(!strcmp("SSD",the_atoms[i]->getType())){
283     dAtom->setI( waterI );
284     dAtom->setSSD( 1 );
285     }
286     else if(!strcmp("HEAD",the_atoms[i]->getType())){
287     dAtom->setI( headI );
288     dAtom->setSSD( 0 );
289     }
290     else{
291     fprintf(stderr,
292     "AtmType error, %s does not have a moment of inertia set.\n",
293     the_atoms[i]->getType() );
294     exit(8);
295     }
296     entry_plug->n_dipoles++;
297     }
298     else{
299     std::cerr
300     << "TraPPE_ExFF error: Atom \""
301     << currentAtomType->name << "\" is a dipole, yet no standard"
302     << " orientation was specifed in the BASS file.\n";
303     exit(8);
304     }
305     }
306     else{
307     if( the_atoms[i]->isDirectional() ){
308     std::cerr
309     << "TraPPE_ExFF error: Atom \""
310     << currentAtomType->name << "\" was given a standard orientation"
311     << " in the BASS file, yet it is not a dipole.\n";
312     exit(8);
313     }
314     }
315     }
316    
317    
318     // clean up the memory
319    
320     delete headAtomType;
321     }
322    
323     void TraPPE_ExFF::initializeBonds( bond_pair* the_bonds ){
324    
325     class LinkedType {
326     public:
327     LinkedType(){
328     next = NULL;
329     nameA[0] = '\0';
330     nameB[0] = '\0';
331     type[0] = '\0';
332     }
333     ~LinkedType(){ if( next != NULL ) delete next; }
334    
335     LinkedType* find(char* key1, char* key2){
336     if( !strcmp(nameA, key1 ) && !strcmp( nameB, key2 ) ) return this;
337     if( !strcmp(nameA, key2 ) && !strcmp( nameB, key1 ) ) return this;
338     if( next != NULL ) return next->find(key1, key2);
339     return NULL;
340     }
341    
342     char nameA[15];
343     char nameB[15];
344     char type[30];
345     double d0;
346    
347     LinkedType* next;
348     };
349    
350    
351    
352     LinkedType* headBondType;
353     LinkedType* currentBondType;
354     LinkedType* tempBondType;
355    
356     char readLine[500];
357     char* the_token;
358     char* eof_test;
359     int foundBond = 0;
360     int lineNum = 0;
361     int i, a, b;
362     char* atomA;
363     char* atomB;
364    
365     SRI **the_sris;
366     Atom** the_atoms;
367     int nBonds;
368     the_sris = entry_plug->sr_interactions;
369     the_atoms = entry_plug->atoms;
370     nBonds = entry_plug->n_bonds;
371    
372     // read in the bond types.
373    
374     rewind( frcFile );
375     headBondType = new LinkedType;
376     currentBondType = headBondType;
377    
378     eof_test = fgets( readLine, sizeof(readLine), frcFile );
379     lineNum++;
380     if( eof_test == NULL ){
381     fprintf( stderr, "Error in reading Bonds from force file.\n" );
382     exit(8);
383     }
384    
385    
386     while( !foundBond ){
387     while( eof_test != NULL && readLine[0] != '#' ){
388     eof_test = fgets( readLine, sizeof(readLine), frcFile );
389     lineNum++;
390     }
391     if( eof_test == NULL ){
392     fprintf( stderr,
393     "Error in reading Bonds from force file at line %d.\n",
394     lineNum );
395     exit(8);
396     }
397    
398     the_token = strtok( readLine, " ,;\t#\n" );
399     foundBond = !strcmp( "BondTypes", the_token );
400    
401     if( !foundBond ){
402     eof_test = fgets( readLine, sizeof(readLine), frcFile );
403     lineNum++;
404    
405     if( eof_test == NULL ){
406     fprintf( stderr,
407     "Error in reading Bonds from force file at line %d.\n",
408     lineNum );
409     exit(8);
410     }
411     }
412     }
413    
414     // we are now at the BondTypes section.
415    
416     eof_test = fgets( readLine, sizeof(readLine), frcFile );
417     lineNum++;
418    
419     if( eof_test == NULL ){
420     fprintf( stderr,
421     "Error in reading Bonds from force file at line %d.\n",
422     lineNum );
423     exit(8);
424     }
425    
426     while( readLine[0] != '#' && eof_test != NULL ){
427    
428     if( readLine[0] != '!' ){
429    
430     the_token = strtok( readLine, " \n\t,;" );
431     if( the_token != NULL ){
432    
433     strcpy( currentBondType->nameA, the_token );
434    
435     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
436     fprintf( stderr, "Error parseing BondTypes: line %d\n", lineNum );
437     exit(8);
438     }
439    
440     strcpy( currentBondType->nameB, the_token );
441    
442     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
443     fprintf( stderr, "Error parseing BondTypes: line %d\n", lineNum );
444     exit(8);
445     }
446    
447     strcpy( currentBondType->type, the_token );
448    
449     if( !strcmp( currentBondType->type, "fixed" ) ){
450     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
451     fprintf( stderr, "Error parseing BondTypes: line %d\n", lineNum );
452     exit(8);
453     }
454    
455     sscanf( the_token, "%lf", &currentBondType->d0 );
456     }
457    
458     else{
459     fprintf(stderr,
460     "Unknown TraPPE_Ex bond type \"%s\" at line %d\n",
461     currentBondType->type,
462     lineNum );
463     exit(8);
464     }
465     }
466     }
467    
468     tempBondType = new LinkedType;
469     currentBondType->next = tempBondType;
470     currentBondType = tempBondType;
471    
472     eof_test = fgets( readLine, sizeof(readLine), frcFile );
473     lineNum++;
474     }
475    
476    
477     // initialize the Bonds
478    
479    
480     for( i=0; i<nBonds; i++ ){
481    
482     a = the_bonds[i].a;
483     b = the_bonds[i].b;
484    
485     atomA = the_atoms[a]->getType();
486     atomB = the_atoms[b]->getType();
487     currentBondType = headBondType->find( atomA, atomB );
488     if( currentBondType == NULL ){
489     fprintf( stderr, "BondType error, %s - %s not found in force file.\n",
490     atomA, atomB );
491     exit(8);
492     }
493    
494     if( !strcmp( currentBondType->type, "fixed" ) ){
495    
496     the_sris[i] = new ConstrainedBond( *the_atoms[a],
497     *the_atoms[b],
498     currentBondType->d0 );
499     entry_plug->n_constraints++;
500     }
501     }
502    
503    
504     // clean up the memory
505    
506     delete headBondType;
507    
508     }
509    
510     void TraPPE_ExFF::initializeBends( bend_set* the_bends ){
511    
512     class LinkedType {
513     public:
514     LinkedType(){
515     next = NULL;
516     nameA[0] = '\0';
517     nameB[0] = '\0';
518     nameC[0] = '\0';
519     type[0] = '\0';
520     }
521     ~LinkedType(){ if( next != NULL ) delete next; }
522    
523     LinkedType* find( char* key1, char* key2, char* key3 ){
524     if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 )
525     && !strcmp( nameC, key3 ) ) return this;
526     if( !strcmp( nameA, key3 ) && !strcmp( nameB, key2 )
527     && !strcmp( nameC, key1 ) ) return this;
528     if( next != NULL ) return next->find(key1, key2, key3);
529     return NULL;
530     }
531    
532     char nameA[15];
533     char nameB[15];
534     char nameC[15];
535     char type[30];
536     double k1, k2, k3, t0;
537    
538     LinkedType* next;
539     };
540    
541     LinkedType* headBendType;
542     LinkedType* currentBendType;
543     LinkedType* tempBendType;
544    
545     char readLine[500];
546     char* the_token;
547     char* eof_test;
548     int foundBend = 0;
549     int lineNum = 0;
550     int i, a, b, c, index;
551     char* atomA;
552     char* atomB;
553     char* atomC;
554     QuadraticBend* qBend;
555    
556     SRI **the_sris;
557     Atom** the_atoms;
558     int nBends;
559     the_sris = entry_plug->sr_interactions;
560     the_atoms = entry_plug->atoms;
561     nBends = entry_plug->n_bends;
562    
563     // read in the bend types.
564    
565     rewind( frcFile );
566     headBendType = new LinkedType;
567     currentBendType = headBendType;
568    
569     eof_test = fgets( readLine, sizeof(readLine), frcFile );
570     lineNum++;
571     if( eof_test == NULL ){
572     fprintf( stderr, "Error in reading Bends from force file.\n" );
573     exit(8);
574     }
575    
576    
577     while( !foundBend ){
578     while( eof_test != NULL && readLine[0] != '#' ){
579     eof_test = fgets( readLine, sizeof(readLine), frcFile );
580     lineNum++;
581     }
582     if( eof_test == NULL ){
583     fprintf( stderr,
584     "Error in reading Bends from force file at line %d.\n",
585     lineNum );
586     exit(8);
587     }
588    
589     the_token = strtok( readLine, " ,;\t#\n" );
590     foundBend = !strcmp( "BendTypes", the_token );
591    
592     if( !foundBend ){
593     eof_test = fgets( readLine, sizeof(readLine), frcFile );
594     lineNum++;
595    
596     if( eof_test == NULL ){
597     fprintf( stderr,
598     "Error in reading Bends from force file at line %d.\n",
599     lineNum );
600     exit(8);
601     }
602     }
603     }
604    
605     // we are now at the BendTypes section.
606    
607     eof_test = fgets( readLine, sizeof(readLine), frcFile );
608     lineNum++;
609    
610     if( eof_test == NULL ){
611     fprintf( stderr,
612     "Error in reading Bends from force file at line %d.\n",
613     lineNum );
614     exit(8);
615     }
616    
617     while( readLine[0] != '#' && eof_test != NULL ){
618    
619     if( readLine[0] != '!' ){
620    
621     the_token = strtok( readLine, " \n\t,;" );
622     if( the_token != NULL ){
623    
624     strcpy( currentBendType->nameA, the_token );
625    
626     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
627     fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
628     exit(8);
629     }
630    
631     strcpy( currentBendType->nameB, the_token );
632    
633     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
634     fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
635     exit(8);
636     }
637    
638     strcpy( currentBendType->nameC, the_token );
639    
640     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
641     fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
642     exit(8);
643     }
644    
645     strcpy( currentBendType->type, the_token );
646    
647     if( !strcmp( currentBendType->type, "quadratic" ) ){
648     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
649     fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
650     exit(8);
651     }
652    
653     sscanf( the_token, "%lf", &currentBendType->k1 );
654    
655     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
656     fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
657     exit(8);
658     }
659    
660     sscanf( the_token, "%lf", &currentBendType->k2 );
661    
662     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
663     fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
664     exit(8);
665     }
666    
667     sscanf( the_token, "%lf", &currentBendType->k3 );
668    
669     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
670     fprintf( stderr, "Error parseing BendTypes: line %d\n", lineNum );
671     exit(8);
672     }
673    
674     sscanf( the_token, "%lf", &currentBendType->t0 );
675     }
676    
677     else{
678     fprintf(stderr,
679     "Unknown TraPPE_Ex bend type \"%s\" at line %d\n",
680     currentBendType->type,
681     lineNum );
682     exit(8);
683     }
684     }
685     }
686    
687     tempBendType = new LinkedType;
688     currentBendType->next = tempBendType;
689     currentBendType = tempBendType;
690    
691     eof_test = fgets( readLine, sizeof(readLine), frcFile );
692     lineNum++;
693     }
694    
695    
696     // initialize the Bends
697    
698     for( i=0; i<nBends; i++ ){
699    
700     a = the_bends[i].a;
701     b = the_bends[i].b;
702     c = the_bends[i].c;
703    
704     atomA = the_atoms[a]->getType();
705     atomB = the_atoms[b]->getType();
706     atomC = the_atoms[c]->getType();
707     currentBendType = headBendType->find( atomA, atomB, atomC );
708     if( currentBendType == NULL ){
709     fprintf( stderr, "BendType error, %s - %s - %s not found"
710     " in force file.\n",
711     atomA, atomB, atomC );
712     exit(8);
713     }
714    
715     if( !strcmp( currentBendType->type, "quadratic" ) ){
716    
717     index = i + entry_plug->n_bonds;
718     qBend = new QuadraticBend( *the_atoms[a],
719     *the_atoms[b],
720     *the_atoms[c] );
721     qBend->setConstants( currentBendType->k1,
722     currentBendType->k2,
723     currentBendType->k3,
724     currentBendType->t0 );
725     the_sris[index] = qBend;
726     }
727     }
728    
729    
730     // clean up the memory
731    
732     delete headBendType;
733    
734     }
735    
736     void TraPPE_ExFF::initializeTorsions( torsion_set* the_torsions ){
737    
738     class LinkedType {
739     public:
740     LinkedType(){
741     next = NULL;
742     nameA[0] = '\0';
743     nameB[0] = '\0';
744     nameC[0] = '\0';
745     type[0] = '\0';
746     }
747     ~LinkedType(){ if( next != NULL ) delete next; }
748    
749     LinkedType* find( char* key1, char* key2, char* key3, char* key4 ){
750     if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 ) &&
751     !strcmp( nameC, key3 ) && !strcmp( nameD, key4 ) ) return this;
752    
753     if( !strcmp( nameA, key4 ) && !strcmp( nameB, key3 ) &&
754     !strcmp( nameC, key2 ) && !strcmp( nameD, key4 ) ) return this;
755    
756     if( next != NULL ) return next->find(key1, key2, key3, key4);
757     return NULL;
758     }
759    
760     char nameA[15];
761     char nameB[15];
762     char nameC[15];
763     char nameD[15];
764     char type[30];
765     double k1, k2, k3, k4;
766    
767     LinkedType* next;
768     };
769    
770     LinkedType* headTorsionType;
771     LinkedType* currentTorsionType;
772     LinkedType* tempTorsionType;
773    
774     char readLine[500];
775     char* the_token;
776     char* eof_test;
777     int foundTorsion = 0;
778     int lineNum = 0;
779     int i, a, b, c, d, index;
780     char* atomA;
781     char* atomB;
782     char* atomC;
783     char* atomD;
784     CubicTorsion* cTors;
785    
786     SRI **the_sris;
787     Atom** the_atoms;
788     int nTorsions;
789     the_sris = entry_plug->sr_interactions;
790     the_atoms = entry_plug->atoms;
791     nTorsions = entry_plug->n_torsions;
792    
793     // read in the torsion types.
794    
795     rewind( frcFile );
796     headTorsionType = new LinkedType;
797     currentTorsionType = headTorsionType;
798    
799     eof_test = fgets( readLine, sizeof(readLine), frcFile );
800     lineNum++;
801     if( eof_test == NULL ){
802     fprintf( stderr, "Error in reading Torsions from force file.\n" );
803     exit(8);
804     }
805    
806    
807     while( !foundTorsion ){
808     while( eof_test != NULL && readLine[0] != '#' ){
809     eof_test = fgets( readLine, sizeof(readLine), frcFile );
810     lineNum++;
811     }
812     if( eof_test == NULL ){
813     fprintf( stderr,
814     "Error in reading Torsions from force file at line %d.\n",
815     lineNum );
816     exit(8);
817     }
818    
819     the_token = strtok( readLine, " ,;\t#\n" );
820     foundTorsion = !strcmp( "TorsionTypes", the_token );
821    
822     if( !foundTorsion ){
823     eof_test = fgets( readLine, sizeof(readLine), frcFile );
824     lineNum++;
825    
826     if( eof_test == NULL ){
827     fprintf( stderr,
828     "Error in reading Torsions from force file at line %d.\n",
829     lineNum );
830     exit(8);
831     }
832     }
833     }
834    
835     // we are now at the TorsionTypes section.
836    
837     eof_test = fgets( readLine, sizeof(readLine), frcFile );
838     lineNum++;
839    
840     if( eof_test == NULL ){
841     fprintf( stderr,
842     "Error in reading Torsions from force file at line %d.\n",
843     lineNum );
844     exit(8);
845     }
846    
847     while( readLine[0] != '#' && eof_test != NULL ){
848    
849     if( readLine[0] != '!' ){
850    
851     the_token = strtok( readLine, " \n\t,;" );
852     if( the_token != NULL ){
853    
854     strcpy( currentTorsionType->nameA, the_token );
855    
856     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
857     fprintf( stderr, "Error parseing TorsionTypes: line %d\n", lineNum );
858     exit(8);
859     }
860    
861     strcpy( currentTorsionType->nameB, the_token );
862    
863     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
864     fprintf( stderr, "Error parseing TorsionTypes: line %d\n", lineNum );
865     exit(8);
866     }
867    
868     strcpy( currentTorsionType->nameC, the_token );
869    
870     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
871     fprintf( stderr, "Error parseing TorsionTypes: line %d\n", lineNum );
872     exit(8);
873     }
874    
875     strcpy( currentTorsionType->nameD, the_token );
876    
877     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
878     fprintf( stderr, "Error parseing TorsionTypes: line %d\n", lineNum );
879     exit(8);
880     }
881    
882     strcpy( currentTorsionType->type, the_token );
883    
884     if( !strcmp( currentTorsionType->type, "cubic" ) ){
885     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
886     fprintf(stderr,"Error parseing TorsionTypes: line %d\n", lineNum );
887     exit(8);
888     }
889    
890     sscanf( the_token, "%lf", &currentTorsionType->k1 );
891    
892     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
893     fprintf(stderr,"Error parseing TorsionTypes: line %d\n", lineNum );
894     exit(8);
895     }
896    
897     sscanf( the_token, "%lf", &currentTorsionType->k2 );
898    
899     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
900     fprintf(stderr,"Error parseing TorsionTypes: line %d\n", lineNum );
901     exit(8);
902     }
903    
904     sscanf( the_token, "%lf", &currentTorsionType->k3 );
905    
906     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
907     fprintf(stderr,"Error parseing TorsionTypes: line %d\n", lineNum );
908     exit(8);
909     }
910    
911     sscanf( the_token, "%lf", &currentTorsionType->k4 );
912     }
913    
914     else{
915     fprintf(stderr,
916     "Unknown TraPPE_Ex torsion type \"%s\" at line %d\n",
917     currentTorsionType->type,
918     lineNum );
919     exit(8);
920     }
921     }
922     }
923    
924     tempTorsionType = new LinkedType;
925     currentTorsionType->next = tempTorsionType;
926     currentTorsionType = tempTorsionType;
927    
928     eof_test = fgets( readLine, sizeof(readLine), frcFile );
929     lineNum++;
930     }
931    
932    
933     // initialize the Torsions
934    
935     for( i=0; i<nTorsions; i++ ){
936    
937     a = the_torsions[i].a;
938     b = the_torsions[i].b;
939     c = the_torsions[i].c;
940     d = the_torsions[i].d;
941    
942     atomA = the_atoms[a]->getType();
943     atomB = the_atoms[b]->getType();
944     atomC = the_atoms[c]->getType();
945     atomD = the_atoms[d]->getType();
946     currentTorsionType = headTorsionType->find( atomA, atomB, atomC, atomD );
947     if( currentTorsionType == NULL ){
948     fprintf( stderr, "TorsionType error, %s - %s - %s - %s not found"
949     " in force file.\n",
950     atomA, atomB, atomC, atomD );
951     exit(8);
952     }
953    
954     if( !strcmp( currentTorsionType->type, "cubic" ) ){
955     index = i + entry_plug->n_bonds + entry_plug->n_bends;
956    
957     cTors = new CubicTorsion( *the_atoms[a], *the_atoms[b],
958     *the_atoms[c], *the_atoms[d] );
959     cTors->setConstants( currentTorsionType->k1, currentTorsionType->k2,
960     currentTorsionType->k3, currentTorsionType->k4 );
961     the_sris[index] = cTors;
962     }
963     }
964    
965    
966     // clean up the memory
967    
968     delete headTorsionType;
969    
970     }