ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DUFF.cpp
Revision: 941
Committed: Tue Jan 13 23:01:43 2004 UTC (20 years, 5 months ago) by gezelter
File size: 43526 byte(s)
Log Message: 
Changes for adding direct charge-charge interactions (with switching function)

File Contents

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