ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/UseTheForce/DUFF.cpp
Revision: 1670
Committed: Thu Oct 28 16:56:20 2004 UTC (19 years, 8 months ago) by gezelter
File size: 44243 byte(s)
Log Message:
fixed duplicate declaration foo

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     entry_plug->useSticky = 1;
777     }
778    
779     if (currentAtomType->isDipole) {
780     ((DirectionalAtomType*)at)->setDipole();
781     entry_plug->useDipoles = 1;
782     }
783    
784     at->setIdent(currentAtomType->ident);
785     at->setName(currentAtomType->name);
786     at->setLennardJones();
787     at->complete();
788     }
789     currentAtomType = currentAtomType->next;
790     }
791 gezelter 1490
792     currentAtomType = headAtomType->next;;
793 gezelter 1634 while( currentAtomType != NULL ){
794 gezelter 1490
795     if( currentAtomType->name[0] != '\0' ){
796     isError = 0;
797 gezelter 1634 newLJtype( &(currentAtomType->ident), &(currentAtomType->sigma),
798     &(currentAtomType->epslon), &isError);
799 gezelter 1490 if( isError ){
800 gezelter 1634 sprintf( painCave.errMsg,
801     "Error initializing the \"%s\" LJ type in fortran\n",
802     currentAtomType->name );
803     painCave.isFatal = 1;
804     simError();
805 gezelter 1490 }
806 gezelter 1634
807     if (currentAtomType->isDipole) {
808     newDipoleType(&(currentAtomType->ident), &(currentAtomType->dipole),
809     &isError);
810     if( isError ){
811     sprintf( painCave.errMsg,
812     "Error initializing the \"%s\" dipole type in fortran\n",
813     currentAtomType->name );
814     painCave.isFatal = 1;
815     simError();
816     }
817     }
818    
819     if(currentAtomType->isSSD) {
820     makeStickyType( &(currentAtomType->w0), &(currentAtomType->v0),
821     &(currentAtomType->v0p),
822     &(currentAtomType->rl), &(currentAtomType->ru),
823     &(currentAtomType->rlp), &(currentAtomType->rup));
824     }
825    
826 gezelter 1490 }
827     currentAtomType = currentAtomType->next;
828     }
829    
830     #ifdef IS_MPI
831     sprintf( checkPointMsg,
832     "DUFF atom structures successfully sent to fortran\n" );
833     MPIcheckPoint();
834     #endif // is_mpi
835    
836    
837    
838     // read in the bonds
839    
840     #ifdef IS_MPI
841     if( worldRank == 0 ){
842     #endif
843    
844     // read in the bond types.
845    
846     headBondType = new LinkedBondType;
847    
848     fastForward( "BondTypes", "initializeBonds" );
849    
850     // we are now at the bondTypes section
851    
852     eof_test = fgets( readLine, sizeof(readLine), frcFile );
853     lineNum++;
854    
855    
856     // read a line, and start parseing out the atom types
857    
858     if( eof_test == NULL ){
859     sprintf( painCave.errMsg,
860     "Error in reading bonds from force file at line %d.\n",
861     lineNum );
862     painCave.isFatal = 1;
863     simError();
864     }
865    
866     // stop reading at end of file, or at next section
867     while( readLine[0] != '#' && eof_test != NULL ){
868    
869     // toss comment lines
870     if( readLine[0] != '!' ){
871    
872     // the parser returns 0 if the line was blank
873     if( parseBond( readLine, lineNum, bondInfo ) ){
874     headBondType->add( bondInfo );
875     }
876     }
877     eof_test = fgets( readLine, sizeof(readLine), frcFile );
878     lineNum++;
879     }
880    
881     #ifdef IS_MPI
882    
883     // send out the linked list to all the other processes
884    
885     sprintf( checkPointMsg,
886     "DUFF bond structures read successfully." );
887     MPIcheckPoint();
888    
889     currentBondType = headBondType->next;
890     while( currentBondType != NULL ){
891     currentBondType->duplicate( bondInfo );
892     sendFrcStruct( &bondInfo, mpiBondStructType );
893     currentBondType = currentBondType->next;
894     }
895     bondInfo.last = 1;
896     sendFrcStruct( &bondInfo, mpiBondStructType );
897    
898     }
899    
900     else{
901    
902     // listen for node 0 to send out the force params
903    
904     MPIcheckPoint();
905    
906     headBondType = new LinkedBondType;
907     receiveFrcStruct( &bondInfo, mpiBondStructType );
908     while( !bondInfo.last ){
909    
910     headBondType->add( bondInfo );
911     receiveFrcStruct( &bondInfo, mpiBondStructType );
912     }
913     }
914    
915     sprintf( checkPointMsg,
916     "DUFF bond structures broadcast successfully." );
917     MPIcheckPoint();
918    
919     #endif // is_mpi
920    
921    
922     // read in the bends
923    
924     #ifdef IS_MPI
925     if( worldRank == 0 ){
926     #endif
927    
928     // read in the bend types.
929    
930     headBendType = new LinkedBendType;
931    
932     fastForward( "BendTypes", "initializeBends" );
933    
934     // we are now at the bendTypes section
935    
936     eof_test = fgets( readLine, sizeof(readLine), frcFile );
937     lineNum++;
938    
939     // read a line, and start parseing out the bend types
940    
941     if( eof_test == NULL ){
942     sprintf( painCave.errMsg,
943     "Error in reading bends from force file at line %d.\n",
944     lineNum );
945     painCave.isFatal = 1;
946     simError();
947     }
948    
949     // stop reading at end of file, or at next section
950     while( readLine[0] != '#' && eof_test != NULL ){
951    
952     // toss comment lines
953     if( readLine[0] != '!' ){
954    
955     // the parser returns 0 if the line was blank
956     if( parseBend( readLine, lineNum, bendInfo ) ){
957     headBendType->add( bendInfo );
958     }
959     }
960     eof_test = fgets( readLine, sizeof(readLine), frcFile );
961     lineNum++;
962     }
963    
964     #ifdef IS_MPI
965    
966     // send out the linked list to all the other processes
967    
968     sprintf( checkPointMsg,
969     "DUFF bend structures read successfully." );
970     MPIcheckPoint();
971    
972     currentBendType = headBendType->next;
973     while( currentBendType != NULL ){
974     currentBendType->duplicate( bendInfo );
975     sendFrcStruct( &bendInfo, mpiBendStructType );
976     currentBendType = currentBendType->next;
977     }
978     bendInfo.last = 1;
979     sendFrcStruct( &bendInfo, mpiBendStructType );
980    
981     }
982    
983     else{
984    
985     // listen for node 0 to send out the force params
986    
987     MPIcheckPoint();
988    
989     headBendType = new LinkedBendType;
990     receiveFrcStruct( &bendInfo, mpiBendStructType );
991     while( !bendInfo.last ){
992    
993     headBendType->add( bendInfo );
994     receiveFrcStruct( &bendInfo, mpiBendStructType );
995     }
996     }
997    
998     sprintf( checkPointMsg,
999     "DUFF bend structures broadcast successfully." );
1000     MPIcheckPoint();
1001    
1002     #endif // is_mpi
1003    
1004    
1005     // read in the torsions
1006    
1007     #ifdef IS_MPI
1008     if( worldRank == 0 ){
1009     #endif
1010    
1011     // read in the torsion types.
1012    
1013     headTorsionType = new LinkedTorsionType;
1014    
1015     fastForward( "TorsionTypes", "initializeTorsions" );
1016    
1017     // we are now at the torsionTypes section
1018    
1019     eof_test = fgets( readLine, sizeof(readLine), frcFile );
1020     lineNum++;
1021    
1022    
1023     // read a line, and start parseing out the atom types
1024    
1025     if( eof_test == NULL ){
1026     sprintf( painCave.errMsg,
1027     "Error in reading torsions from force file at line %d.\n",
1028     lineNum );
1029     painCave.isFatal = 1;
1030     simError();
1031     }
1032    
1033     // stop reading at end of file, or at next section
1034     while( readLine[0] != '#' && eof_test != NULL ){
1035    
1036     // toss comment lines
1037     if( readLine[0] != '!' ){
1038    
1039     // the parser returns 0 if the line was blank
1040     if( parseTorsion( readLine, lineNum, torsionInfo ) ){
1041     headTorsionType->add( torsionInfo );
1042    
1043     }
1044     }
1045     eof_test = fgets( readLine, sizeof(readLine), frcFile );
1046     lineNum++;
1047     }
1048    
1049     #ifdef IS_MPI
1050    
1051     // send out the linked list to all the other processes
1052    
1053     sprintf( checkPointMsg,
1054     "DUFF torsion structures read successfully." );
1055     MPIcheckPoint();
1056    
1057     currentTorsionType = headTorsionType->next;
1058     while( currentTorsionType != NULL ){
1059     currentTorsionType->duplicate( torsionInfo );
1060     sendFrcStruct( &torsionInfo, mpiTorsionStructType );
1061     currentTorsionType = currentTorsionType->next;
1062     }
1063     torsionInfo.last = 1;
1064     sendFrcStruct( &torsionInfo, mpiTorsionStructType );
1065    
1066     }
1067    
1068     else{
1069    
1070     // listen for node 0 to send out the force params
1071    
1072     MPIcheckPoint();
1073    
1074     headTorsionType = new LinkedTorsionType;
1075     receiveFrcStruct( &torsionInfo, mpiTorsionStructType );
1076     while( !torsionInfo.last ){
1077    
1078     headTorsionType->add( torsionInfo );
1079     receiveFrcStruct( &torsionInfo, mpiTorsionStructType );
1080     }
1081     }
1082    
1083     sprintf( checkPointMsg,
1084     "DUFF torsion structures broadcast successfully." );
1085     MPIcheckPoint();
1086    
1087     #endif // is_mpi
1088    
1089 chrisfen 1636 entry_plug->useLennardJones = 1;
1090 gezelter 1490 }
1091    
1092    
1093    
1094     void DUFF::initializeAtoms( int nAtoms, Atom** the_atoms ){
1095    
1096    
1097     //////////////////////////////////////////////////
1098     // a quick water fix
1099    
1100     double waterI[3][3];
1101     waterI[0][0] = 1.76958347772500;
1102     waterI[0][1] = 0.0;
1103     waterI[0][2] = 0.0;
1104    
1105     waterI[1][0] = 0.0;
1106     waterI[1][1] = 0.614537057924513;
1107     waterI[1][2] = 0.0;
1108    
1109     waterI[2][0] = 0.0;
1110     waterI[2][1] = 0.0;
1111     waterI[2][2] = 1.15504641980049;
1112    
1113    
1114     double headI[3][3];
1115     headI[0][0] = 1125;
1116     headI[0][1] = 0.0;
1117     headI[0][2] = 0.0;
1118    
1119     headI[1][0] = 0.0;
1120     headI[1][1] = 1125;
1121     headI[1][2] = 0.0;
1122    
1123     headI[2][0] = 0.0;
1124     headI[2][1] = 0.0;
1125     headI[2][2] = 250;
1126    
1127     //////////////////////////////////////////////////
1128    
1129    
1130     // initialize the atoms
1131    
1132     DirectionalAtom* dAtom;
1133     double ji[3];
1134    
1135     for(int i=0; i<nAtoms; i++ ){
1136    
1137     currentAtomType = headAtomType->find( the_atoms[i]->getType() );
1138     if( currentAtomType == NULL ){
1139     sprintf( painCave.errMsg,
1140     "AtomType error, %s not found in force file.\n",
1141     the_atoms[i]->getType() );
1142     painCave.isFatal = 1;
1143     simError();
1144     }
1145    
1146     the_atoms[i]->setMass( currentAtomType->mass );
1147     the_atoms[i]->setIdent( currentAtomType->ident );
1148    
1149     if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
1150    
1151     if( currentAtomType->isDipole ){
1152     if( the_atoms[i]->isDirectional() ){
1153    
1154     dAtom = (DirectionalAtom *) the_atoms[i];
1155     dAtom->setHasDipole( 1 );
1156    
1157     ji[0] = 0.0;
1158     ji[1] = 0.0;
1159     ji[2] = 0.0;
1160    
1161     dAtom->setJ( ji );
1162    
1163     if(!strcmp("SSD",the_atoms[i]->getType())){
1164     dAtom->setI( waterI );
1165     }
1166     else if(!strcmp("HEAD",the_atoms[i]->getType())){
1167     dAtom->setI( headI );
1168     }
1169     else{
1170     sprintf(painCave.errMsg,
1171     "AtmType error, %s does not have a moment of inertia set.\n",
1172     the_atoms[i]->getType() );
1173     painCave.isFatal = 1;
1174     simError();
1175     }
1176     entry_plug->n_dipoles++;
1177     }
1178     else{
1179    
1180     sprintf( painCave.errMsg,
1181     "DUFF error: Atom \"%s\" is a dipole, yet no standard"
1182     " orientation was specifed in the meta-data file.\n",
1183     currentAtomType->name );
1184     painCave.isFatal = 1;
1185     simError();
1186     }
1187     }
1188     else{
1189     if( the_atoms[i]->isDirectional() ){
1190     sprintf( painCave.errMsg,
1191     "DUFF error: Atom \"%s\" was given a standard "
1192     "orientation in the meta-data file, yet it is not a dipole.\n",
1193     currentAtomType->name);
1194     painCave.isFatal = 1;
1195     simError();
1196     }
1197     }
1198     }
1199     }
1200    
1201     void DUFF::initializeBonds( int nBonds, Bond** bondArray,
1202     bond_pair* the_bonds ){
1203     int i,a,b;
1204     char* atomA;
1205     char* atomB;
1206    
1207     Atom** the_atoms;
1208     the_atoms = entry_plug->atoms;
1209    
1210    
1211     // initialize the Bonds
1212    
1213     for( i=0; i<nBonds; i++ ){
1214    
1215     a = the_bonds[i].a;
1216     b = the_bonds[i].b;
1217    
1218     atomA = the_atoms[a]->getType();
1219     atomB = the_atoms[b]->getType();
1220     currentBondType = headBondType->find( atomA, atomB );
1221     if( currentBondType == NULL ){
1222     sprintf( painCave.errMsg,
1223     "BondType error, %s - %s not found in force file.\n",
1224     atomA, atomB );
1225     painCave.isFatal = 1;
1226     simError();
1227     }
1228    
1229     switch( currentBondType->type ){
1230    
1231     case FIXED_BOND:
1232    
1233     bondArray[i] = new ConstrainedBond( *the_atoms[a],
1234     *the_atoms[b],
1235     currentBondType->d0 );
1236     entry_plug->n_constraints++;
1237     break;
1238    
1239     case HARMONIC_BOND:
1240    
1241     bondArray[i] = new HarmonicBond( *the_atoms[a],
1242     *the_atoms[b],
1243     currentBondType->d0,
1244     currentBondType->k0 );
1245     break;
1246    
1247     default:
1248    
1249     break;
1250     // do nothing
1251     }
1252     }
1253     }
1254    
1255     void DUFF::initializeBends( int nBends, Bend** bendArray,
1256     bend_set* the_bends ){
1257    
1258     QuadraticBend* qBend;
1259     GhostBend* gBend;
1260     Atom** the_atoms;
1261     the_atoms = entry_plug->atoms;
1262    
1263     int i, a, b, c;
1264     char* atomA;
1265     char* atomB;
1266     char* atomC;
1267    
1268     // initialize the Bends
1269    
1270     for( i=0; i<nBends; i++ ){
1271    
1272     a = the_bends[i].a;
1273     b = the_bends[i].b;
1274     c = the_bends[i].c;
1275    
1276     atomA = the_atoms[a]->getType();
1277     atomB = the_atoms[b]->getType();
1278    
1279     if( the_bends[i].isGhost ) atomC = "GHOST";
1280     else atomC = the_atoms[c]->getType();
1281    
1282     currentBendType = headBendType->find( atomA, atomB, atomC );
1283     if( currentBendType == NULL ){
1284     sprintf( painCave.errMsg, "BendType error, %s - %s - %s not found"
1285     " in force file.\n",
1286     atomA, atomB, atomC );
1287     painCave.isFatal = 1;
1288     simError();
1289     }
1290    
1291     if( !strcmp( currentBendType->type, "quadratic" ) ){
1292    
1293     if( the_bends[i].isGhost){
1294    
1295     if( the_bends[i].ghost == b ){
1296     // do nothing
1297     }
1298     else if( the_bends[i].ghost == a ){
1299     c = a;
1300     a = b;
1301     b = c;
1302     }
1303     else{
1304     sprintf( painCave.errMsg,
1305     "BendType error, %s - %s - %s,\n"
1306     " --> central atom is not "
1307     "correctly identified with the "
1308     "\"ghostVectorSource = \" tag.\n",
1309     atomA, atomB, atomC );
1310     painCave.isFatal = 1;
1311     simError();
1312     }
1313    
1314     gBend = new GhostBend( *the_atoms[a],
1315     *the_atoms[b]);
1316    
1317     gBend->setConstants( currentBendType->k1,
1318     currentBendType->k2,
1319     currentBendType->k3,
1320     currentBendType->t0 );
1321     bendArray[i] = gBend;
1322     }
1323     else{
1324     qBend = new QuadraticBend( *the_atoms[a],
1325     *the_atoms[b],
1326     *the_atoms[c] );
1327     qBend->setConstants( currentBendType->k1,
1328     currentBendType->k2,
1329     currentBendType->k3,
1330     currentBendType->t0 );
1331     bendArray[i] = qBend;
1332     }
1333     }
1334     }
1335     }
1336    
1337     void DUFF::initializeTorsions( int nTorsions, Torsion** torsionArray,
1338     torsion_set* the_torsions ){
1339    
1340     int i, a, b, c, d;
1341     char* atomA;
1342     char* atomB;
1343     char* atomC;
1344     char* atomD;
1345    
1346     CubicTorsion* cTors;
1347     Atom** the_atoms;
1348     the_atoms = entry_plug->atoms;
1349    
1350     // initialize the Torsions
1351    
1352     for( i=0; i<nTorsions; i++ ){
1353    
1354     a = the_torsions[i].a;
1355     b = the_torsions[i].b;
1356     c = the_torsions[i].c;
1357     d = the_torsions[i].d;
1358    
1359     atomA = the_atoms[a]->getType();
1360     atomB = the_atoms[b]->getType();
1361     atomC = the_atoms[c]->getType();
1362     atomD = the_atoms[d]->getType();
1363     currentTorsionType = headTorsionType->find( atomA, atomB, atomC, atomD );
1364     if( currentTorsionType == NULL ){
1365     sprintf( painCave.errMsg,
1366     "TorsionType error, %s - %s - %s - %s not found"
1367     " in force file.\n",
1368     atomA, atomB, atomC, atomD );
1369     painCave.isFatal = 1;
1370     simError();
1371     }
1372    
1373     if( !strcmp( currentTorsionType->type, "cubic" ) ){
1374    
1375     cTors = new CubicTorsion( *the_atoms[a], *the_atoms[b],
1376     *the_atoms[c], *the_atoms[d] );
1377     cTors->setConstants( currentTorsionType->k1, currentTorsionType->k2,
1378     currentTorsionType->k3, currentTorsionType->k4 );
1379     torsionArray[i] = cTors;
1380     }
1381     }
1382     }
1383    
1384     void DUFF::fastForward( char* stopText, char* searchOwner ){
1385    
1386     int foundText = 0;
1387     char* the_token;
1388    
1389     rewind( frcFile );
1390     lineNum = 0;
1391    
1392     eof_test = fgets( readLine, sizeof(readLine), frcFile );
1393     lineNum++;
1394     if( eof_test == NULL ){
1395     sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
1396     " file is empty.\n",
1397     searchOwner );
1398     painCave.isFatal = 1;
1399     simError();
1400     }
1401    
1402    
1403     while( !foundText ){
1404     while( eof_test != NULL && readLine[0] != '#' ){
1405     eof_test = fgets( readLine, sizeof(readLine), frcFile );
1406     lineNum++;
1407     }
1408     if( eof_test == NULL ){
1409     sprintf( painCave.errMsg,
1410     "Error fast forwarding force file for %s at "
1411     "line %d: file ended unexpectedly.\n",
1412     searchOwner,
1413     lineNum );
1414     painCave.isFatal = 1;
1415     simError();
1416     }
1417    
1418     the_token = strtok( readLine, " ,;\t#\n" );
1419     foundText = !strcmp( stopText, the_token );
1420    
1421     if( !foundText ){
1422     eof_test = fgets( readLine, sizeof(readLine), frcFile );
1423     lineNum++;
1424    
1425     if( eof_test == NULL ){
1426     sprintf( painCave.errMsg,
1427     "Error fast forwarding force file for %s at "
1428     "line %d: file ended unexpectedly.\n",
1429     searchOwner,
1430     lineNum );
1431     painCave.isFatal = 1;
1432     simError();
1433     }
1434     }
1435     }
1436     }
1437    
1438    
1439     int DUFF_NS::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
1440    
1441     char* the_token;
1442    
1443     the_token = strtok( lineBuffer, " \n\t,;" );
1444     if( the_token != NULL ){
1445    
1446     strcpy( info.name, the_token );
1447    
1448     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1449     sprintf( painCave.errMsg,
1450     "Error parseing AtomTypes: line %d\n", lineNum );
1451     painCave.isFatal = 1;
1452     simError();
1453     }
1454    
1455     info.isDipole = atoi( the_token );
1456    
1457     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1458     sprintf( painCave.errMsg,
1459     "Error parseing AtomTypes: line %d\n", lineNum );
1460     painCave.isFatal = 1;
1461     simError();
1462     }
1463    
1464     info.isSSD = atoi( the_token );
1465    
1466     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1467     sprintf( painCave.errMsg,
1468     "Error parseing AtomTypes: line %d\n", lineNum );
1469     painCave.isFatal = 1;
1470     simError();
1471     }
1472    
1473     info.mass = atof( the_token );
1474    
1475     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1476     sprintf( painCave.errMsg,
1477     "Error parseing AtomTypes: line %d\n", lineNum );
1478     painCave.isFatal = 1;
1479     simError();
1480     }
1481    
1482     info.epslon = atof( the_token );
1483    
1484     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1485     sprintf( painCave.errMsg,
1486     "Error parseing AtomTypes: line %d\n", lineNum );
1487     painCave.isFatal = 1;
1488     simError();
1489     }
1490    
1491     info.sigma = atof( the_token );
1492    
1493     if( info.isDipole ){
1494    
1495     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1496     sprintf( painCave.errMsg,
1497     "Error parseing AtomTypes: line %d\n", lineNum );
1498     painCave.isFatal = 1;
1499     simError();
1500     }
1501    
1502     info.dipole = atof( the_token );
1503     }
1504     else info.dipole = 0.0;
1505    
1506     if( info.isSSD ){
1507    
1508     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1509     sprintf( painCave.errMsg,
1510     "Error parseing AtomTypes: line %d\n", lineNum );
1511     painCave.isFatal = 1;
1512     simError();
1513     }
1514    
1515     info.w0 = atof( the_token );
1516    
1517     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1518     sprintf( painCave.errMsg,
1519     "Error parseing AtomTypes: line %d\n", lineNum );
1520     painCave.isFatal = 1;
1521     simError();
1522     }
1523    
1524     info.v0 = atof( the_token );
1525     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1526     sprintf( painCave.errMsg,
1527     "Error parseing AtomTypes: line %d\n", lineNum );
1528     painCave.isFatal = 1;
1529     simError();
1530     }
1531    
1532     info.v0p = atof( the_token );
1533    
1534     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1535     sprintf( painCave.errMsg,
1536     "Error parseing AtomTypes: line %d\n", lineNum );
1537     painCave.isFatal = 1;
1538     simError();
1539     }
1540    
1541     info.rl = atof( the_token );
1542    
1543     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1544     sprintf( painCave.errMsg,
1545     "Error parseing AtomTypes: line %d\n", lineNum );
1546     painCave.isFatal = 1;
1547     simError();
1548     }
1549    
1550     info.ru = atof( the_token );
1551    
1552     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1553     sprintf( painCave.errMsg,
1554     "Error parseing AtomTypes: line %d\n", lineNum );
1555     painCave.isFatal = 1;
1556     simError();
1557     }
1558    
1559     info.rlp = atof( the_token );
1560    
1561     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1562     sprintf( painCave.errMsg,
1563     "Error parseing AtomTypes: line %d\n", lineNum );
1564     painCave.isFatal = 1;
1565     simError();
1566     }
1567    
1568     info.rup = atof( the_token );
1569     }
1570     else info.v0 = info.w0 = info.v0p = info.rl = info.ru = info.rlp = info.rup = 0.0;
1571    
1572     return 1;
1573     }
1574     else return 0;
1575     }
1576    
1577     int DUFF_NS::parseBond( char *lineBuffer, int lineNum, bondStruct &info ){
1578    
1579     char* the_token;
1580     char bondType[30];
1581    
1582     the_token = strtok( lineBuffer, " \n\t,;" );
1583     if( the_token != NULL ){
1584    
1585     strcpy( info.nameA, the_token );
1586    
1587     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1588     sprintf( painCave.errMsg,
1589     "Error parseing BondTypes: line %d\n", lineNum );
1590     painCave.isFatal = 1;
1591     simError();
1592     }
1593    
1594     strcpy( info.nameB, the_token );
1595    
1596     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1597     sprintf( painCave.errMsg,
1598     "Error parseing BondTypes: line %d\n", lineNum );
1599     painCave.isFatal = 1;
1600     simError();
1601     }
1602    
1603     strcpy( bondType, the_token );
1604    
1605     if( !strcmp( bondType, "fixed" ) ){
1606     info.type = FIXED_BOND;
1607    
1608     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1609     sprintf( painCave.errMsg,
1610     "Error parseing BondTypes: line %d\n", lineNum );
1611     painCave.isFatal = 1;
1612     simError();
1613     }
1614    
1615     info.d0 = atof( the_token );
1616    
1617     info.k0=0.0;
1618     }
1619     else if( !strcmp( bondType, "harmonic" ) ){
1620     info.type = HARMONIC_BOND;
1621    
1622     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1623     sprintf( painCave.errMsg,
1624     "Error parseing BondTypes: line %d\n", lineNum );
1625     painCave.isFatal = 1;
1626     simError();
1627     }
1628    
1629     info.d0 = atof( the_token );
1630    
1631     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1632     sprintf( painCave.errMsg,
1633     "Error parseing BondTypes: line %d\n", lineNum );
1634     painCave.isFatal = 1;
1635     simError();
1636     }
1637    
1638     info.k0 = atof( the_token );
1639     }
1640    
1641     else{
1642     sprintf( painCave.errMsg,
1643     "Unknown DUFF bond type \"%s\" at line %d\n",
1644     bondType,
1645     lineNum );
1646     painCave.isFatal = 1;
1647     simError();
1648     }
1649    
1650     return 1;
1651     }
1652     else return 0;
1653     }
1654    
1655    
1656     int DUFF_NS::parseBend( char *lineBuffer, int lineNum, bendStruct &info ){
1657    
1658     char* the_token;
1659    
1660     the_token = strtok( lineBuffer, " \n\t,;" );
1661     if( the_token != NULL ){
1662    
1663     strcpy( info.nameA, the_token );
1664    
1665     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1666     sprintf( painCave.errMsg,
1667     "Error parseing BendTypes: line %d\n", lineNum );
1668     painCave.isFatal = 1;
1669     simError();
1670     }
1671    
1672     strcpy( info.nameB, the_token );
1673    
1674     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1675     sprintf( painCave.errMsg,
1676     "Error parseing BendTypes: line %d\n", lineNum );
1677     painCave.isFatal = 1;
1678     simError();
1679     }
1680    
1681     strcpy( info.nameC, the_token );
1682    
1683     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1684     sprintf( painCave.errMsg,
1685     "Error parseing BendTypes: line %d\n", lineNum );
1686     painCave.isFatal = 1;
1687     simError();
1688     }
1689    
1690     strcpy( info.type, the_token );
1691    
1692     if( !strcmp( info.type, "quadratic" ) ){
1693     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1694     sprintf( painCave.errMsg,
1695     "Error parseing BendTypes: line %d\n", lineNum );
1696     painCave.isFatal = 1;
1697     simError();
1698     }
1699    
1700     info.k1 = atof( the_token );
1701    
1702     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1703     sprintf( painCave.errMsg,
1704     "Error parseing BendTypes: line %d\n", lineNum );
1705     painCave.isFatal = 1;
1706     simError();
1707     }
1708    
1709     info.k2 = atof( the_token );
1710    
1711     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1712     sprintf( painCave.errMsg,
1713     "Error parseing BendTypes: line %d\n", lineNum );
1714     painCave.isFatal = 1;
1715     simError();
1716     }
1717    
1718     info.k3 = atof( the_token );
1719    
1720     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1721     sprintf( painCave.errMsg,
1722     "Error parseing BendTypes: line %d\n", lineNum );
1723     painCave.isFatal = 1;
1724     simError();
1725     }
1726    
1727     info.t0 = atof( the_token );
1728     }
1729    
1730     else{
1731     sprintf( painCave.errMsg,
1732     "Unknown DUFF bend type \"%s\" at line %d\n",
1733     info.type,
1734     lineNum );
1735     painCave.isFatal = 1;
1736     simError();
1737     }
1738    
1739     return 1;
1740     }
1741     else return 0;
1742     }
1743    
1744     int DUFF_NS::parseTorsion( char *lineBuffer, int lineNum, torsionStruct &info ){
1745    
1746     char* the_token;
1747    
1748     the_token = strtok( lineBuffer, " \n\t,;" );
1749     if( the_token != NULL ){
1750    
1751     strcpy( info.nameA, the_token );
1752    
1753     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1754     sprintf( painCave.errMsg,
1755     "Error parseing TorsionTypes: line %d\n", lineNum );
1756     painCave.isFatal = 1;
1757     simError();
1758     }
1759    
1760     strcpy( info.nameB, the_token );
1761    
1762     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1763     sprintf( painCave.errMsg,
1764     "Error parseing TorsionTypes: line %d\n", lineNum );
1765     painCave.isFatal = 1;
1766     simError();
1767     }
1768    
1769     strcpy( info.nameC, the_token );
1770    
1771     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1772     sprintf( painCave.errMsg,
1773     "Error parseing TorsionTypes: line %d\n", lineNum );
1774     painCave.isFatal = 1;
1775     simError();
1776     }
1777    
1778     strcpy( info.nameD, the_token );
1779    
1780     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1781     sprintf( painCave.errMsg,
1782     "Error parseing TorsionTypes: line %d\n", lineNum );
1783     painCave.isFatal = 1;
1784     simError();
1785     }
1786    
1787     strcpy( info.type, the_token );
1788    
1789     if( !strcmp( info.type, "cubic" ) ){
1790     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1791     sprintf( painCave.errMsg,
1792     "Error parseing TorsionTypes: line %d\n", lineNum );
1793     painCave.isFatal = 1;
1794     simError();
1795     }
1796    
1797     info.k1 = atof( the_token );
1798    
1799     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1800     sprintf( painCave.errMsg,
1801     "Error parseing TorsionTypes: line %d\n", lineNum );
1802     painCave.isFatal = 1;
1803     simError();
1804     }
1805    
1806     info.k2 = atof( the_token );
1807    
1808     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1809     sprintf( painCave.errMsg,
1810     "Error parseing TorsionTypes: line %d\n", lineNum );
1811     painCave.isFatal = 1;
1812     simError();
1813     }
1814    
1815     info.k3 = atof( the_token );
1816    
1817     if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1818     sprintf( painCave.errMsg,
1819     "Error parseing TorsionTypes: line %d\n", lineNum );
1820     painCave.isFatal = 1;
1821     simError();
1822     }
1823    
1824     info.k4 = atof( the_token );
1825    
1826     }
1827    
1828     else{
1829     sprintf( painCave.errMsg,
1830     "Unknown DUFF torsion type \"%s\" at line %d\n",
1831     info.type,
1832     lineNum );
1833     painCave.isFatal = 1;
1834     simError();
1835     }
1836    
1837     return 1;
1838     }
1839    
1840     else return 0;
1841     }