ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/TraPPE_ExFF.cpp
Revision: 374
Committed: Thu Mar 20 19:50:42 2003 UTC (21 years, 3 months ago) by chuckv
File size: 41043 byte(s)
Log Message:
*** empty log message ***

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