ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DUFF.cpp
Revision: 829
Committed: Tue Oct 28 16:03:37 2003 UTC (20 years, 8 months ago) by gezelter
File size: 43421 byte(s)
Log Message:
replace c++ header stuff with more portable c header stuff
Also, mod file fixes and portability changes
Some fortran changes will need to be reversed.

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     double dipole;
36     double w0;
37     double v0;
38 gezelter 635 double v0p;
39     double rl;
40     double ru;
41     double rlp;
42     double rup;
43 mmeineke 559 int isSSD;
44     int isDipole;
45     int ident;
46     int last; // 0 -> default
47     // 1 -> tells nodes to stop listening
48     } atomStruct;
49    
50    
51     typedef struct{
52     char nameA[15];
53     char nameB[15];
54     double d0;
55 mmeineke 564 double k0;
56 mmeineke 559 int last; // 0 -> default
57     // 1 -> tells nodes to stop listening
58 mmeineke 564 int type;
59 mmeineke 559 } bondStruct;
60    
61    
62     typedef struct{
63     char nameA[15];
64     char nameB[15];
65     char nameC[15];
66     char type[30];
67     double k1, k2, k3, t0;
68     int last; // 0 -> default
69     // 1 -> tells nodes to stop listening
70     } bendStruct;
71    
72    
73     typedef struct{
74     char nameA[15];
75     char nameB[15];
76     char nameC[15];
77     char nameD[15];
78     char type[30];
79     double k1, k2, k3, k4;
80     int last; // 0 -> default
81     // 1 -> tells nodes to stop listening
82     } torsionStruct;
83    
84    
85     int parseAtom( char *lineBuffer, int lineNum, atomStruct &info );
86     int parseBond( char *lineBuffer, int lineNum, bondStruct &info );
87     int parseBend( char *lineBuffer, int lineNum, bendStruct &info );
88     int parseTorsion( char *lineBuffer, int lineNum, torsionStruct &info );
89    
90    
91     #ifdef IS_MPI
92    
93     MPI_Datatype mpiAtomStructType;
94     MPI_Datatype mpiBondStructType;
95     MPI_Datatype mpiBendStructType;
96     MPI_Datatype mpiTorsionStructType;
97    
98     #endif
99    
100     class LinkedAtomType {
101     public:
102     LinkedAtomType(){
103     next = NULL;
104     name[0] = '\0';
105     }
106     ~LinkedAtomType(){ if( next != NULL ) delete next; }
107    
108     LinkedAtomType* find(char* key){
109     if( !strcmp(name, key) ) return this;
110     if( next != NULL ) return next->find(key);
111     return NULL;
112     }
113    
114     void printMe( void ){
115    
116     std::cerr << "LinkedAtype " << name << ": ident = " << ident << "\n";
117     if( next != NULL ) next->printMe();
118    
119     }
120    
121     void add( atomStruct &info ){
122    
123     // check for duplicates
124    
125     if( !strcmp( info.name, name ) ){
126     sprintf( painCave.errMsg,
127 mmeineke 561 "Duplicate DUFF atom type \"%s\" found in "
128     "the DUFF param file./n",
129 mmeineke 559 name );
130     painCave.isFatal = 1;
131     simError();
132     }
133    
134     if( next != NULL ) next->add(info);
135     else{
136     next = new LinkedAtomType();
137     strcpy(next->name, info.name);
138     next->isDipole = info.isDipole;
139     next->isSSD = info.isSSD;
140     next->mass = info.mass;
141     next->epslon = info.epslon;
142     next->sigma = info.sigma;
143     next->dipole = info.dipole;
144     next->w0 = info.w0;
145     next->v0 = info.v0;
146 gezelter 635 next->v0p = info.v0p;
147     next->rl = info.rl;
148     next->ru = info.ru;
149     next->rlp = info.rlp;
150     next->rup = info.rup;
151 mmeineke 559 next->ident = info.ident;
152     }
153     }
154    
155     #ifdef IS_MPI
156    
157     void duplicate( atomStruct &info ){
158     strcpy(info.name, name);
159     info.isDipole = isDipole;
160     info.isSSD = isSSD;
161     info.mass = mass;
162     info.epslon = epslon;
163     info.sigma = sigma;
164     info.dipole = dipole;
165     info.w0 = w0;
166     info.v0 = v0;
167 gezelter 635 info.v0p = v0p;
168     info.rl = rl;
169     info.ru = ru;
170     info.rlp = rlp;
171     info.rup = rup;
172 mmeineke 559 info.ident = ident;
173     info.last = 0;
174     }
175    
176    
177     #endif
178    
179     char name[15];
180     int isDipole;
181     int isSSD;
182     double mass;
183     double epslon;
184     double sigma;
185     double dipole;
186     double w0;
187     double v0;
188 gezelter 635 double v0p;
189     double rl;
190     double ru;
191     double rlp;
192     double rup;
193 mmeineke 559 int ident;
194     LinkedAtomType* next;
195     };
196    
197     class LinkedBondType {
198     public:
199     LinkedBondType(){
200     next = NULL;
201     nameA[0] = '\0';
202     nameB[0] = '\0';
203     }
204     ~LinkedBondType(){ if( next != NULL ) delete next; }
205    
206     LinkedBondType* find(char* key1, char* key2){
207     if( !strcmp(nameA, key1 ) && !strcmp( nameB, key2 ) ) return this;
208     if( !strcmp(nameA, key2 ) && !strcmp( nameB, key1 ) ) return this;
209     if( next != NULL ) return next->find(key1, key2);
210     return NULL;
211     }
212    
213    
214     void add( bondStruct &info ){
215    
216     // check for duplicates
217     int dup = 0;
218    
219     if( !strcmp(nameA, info.nameA ) && !strcmp( nameB, info.nameB ) ) dup = 1;
220     if( !strcmp(nameA, info.nameB ) && !strcmp( nameB, info.nameA ) ) dup = 1;
221    
222     if(dup){
223     sprintf( painCave.errMsg,
224 mmeineke 561 "Duplicate DUFF bond type \"%s - %s\" found in "
225     "the DUFF param file./n",
226 mmeineke 559 nameA, nameB );
227     painCave.isFatal = 1;
228     simError();
229     }
230    
231    
232     if( next != NULL ) next->add(info);
233     else{
234     next = new LinkedBondType();
235     strcpy(next->nameA, info.nameA);
236     strcpy(next->nameB, info.nameB);
237 mmeineke 564 next->type = info.type;
238 mmeineke 559 next->d0 = info.d0;
239 mmeineke 564 next->k0 = info.k0;
240 mmeineke 559 }
241     }
242    
243     #ifdef IS_MPI
244     void duplicate( bondStruct &info ){
245     strcpy(info.nameA, nameA);
246     strcpy(info.nameB, nameB);
247 mmeineke 564 info.type = type;
248 mmeineke 559 info.d0 = d0;
249 mmeineke 564 info.k0 = k0;
250 mmeineke 559 info.last = 0;
251     }
252    
253    
254     #endif
255    
256     char nameA[15];
257     char nameB[15];
258 mmeineke 564 int type;
259 mmeineke 559 double d0;
260 mmeineke 564 double k0;
261 mmeineke 559
262     LinkedBondType* next;
263     };
264    
265     class LinkedBendType {
266     public:
267     LinkedBendType(){
268     next = NULL;
269     nameA[0] = '\0';
270     nameB[0] = '\0';
271     nameC[0] = '\0';
272     type[0] = '\0';
273     }
274     ~LinkedBendType(){ if( next != NULL ) delete next; }
275    
276     LinkedBendType* find( char* key1, char* key2, char* key3 ){
277     if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 )
278     && !strcmp( nameC, key3 ) ) return this;
279     if( !strcmp( nameA, key3 ) && !strcmp( nameB, key2 )
280     && !strcmp( nameC, key1 ) ) return this;
281     if( next != NULL ) return next->find(key1, key2, key3);
282     return NULL;
283     }
284    
285     void add( bendStruct &info ){
286    
287     // check for duplicates
288     int dup = 0;
289    
290     if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB )
291     && !strcmp( nameC, info.nameC ) ) dup = 1;
292     if( !strcmp( nameA, info.nameC ) && !strcmp( nameB, info.nameB )
293     && !strcmp( nameC, info.nameA ) ) dup = 1;
294    
295     if(dup){
296     sprintf( painCave.errMsg,
297 mmeineke 561 "Duplicate DUFF bend type \"%s - %s - %s\" found in "
298     "the DUFF param file./n",
299 mmeineke 559 nameA, nameB, nameC );
300     painCave.isFatal = 1;
301     simError();
302     }
303    
304     if( next != NULL ) next->add(info);
305     else{
306     next = new LinkedBendType();
307     strcpy(next->nameA, info.nameA);
308     strcpy(next->nameB, info.nameB);
309     strcpy(next->nameC, info.nameC);
310     strcpy(next->type, info.type);
311     next->k1 = info.k1;
312     next->k2 = info.k2;
313     next->k3 = info.k3;
314     next->t0 = info.t0;
315     }
316     }
317    
318     #ifdef IS_MPI
319    
320     void duplicate( bendStruct &info ){
321     strcpy(info.nameA, nameA);
322     strcpy(info.nameB, nameB);
323     strcpy(info.nameC, nameC);
324     strcpy(info.type, type);
325     info.k1 = k1;
326     info.k2 = k2;
327     info.k3 = k3;
328     info.t0 = t0;
329     info.last = 0;
330     }
331    
332     #endif // is_mpi
333    
334     char nameA[15];
335     char nameB[15];
336     char nameC[15];
337     char type[30];
338     double k1, k2, k3, t0;
339    
340     LinkedBendType* next;
341     };
342    
343     class LinkedTorsionType {
344     public:
345     LinkedTorsionType(){
346     next = NULL;
347     nameA[0] = '\0';
348     nameB[0] = '\0';
349     nameC[0] = '\0';
350     type[0] = '\0';
351     }
352     ~LinkedTorsionType(){ if( next != NULL ) delete next; }
353    
354     LinkedTorsionType* find( char* key1, char* key2, char* key3, char* key4 ){
355    
356    
357    
358    
359     if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 ) &&
360     !strcmp( nameC, key3 ) && !strcmp( nameD, key4 ) ) return this;
361    
362     if( !strcmp( nameA, key4 ) && !strcmp( nameB, key3 ) &&
363     !strcmp( nameC, key2 ) && !strcmp( nameD, key1 ) ) return this;
364    
365     if( next != NULL ) return next->find(key1, key2, key3, key4);
366     return NULL;
367     }
368    
369     void add( torsionStruct &info ){
370    
371     // check for duplicates
372     int dup = 0;
373    
374     if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB ) &&
375     !strcmp( nameC, info.nameC ) && !strcmp( nameD, info.nameD ) ) dup = 1;
376    
377     if( !strcmp( nameA, info.nameD ) && !strcmp( nameB, info.nameC ) &&
378     !strcmp( nameC, info.nameB ) && !strcmp( nameD, info.nameA ) ) dup = 1;
379    
380     if(dup){
381     sprintf( painCave.errMsg,
382 mmeineke 561 "Duplicate DUFF torsion type \"%s - %s - %s - %s\" found in "
383     "the DUFF param file./n", nameA, nameB, nameC, nameD );
384 mmeineke 559 painCave.isFatal = 1;
385     simError();
386     }
387    
388     if( next != NULL ) next->add(info);
389     else{
390     next = new LinkedTorsionType();
391     strcpy(next->nameA, info.nameA);
392     strcpy(next->nameB, info.nameB);
393     strcpy(next->nameC, info.nameC);
394     strcpy(next->nameD, info.nameD);
395     strcpy(next->type, info.type);
396     next->k1 = info.k1;
397     next->k2 = info.k2;
398     next->k3 = info.k3;
399     next->k4 = info.k4;
400    
401     }
402     }
403    
404     #ifdef IS_MPI
405    
406     void duplicate( torsionStruct &info ){
407     strcpy(info.nameA, nameA);
408     strcpy(info.nameB, nameB);
409     strcpy(info.nameC, nameC);
410     strcpy(info.nameD, nameD);
411     strcpy(info.type, type);
412     info.k1 = k1;
413     info.k2 = k2;
414     info.k3 = k3;
415     info.k4 = k4;
416     info.last = 0;
417     }
418    
419     #endif
420    
421     char nameA[15];
422     char nameB[15];
423     char nameC[15];
424     char nameD[15];
425     char type[30];
426     double k1, k2, k3, k4;
427    
428     LinkedTorsionType* next;
429     };
430    
431    
432     LinkedAtomType* headAtomType;
433     LinkedAtomType* currentAtomType;
434     LinkedBondType* headBondType;
435     LinkedBondType* currentBondType;
436     LinkedBendType* headBendType;
437     LinkedBendType* currentBendType;
438     LinkedTorsionType* headTorsionType;
439     LinkedTorsionType* currentTorsionType;
440    
441     } // namespace
442    
443 mmeineke 594 using namespace DUFF_NS;
444 mmeineke 559
445    
446     //****************************************************************
447     // begins the actual forcefield stuff.
448     //****************************************************************
449    
450    
451 mmeineke 561 DUFF::DUFF(){
452 mmeineke 559
453     char fileName[200];
454     char* ffPath_env = "FORCE_PARAM_PATH";
455     char* ffPath;
456     char temp[200];
457    
458     headAtomType = NULL;
459     currentAtomType = NULL;
460     headBondType = NULL;
461     currentBondType = NULL;
462     headBendType = NULL;
463     currentBendType = NULL;
464     headTorsionType = NULL;
465     currentTorsionType = NULL;
466    
467     // do the funtion wrapping
468     wrapMeFF( this );
469    
470    
471     #ifdef IS_MPI
472     int i;
473    
474     // **********************************************************************
475     // Init the atomStruct mpi type
476    
477     atomStruct atomProto; // mpiPrototype
478 mmeineke 721 int atomBC[3] = {15,11,4}; // block counts
479 mmeineke 559 MPI_Aint atomDspls[3]; // displacements
480     MPI_Datatype atomMbrTypes[3]; // member mpi types
481    
482     MPI_Address(&atomProto.name, &atomDspls[0]);
483     MPI_Address(&atomProto.mass, &atomDspls[1]);
484     MPI_Address(&atomProto.isSSD, &atomDspls[2]);
485    
486     atomMbrTypes[0] = MPI_CHAR;
487     atomMbrTypes[1] = MPI_DOUBLE;
488     atomMbrTypes[2] = MPI_INT;
489    
490     for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
491    
492     MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
493     MPI_Type_commit(&mpiAtomStructType);
494    
495    
496     // **********************************************************************
497     // Init the bondStruct mpi type
498    
499     bondStruct bondProto; // mpiPrototype
500 mmeineke 564 int bondBC[3] = {30,2,2}; // block counts
501 mmeineke 559 MPI_Aint bondDspls[3]; // displacements
502     MPI_Datatype bondMbrTypes[3]; // member mpi types
503    
504     MPI_Address(&bondProto.nameA, &bondDspls[0]);
505     MPI_Address(&bondProto.d0, &bondDspls[1]);
506     MPI_Address(&bondProto.last, &bondDspls[2]);
507    
508     bondMbrTypes[0] = MPI_CHAR;
509     bondMbrTypes[1] = MPI_DOUBLE;
510     bondMbrTypes[2] = MPI_INT;
511    
512     for (i=2; i >= 0; i--) bondDspls[i] -= bondDspls[0];
513    
514     MPI_Type_struct(3, bondBC, bondDspls, bondMbrTypes, &mpiBondStructType);
515     MPI_Type_commit(&mpiBondStructType);
516    
517    
518     // **********************************************************************
519     // Init the bendStruct mpi type
520    
521     bendStruct bendProto; // mpiPrototype
522     int bendBC[3] = {75,4,1}; // block counts
523     MPI_Aint bendDspls[3]; // displacements
524     MPI_Datatype bendMbrTypes[3]; // member mpi types
525    
526     MPI_Address(&bendProto.nameA, &bendDspls[0]);
527     MPI_Address(&bendProto.k1, &bendDspls[1]);
528     MPI_Address(&bendProto.last, &bendDspls[2]);
529    
530     bendMbrTypes[0] = MPI_CHAR;
531     bendMbrTypes[1] = MPI_DOUBLE;
532     bendMbrTypes[2] = MPI_INT;
533    
534     for (i=2; i >= 0; i--) bendDspls[i] -= bendDspls[0];
535    
536     MPI_Type_struct(3, bendBC, bendDspls, bendMbrTypes, &mpiBendStructType);
537     MPI_Type_commit(&mpiBendStructType);
538    
539    
540     // **********************************************************************
541     // Init the torsionStruct mpi type
542    
543     torsionStruct torsionProto; // mpiPrototype
544     int torsionBC[3] = {90,4,1}; // block counts
545     MPI_Aint torsionDspls[3]; // displacements
546     MPI_Datatype torsionMbrTypes[3]; // member mpi types
547    
548     MPI_Address(&torsionProto.nameA, &torsionDspls[0]);
549     MPI_Address(&torsionProto.k1, &torsionDspls[1]);
550     MPI_Address(&torsionProto.last, &torsionDspls[2]);
551    
552     torsionMbrTypes[0] = MPI_CHAR;
553     torsionMbrTypes[1] = MPI_DOUBLE;
554     torsionMbrTypes[2] = MPI_INT;
555    
556     for (i=2; i >= 0; i--) torsionDspls[i] -= torsionDspls[0];
557    
558     MPI_Type_struct(3, torsionBC, torsionDspls, torsionMbrTypes,
559     &mpiTorsionStructType);
560     MPI_Type_commit(&mpiTorsionStructType);
561    
562     // ***********************************************************************
563    
564     if( worldRank == 0 ){
565     #endif
566    
567     // generate the force file name
568    
569 mmeineke 561 strcpy( fileName, "DUFF.frc" );
570 mmeineke 559 // fprintf( stderr,"Trying to open %s\n", fileName );
571    
572     // attempt to open the file in the current directory first.
573    
574     frcFile = fopen( fileName, "r" );
575    
576     if( frcFile == NULL ){
577    
578     // next see if the force path enviorment variable is set
579    
580     ffPath = getenv( ffPath_env );
581     if( ffPath == NULL ) {
582     STR_DEFINE(ffPath, FRC_PATH );
583     }
584    
585    
586     strcpy( temp, ffPath );
587     strcat( temp, "/" );
588     strcat( temp, fileName );
589     strcpy( fileName, temp );
590    
591     frcFile = fopen( fileName, "r" );
592    
593     if( frcFile == NULL ){
594    
595     sprintf( painCave.errMsg,
596     "Error opening the force field parameter file: %s\n"
597     "Have you tried setting the FORCE_PARAM_PATH environment "
598     "vairable?\n",
599     fileName );
600     painCave.isFatal = 1;
601     simError();
602     }
603     }
604    
605     #ifdef IS_MPI
606     }
607    
608 mmeineke 561 sprintf( checkPointMsg, "DUFF file opened sucessfully." );
609 mmeineke 559 MPIcheckPoint();
610    
611     #endif // is_mpi
612     }
613    
614    
615 mmeineke 561 DUFF::~DUFF(){
616 mmeineke 559
617     if( headAtomType != NULL ) delete headAtomType;
618     if( headBondType != NULL ) delete headBondType;
619     if( headBendType != NULL ) delete headBendType;
620     if( headTorsionType != NULL ) delete headTorsionType;
621    
622     #ifdef IS_MPI
623     if( worldRank == 0 ){
624     #endif // is_mpi
625    
626     fclose( frcFile );
627    
628     #ifdef IS_MPI
629     }
630     #endif // is_mpi
631     }
632    
633 mmeineke 561 void DUFF::cleanMe( void ){
634 mmeineke 559
635     #ifdef IS_MPI
636    
637     // keep the linked lists in the mpi version
638    
639     #else // is_mpi
640    
641     // delete the linked lists in the single processor version
642    
643     if( headAtomType != NULL ) delete headAtomType;
644     if( headBondType != NULL ) delete headBondType;
645     if( headBendType != NULL ) delete headBendType;
646     if( headTorsionType != NULL ) delete headTorsionType;
647    
648     #endif // is_mpi
649     }
650    
651    
652 mmeineke 561 void DUFF::initForceField( int ljMixRule ){
653 mmeineke 559
654     initFortran( ljMixRule, entry_plug->useReactionField );
655     }
656    
657    
658 mmeineke 561 void DUFF::readParams( void ){
659 mmeineke 559
660     int identNum;
661    
662     atomStruct atomInfo;
663     bondStruct bondInfo;
664     bendStruct bendInfo;
665     torsionStruct torsionInfo;
666    
667     bigSigma = 0.0;
668    
669     atomInfo.last = 1;
670     bondInfo.last = 1;
671     bendInfo.last = 1;
672     torsionInfo.last = 1;
673    
674     // read in the atom info
675    
676     #ifdef IS_MPI
677     if( worldRank == 0 ){
678     #endif
679    
680     // read in the atom types.
681    
682     headAtomType = new LinkedAtomType;
683    
684     fastForward( "AtomTypes", "initializeAtoms" );
685    
686     // we are now at the AtomTypes section.
687    
688     eof_test = fgets( readLine, sizeof(readLine), frcFile );
689     lineNum++;
690    
691    
692     // read a line, and start parseing out the atom types
693    
694     if( eof_test == NULL ){
695     sprintf( painCave.errMsg,
696     "Error in reading Atoms from force file at line %d.\n",
697     lineNum );
698     painCave.isFatal = 1;
699     simError();
700     }
701    
702     identNum = 1;
703     // stop reading at end of file, or at next section
704     while( readLine[0] != '#' && eof_test != NULL ){
705    
706     // toss comment lines
707     if( readLine[0] != '!' ){
708    
709     // the parser returns 0 if the line was blank
710     if( parseAtom( readLine, lineNum, atomInfo ) ){
711     atomInfo.ident = identNum;
712     headAtomType->add( atomInfo );;
713     identNum++;
714     }
715     }
716     eof_test = fgets( readLine, sizeof(readLine), frcFile );
717     lineNum++;
718     }
719    
720     #ifdef IS_MPI
721    
722     // send out the linked list to all the other processes
723    
724     sprintf( checkPointMsg,
725 mmeineke 561 "DUFF atom structures read successfully." );
726 mmeineke 559 MPIcheckPoint();
727    
728     currentAtomType = headAtomType->next; //skip the first element who is a place holder.
729     while( currentAtomType != NULL ){
730     currentAtomType->duplicate( atomInfo );
731    
732    
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    
752     MPIcheckPoint();
753    
754     headAtomType = new LinkedAtomType;
755     recieveFrcStruct( &atomInfo, mpiAtomStructType );
756    
757     while( !atomInfo.last ){
758    
759    
760    
761     headAtomType->add( atomInfo );
762    
763     MPIcheckPoint();
764    
765     recieveFrcStruct( &atomInfo, mpiAtomStructType );
766     }
767     }
768    
769     #endif // is_mpi
770    
771    
772    
773     // call new A_types in fortran
774    
775     int isError;
776    
777     // dummy variables
778    
779     int isGB = 0;
780     int isLJ = 1;
781 chuckv 631 int isEAM =0;
782 tim 727
783 mmeineke 559 currentAtomType = headAtomType->next;;
784     while( currentAtomType != NULL ){
785    
786     if(currentAtomType->isDipole) entry_plug->useDipole = 1;
787     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 mmeineke 559 &(currentAtomType->epslon),
804     &(currentAtomType->sigma),
805     &(currentAtomType->dipole),
806     &isError );
807     if( isError ){
808     sprintf( painCave.errMsg,
809     "Error initializing the \"%s\" atom type in fortran\n",
810     currentAtomType->name );
811     painCave.isFatal = 1;
812     simError();
813     }
814     }
815     currentAtomType = currentAtomType->next;
816     }
817    
818     #ifdef IS_MPI
819     sprintf( checkPointMsg,
820 mmeineke 561 "DUFF atom structures successfully sent to fortran\n" );
821 mmeineke 559 MPIcheckPoint();
822     #endif // is_mpi
823    
824    
825    
826     // read in the bonds
827    
828     #ifdef IS_MPI
829     if( worldRank == 0 ){
830     #endif
831    
832     // read in the bond types.
833    
834     headBondType = new LinkedBondType;
835    
836     fastForward( "BondTypes", "initializeBonds" );
837    
838     // we are now at the bondTypes section
839    
840     eof_test = fgets( readLine, sizeof(readLine), frcFile );
841     lineNum++;
842    
843    
844     // read a line, and start parseing out the atom types
845    
846     if( eof_test == NULL ){
847     sprintf( painCave.errMsg,
848     "Error in reading bonds from force file at line %d.\n",
849     lineNum );
850     painCave.isFatal = 1;
851     simError();
852     }
853    
854     // stop reading at end of file, or at next section
855     while( readLine[0] != '#' && eof_test != NULL ){
856    
857     // toss comment lines
858     if( readLine[0] != '!' ){
859    
860     // the parser returns 0 if the line was blank
861     if( parseBond( readLine, lineNum, bondInfo ) ){
862     headBondType->add( bondInfo );
863     }
864     }
865     eof_test = fgets( readLine, sizeof(readLine), frcFile );
866     lineNum++;
867     }
868    
869     #ifdef IS_MPI
870    
871     // send out the linked list to all the other processes
872    
873     sprintf( checkPointMsg,
874 mmeineke 561 "DUFF bond structures read successfully." );
875 mmeineke 559 MPIcheckPoint();
876    
877     currentBondType = headBondType->next;
878     while( currentBondType != NULL ){
879     currentBondType->duplicate( bondInfo );
880     sendFrcStruct( &bondInfo, mpiBondStructType );
881     currentBondType = currentBondType->next;
882     }
883     bondInfo.last = 1;
884     sendFrcStruct( &bondInfo, mpiBondStructType );
885    
886     }
887    
888     else{
889    
890     // listen for node 0 to send out the force params
891    
892     MPIcheckPoint();
893    
894     headBondType = new LinkedBondType;
895     recieveFrcStruct( &bondInfo, mpiBondStructType );
896     while( !bondInfo.last ){
897    
898     headBondType->add( bondInfo );
899     recieveFrcStruct( &bondInfo, mpiBondStructType );
900     }
901     }
902    
903     sprintf( checkPointMsg,
904 mmeineke 561 "DUFF bond structures broadcast successfully." );
905 mmeineke 559 MPIcheckPoint();
906    
907     #endif // is_mpi
908    
909    
910     // read in the bends
911    
912     #ifdef IS_MPI
913     if( worldRank == 0 ){
914     #endif
915    
916     // read in the bend types.
917    
918     headBendType = new LinkedBendType;
919    
920     fastForward( "BendTypes", "initializeBends" );
921    
922     // we are now at the bendTypes section
923    
924     eof_test = fgets( readLine, sizeof(readLine), frcFile );
925     lineNum++;
926    
927     // read a line, and start parseing out the bend types
928    
929     if( eof_test == NULL ){
930     sprintf( painCave.errMsg,
931     "Error in reading bends from force file at line %d.\n",
932     lineNum );
933     painCave.isFatal = 1;
934     simError();
935     }
936    
937     // stop reading at end of file, or at next section
938     while( readLine[0] != '#' && eof_test != NULL ){
939    
940     // toss comment lines
941     if( readLine[0] != '!' ){
942    
943     // the parser returns 0 if the line was blank
944     if( parseBend( readLine, lineNum, bendInfo ) ){
945     headBendType->add( bendInfo );
946     }
947     }
948     eof_test = fgets( readLine, sizeof(readLine), frcFile );
949     lineNum++;
950     }
951    
952     #ifdef IS_MPI
953    
954     // send out the linked list to all the other processes
955    
956     sprintf( checkPointMsg,
957 mmeineke 561 "DUFF bend structures read successfully." );
958 mmeineke 559 MPIcheckPoint();
959    
960     currentBendType = headBendType->next;
961     while( currentBendType != NULL ){
962     currentBendType->duplicate( bendInfo );
963     sendFrcStruct( &bendInfo, mpiBendStructType );
964     currentBendType = currentBendType->next;
965     }
966     bendInfo.last = 1;
967     sendFrcStruct( &bendInfo, mpiBendStructType );
968    
969     }
970    
971     else{
972    
973     // listen for node 0 to send out the force params
974    
975     MPIcheckPoint();
976    
977     headBendType = new LinkedBendType;
978     recieveFrcStruct( &bendInfo, mpiBendStructType );
979     while( !bendInfo.last ){
980    
981     headBendType->add( bendInfo );
982     recieveFrcStruct( &bendInfo, mpiBendStructType );
983     }
984     }
985    
986     sprintf( checkPointMsg,
987 mmeineke 561 "DUFF bend structures broadcast successfully." );
988 mmeineke 559 MPIcheckPoint();
989    
990     #endif // is_mpi
991    
992    
993     // read in the torsions
994    
995     #ifdef IS_MPI
996     if( worldRank == 0 ){
997     #endif
998    
999     // read in the torsion types.
1000    
1001     headTorsionType = new LinkedTorsionType;
1002    
1003     fastForward( "TorsionTypes", "initializeTorsions" );
1004    
1005     // we are now at the torsionTypes section
1006    
1007     eof_test = fgets( readLine, sizeof(readLine), frcFile );
1008     lineNum++;
1009    
1010    
1011     // read a line, and start parseing out the atom types
1012    
1013     if( eof_test == NULL ){
1014     sprintf( painCave.errMsg,
1015     "Error in reading torsions from force file at line %d.\n",
1016     lineNum );
1017     painCave.isFatal = 1;
1018     simError();
1019     }
1020    
1021     // stop reading at end of file, or at next section
1022     while( readLine[0] != '#' && eof_test != NULL ){
1023    
1024     // toss comment lines
1025     if( readLine[0] != '!' ){
1026    
1027     // the parser returns 0 if the line was blank
1028     if( parseTorsion( readLine, lineNum, torsionInfo ) ){
1029     headTorsionType->add( torsionInfo );
1030    
1031     }
1032     }
1033     eof_test = fgets( readLine, sizeof(readLine), frcFile );
1034     lineNum++;
1035     }
1036    
1037     #ifdef IS_MPI
1038    
1039     // send out the linked list to all the other processes
1040    
1041     sprintf( checkPointMsg,
1042 mmeineke 561 "DUFF torsion structures read successfully." );
1043 mmeineke 559 MPIcheckPoint();
1044    
1045     currentTorsionType = headTorsionType->next;
1046     while( currentTorsionType != NULL ){
1047     currentTorsionType->duplicate( torsionInfo );
1048     sendFrcStruct( &torsionInfo, mpiTorsionStructType );
1049     currentTorsionType = currentTorsionType->next;
1050     }
1051     torsionInfo.last = 1;
1052     sendFrcStruct( &torsionInfo, mpiTorsionStructType );
1053    
1054     }
1055    
1056     else{
1057    
1058     // listen for node 0 to send out the force params
1059    
1060     MPIcheckPoint();
1061    
1062     headTorsionType = new LinkedTorsionType;
1063     recieveFrcStruct( &torsionInfo, mpiTorsionStructType );
1064     while( !torsionInfo.last ){
1065    
1066     headTorsionType->add( torsionInfo );
1067     recieveFrcStruct( &torsionInfo, mpiTorsionStructType );
1068     }
1069     }
1070    
1071     sprintf( checkPointMsg,
1072 mmeineke 561 "DUFF torsion structures broadcast successfully." );
1073 mmeineke 559 MPIcheckPoint();
1074    
1075     #endif // is_mpi
1076    
1077     entry_plug->useLJ = 1;
1078     }
1079    
1080    
1081    
1082 mmeineke 561 void DUFF::initializeAtoms( int nAtoms, Atom** the_atoms ){
1083 mmeineke 559
1084    
1085     //////////////////////////////////////////////////
1086     // a quick water fix
1087    
1088     double waterI[3][3];
1089     waterI[0][0] = 1.76958347772500;
1090     waterI[0][1] = 0.0;
1091     waterI[0][2] = 0.0;
1092    
1093     waterI[1][0] = 0.0;
1094     waterI[1][1] = 0.614537057924513;
1095     waterI[1][2] = 0.0;
1096    
1097     waterI[2][0] = 0.0;
1098     waterI[2][1] = 0.0;
1099     waterI[2][2] = 1.15504641980049;
1100    
1101    
1102     double headI[3][3];
1103     headI[0][0] = 1125;
1104     headI[0][1] = 0.0;
1105     headI[0][2] = 0.0;
1106    
1107     headI[1][0] = 0.0;
1108     headI[1][1] = 1125;
1109     headI[1][2] = 0.0;
1110    
1111     headI[2][0] = 0.0;
1112     headI[2][1] = 0.0;
1113     headI[2][2] = 250;
1114    
1115     //////////////////////////////////////////////////
1116    
1117    
1118     // initialize the atoms
1119    
1120     DirectionalAtom* dAtom;
1121    
1122     for(int i=0; i<nAtoms; i++ ){
1123    
1124     currentAtomType = headAtomType->find( the_atoms[i]->getType() );
1125     if( currentAtomType == NULL ){
1126     sprintf( painCave.errMsg,
1127     "AtomType error, %s not found in force file.\n",
1128     the_atoms[i]->getType() );
1129     painCave.isFatal = 1;
1130     simError();
1131     }
1132    
1133     the_atoms[i]->setMass( currentAtomType->mass );
1134     the_atoms[i]->setEpslon( currentAtomType->epslon );
1135     the_atoms[i]->setSigma( currentAtomType->sigma );
1136     the_atoms[i]->setIdent( currentAtomType->ident );
1137     the_atoms[i]->setLJ();
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->setMu( currentAtomType->dipole );
1146     dAtom->setHasDipole( 1 );
1147     dAtom->setJx( 0.0 );
1148     dAtom->setJy( 0.0 );
1149     dAtom->setJz( 0.0 );
1150    
1151     if(!strcmp("SSD",the_atoms[i]->getType())){
1152     dAtom->setI( waterI );
1153     dAtom->setSSD( 1 );
1154     }
1155     else if(!strcmp("HEAD",the_atoms[i]->getType())){
1156     dAtom->setI( headI );
1157     dAtom->setSSD( 0 );
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 mmeineke 561 "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     }