ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/UseTheForce/DUFF.cpp
Revision: 1706
Committed: Thu Nov 4 16:20:28 2004 UTC (19 years, 8 months ago) by gezelter
File size: 44344 byte(s)
Log Message:
fixed useXXX in the entry_plug so that it only is
set if the atoms really are in the simulation

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