ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DUFF.cpp
Revision: 787
Committed: Thu Sep 25 19:27:15 2003 UTC (20 years, 9 months ago) by mmeineke
File size: 43418 byte(s)
Log Message:
cleaned things with gcc -Wall and g++ -Wall

File Contents

# Content
1 #include <cstdlib>
2 #include <cstdio>
3 #include <cstring>
4
5 #include <iostream>
6 using namespace std;
7
8 #include "ForceFields.hpp"
9 #include "SRI.hpp"
10 #include "simError.h"
11
12 #include "fortranWrappers.hpp"
13
14 #ifdef IS_MPI
15 #include "mpiForceField.h"
16 #endif // is_mpi
17
18
19 // define some bond Types
20
21 #define FIXED_BOND 0
22 #define HARMONIC_BOND 1
23
24
25 namespace DUFF_NS { // restrict the access of the folowing to this file only.
26
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 double v0p;
39 double rl;
40 double ru;
41 double rlp;
42 double rup;
43 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 double k0;
56 int last; // 0 -> default
57 // 1 -> tells nodes to stop listening
58 int type;
59 } 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 "Duplicate DUFF atom type \"%s\" found in "
128 "the DUFF param file./n",
129 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 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 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 info.v0p = v0p;
168 info.rl = rl;
169 info.ru = ru;
170 info.rlp = rlp;
171 info.rup = rup;
172 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 double v0p;
189 double rl;
190 double ru;
191 double rlp;
192 double rup;
193 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 "Duplicate DUFF bond type \"%s - %s\" found in "
225 "the DUFF param file./n",
226 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 next->type = info.type;
238 next->d0 = info.d0;
239 next->k0 = info.k0;
240 }
241 }
242
243 #ifdef IS_MPI
244 void duplicate( bondStruct &info ){
245 strcpy(info.nameA, nameA);
246 strcpy(info.nameB, nameB);
247 info.type = type;
248 info.d0 = d0;
249 info.k0 = k0;
250 info.last = 0;
251 }
252
253
254 #endif
255
256 char nameA[15];
257 char nameB[15];
258 int type;
259 double d0;
260 double k0;
261
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 "Duplicate DUFF bend type \"%s - %s - %s\" found in "
298 "the DUFF param file./n",
299 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 "Duplicate DUFF torsion type \"%s - %s - %s - %s\" found in "
383 "the DUFF param file./n", nameA, nameB, nameC, nameD );
384 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 using namespace DUFF_NS;
444
445
446 //****************************************************************
447 // begins the actual forcefield stuff.
448 //****************************************************************
449
450
451 DUFF::DUFF(){
452
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 int atomBC[3] = {15,11,4}; // block counts
479 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 int bondBC[3] = {30,2,2}; // block counts
501 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 strcpy( fileName, "DUFF.frc" );
570 // 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 sprintf( checkPointMsg, "DUFF file opened sucessfully." );
609 MPIcheckPoint();
610
611 #endif // is_mpi
612 }
613
614
615 DUFF::~DUFF(){
616
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 void DUFF::cleanMe( void ){
634
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 void DUFF::initForceField( int ljMixRule ){
653
654 initFortran( ljMixRule, entry_plug->useReactionField );
655 }
656
657
658 void DUFF::readParams( void ){
659
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 "DUFF atom structures read successfully." );
726 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 "successfully sent DUFF force type: \"%s\"\n",
738 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 int isEAM =0;
782
783 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 set_sticky_params( &(currentAtomType->w0), &(currentAtomType->v0),
790 &(currentAtomType->v0p),
791 &(currentAtomType->rl), &(currentAtomType->ru),
792 &(currentAtomType->rlp), &(currentAtomType->rup));
793 }
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 &isEAM,
803 &(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 "DUFF atom structures successfully sent to fortran\n" );
821 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 "DUFF bond structures read successfully." );
875 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 "DUFF bond structures broadcast successfully." );
905 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 "DUFF bend structures read successfully." );
958 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 "DUFF bend structures broadcast successfully." );
988 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 "DUFF torsion structures read successfully." );
1043 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 "DUFF torsion structures broadcast successfully." );
1073 MPIcheckPoint();
1074
1075 #endif // is_mpi
1076
1077 entry_plug->useLJ = 1;
1078 }
1079
1080
1081
1082 void DUFF::initializeAtoms( int nAtoms, Atom** the_atoms ){
1083
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 "DUFF error: Atom \"%s\" is a dipole, yet no standard"
1172 " 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 "DUFF error: Atom \"%s\" was given a standard"
1182 "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 void DUFF::initializeBonds( int nBonds, Bond** bondArray,
1192 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 switch( currentBondType->type ){
1220
1221 case FIXED_BOND:
1222
1223 bondArray[i] = new ConstrainedBond( *the_atoms[a],
1224 *the_atoms[b],
1225 currentBondType->d0 );
1226 entry_plug->n_constraints++;
1227 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 }
1242 }
1243 }
1244
1245 void DUFF::initializeBends( int nBends, Bend** bendArray,
1246 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 *the_atoms[b]);
1306
1307 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 }
1323 }
1324 }
1325 }
1326
1327 void DUFF::initializeTorsions( int nTorsions, Torsion** torsionArray,
1328 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 void DUFF::fastForward( char* stopText, char* searchOwner ){
1375
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 int DUFF_NS::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
1430
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 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 }
1560 else info.v0 = info.w0 = info.v0p = info.rl = info.ru = info.rlp = info.rup = 0.0;
1561
1562 return 1;
1563 }
1564 else return 0;
1565 }
1566
1567 int DUFF_NS::parseBond( char *lineBuffer, int lineNum, bondStruct &info ){
1568
1569 char* the_token;
1570 char bondType[30];
1571
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 strcpy( bondType, the_token );
1594
1595 if( !strcmp( bondType, "fixed" ) ){
1596 info.type = FIXED_BOND;
1597
1598 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1599 sprintf( painCave.errMsg,
1600 "Error parseing BondTypes: line %d\n", lineNum );
1601 painCave.isFatal = 1;
1602 simError();
1603 }
1604
1605 info.d0 = atof( the_token );
1606
1607 info.k0=0.0;
1608 }
1609 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 else{
1632 sprintf( painCave.errMsg,
1633 "Unknown DUFF bond type \"%s\" at line %d\n",
1634 bondType,
1635 lineNum );
1636 painCave.isFatal = 1;
1637 simError();
1638 }
1639
1640 return 1;
1641 }
1642 else return 0;
1643 }
1644
1645
1646 int DUFF_NS::parseBend( char *lineBuffer, int lineNum, bendStruct &info ){
1647
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 "Unknown DUFF bend type \"%s\" at line %d\n",
1723 info.type,
1724 lineNum );
1725 painCave.isFatal = 1;
1726 simError();
1727 }
1728
1729 return 1;
1730 }
1731 else return 0;
1732 }
1733
1734 int DUFF_NS::parseTorsion( char *lineBuffer, int lineNum, torsionStruct &info ){
1735
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 "Unknown DUFF torsion type \"%s\" at line %d\n",
1821 info.type,
1822 lineNum );
1823 painCave.isFatal = 1;
1824 simError();
1825 }
1826
1827 return 1;
1828 }
1829
1830 else return 0;
1831 }