ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DUFF.cpp
Revision: 594
Committed: Fri Jul 11 22:34:48 2003 UTC (20 years, 11 months ago) by mmeineke
File size: 41631 byte(s)
Log Message:
working on som integrator bugs

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