ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/TraPPE_ExFF.cpp
Revision: 433
Committed: Fri Mar 28 15:28:53 2003 UTC (21 years, 3 months ago) by mmeineke
File size: 40878 byte(s)
Log Message:
fixed long range interactions in Trappe

File Contents

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