ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/UseTheForce/DUFF.cpp
Revision: 1706
Committed: Thu Nov 4 16:20:28 2004 UTC (19 years, 7 months ago) by gezelter
File size: 44344 byte(s)
Log Message:
fixed useXXX in the entry_plug so that it only is
set if the atoms really are in the simulation

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