ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DUFF.cpp
Revision: 635
Committed: Thu Jul 17 20:32:24 2003 UTC (20 years, 11 months ago) by gezelter
File size: 43496 byte(s)
Log Message: 
Changes for SSD/E

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 char errMsg[1000];
458
459 headAtomType = NULL;
460 currentAtomType = NULL;
461 headBondType = NULL;
462 currentBondType = NULL;
463 headBendType = NULL;
464 currentBendType = NULL;
465 headTorsionType = NULL;
466 currentTorsionType = NULL;
467
468 // do the funtion wrapping
469 wrapMeFF( this );
470
471
472 #ifdef IS_MPI
473 int i;
474
475 // **********************************************************************
476 // Init the atomStruct mpi type
477
478 atomStruct atomProto; // mpiPrototype
479 int atomBC[3] = {15,6,4}; // block counts
480 MPI_Aint atomDspls[3]; // displacements
481 MPI_Datatype atomMbrTypes[3]; // member mpi types
482
483 MPI_Address(&atomProto.name, &atomDspls[0]);
484 MPI_Address(&atomProto.mass, &atomDspls[1]);
485 MPI_Address(&atomProto.isSSD, &atomDspls[2]);
486
487 atomMbrTypes[0] = MPI_CHAR;
488 atomMbrTypes[1] = MPI_DOUBLE;
489 atomMbrTypes[2] = MPI_INT;
490
491 for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
492
493 MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
494 MPI_Type_commit(&mpiAtomStructType);
495
496
497 // **********************************************************************
498 // Init the bondStruct mpi type
499
500 bondStruct bondProto; // mpiPrototype
501 int bondBC[3] = {30,2,2}; // block counts
502 MPI_Aint bondDspls[3]; // displacements
503 MPI_Datatype bondMbrTypes[3]; // member mpi types
504
505 MPI_Address(&bondProto.nameA, &bondDspls[0]);
506 MPI_Address(&bondProto.d0, &bondDspls[1]);
507 MPI_Address(&bondProto.last, &bondDspls[2]);
508
509 bondMbrTypes[0] = MPI_CHAR;
510 bondMbrTypes[1] = MPI_DOUBLE;
511 bondMbrTypes[2] = MPI_INT;
512
513 for (i=2; i >= 0; i--) bondDspls[i] -= bondDspls[0];
514
515 MPI_Type_struct(3, bondBC, bondDspls, bondMbrTypes, &mpiBondStructType);
516 MPI_Type_commit(&mpiBondStructType);
517
518
519 // **********************************************************************
520 // Init the bendStruct mpi type
521
522 bendStruct bendProto; // mpiPrototype
523 int bendBC[3] = {75,4,1}; // block counts
524 MPI_Aint bendDspls[3]; // displacements
525 MPI_Datatype bendMbrTypes[3]; // member mpi types
526
527 MPI_Address(&bendProto.nameA, &bendDspls[0]);
528 MPI_Address(&bendProto.k1, &bendDspls[1]);
529 MPI_Address(&bendProto.last, &bendDspls[2]);
530
531 bendMbrTypes[0] = MPI_CHAR;
532 bendMbrTypes[1] = MPI_DOUBLE;
533 bendMbrTypes[2] = MPI_INT;
534
535 for (i=2; i >= 0; i--) bendDspls[i] -= bendDspls[0];
536
537 MPI_Type_struct(3, bendBC, bendDspls, bendMbrTypes, &mpiBendStructType);
538 MPI_Type_commit(&mpiBendStructType);
539
540
541 // **********************************************************************
542 // Init the torsionStruct mpi type
543
544 torsionStruct torsionProto; // mpiPrototype
545 int torsionBC[3] = {90,4,1}; // block counts
546 MPI_Aint torsionDspls[3]; // displacements
547 MPI_Datatype torsionMbrTypes[3]; // member mpi types
548
549 MPI_Address(&torsionProto.nameA, &torsionDspls[0]);
550 MPI_Address(&torsionProto.k1, &torsionDspls[1]);
551 MPI_Address(&torsionProto.last, &torsionDspls[2]);
552
553 torsionMbrTypes[0] = MPI_CHAR;
554 torsionMbrTypes[1] = MPI_DOUBLE;
555 torsionMbrTypes[2] = MPI_INT;
556
557 for (i=2; i >= 0; i--) torsionDspls[i] -= torsionDspls[0];
558
559 MPI_Type_struct(3, torsionBC, torsionDspls, torsionMbrTypes,
560 &mpiTorsionStructType);
561 MPI_Type_commit(&mpiTorsionStructType);
562
563 // ***********************************************************************
564
565 if( worldRank == 0 ){
566 #endif
567
568 // generate the force file name
569
570 strcpy( fileName, "DUFF.frc" );
571 // fprintf( stderr,"Trying to open %s\n", fileName );
572
573 // attempt to open the file in the current directory first.
574
575 frcFile = fopen( fileName, "r" );
576
577 if( frcFile == NULL ){
578
579 // next see if the force path enviorment variable is set
580
581 ffPath = getenv( ffPath_env );
582 if( ffPath == NULL ) {
583 STR_DEFINE(ffPath, FRC_PATH );
584 }
585
586
587 strcpy( temp, ffPath );
588 strcat( temp, "/" );
589 strcat( temp, fileName );
590 strcpy( fileName, temp );
591
592 frcFile = fopen( fileName, "r" );
593
594 if( frcFile == NULL ){
595
596 sprintf( painCave.errMsg,
597 "Error opening the force field parameter file: %s\n"
598 "Have you tried setting the FORCE_PARAM_PATH environment "
599 "vairable?\n",
600 fileName );
601 painCave.isFatal = 1;
602 simError();
603 }
604 }
605
606 #ifdef IS_MPI
607 }
608
609 sprintf( checkPointMsg, "DUFF file opened sucessfully." );
610 MPIcheckPoint();
611
612 #endif // is_mpi
613 }
614
615
616 DUFF::~DUFF(){
617
618 if( headAtomType != NULL ) delete headAtomType;
619 if( headBondType != NULL ) delete headBondType;
620 if( headBendType != NULL ) delete headBendType;
621 if( headTorsionType != NULL ) delete headTorsionType;
622
623 #ifdef IS_MPI
624 if( worldRank == 0 ){
625 #endif // is_mpi
626
627 fclose( frcFile );
628
629 #ifdef IS_MPI
630 }
631 #endif // is_mpi
632 }
633
634 void DUFF::cleanMe( void ){
635
636 #ifdef IS_MPI
637
638 // keep the linked lists in the mpi version
639
640 #else // is_mpi
641
642 // delete the linked lists in the single processor version
643
644 if( headAtomType != NULL ) delete headAtomType;
645 if( headBondType != NULL ) delete headBondType;
646 if( headBendType != NULL ) delete headBendType;
647 if( headTorsionType != NULL ) delete headTorsionType;
648
649 #endif // is_mpi
650 }
651
652
653 void DUFF::initForceField( int ljMixRule ){
654
655 initFortran( ljMixRule, entry_plug->useReactionField );
656 }
657
658
659 void DUFF::readParams( void ){
660
661 int i, a, b, c, d;
662 int identNum;
663 char* atomA;
664 char* atomB;
665 char* atomC;
666 char* atomD;
667
668 atomStruct atomInfo;
669 bondStruct bondInfo;
670 bendStruct bendInfo;
671 torsionStruct torsionInfo;
672
673 bigSigma = 0.0;
674
675 atomInfo.last = 1;
676 bondInfo.last = 1;
677 bendInfo.last = 1;
678 torsionInfo.last = 1;
679
680 // read in the atom info
681
682 #ifdef IS_MPI
683 if( worldRank == 0 ){
684 #endif
685
686 // read in the atom types.
687
688 headAtomType = new LinkedAtomType;
689
690 fastForward( "AtomTypes", "initializeAtoms" );
691
692 // we are now at the AtomTypes section.
693
694 eof_test = fgets( readLine, sizeof(readLine), frcFile );
695 lineNum++;
696
697
698 // read a line, and start parseing out the atom types
699
700 if( eof_test == NULL ){
701 sprintf( painCave.errMsg,
702 "Error in reading Atoms from force file at line %d.\n",
703 lineNum );
704 painCave.isFatal = 1;
705 simError();
706 }
707
708 identNum = 1;
709 // stop reading at end of file, or at next section
710 while( readLine[0] != '#' && eof_test != NULL ){
711
712 // toss comment lines
713 if( readLine[0] != '!' ){
714
715 // the parser returns 0 if the line was blank
716 if( parseAtom( readLine, lineNum, atomInfo ) ){
717 atomInfo.ident = identNum;
718 headAtomType->add( atomInfo );;
719 identNum++;
720 }
721 }
722 eof_test = fgets( readLine, sizeof(readLine), frcFile );
723 lineNum++;
724 }
725
726 #ifdef IS_MPI
727
728 // send out the linked list to all the other processes
729
730 sprintf( checkPointMsg,
731 "DUFF atom structures read successfully." );
732 MPIcheckPoint();
733
734 currentAtomType = headAtomType->next; //skip the first element who is a place holder.
735 while( currentAtomType != NULL ){
736 currentAtomType->duplicate( atomInfo );
737
738
739
740 sendFrcStruct( &atomInfo, mpiAtomStructType );
741
742 sprintf( checkPointMsg,
743 "successfully sent DUFF force type: \"%s\"\n",
744 atomInfo.name );
745 MPIcheckPoint();
746
747 currentAtomType = currentAtomType->next;
748 }
749 atomInfo.last = 1;
750 sendFrcStruct( &atomInfo, mpiAtomStructType );
751
752 }
753
754 else{
755
756 // listen for node 0 to send out the force params
757
758 MPIcheckPoint();
759
760 headAtomType = new LinkedAtomType;
761 recieveFrcStruct( &atomInfo, mpiAtomStructType );
762
763 while( !atomInfo.last ){
764
765
766
767 headAtomType->add( atomInfo );
768
769 MPIcheckPoint();
770
771 recieveFrcStruct( &atomInfo, mpiAtomStructType );
772 }
773 }
774
775 #endif // is_mpi
776
777
778
779 // call new A_types in fortran
780
781 int isError;
782
783 // dummy variables
784
785 int isGB = 0;
786 int isLJ = 1;
787 int isEAM =0;
788 double GB_dummy = 0.0;
789
790
791 currentAtomType = headAtomType->next;;
792 while( currentAtomType != NULL ){
793
794 if(currentAtomType->isDipole) entry_plug->useDipole = 1;
795 if(currentAtomType->isSSD) {
796 entry_plug->useSticky = 1;
797 set_sticky_params( &(currentAtomType->w0), &(currentAtomType->v0),
798 &(currentAtomType->v0p),
799 &(currentAtomType->rl), &(currentAtomType->ru),
800 &(currentAtomType->rlp), &(currentAtomType->rup));
801 }
802
803 if( currentAtomType->name[0] != '\0' ){
804 isError = 0;
805 makeAtype( &(currentAtomType->ident),
806 &isLJ,
807 &(currentAtomType->isSSD),
808 &(currentAtomType->isDipole),
809 &isGB,
810 &isEAM,
811 &(currentAtomType->epslon),
812 &(currentAtomType->sigma),
813 &(currentAtomType->dipole),
814 &isError );
815 if( isError ){
816 sprintf( painCave.errMsg,
817 "Error initializing the \"%s\" atom type in fortran\n",
818 currentAtomType->name );
819 painCave.isFatal = 1;
820 simError();
821 }
822 }
823 currentAtomType = currentAtomType->next;
824 }
825
826 #ifdef IS_MPI
827 sprintf( checkPointMsg,
828 "DUFF atom structures successfully sent to fortran\n" );
829 MPIcheckPoint();
830 #endif // is_mpi
831
832
833
834 // read in the bonds
835
836 #ifdef IS_MPI
837 if( worldRank == 0 ){
838 #endif
839
840 // read in the bond types.
841
842 headBondType = new LinkedBondType;
843
844 fastForward( "BondTypes", "initializeBonds" );
845
846 // we are now at the bondTypes section
847
848 eof_test = fgets( readLine, sizeof(readLine), frcFile );
849 lineNum++;
850
851
852 // read a line, and start parseing out the atom types
853
854 if( eof_test == NULL ){
855 sprintf( painCave.errMsg,
856 "Error in reading bonds from force file at line %d.\n",
857 lineNum );
858 painCave.isFatal = 1;
859 simError();
860 }
861
862 // stop reading at end of file, or at next section
863 while( readLine[0] != '#' && eof_test != NULL ){
864
865 // toss comment lines
866 if( readLine[0] != '!' ){
867
868 // the parser returns 0 if the line was blank
869 if( parseBond( readLine, lineNum, bondInfo ) ){
870 headBondType->add( bondInfo );
871 }
872 }
873 eof_test = fgets( readLine, sizeof(readLine), frcFile );
874 lineNum++;
875 }
876
877 #ifdef IS_MPI
878
879 // send out the linked list to all the other processes
880
881 sprintf( checkPointMsg,
882 "DUFF bond structures read successfully." );
883 MPIcheckPoint();
884
885 currentBondType = headBondType->next;
886 while( currentBondType != NULL ){
887 currentBondType->duplicate( bondInfo );
888 sendFrcStruct( &bondInfo, mpiBondStructType );
889 currentBondType = currentBondType->next;
890 }
891 bondInfo.last = 1;
892 sendFrcStruct( &bondInfo, mpiBondStructType );
893
894 }
895
896 else{
897
898 // listen for node 0 to send out the force params
899
900 MPIcheckPoint();
901
902 headBondType = new LinkedBondType;
903 recieveFrcStruct( &bondInfo, mpiBondStructType );
904 while( !bondInfo.last ){
905
906 headBondType->add( bondInfo );
907 recieveFrcStruct( &bondInfo, mpiBondStructType );
908 }
909 }
910
911 sprintf( checkPointMsg,
912 "DUFF bond structures broadcast successfully." );
913 MPIcheckPoint();
914
915 #endif // is_mpi
916
917
918 // read in the bends
919
920 #ifdef IS_MPI
921 if( worldRank == 0 ){
922 #endif
923
924 // read in the bend types.
925
926 headBendType = new LinkedBendType;
927
928 fastForward( "BendTypes", "initializeBends" );
929
930 // we are now at the bendTypes section
931
932 eof_test = fgets( readLine, sizeof(readLine), frcFile );
933 lineNum++;
934
935 // read a line, and start parseing out the bend types
936
937 if( eof_test == NULL ){
938 sprintf( painCave.errMsg,
939 "Error in reading bends from force file at line %d.\n",
940 lineNum );
941 painCave.isFatal = 1;
942 simError();
943 }
944
945 // stop reading at end of file, or at next section
946 while( readLine[0] != '#' && eof_test != NULL ){
947
948 // toss comment lines
949 if( readLine[0] != '!' ){
950
951 // the parser returns 0 if the line was blank
952 if( parseBend( readLine, lineNum, bendInfo ) ){
953 headBendType->add( bendInfo );
954 }
955 }
956 eof_test = fgets( readLine, sizeof(readLine), frcFile );
957 lineNum++;
958 }
959
960 #ifdef IS_MPI
961
962 // send out the linked list to all the other processes
963
964 sprintf( checkPointMsg,
965 "DUFF bend structures read successfully." );
966 MPIcheckPoint();
967
968 currentBendType = headBendType->next;
969 while( currentBendType != NULL ){
970 currentBendType->duplicate( bendInfo );
971 sendFrcStruct( &bendInfo, mpiBendStructType );
972 currentBendType = currentBendType->next;
973 }
974 bendInfo.last = 1;
975 sendFrcStruct( &bendInfo, mpiBendStructType );
976
977 }
978
979 else{
980
981 // listen for node 0 to send out the force params
982
983 MPIcheckPoint();
984
985 headBendType = new LinkedBendType;
986 recieveFrcStruct( &bendInfo, mpiBendStructType );
987 while( !bendInfo.last ){
988
989 headBendType->add( bendInfo );
990 recieveFrcStruct( &bendInfo, mpiBendStructType );
991 }
992 }
993
994 sprintf( checkPointMsg,
995 "DUFF bend structures broadcast successfully." );
996 MPIcheckPoint();
997
998 #endif // is_mpi
999
1000
1001 // read in the torsions
1002
1003 #ifdef IS_MPI
1004 if( worldRank == 0 ){
1005 #endif
1006
1007 // read in the torsion types.
1008
1009 headTorsionType = new LinkedTorsionType;
1010
1011 fastForward( "TorsionTypes", "initializeTorsions" );
1012
1013 // we are now at the torsionTypes section
1014
1015 eof_test = fgets( readLine, sizeof(readLine), frcFile );
1016 lineNum++;
1017
1018
1019 // read a line, and start parseing out the atom types
1020
1021 if( eof_test == NULL ){
1022 sprintf( painCave.errMsg,
1023 "Error in reading torsions from force file at line %d.\n",
1024 lineNum );
1025 painCave.isFatal = 1;
1026 simError();
1027 }
1028
1029 // stop reading at end of file, or at next section
1030 while( readLine[0] != '#' && eof_test != NULL ){
1031
1032 // toss comment lines
1033 if( readLine[0] != '!' ){
1034
1035 // the parser returns 0 if the line was blank
1036 if( parseTorsion( readLine, lineNum, torsionInfo ) ){
1037 headTorsionType->add( torsionInfo );
1038
1039 }
1040 }
1041 eof_test = fgets( readLine, sizeof(readLine), frcFile );
1042 lineNum++;
1043 }
1044
1045 #ifdef IS_MPI
1046
1047 // send out the linked list to all the other processes
1048
1049 sprintf( checkPointMsg,
1050 "DUFF torsion structures read successfully." );
1051 MPIcheckPoint();
1052
1053 currentTorsionType = headTorsionType->next;
1054 while( currentTorsionType != NULL ){
1055 currentTorsionType->duplicate( torsionInfo );
1056 sendFrcStruct( &torsionInfo, mpiTorsionStructType );
1057 currentTorsionType = currentTorsionType->next;
1058 }
1059 torsionInfo.last = 1;
1060 sendFrcStruct( &torsionInfo, mpiTorsionStructType );
1061
1062 }
1063
1064 else{
1065
1066 // listen for node 0 to send out the force params
1067
1068 MPIcheckPoint();
1069
1070 headTorsionType = new LinkedTorsionType;
1071 recieveFrcStruct( &torsionInfo, mpiTorsionStructType );
1072 while( !torsionInfo.last ){
1073
1074 headTorsionType->add( torsionInfo );
1075 recieveFrcStruct( &torsionInfo, mpiTorsionStructType );
1076 }
1077 }
1078
1079 sprintf( checkPointMsg,
1080 "DUFF torsion structures broadcast successfully." );
1081 MPIcheckPoint();
1082
1083 #endif // is_mpi
1084
1085 entry_plug->useLJ = 1;
1086 }
1087
1088
1089
1090 void DUFF::initializeAtoms( int nAtoms, Atom** the_atoms ){
1091
1092
1093 //////////////////////////////////////////////////
1094 // a quick water fix
1095
1096 double waterI[3][3];
1097 waterI[0][0] = 1.76958347772500;
1098 waterI[0][1] = 0.0;
1099 waterI[0][2] = 0.0;
1100
1101 waterI[1][0] = 0.0;
1102 waterI[1][1] = 0.614537057924513;
1103 waterI[1][2] = 0.0;
1104
1105 waterI[2][0] = 0.0;
1106 waterI[2][1] = 0.0;
1107 waterI[2][2] = 1.15504641980049;
1108
1109
1110 double headI[3][3];
1111 headI[0][0] = 1125;
1112 headI[0][1] = 0.0;
1113 headI[0][2] = 0.0;
1114
1115 headI[1][0] = 0.0;
1116 headI[1][1] = 1125;
1117 headI[1][2] = 0.0;
1118
1119 headI[2][0] = 0.0;
1120 headI[2][1] = 0.0;
1121 headI[2][2] = 250;
1122
1123 //////////////////////////////////////////////////
1124
1125
1126 // initialize the atoms
1127
1128 DirectionalAtom* dAtom;
1129
1130 for(int i=0; i<nAtoms; i++ ){
1131
1132 currentAtomType = headAtomType->find( the_atoms[i]->getType() );
1133 if( currentAtomType == NULL ){
1134 sprintf( painCave.errMsg,
1135 "AtomType error, %s not found in force file.\n",
1136 the_atoms[i]->getType() );
1137 painCave.isFatal = 1;
1138 simError();
1139 }
1140
1141 the_atoms[i]->setMass( currentAtomType->mass );
1142 the_atoms[i]->setEpslon( currentAtomType->epslon );
1143 the_atoms[i]->setSigma( currentAtomType->sigma );
1144 the_atoms[i]->setIdent( currentAtomType->ident );
1145 the_atoms[i]->setLJ();
1146
1147 if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
1148
1149 if( currentAtomType->isDipole ){
1150 if( the_atoms[i]->isDirectional() ){
1151
1152 dAtom = (DirectionalAtom *) the_atoms[i];
1153 dAtom->setMu( currentAtomType->dipole );
1154 dAtom->setHasDipole( 1 );
1155 dAtom->setJx( 0.0 );
1156 dAtom->setJy( 0.0 );
1157 dAtom->setJz( 0.0 );
1158
1159 if(!strcmp("SSD",the_atoms[i]->getType())){
1160 dAtom->setI( waterI );
1161 dAtom->setSSD( 1 );
1162 }
1163 else if(!strcmp("HEAD",the_atoms[i]->getType())){
1164 dAtom->setI( headI );
1165 dAtom->setSSD( 0 );
1166 }
1167 else{
1168 sprintf(painCave.errMsg,
1169 "AtmType error, %s does not have a moment of inertia set.\n",
1170 the_atoms[i]->getType() );
1171 painCave.isFatal = 1;
1172 simError();
1173 }
1174 entry_plug->n_dipoles++;
1175 }
1176 else{
1177
1178 sprintf( painCave.errMsg,
1179 "DUFF error: Atom \"%s\" is a dipole, yet no standard"
1180 " orientation was specifed in the BASS file.\n",
1181 currentAtomType->name );
1182 painCave.isFatal = 1;
1183 simError();
1184 }
1185 }
1186 else{
1187 if( the_atoms[i]->isDirectional() ){
1188 sprintf( painCave.errMsg,
1189 "DUFF error: Atom \"%s\" was given a standard"
1190 "orientation in the BASS file, yet it is not a dipole.\n",
1191 currentAtomType->name);
1192 painCave.isFatal = 1;
1193 simError();
1194 }
1195 }
1196 }
1197 }
1198
1199 void DUFF::initializeBonds( int nBonds, Bond** bondArray,
1200 bond_pair* the_bonds ){
1201 int i,a,b;
1202 char* atomA;
1203 char* atomB;
1204
1205 Atom** the_atoms;
1206 the_atoms = entry_plug->atoms;
1207
1208
1209 // initialize the Bonds
1210
1211 for( i=0; i<nBonds; i++ ){
1212
1213 a = the_bonds[i].a;
1214 b = the_bonds[i].b;
1215
1216 atomA = the_atoms[a]->getType();
1217 atomB = the_atoms[b]->getType();
1218 currentBondType = headBondType->find( atomA, atomB );
1219 if( currentBondType == NULL ){
1220 sprintf( painCave.errMsg,
1221 "BondType error, %s - %s not found in force file.\n",
1222 atomA, atomB );
1223 painCave.isFatal = 1;
1224 simError();
1225 }
1226
1227 switch( currentBondType->type ){
1228
1229 case FIXED_BOND:
1230
1231 bondArray[i] = new ConstrainedBond( *the_atoms[a],
1232 *the_atoms[b],
1233 currentBondType->d0 );
1234 entry_plug->n_constraints++;
1235 break;
1236
1237 case HARMONIC_BOND:
1238
1239 bondArray[i] = new HarmonicBond( *the_atoms[a],
1240 *the_atoms[b],
1241 currentBondType->d0,
1242 currentBondType->k0 );
1243 break;
1244
1245 default:
1246
1247 break;
1248 // do nothing
1249 }
1250 }
1251 }
1252
1253 void DUFF::initializeBends( int nBends, Bend** bendArray,
1254 bend_set* the_bends ){
1255
1256 QuadraticBend* qBend;
1257 GhostBend* gBend;
1258 Atom** the_atoms;
1259 the_atoms = entry_plug->atoms;
1260
1261 int i, a, b, c;
1262 char* atomA;
1263 char* atomB;
1264 char* atomC;
1265
1266 // initialize the Bends
1267
1268 for( i=0; i<nBends; i++ ){
1269
1270 a = the_bends[i].a;
1271 b = the_bends[i].b;
1272 c = the_bends[i].c;
1273
1274 atomA = the_atoms[a]->getType();
1275 atomB = the_atoms[b]->getType();
1276
1277 if( the_bends[i].isGhost ) atomC = "GHOST";
1278 else atomC = the_atoms[c]->getType();
1279
1280 currentBendType = headBendType->find( atomA, atomB, atomC );
1281 if( currentBendType == NULL ){
1282 sprintf( painCave.errMsg, "BendType error, %s - %s - %s not found"
1283 " in force file.\n",
1284 atomA, atomB, atomC );
1285 painCave.isFatal = 1;
1286 simError();
1287 }
1288
1289 if( !strcmp( currentBendType->type, "quadratic" ) ){
1290
1291 if( the_bends[i].isGhost){
1292
1293 if( the_bends[i].ghost == b ){
1294 // do nothing
1295 }
1296 else if( the_bends[i].ghost == a ){
1297 c = a;
1298 a = b;
1299 b = c;
1300 }
1301 else{
1302 sprintf( painCave.errMsg,
1303 "BendType error, %s - %s - %s,\n"
1304 " --> central atom is not "
1305 "correctly identified with the "
1306 "\"ghostVectorSource = \" tag.\n",
1307 atomA, atomB, atomC );
1308 painCave.isFatal = 1;
1309 simError();
1310 }
1311
1312 gBend = new GhostBend( *the_atoms[a],
1313 *the_atoms[b] );
1314 gBend->setConstants( currentBendType->k1,
1315 currentBendType->k2,
1316 currentBendType->k3,
1317 currentBendType->t0 );
1318 bendArray[i] = gBend;
1319 }
1320 else{
1321 qBend = new QuadraticBend( *the_atoms[a],
1322 *the_atoms[b],
1323 *the_atoms[c] );
1324 qBend->setConstants( currentBendType->k1,
1325 currentBendType->k2,
1326 currentBendType->k3,
1327 currentBendType->t0 );
1328 bendArray[i] = qBend;
1329 }
1330 }
1331 }
1332 }
1333
1334 void DUFF::initializeTorsions( int nTorsions, Torsion** torsionArray,
1335 torsion_set* the_torsions ){
1336
1337 int i, a, b, c, d;
1338 char* atomA;
1339 char* atomB;
1340 char* atomC;
1341 char* atomD;
1342
1343 CubicTorsion* cTors;
1344 Atom** the_atoms;
1345 the_atoms = entry_plug->atoms;
1346
1347 // initialize the Torsions
1348
1349 for( i=0; i<nTorsions; i++ ){
1350
1351 a = the_torsions[i].a;
1352 b = the_torsions[i].b;
1353 c = the_torsions[i].c;
1354 d = the_torsions[i].d;
1355
1356 atomA = the_atoms[a]->getType();
1357 atomB = the_atoms[b]->getType();
1358 atomC = the_atoms[c]->getType();
1359 atomD = the_atoms[d]->getType();
1360 currentTorsionType = headTorsionType->find( atomA, atomB, atomC, atomD );
1361 if( currentTorsionType == NULL ){
1362 sprintf( painCave.errMsg,
1363 "TorsionType error, %s - %s - %s - %s not found"
1364 " in force file.\n",
1365 atomA, atomB, atomC, atomD );
1366 painCave.isFatal = 1;
1367 simError();
1368 }
1369
1370 if( !strcmp( currentTorsionType->type, "cubic" ) ){
1371
1372 cTors = new CubicTorsion( *the_atoms[a], *the_atoms[b],
1373 *the_atoms[c], *the_atoms[d] );
1374 cTors->setConstants( currentTorsionType->k1, currentTorsionType->k2,
1375 currentTorsionType->k3, currentTorsionType->k4 );
1376 torsionArray[i] = cTors;
1377 }
1378 }
1379 }
1380
1381 void DUFF::fastForward( char* stopText, char* searchOwner ){
1382
1383 int foundText = 0;
1384 char* the_token;
1385
1386 rewind( frcFile );
1387 lineNum = 0;
1388
1389 eof_test = fgets( readLine, sizeof(readLine), frcFile );
1390 lineNum++;
1391 if( eof_test == NULL ){
1392 sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
1393 " file is empty.\n",
1394 searchOwner );
1395 painCave.isFatal = 1;
1396 simError();
1397 }
1398
1399
1400 while( !foundText ){
1401 while( eof_test != NULL && readLine[0] != '#' ){
1402 eof_test = fgets( readLine, sizeof(readLine), frcFile );
1403 lineNum++;
1404 }
1405 if( eof_test == NULL ){
1406 sprintf( painCave.errMsg,
1407 "Error fast forwarding force file for %s at "
1408 "line %d: file ended unexpectedly.\n",
1409 searchOwner,
1410 lineNum );
1411 painCave.isFatal = 1;
1412 simError();
1413 }
1414
1415 the_token = strtok( readLine, " ,;\t#\n" );
1416 foundText = !strcmp( stopText, the_token );
1417
1418 if( !foundText ){
1419 eof_test = fgets( readLine, sizeof(readLine), frcFile );
1420 lineNum++;
1421
1422 if( eof_test == NULL ){
1423 sprintf( painCave.errMsg,
1424 "Error fast forwarding force file for %s at "
1425 "line %d: file ended unexpectedly.\n",
1426 searchOwner,
1427 lineNum );
1428 painCave.isFatal = 1;
1429 simError();
1430 }
1431 }
1432 }
1433 }
1434
1435
1436 int DUFF_NS::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
1437
1438 char* the_token;
1439
1440 the_token = strtok( lineBuffer, " \n\t,;" );
1441 if( the_token != NULL ){
1442
1443 strcpy( info.name, the_token );
1444
1445 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1446 sprintf( painCave.errMsg,
1447 "Error parseing AtomTypes: line %d\n", lineNum );
1448 painCave.isFatal = 1;
1449 simError();
1450 }
1451
1452 info.isDipole = atoi( the_token );
1453
1454 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1455 sprintf( painCave.errMsg,
1456 "Error parseing AtomTypes: line %d\n", lineNum );
1457 painCave.isFatal = 1;
1458 simError();
1459 }
1460
1461 info.isSSD = atoi( the_token );
1462
1463 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1464 sprintf( painCave.errMsg,
1465 "Error parseing AtomTypes: line %d\n", lineNum );
1466 painCave.isFatal = 1;
1467 simError();
1468 }
1469
1470 info.mass = atof( the_token );
1471
1472 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1473 sprintf( painCave.errMsg,
1474 "Error parseing AtomTypes: line %d\n", lineNum );
1475 painCave.isFatal = 1;
1476 simError();
1477 }
1478
1479 info.epslon = atof( the_token );
1480
1481 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1482 sprintf( painCave.errMsg,
1483 "Error parseing AtomTypes: line %d\n", lineNum );
1484 painCave.isFatal = 1;
1485 simError();
1486 }
1487
1488 info.sigma = atof( the_token );
1489
1490 if( info.isDipole ){
1491
1492 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1493 sprintf( painCave.errMsg,
1494 "Error parseing AtomTypes: line %d\n", lineNum );
1495 painCave.isFatal = 1;
1496 simError();
1497 }
1498
1499 info.dipole = atof( the_token );
1500 }
1501 else info.dipole = 0.0;
1502
1503 if( info.isSSD ){
1504
1505 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1506 sprintf( painCave.errMsg,
1507 "Error parseing AtomTypes: line %d\n", lineNum );
1508 painCave.isFatal = 1;
1509 simError();
1510 }
1511
1512 info.w0 = atof( the_token );
1513
1514 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1515 sprintf( painCave.errMsg,
1516 "Error parseing AtomTypes: line %d\n", lineNum );
1517 painCave.isFatal = 1;
1518 simError();
1519 }
1520
1521 info.v0 = atof( the_token );
1522 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1523 sprintf( painCave.errMsg,
1524 "Error parseing AtomTypes: line %d\n", lineNum );
1525 painCave.isFatal = 1;
1526 simError();
1527 }
1528
1529 info.v0p = atof( the_token );
1530
1531 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1532 sprintf( painCave.errMsg,
1533 "Error parseing AtomTypes: line %d\n", lineNum );
1534 painCave.isFatal = 1;
1535 simError();
1536 }
1537
1538 info.rl = atof( the_token );
1539
1540 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1541 sprintf( painCave.errMsg,
1542 "Error parseing AtomTypes: line %d\n", lineNum );
1543 painCave.isFatal = 1;
1544 simError();
1545 }
1546
1547 info.ru = atof( the_token );
1548
1549 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1550 sprintf( painCave.errMsg,
1551 "Error parseing AtomTypes: line %d\n", lineNum );
1552 painCave.isFatal = 1;
1553 simError();
1554 }
1555
1556 info.rlp = atof( the_token );
1557
1558 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1559 sprintf( painCave.errMsg,
1560 "Error parseing AtomTypes: line %d\n", lineNum );
1561 painCave.isFatal = 1;
1562 simError();
1563 }
1564
1565 info.rup = atof( the_token );
1566 }
1567 else info.v0 = info.w0 = info.v0p = info.rl = info.ru = info.rlp = info.rup = 0.0;
1568
1569 return 1;
1570 }
1571 else return 0;
1572 }
1573
1574 int DUFF_NS::parseBond( char *lineBuffer, int lineNum, bondStruct &info ){
1575
1576 char* the_token;
1577 char bondType[30];
1578
1579 the_token = strtok( lineBuffer, " \n\t,;" );
1580 if( the_token != NULL ){
1581
1582 strcpy( info.nameA, the_token );
1583
1584 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1585 sprintf( painCave.errMsg,
1586 "Error parseing BondTypes: line %d\n", lineNum );
1587 painCave.isFatal = 1;
1588 simError();
1589 }
1590
1591 strcpy( info.nameB, the_token );
1592
1593 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1594 sprintf( painCave.errMsg,
1595 "Error parseing BondTypes: line %d\n", lineNum );
1596 painCave.isFatal = 1;
1597 simError();
1598 }
1599
1600 strcpy( bondType, the_token );
1601
1602 if( !strcmp( bondType, "fixed" ) ){
1603 info.type = FIXED_BOND;
1604
1605 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1606 sprintf( painCave.errMsg,
1607 "Error parseing BondTypes: line %d\n", lineNum );
1608 painCave.isFatal = 1;
1609 simError();
1610 }
1611
1612 info.d0 = atof( the_token );
1613 }
1614 else if( !strcmp( bondType, "harmonic" ) ){
1615 info.type = HARMONIC_BOND;
1616
1617 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1618 sprintf( painCave.errMsg,
1619 "Error parseing BondTypes: line %d\n", lineNum );
1620 painCave.isFatal = 1;
1621 simError();
1622 }
1623
1624 info.d0 = atof( the_token );
1625
1626 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1627 sprintf( painCave.errMsg,
1628 "Error parseing BondTypes: line %d\n", lineNum );
1629 painCave.isFatal = 1;
1630 simError();
1631 }
1632
1633 info.k0 = atof( the_token );
1634 }
1635
1636 else{
1637 sprintf( painCave.errMsg,
1638 "Unknown DUFF bond type \"%s\" at line %d\n",
1639 info.type,
1640 lineNum );
1641 painCave.isFatal = 1;
1642 simError();
1643 }
1644
1645 return 1;
1646 }
1647 else return 0;
1648 }
1649
1650
1651 int DUFF_NS::parseBend( char *lineBuffer, int lineNum, bendStruct &info ){
1652
1653 char* the_token;
1654
1655 the_token = strtok( lineBuffer, " \n\t,;" );
1656 if( the_token != NULL ){
1657
1658 strcpy( info.nameA, the_token );
1659
1660 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1661 sprintf( painCave.errMsg,
1662 "Error parseing BendTypes: line %d\n", lineNum );
1663 painCave.isFatal = 1;
1664 simError();
1665 }
1666
1667 strcpy( info.nameB, the_token );
1668
1669 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1670 sprintf( painCave.errMsg,
1671 "Error parseing BendTypes: line %d\n", lineNum );
1672 painCave.isFatal = 1;
1673 simError();
1674 }
1675
1676 strcpy( info.nameC, the_token );
1677
1678 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1679 sprintf( painCave.errMsg,
1680 "Error parseing BendTypes: line %d\n", lineNum );
1681 painCave.isFatal = 1;
1682 simError();
1683 }
1684
1685 strcpy( info.type, the_token );
1686
1687 if( !strcmp( info.type, "quadratic" ) ){
1688 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1689 sprintf( painCave.errMsg,
1690 "Error parseing BendTypes: line %d\n", lineNum );
1691 painCave.isFatal = 1;
1692 simError();
1693 }
1694
1695 info.k1 = atof( the_token );
1696
1697 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1698 sprintf( painCave.errMsg,
1699 "Error parseing BendTypes: line %d\n", lineNum );
1700 painCave.isFatal = 1;
1701 simError();
1702 }
1703
1704 info.k2 = atof( the_token );
1705
1706 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1707 sprintf( painCave.errMsg,
1708 "Error parseing BendTypes: line %d\n", lineNum );
1709 painCave.isFatal = 1;
1710 simError();
1711 }
1712
1713 info.k3 = atof( the_token );
1714
1715 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1716 sprintf( painCave.errMsg,
1717 "Error parseing BendTypes: line %d\n", lineNum );
1718 painCave.isFatal = 1;
1719 simError();
1720 }
1721
1722 info.t0 = atof( the_token );
1723 }
1724
1725 else{
1726 sprintf( painCave.errMsg,
1727 "Unknown DUFF bend type \"%s\" at line %d\n",
1728 info.type,
1729 lineNum );
1730 painCave.isFatal = 1;
1731 simError();
1732 }
1733
1734 return 1;
1735 }
1736 else return 0;
1737 }
1738
1739 int DUFF_NS::parseTorsion( char *lineBuffer, int lineNum, torsionStruct &info ){
1740
1741 char* the_token;
1742
1743 the_token = strtok( lineBuffer, " \n\t,;" );
1744 if( the_token != NULL ){
1745
1746 strcpy( info.nameA, the_token );
1747
1748 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1749 sprintf( painCave.errMsg,
1750 "Error parseing TorsionTypes: line %d\n", lineNum );
1751 painCave.isFatal = 1;
1752 simError();
1753 }
1754
1755 strcpy( info.nameB, the_token );
1756
1757 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1758 sprintf( painCave.errMsg,
1759 "Error parseing TorsionTypes: line %d\n", lineNum );
1760 painCave.isFatal = 1;
1761 simError();
1762 }
1763
1764 strcpy( info.nameC, the_token );
1765
1766 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1767 sprintf( painCave.errMsg,
1768 "Error parseing TorsionTypes: line %d\n", lineNum );
1769 painCave.isFatal = 1;
1770 simError();
1771 }
1772
1773 strcpy( info.nameD, the_token );
1774
1775 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1776 sprintf( painCave.errMsg,
1777 "Error parseing TorsionTypes: line %d\n", lineNum );
1778 painCave.isFatal = 1;
1779 simError();
1780 }
1781
1782 strcpy( info.type, the_token );
1783
1784 if( !strcmp( info.type, "cubic" ) ){
1785 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1786 sprintf( painCave.errMsg,
1787 "Error parseing TorsionTypes: line %d\n", lineNum );
1788 painCave.isFatal = 1;
1789 simError();
1790 }
1791
1792 info.k1 = atof( the_token );
1793
1794 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1795 sprintf( painCave.errMsg,
1796 "Error parseing TorsionTypes: line %d\n", lineNum );
1797 painCave.isFatal = 1;
1798 simError();
1799 }
1800
1801 info.k2 = atof( the_token );
1802
1803 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1804 sprintf( painCave.errMsg,
1805 "Error parseing TorsionTypes: line %d\n", lineNum );
1806 painCave.isFatal = 1;
1807 simError();
1808 }
1809
1810 info.k3 = atof( the_token );
1811
1812 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1813 sprintf( painCave.errMsg,
1814 "Error parseing TorsionTypes: line %d\n", lineNum );
1815 painCave.isFatal = 1;
1816 simError();
1817 }
1818
1819 info.k4 = atof( the_token );
1820
1821 }
1822
1823 else{
1824 sprintf( painCave.errMsg,
1825 "Unknown DUFF torsion type \"%s\" at line %d\n",
1826 info.type,
1827 lineNum );
1828 painCave.isFatal = 1;
1829 simError();
1830 }
1831
1832 return 1;
1833 }
1834
1835 else return 0;
1836 }