ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/TraPPE_ExFF.cpp
Revision: 365
Committed: Tue Mar 18 22:17:31 2003 UTC (21 years, 3 months ago) by mmeineke
File size: 40957 byte(s)
Log Message:
some initial bug fixes

File Contents

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