ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DUFF.cpp
Revision: 1097
Committed: Mon Apr 12 20:32:20 2004 UTC (20 years, 2 months ago) by gezelter
File size: 43347 byte(s)
Log Message:
Changes for RigidBody dynamics (Somewhat extensive)

File Contents

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