ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/UseTheForce/DUFF.cpp
Revision: 1628
Committed: Thu Oct 21 20:15:31 2004 UTC (19 years, 8 months ago) by gezelter
File size: 43823 byte(s)
Log Message: 
Breaky Breaky.   Fixey Fixey.

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 "UseTheForce/ForceFields.hpp"
9 #include "primitives/SRI.hpp"
10 #include "utils/simError.h"
11 #include "UseTheForce/DarkSide/sticky_interface.h"
12 #include "UseTheForce/DarkSide/atype_interface.h"
13
14 //#include "UseTheForce/fortranWrappers.hpp"
15
16
17 #ifdef IS_MPI
18 #include "UseTheForce/mpiForceField.h"
19 #endif // is_mpi
20
21
22 // define some bond Types
23
24 #define FIXED_BOND 0
25 #define HARMONIC_BOND 1
26
27
28 namespace DUFF_NS { // restrict the access of the folowing to this file only.
29
30
31 // Declare the structures that will be passed by MPI
32
33 typedef struct{
34 char name[15];
35 double mass;
36 double epslon;
37 double sigma;
38 double charge;
39 double dipole;
40 double w0;
41 double v0;
42 double v0p;
43 double rl;
44 double ru;
45 double rlp;
46 double rup;
47 int isSSD;
48 int isCharge;
49 int isDipole;
50 int ident;
51 int last; // 0 -> default
52 // 1 -> tells nodes to stop listening
53 } atomStruct;
54
55
56 typedef struct{
57 char nameA[15];
58 char nameB[15];
59 double d0;
60 double k0;
61 int last; // 0 -> default
62 // 1 -> tells nodes to stop listening
63 int type;
64 } bondStruct;
65
66
67 typedef struct{
68 char nameA[15];
69 char nameB[15];
70 char nameC[15];
71 char type[30];
72 double k1, k2, k3, t0;
73 int last; // 0 -> default
74 // 1 -> tells nodes to stop listening
75 } bendStruct;
76
77
78 typedef struct{
79 char nameA[15];
80 char nameB[15];
81 char nameC[15];
82 char nameD[15];
83 char type[30];
84 double k1, k2, k3, k4;
85 int last; // 0 -> default
86 // 1 -> tells nodes to stop listening
87 } torsionStruct;
88
89
90 int parseAtom( char *lineBuffer, int lineNum, atomStruct &info );
91 int parseBond( char *lineBuffer, int lineNum, bondStruct &info );
92 int parseBend( char *lineBuffer, int lineNum, bendStruct &info );
93 int parseTorsion( char *lineBuffer, int lineNum, torsionStruct &info );
94
95
96 #ifdef IS_MPI
97
98 MPI_Datatype mpiAtomStructType;
99 MPI_Datatype mpiBondStructType;
100 MPI_Datatype mpiBendStructType;
101 MPI_Datatype mpiTorsionStructType;
102
103 #endif
104
105 class LinkedAtomType {
106 public:
107 LinkedAtomType(){
108 next = NULL;
109 name[0] = '\0';
110 }
111 ~LinkedAtomType(){ if( next != NULL ) delete next; }
112
113 LinkedAtomType* find(char* key){
114 if( !strcmp(name, key) ) return this;
115 if( next != NULL ) return next->find(key);
116 return NULL;
117 }
118
119 void printMe( void ){
120
121 std::cerr << "LinkedAtype " << name << ": ident = " << ident << "\n";
122 // if( next != NULL ) next->printMe();
123
124 }
125
126 void add( atomStruct &info ){
127
128 // check for duplicates
129
130 if( !strcmp( info.name, name ) ){
131 sprintf( painCave.errMsg,
132 "Duplicate DUFF atom type \"%s\" found in "
133 "the DUFF param file./n",
134 name );
135 painCave.isFatal = 1;
136 simError();
137 }
138
139 if( next != NULL ) next->add(info);
140 else{
141 next = new LinkedAtomType();
142 strcpy(next->name, info.name);
143 next->isDipole = info.isDipole;
144 next->isSSD = info.isSSD;
145 next->mass = info.mass;
146 next->epslon = info.epslon;
147 next->sigma = info.sigma;
148 next->dipole = info.dipole;
149 next->w0 = info.w0;
150 next->v0 = info.v0;
151 next->v0p = info.v0p;
152 next->rl = info.rl;
153 next->ru = info.ru;
154 next->rlp = info.rlp;
155 next->rup = info.rup;
156 next->ident = info.ident;
157 }
158 }
159
160 #ifdef IS_MPI
161
162 void duplicate( atomStruct &info ){
163 strcpy(info.name, name);
164 info.isDipole = isDipole;
165 info.isSSD = isSSD;
166 info.mass = mass;
167 info.epslon = epslon;
168 info.sigma = sigma;
169 info.dipole = dipole;
170 info.w0 = w0;
171 info.v0 = v0;
172 info.v0p = v0p;
173 info.rl = rl;
174 info.ru = ru;
175 info.rlp = rlp;
176 info.rup = rup;
177 info.ident = ident;
178 info.last = 0;
179 }
180
181
182 #endif
183
184 char name[15];
185 int isDipole;
186 int isSSD;
187 double mass;
188 double epslon;
189 double sigma;
190 double dipole;
191 double w0;
192 double v0;
193 double v0p;
194 double rl;
195 double ru;
196 double rlp;
197 double rup;
198 int ident;
199 LinkedAtomType* next;
200 };
201
202 class LinkedBondType {
203 public:
204 LinkedBondType(){
205 next = NULL;
206 nameA[0] = '\0';
207 nameB[0] = '\0';
208 }
209 ~LinkedBondType(){ if( next != NULL ) delete next; }
210
211 LinkedBondType* find(char* key1, char* key2){
212 if( !strcmp(nameA, key1 ) && !strcmp( nameB, key2 ) ) return this;
213 if( !strcmp(nameA, key2 ) && !strcmp( nameB, key1 ) ) return this;
214 if( next != NULL ) return next->find(key1, key2);
215 return NULL;
216 }
217
218
219 void add( bondStruct &info ){
220
221 // check for duplicates
222 int dup = 0;
223
224 if( !strcmp(nameA, info.nameA ) && !strcmp( nameB, info.nameB ) ) dup = 1;
225 if( !strcmp(nameA, info.nameB ) && !strcmp( nameB, info.nameA ) ) dup = 1;
226
227 if(dup){
228 sprintf( painCave.errMsg,
229 "Duplicate DUFF bond type \"%s - %s\" found in "
230 "the DUFF param file./n",
231 nameA, nameB );
232 painCave.isFatal = 1;
233 simError();
234 }
235
236
237 if( next != NULL ) next->add(info);
238 else{
239 next = new LinkedBondType();
240 strcpy(next->nameA, info.nameA);
241 strcpy(next->nameB, info.nameB);
242 next->type = info.type;
243 next->d0 = info.d0;
244 next->k0 = info.k0;
245 }
246 }
247
248 #ifdef IS_MPI
249 void duplicate( bondStruct &info ){
250 strcpy(info.nameA, nameA);
251 strcpy(info.nameB, nameB);
252 info.type = type;
253 info.d0 = d0;
254 info.k0 = k0;
255 info.last = 0;
256 }
257
258
259 #endif
260
261 char nameA[15];
262 char nameB[15];
263 int type;
264 double d0;
265 double k0;
266
267 LinkedBondType* next;
268 };
269
270 class LinkedBendType {
271 public:
272 LinkedBendType(){
273 next = NULL;
274 nameA[0] = '\0';
275 nameB[0] = '\0';
276 nameC[0] = '\0';
277 type[0] = '\0';
278 }
279 ~LinkedBendType(){ if( next != NULL ) delete next; }
280
281 LinkedBendType* find( char* key1, char* key2, char* key3 ){
282 if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 )
283 && !strcmp( nameC, key3 ) ) return this;
284 if( !strcmp( nameA, key3 ) && !strcmp( nameB, key2 )
285 && !strcmp( nameC, key1 ) ) return this;
286 if( next != NULL ) return next->find(key1, key2, key3);
287 return NULL;
288 }
289
290 void add( bendStruct &info ){
291
292 // check for duplicates
293 int dup = 0;
294
295 if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB )
296 && !strcmp( nameC, info.nameC ) ) dup = 1;
297 if( !strcmp( nameA, info.nameC ) && !strcmp( nameB, info.nameB )
298 && !strcmp( nameC, info.nameA ) ) dup = 1;
299
300 if(dup){
301 sprintf( painCave.errMsg,
302 "Duplicate DUFF bend type \"%s - %s - %s\" found in "
303 "the DUFF param file./n",
304 nameA, nameB, nameC );
305 painCave.isFatal = 1;
306 simError();
307 }
308
309 if( next != NULL ) next->add(info);
310 else{
311 next = new LinkedBendType();
312 strcpy(next->nameA, info.nameA);
313 strcpy(next->nameB, info.nameB);
314 strcpy(next->nameC, info.nameC);
315 strcpy(next->type, info.type);
316 next->k1 = info.k1;
317 next->k2 = info.k2;
318 next->k3 = info.k3;
319 next->t0 = info.t0;
320 }
321 }
322
323 #ifdef IS_MPI
324
325 void duplicate( bendStruct &info ){
326 strcpy(info.nameA, nameA);
327 strcpy(info.nameB, nameB);
328 strcpy(info.nameC, nameC);
329 strcpy(info.type, type);
330 info.k1 = k1;
331 info.k2 = k2;
332 info.k3 = k3;
333 info.t0 = t0;
334 info.last = 0;
335 }
336
337 #endif // is_mpi
338
339 char nameA[15];
340 char nameB[15];
341 char nameC[15];
342 char type[30];
343 double k1, k2, k3, t0;
344
345 LinkedBendType* next;
346 };
347
348 class LinkedTorsionType {
349 public:
350 LinkedTorsionType(){
351 next = NULL;
352 nameA[0] = '\0';
353 nameB[0] = '\0';
354 nameC[0] = '\0';
355 type[0] = '\0';
356 }
357 ~LinkedTorsionType(){ if( next != NULL ) delete next; }
358
359 LinkedTorsionType* find( char* key1, char* key2, char* key3, char* key4 ){
360
361
362
363
364 if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 ) &&
365 !strcmp( nameC, key3 ) && !strcmp( nameD, key4 ) ) return this;
366
367 if( !strcmp( nameA, key4 ) && !strcmp( nameB, key3 ) &&
368 !strcmp( nameC, key2 ) && !strcmp( nameD, key1 ) ) return this;
369
370 if( next != NULL ) return next->find(key1, key2, key3, key4);
371 return NULL;
372 }
373
374 void add( torsionStruct &info ){
375
376 // check for duplicates
377 int dup = 0;
378
379 if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB ) &&
380 !strcmp( nameC, info.nameC ) && !strcmp( nameD, info.nameD ) ) dup = 1;
381
382 if( !strcmp( nameA, info.nameD ) && !strcmp( nameB, info.nameC ) &&
383 !strcmp( nameC, info.nameB ) && !strcmp( nameD, info.nameA ) ) dup = 1;
384
385 if(dup){
386 sprintf( painCave.errMsg,
387 "Duplicate DUFF torsion type \"%s - %s - %s - %s\" found in "
388 "the DUFF param file./n", nameA, nameB, nameC, nameD );
389 painCave.isFatal = 1;
390 simError();
391 }
392
393 if( next != NULL ) next->add(info);
394 else{
395 next = new LinkedTorsionType();
396 strcpy(next->nameA, info.nameA);
397 strcpy(next->nameB, info.nameB);
398 strcpy(next->nameC, info.nameC);
399 strcpy(next->nameD, info.nameD);
400 strcpy(next->type, info.type);
401 next->k1 = info.k1;
402 next->k2 = info.k2;
403 next->k3 = info.k3;
404 next->k4 = info.k4;
405
406 }
407 }
408
409 #ifdef IS_MPI
410
411 void duplicate( torsionStruct &info ){
412 strcpy(info.nameA, nameA);
413 strcpy(info.nameB, nameB);
414 strcpy(info.nameC, nameC);
415 strcpy(info.nameD, nameD);
416 strcpy(info.type, type);
417 info.k1 = k1;
418 info.k2 = k2;
419 info.k3 = k3;
420 info.k4 = k4;
421 info.last = 0;
422 }
423
424 #endif
425
426 char nameA[15];
427 char nameB[15];
428 char nameC[15];
429 char nameD[15];
430 char type[30];
431 double k1, k2, k3, k4;
432
433 LinkedTorsionType* next;
434 };
435
436
437 LinkedAtomType* headAtomType;
438 LinkedAtomType* currentAtomType;
439 LinkedBondType* headBondType;
440 LinkedBondType* currentBondType;
441 LinkedBendType* headBendType;
442 LinkedBendType* currentBendType;
443 LinkedTorsionType* headTorsionType;
444 LinkedTorsionType* currentTorsionType;
445
446 } // namespace
447
448 using namespace DUFF_NS;
449
450
451 //****************************************************************
452 // begins the actual forcefield stuff.
453 //****************************************************************
454
455
456 DUFF::DUFF(){
457
458 char fileName[200];
459 char* ffPath_env = "FORCE_PARAM_PATH";
460 char* ffPath;
461 char temp[200];
462
463 headAtomType = NULL;
464 currentAtomType = NULL;
465 headBondType = NULL;
466 currentBondType = NULL;
467 headBendType = NULL;
468 currentBendType = NULL;
469 headTorsionType = NULL;
470 currentTorsionType = NULL;
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,12,5}; // 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:\n"
599 "\t%s\n"
600 "\tHave you tried setting the FORCE_PARAM_PATH environment "
601 "variable?\n",
602 fileName );
603 painCave.severity = OOPSE_ERROR;
604 painCave.isFatal = 1;
605 simError();
606 }
607 }
608
609 #ifdef IS_MPI
610 }
611
612 sprintf( checkPointMsg, "DUFF file opened sucessfully." );
613 MPIcheckPoint();
614
615 #endif // is_mpi
616 }
617
618
619 DUFF::~DUFF(){
620
621 if( headAtomType != NULL ) delete headAtomType;
622 if( headBondType != NULL ) delete headBondType;
623 if( headBendType != NULL ) delete headBendType;
624 if( headTorsionType != NULL ) delete headTorsionType;
625
626 #ifdef IS_MPI
627 if( worldRank == 0 ){
628 #endif // is_mpi
629
630 fclose( frcFile );
631
632 #ifdef IS_MPI
633 }
634 #endif // is_mpi
635 }
636
637 void DUFF::cleanMe( void ){
638
639 #ifdef IS_MPI
640
641 // keep the linked lists in the mpi version
642
643 #else // is_mpi
644
645 // delete the linked lists in the single processor version
646
647 if( headAtomType != NULL ) delete headAtomType;
648 if( headBondType != NULL ) delete headBondType;
649 if( headBendType != NULL ) delete headBendType;
650 if( headTorsionType != NULL ) delete headTorsionType;
651
652 #endif // is_mpi
653 }
654
655
656 void DUFF::initForceField(){
657
658 initFortran( entry_plug->useReactionField );
659 }
660
661 double DUFF::getAtomTypeMass (char* atomType) {
662
663 currentAtomType = headAtomType->find( atomType );
664 if( currentAtomType == NULL ){
665 sprintf( painCave.errMsg,
666 "AtomType error, %s not found in force file.\n",
667 atomType );
668 painCave.isFatal = 1;
669 simError();
670 }
671
672 return currentAtomType->mass;
673 }
674
675 void DUFF::readParams( void ){
676
677 int identNum;
678
679 atomStruct atomInfo;
680 bondStruct bondInfo;
681 bendStruct bendInfo;
682 torsionStruct torsionInfo;
683
684 bigSigma = 0.0;
685
686 atomInfo.last = 1;
687 bondInfo.last = 1;
688 bendInfo.last = 1;
689 torsionInfo.last = 1;
690
691 // read in the atom info
692
693 #ifdef IS_MPI
694 if( worldRank == 0 ){
695 #endif
696
697 // read in the atom types.
698
699 headAtomType = new LinkedAtomType;
700
701 fastForward( "AtomTypes", "initializeAtoms" );
702
703 // we are now at the AtomTypes section.
704
705 eof_test = fgets( readLine, sizeof(readLine), frcFile );
706 lineNum++;
707
708
709 // read a line, and start parseing out the atom types
710
711 if( eof_test == NULL ){
712 sprintf( painCave.errMsg,
713 "Error in reading Atoms from force file at line %d.\n",
714 lineNum );
715 painCave.isFatal = 1;
716 simError();
717 }
718
719 identNum = 1;
720 // stop reading at end of file, or at next section
721 while( readLine[0] != '#' && eof_test != NULL ){
722
723 // toss comment lines
724 if( readLine[0] != '!' ){
725
726 // the parser returns 0 if the line was blank
727 if( parseAtom( readLine, lineNum, atomInfo ) ){
728 atomInfo.ident = identNum;
729 headAtomType->add( atomInfo );;
730 identNum++;
731 }
732 }
733 eof_test = fgets( readLine, sizeof(readLine), frcFile );
734 lineNum++;
735 }
736
737 #ifdef IS_MPI
738
739 // send out the linked list to all the other processes
740
741 sprintf( checkPointMsg,
742 "DUFF atom structures read successfully." );
743 MPIcheckPoint();
744
745 currentAtomType = headAtomType->next; //skip the first element who is a place holder.
746 while( currentAtomType != NULL ){
747 currentAtomType->duplicate( atomInfo );
748
749 sendFrcStruct( &atomInfo, mpiAtomStructType );
750
751 sprintf( checkPointMsg,
752 "successfully sent DUFF force type: \"%s\"\n",
753 atomInfo.name );
754 MPIcheckPoint();
755
756 currentAtomType = currentAtomType->next;
757 }
758 atomInfo.last = 1;
759 sendFrcStruct( &atomInfo, mpiAtomStructType );
760
761 }
762
763 else{
764
765 // listen for node 0 to send out the force params
766
767 MPIcheckPoint();
768
769 headAtomType = new LinkedAtomType;
770 receiveFrcStruct( &atomInfo, mpiAtomStructType );
771
772 while( !atomInfo.last ){
773
774 headAtomType->add( atomInfo );
775
776 MPIcheckPoint();
777
778 receiveFrcStruct( &atomInfo, mpiAtomStructType );
779 }
780 }
781
782 #endif // is_mpi
783
784
785
786 // call new A_types in fortran
787
788 int isError;
789
790 // dummy variables
791
792 int isGB = 0;
793 int isLJ = 1;
794 int isEAM =0;
795 int isCharge = 0;
796 double charge=0.0;
797
798 currentAtomType = headAtomType->next;;
799 while( currentAtomType != NULL ){
800
801 if(currentAtomType->isDipole) entry_plug->useDipoles = 1;
802 if(currentAtomType->isSSD) {
803 entry_plug->useSticky = 1;
804 makeStickyType( &(currentAtomType->w0), &(currentAtomType->v0),
805 &(currentAtomType->v0p),
806 &(currentAtomType->rl), &(currentAtomType->ru),
807 &(currentAtomType->rlp), &(currentAtomType->rup));
808 }
809
810 if( currentAtomType->name[0] != '\0' ){
811 isError = 0;
812 makeAtype( &(currentAtomType->ident),
813 &isLJ,
814 &(currentAtomType->isSSD),
815 &(currentAtomType->isDipole),
816 &isGB,
817 &isEAM,
818 &isCharge,
819 &(currentAtomType->epslon),
820 &(currentAtomType->sigma),
821 &charge,
822 &(currentAtomType->dipole),
823 &isError );
824 if( isError ){
825 sprintf( painCave.errMsg,
826 "Error initializing the \"%s\" atom type in fortran\n",
827 currentAtomType->name );
828 painCave.isFatal = 1;
829 simError();
830 }
831 }
832 currentAtomType = currentAtomType->next;
833 }
834
835 #ifdef IS_MPI
836 sprintf( checkPointMsg,
837 "DUFF atom structures successfully sent to fortran\n" );
838 MPIcheckPoint();
839 #endif // is_mpi
840
841
842
843 // read in the bonds
844
845 #ifdef IS_MPI
846 if( worldRank == 0 ){
847 #endif
848
849 // read in the bond types.
850
851 headBondType = new LinkedBondType;
852
853 fastForward( "BondTypes", "initializeBonds" );
854
855 // we are now at the bondTypes section
856
857 eof_test = fgets( readLine, sizeof(readLine), frcFile );
858 lineNum++;
859
860
861 // read a line, and start parseing out the atom types
862
863 if( eof_test == NULL ){
864 sprintf( painCave.errMsg,
865 "Error in reading bonds from force file at line %d.\n",
866 lineNum );
867 painCave.isFatal = 1;
868 simError();
869 }
870
871 // stop reading at end of file, or at next section
872 while( readLine[0] != '#' && eof_test != NULL ){
873
874 // toss comment lines
875 if( readLine[0] != '!' ){
876
877 // the parser returns 0 if the line was blank
878 if( parseBond( readLine, lineNum, bondInfo ) ){
879 headBondType->add( bondInfo );
880 }
881 }
882 eof_test = fgets( readLine, sizeof(readLine), frcFile );
883 lineNum++;
884 }
885
886 #ifdef IS_MPI
887
888 // send out the linked list to all the other processes
889
890 sprintf( checkPointMsg,
891 "DUFF bond structures read successfully." );
892 MPIcheckPoint();
893
894 currentBondType = headBondType->next;
895 while( currentBondType != NULL ){
896 currentBondType->duplicate( bondInfo );
897 sendFrcStruct( &bondInfo, mpiBondStructType );
898 currentBondType = currentBondType->next;
899 }
900 bondInfo.last = 1;
901 sendFrcStruct( &bondInfo, mpiBondStructType );
902
903 }
904
905 else{
906
907 // listen for node 0 to send out the force params
908
909 MPIcheckPoint();
910
911 headBondType = new LinkedBondType;
912 receiveFrcStruct( &bondInfo, mpiBondStructType );
913 while( !bondInfo.last ){
914
915 headBondType->add( bondInfo );
916 receiveFrcStruct( &bondInfo, mpiBondStructType );
917 }
918 }
919
920 sprintf( checkPointMsg,
921 "DUFF bond structures broadcast successfully." );
922 MPIcheckPoint();
923
924 #endif // is_mpi
925
926
927 // read in the bends
928
929 #ifdef IS_MPI
930 if( worldRank == 0 ){
931 #endif
932
933 // read in the bend types.
934
935 headBendType = new LinkedBendType;
936
937 fastForward( "BendTypes", "initializeBends" );
938
939 // we are now at the bendTypes section
940
941 eof_test = fgets( readLine, sizeof(readLine), frcFile );
942 lineNum++;
943
944 // read a line, and start parseing out the bend types
945
946 if( eof_test == NULL ){
947 sprintf( painCave.errMsg,
948 "Error in reading bends from force file at line %d.\n",
949 lineNum );
950 painCave.isFatal = 1;
951 simError();
952 }
953
954 // stop reading at end of file, or at next section
955 while( readLine[0] != '#' && eof_test != NULL ){
956
957 // toss comment lines
958 if( readLine[0] != '!' ){
959
960 // the parser returns 0 if the line was blank
961 if( parseBend( readLine, lineNum, bendInfo ) ){
962 headBendType->add( bendInfo );
963 }
964 }
965 eof_test = fgets( readLine, sizeof(readLine), frcFile );
966 lineNum++;
967 }
968
969 #ifdef IS_MPI
970
971 // send out the linked list to all the other processes
972
973 sprintf( checkPointMsg,
974 "DUFF bend structures read successfully." );
975 MPIcheckPoint();
976
977 currentBendType = headBendType->next;
978 while( currentBendType != NULL ){
979 currentBendType->duplicate( bendInfo );
980 sendFrcStruct( &bendInfo, mpiBendStructType );
981 currentBendType = currentBendType->next;
982 }
983 bendInfo.last = 1;
984 sendFrcStruct( &bendInfo, mpiBendStructType );
985
986 }
987
988 else{
989
990 // listen for node 0 to send out the force params
991
992 MPIcheckPoint();
993
994 headBendType = new LinkedBendType;
995 receiveFrcStruct( &bendInfo, mpiBendStructType );
996 while( !bendInfo.last ){
997
998 headBendType->add( bendInfo );
999 receiveFrcStruct( &bendInfo, mpiBendStructType );
1000 }
1001 }
1002
1003 sprintf( checkPointMsg,
1004 "DUFF bend structures broadcast successfully." );
1005 MPIcheckPoint();
1006
1007 #endif // is_mpi
1008
1009
1010 // read in the torsions
1011
1012 #ifdef IS_MPI
1013 if( worldRank == 0 ){
1014 #endif
1015
1016 // read in the torsion types.
1017
1018 headTorsionType = new LinkedTorsionType;
1019
1020 fastForward( "TorsionTypes", "initializeTorsions" );
1021
1022 // we are now at the torsionTypes section
1023
1024 eof_test = fgets( readLine, sizeof(readLine), frcFile );
1025 lineNum++;
1026
1027
1028 // read a line, and start parseing out the atom types
1029
1030 if( eof_test == NULL ){
1031 sprintf( painCave.errMsg,
1032 "Error in reading torsions from force file at line %d.\n",
1033 lineNum );
1034 painCave.isFatal = 1;
1035 simError();
1036 }
1037
1038 // stop reading at end of file, or at next section
1039 while( readLine[0] != '#' && eof_test != NULL ){
1040
1041 // toss comment lines
1042 if( readLine[0] != '!' ){
1043
1044 // the parser returns 0 if the line was blank
1045 if( parseTorsion( readLine, lineNum, torsionInfo ) ){
1046 headTorsionType->add( torsionInfo );
1047
1048 }
1049 }
1050 eof_test = fgets( readLine, sizeof(readLine), frcFile );
1051 lineNum++;
1052 }
1053
1054 #ifdef IS_MPI
1055
1056 // send out the linked list to all the other processes
1057
1058 sprintf( checkPointMsg,
1059 "DUFF torsion structures read successfully." );
1060 MPIcheckPoint();
1061
1062 currentTorsionType = headTorsionType->next;
1063 while( currentTorsionType != NULL ){
1064 currentTorsionType->duplicate( torsionInfo );
1065 sendFrcStruct( &torsionInfo, mpiTorsionStructType );
1066 currentTorsionType = currentTorsionType->next;
1067 }
1068 torsionInfo.last = 1;
1069 sendFrcStruct( &torsionInfo, mpiTorsionStructType );
1070
1071 }
1072
1073 else{
1074
1075 // listen for node 0 to send out the force params
1076
1077 MPIcheckPoint();
1078
1079 headTorsionType = new LinkedTorsionType;
1080 receiveFrcStruct( &torsionInfo, mpiTorsionStructType );
1081 while( !torsionInfo.last ){
1082
1083 headTorsionType->add( torsionInfo );
1084 receiveFrcStruct( &torsionInfo, mpiTorsionStructType );
1085 }
1086 }
1087
1088 sprintf( checkPointMsg,
1089 "DUFF torsion structures broadcast successfully." );
1090 MPIcheckPoint();
1091
1092 #endif // is_mpi
1093
1094 entry_plug->useLJ = 1;
1095 }
1096
1097
1098
1099 void DUFF::initializeAtoms( int nAtoms, Atom** the_atoms ){
1100
1101
1102 //////////////////////////////////////////////////
1103 // a quick water fix
1104
1105 double waterI[3][3];
1106 waterI[0][0] = 1.76958347772500;
1107 waterI[0][1] = 0.0;
1108 waterI[0][2] = 0.0;
1109
1110 waterI[1][0] = 0.0;
1111 waterI[1][1] = 0.614537057924513;
1112 waterI[1][2] = 0.0;
1113
1114 waterI[2][0] = 0.0;
1115 waterI[2][1] = 0.0;
1116 waterI[2][2] = 1.15504641980049;
1117
1118
1119 double headI[3][3];
1120 headI[0][0] = 1125;
1121 headI[0][1] = 0.0;
1122 headI[0][2] = 0.0;
1123
1124 headI[1][0] = 0.0;
1125 headI[1][1] = 1125;
1126 headI[1][2] = 0.0;
1127
1128 headI[2][0] = 0.0;
1129 headI[2][1] = 0.0;
1130 headI[2][2] = 250;
1131
1132 //////////////////////////////////////////////////
1133
1134
1135 // initialize the atoms
1136
1137 DirectionalAtom* dAtom;
1138 double ji[3];
1139
1140 for(int i=0; i<nAtoms; i++ ){
1141
1142 currentAtomType = headAtomType->find( the_atoms[i]->getType() );
1143 if( currentAtomType == NULL ){
1144 sprintf( painCave.errMsg,
1145 "AtomType error, %s not found in force file.\n",
1146 the_atoms[i]->getType() );
1147 painCave.isFatal = 1;
1148 simError();
1149 }
1150
1151 the_atoms[i]->setMass( currentAtomType->mass );
1152 the_atoms[i]->setIdent( currentAtomType->ident );
1153
1154 if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
1155
1156 if( currentAtomType->isDipole ){
1157 if( the_atoms[i]->isDirectional() ){
1158
1159 dAtom = (DirectionalAtom *) the_atoms[i];
1160 dAtom->setHasDipole( 1 );
1161
1162 ji[0] = 0.0;
1163 ji[1] = 0.0;
1164 ji[2] = 0.0;
1165
1166 dAtom->setJ( ji );
1167
1168 if(!strcmp("SSD",the_atoms[i]->getType())){
1169 dAtom->setI( waterI );
1170 }
1171 else if(!strcmp("HEAD",the_atoms[i]->getType())){
1172 dAtom->setI( headI );
1173 }
1174 else{
1175 sprintf(painCave.errMsg,
1176 "AtmType error, %s does not have a moment of inertia set.\n",
1177 the_atoms[i]->getType() );
1178 painCave.isFatal = 1;
1179 simError();
1180 }
1181 entry_plug->n_dipoles++;
1182 }
1183 else{
1184
1185 sprintf( painCave.errMsg,
1186 "DUFF error: Atom \"%s\" is a dipole, yet no standard"
1187 " orientation was specifed in the meta-data file.\n",
1188 currentAtomType->name );
1189 painCave.isFatal = 1;
1190 simError();
1191 }
1192 }
1193 else{
1194 if( the_atoms[i]->isDirectional() ){
1195 sprintf( painCave.errMsg,
1196 "DUFF error: Atom \"%s\" was given a standard "
1197 "orientation in the meta-data file, yet it is not a dipole.\n",
1198 currentAtomType->name);
1199 painCave.isFatal = 1;
1200 simError();
1201 }
1202 }
1203 }
1204 }
1205
1206 void DUFF::initializeBonds( int nBonds, Bond** bondArray,
1207 bond_pair* the_bonds ){
1208 int i,a,b;
1209 char* atomA;
1210 char* atomB;
1211
1212 Atom** the_atoms;
1213 the_atoms = entry_plug->atoms;
1214
1215
1216 // initialize the Bonds
1217
1218 for( i=0; i<nBonds; i++ ){
1219
1220 a = the_bonds[i].a;
1221 b = the_bonds[i].b;
1222
1223 atomA = the_atoms[a]->getType();
1224 atomB = the_atoms[b]->getType();
1225 currentBondType = headBondType->find( atomA, atomB );
1226 if( currentBondType == NULL ){
1227 sprintf( painCave.errMsg,
1228 "BondType error, %s - %s not found in force file.\n",
1229 atomA, atomB );
1230 painCave.isFatal = 1;
1231 simError();
1232 }
1233
1234 switch( currentBondType->type ){
1235
1236 case FIXED_BOND:
1237
1238 bondArray[i] = new ConstrainedBond( *the_atoms[a],
1239 *the_atoms[b],
1240 currentBondType->d0 );
1241 entry_plug->n_constraints++;
1242 break;
1243
1244 case HARMONIC_BOND:
1245
1246 bondArray[i] = new HarmonicBond( *the_atoms[a],
1247 *the_atoms[b],
1248 currentBondType->d0,
1249 currentBondType->k0 );
1250 break;
1251
1252 default:
1253
1254 break;
1255 // do nothing
1256 }
1257 }
1258 }
1259
1260 void DUFF::initializeBends( int nBends, Bend** bendArray,
1261 bend_set* the_bends ){
1262
1263 QuadraticBend* qBend;
1264 GhostBend* gBend;
1265 Atom** the_atoms;
1266 the_atoms = entry_plug->atoms;
1267
1268 int i, a, b, c;
1269 char* atomA;
1270 char* atomB;
1271 char* atomC;
1272
1273 // initialize the Bends
1274
1275 for( i=0; i<nBends; i++ ){
1276
1277 a = the_bends[i].a;
1278 b = the_bends[i].b;
1279 c = the_bends[i].c;
1280
1281 atomA = the_atoms[a]->getType();
1282 atomB = the_atoms[b]->getType();
1283
1284 if( the_bends[i].isGhost ) atomC = "GHOST";
1285 else atomC = the_atoms[c]->getType();
1286
1287 currentBendType = headBendType->find( atomA, atomB, atomC );
1288 if( currentBendType == NULL ){
1289 sprintf( painCave.errMsg, "BendType error, %s - %s - %s not found"
1290 " in force file.\n",
1291 atomA, atomB, atomC );
1292 painCave.isFatal = 1;
1293 simError();
1294 }
1295
1296 if( !strcmp( currentBendType->type, "quadratic" ) ){
1297
1298 if( the_bends[i].isGhost){
1299
1300 if( the_bends[i].ghost == b ){
1301 // do nothing
1302 }
1303 else if( the_bends[i].ghost == a ){
1304 c = a;
1305 a = b;
1306 b = c;
1307 }
1308 else{
1309 sprintf( painCave.errMsg,
1310 "BendType error, %s - %s - %s,\n"
1311 " --> central atom is not "
1312 "correctly identified with the "
1313 "\"ghostVectorSource = \" tag.\n",
1314 atomA, atomB, atomC );
1315 painCave.isFatal = 1;
1316 simError();
1317 }
1318
1319 gBend = new GhostBend( *the_atoms[a],
1320 *the_atoms[b]);
1321
1322 gBend->setConstants( currentBendType->k1,
1323 currentBendType->k2,
1324 currentBendType->k3,
1325 currentBendType->t0 );
1326 bendArray[i] = gBend;
1327 }
1328 else{
1329 qBend = new QuadraticBend( *the_atoms[a],
1330 *the_atoms[b],
1331 *the_atoms[c] );
1332 qBend->setConstants( currentBendType->k1,
1333 currentBendType->k2,
1334 currentBendType->k3,
1335 currentBendType->t0 );
1336 bendArray[i] = qBend;
1337 }
1338 }
1339 }
1340 }
1341
1342 void DUFF::initializeTorsions( int nTorsions, Torsion** torsionArray,
1343 torsion_set* the_torsions ){
1344
1345 int i, a, b, c, d;
1346 char* atomA;
1347 char* atomB;
1348 char* atomC;
1349 char* atomD;
1350
1351 CubicTorsion* cTors;
1352 Atom** the_atoms;
1353 the_atoms = entry_plug->atoms;
1354
1355 // initialize the Torsions
1356
1357 for( i=0; i<nTorsions; i++ ){
1358
1359 a = the_torsions[i].a;
1360 b = the_torsions[i].b;
1361 c = the_torsions[i].c;
1362 d = the_torsions[i].d;
1363
1364 atomA = the_atoms[a]->getType();
1365 atomB = the_atoms[b]->getType();
1366 atomC = the_atoms[c]->getType();
1367 atomD = the_atoms[d]->getType();
1368 currentTorsionType = headTorsionType->find( atomA, atomB, atomC, atomD );
1369 if( currentTorsionType == NULL ){
1370 sprintf( painCave.errMsg,
1371 "TorsionType error, %s - %s - %s - %s not found"
1372 " in force file.\n",
1373 atomA, atomB, atomC, atomD );
1374 painCave.isFatal = 1;
1375 simError();
1376 }
1377
1378 if( !strcmp( currentTorsionType->type, "cubic" ) ){
1379
1380 cTors = new CubicTorsion( *the_atoms[a], *the_atoms[b],
1381 *the_atoms[c], *the_atoms[d] );
1382 cTors->setConstants( currentTorsionType->k1, currentTorsionType->k2,
1383 currentTorsionType->k3, currentTorsionType->k4 );
1384 torsionArray[i] = cTors;
1385 }
1386 }
1387 }
1388
1389 void DUFF::fastForward( char* stopText, char* searchOwner ){
1390
1391 int foundText = 0;
1392 char* the_token;
1393
1394 rewind( frcFile );
1395 lineNum = 0;
1396
1397 eof_test = fgets( readLine, sizeof(readLine), frcFile );
1398 lineNum++;
1399 if( eof_test == NULL ){
1400 sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
1401 " file is empty.\n",
1402 searchOwner );
1403 painCave.isFatal = 1;
1404 simError();
1405 }
1406
1407
1408 while( !foundText ){
1409 while( eof_test != NULL && readLine[0] != '#' ){
1410 eof_test = fgets( readLine, sizeof(readLine), frcFile );
1411 lineNum++;
1412 }
1413 if( eof_test == NULL ){
1414 sprintf( painCave.errMsg,
1415 "Error fast forwarding force file for %s at "
1416 "line %d: file ended unexpectedly.\n",
1417 searchOwner,
1418 lineNum );
1419 painCave.isFatal = 1;
1420 simError();
1421 }
1422
1423 the_token = strtok( readLine, " ,;\t#\n" );
1424 foundText = !strcmp( stopText, the_token );
1425
1426 if( !foundText ){
1427 eof_test = fgets( readLine, sizeof(readLine), frcFile );
1428 lineNum++;
1429
1430 if( eof_test == NULL ){
1431 sprintf( painCave.errMsg,
1432 "Error fast forwarding force file for %s at "
1433 "line %d: file ended unexpectedly.\n",
1434 searchOwner,
1435 lineNum );
1436 painCave.isFatal = 1;
1437 simError();
1438 }
1439 }
1440 }
1441 }
1442
1443
1444 int DUFF_NS::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
1445
1446 char* the_token;
1447
1448 the_token = strtok( lineBuffer, " \n\t,;" );
1449 if( the_token != NULL ){
1450
1451 strcpy( info.name, 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.isDipole = 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.isSSD = atoi( 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.mass = 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.epslon = 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.sigma = atof( the_token );
1497
1498 if( info.isDipole ){
1499
1500 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1501 sprintf( painCave.errMsg,
1502 "Error parseing AtomTypes: line %d\n", lineNum );
1503 painCave.isFatal = 1;
1504 simError();
1505 }
1506
1507 info.dipole = atof( the_token );
1508 }
1509 else info.dipole = 0.0;
1510
1511 if( info.isSSD ){
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.w0 = atof( the_token );
1521
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.v0 = atof( the_token );
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.v0p = 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.rl = 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.ru = 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.rlp = atof( the_token );
1565
1566 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1567 sprintf( painCave.errMsg,
1568 "Error parseing AtomTypes: line %d\n", lineNum );
1569 painCave.isFatal = 1;
1570 simError();
1571 }
1572
1573 info.rup = atof( the_token );
1574 }
1575 else info.v0 = info.w0 = info.v0p = info.rl = info.ru = info.rlp = info.rup = 0.0;
1576
1577 return 1;
1578 }
1579 else return 0;
1580 }
1581
1582 int DUFF_NS::parseBond( char *lineBuffer, int lineNum, bondStruct &info ){
1583
1584 char* the_token;
1585 char bondType[30];
1586
1587 the_token = strtok( lineBuffer, " \n\t,;" );
1588 if( the_token != NULL ){
1589
1590 strcpy( info.nameA, 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( info.nameB, the_token );
1600
1601 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1602 sprintf( painCave.errMsg,
1603 "Error parseing BondTypes: line %d\n", lineNum );
1604 painCave.isFatal = 1;
1605 simError();
1606 }
1607
1608 strcpy( bondType, the_token );
1609
1610 if( !strcmp( bondType, "fixed" ) ){
1611 info.type = FIXED_BOND;
1612
1613 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1614 sprintf( painCave.errMsg,
1615 "Error parseing BondTypes: line %d\n", lineNum );
1616 painCave.isFatal = 1;
1617 simError();
1618 }
1619
1620 info.d0 = atof( the_token );
1621
1622 info.k0=0.0;
1623 }
1624 else if( !strcmp( bondType, "harmonic" ) ){
1625 info.type = HARMONIC_BOND;
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.d0 = atof( the_token );
1635
1636 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1637 sprintf( painCave.errMsg,
1638 "Error parseing BondTypes: line %d\n", lineNum );
1639 painCave.isFatal = 1;
1640 simError();
1641 }
1642
1643 info.k0 = atof( the_token );
1644 }
1645
1646 else{
1647 sprintf( painCave.errMsg,
1648 "Unknown DUFF bond type \"%s\" at line %d\n",
1649 bondType,
1650 lineNum );
1651 painCave.isFatal = 1;
1652 simError();
1653 }
1654
1655 return 1;
1656 }
1657 else return 0;
1658 }
1659
1660
1661 int DUFF_NS::parseBend( char *lineBuffer, int lineNum, bendStruct &info ){
1662
1663 char* the_token;
1664
1665 the_token = strtok( lineBuffer, " \n\t,;" );
1666 if( the_token != NULL ){
1667
1668 strcpy( info.nameA, 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.nameB, 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.nameC, the_token );
1687
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 strcpy( info.type, the_token );
1696
1697 if( !strcmp( info.type, "quadratic" ) ){
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.k1 = 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.k2 = 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.k3 = atof( the_token );
1724
1725 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1726 sprintf( painCave.errMsg,
1727 "Error parseing BendTypes: line %d\n", lineNum );
1728 painCave.isFatal = 1;
1729 simError();
1730 }
1731
1732 info.t0 = atof( the_token );
1733 }
1734
1735 else{
1736 sprintf( painCave.errMsg,
1737 "Unknown DUFF bend type \"%s\" at line %d\n",
1738 info.type,
1739 lineNum );
1740 painCave.isFatal = 1;
1741 simError();
1742 }
1743
1744 return 1;
1745 }
1746 else return 0;
1747 }
1748
1749 int DUFF_NS::parseTorsion( char *lineBuffer, int lineNum, torsionStruct &info ){
1750
1751 char* the_token;
1752
1753 the_token = strtok( lineBuffer, " \n\t,;" );
1754 if( the_token != NULL ){
1755
1756 strcpy( info.nameA, 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.nameB, 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.nameC, 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.nameD, the_token );
1784
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 strcpy( info.type, the_token );
1793
1794 if( !strcmp( info.type, "cubic" ) ){
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.k1 = 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.k2 = 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.k3 = atof( the_token );
1821
1822 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1823 sprintf( painCave.errMsg,
1824 "Error parseing TorsionTypes: line %d\n", lineNum );
1825 painCave.isFatal = 1;
1826 simError();
1827 }
1828
1829 info.k4 = atof( the_token );
1830
1831 }
1832
1833 else{
1834 sprintf( painCave.errMsg,
1835 "Unknown DUFF torsion type \"%s\" at line %d\n",
1836 info.type,
1837 lineNum );
1838 painCave.isFatal = 1;
1839 simError();
1840 }
1841
1842 return 1;
1843 }
1844
1845 else return 0;
1846 }