ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/TraPPE_ExFF.cpp
Revision: 367
Committed: Wed Mar 19 17:29:49 2003 UTC (21 years, 3 months ago) by mmeineke
File size: 40996 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 if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 ) &&
1132 !strcmp( nameC, key3 ) && !strcmp( nameD, key4 ) ) return this;
1133
1134 if( !strcmp( nameA, key4 ) && !strcmp( nameB, key3 ) &&
1135 !strcmp( nameC, key2 ) && !strcmp( nameD, key1 ) ) return this;
1136
1137 if( next != NULL ) return next->find(key1, key2, key3, key4);
1138 return NULL;
1139 }
1140
1141 void add( torsionStruct &info ){
1142
1143 // check for duplicates
1144 int dup = 0;
1145
1146 if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB ) &&
1147 !strcmp( nameC, info.nameC ) && !strcmp( nameD, info.nameD ) ) dup = 1;
1148
1149 if( !strcmp( nameA, info.nameD ) && !strcmp( nameB, info.nameC ) &&
1150 !strcmp( nameC, info.nameB ) && !strcmp( nameD, info.nameA ) ) dup = 1;
1151
1152 if(dup){
1153 sprintf( painCave.errMsg,
1154 "Duplicate TraPPE_Ex torsion type \"%s - %s - %s - %s\" found in "
1155 "the TraPPE_ExFF param file./n", nameA, nameB, nameC, nameD );
1156 painCave.isFatal = 1;
1157 simError();
1158 }
1159
1160 if( next != NULL ) next->add(info);
1161 else{
1162 next = new LinkedType();
1163 strcpy(next->nameA, info.nameA);
1164 strcpy(next->nameB, info.nameB);
1165 strcpy(next->nameC, info.nameC);
1166 strcpy(next->type, info.type);
1167 next->k1 = info.k1;
1168 next->k2 = info.k2;
1169 next->k3 = info.k3;
1170 next->k4 = info.k4;
1171 }
1172 }
1173
1174 #ifdef IS_MPI
1175
1176 void duplicate( torsionStruct &info ){
1177 strcpy(info.nameA, nameA);
1178 strcpy(info.nameB, nameB);
1179 strcpy(info.nameC, nameC);
1180 strcpy(info.nameD, nameD);
1181 strcpy(info.type, type);
1182 info.k1 = k1;
1183 info.k2 = k2;
1184 info.k3 = k3;
1185 info.k4 = k4;
1186 info.last = 0;
1187 }
1188
1189 #endif
1190
1191 char nameA[15];
1192 char nameB[15];
1193 char nameC[15];
1194 char nameD[15];
1195 char type[30];
1196 double k1, k2, k3, k4;
1197
1198 LinkedType* next;
1199 };
1200
1201 LinkedType* headTorsionType;
1202 LinkedType* currentTorsionType;
1203 torsionStruct info;
1204 info.last = 1; // initialize last to have the last set.
1205 // if things go well, last will be set to 0
1206
1207 int i, a, b, c, d, index;
1208 char* atomA;
1209 char* atomB;
1210 char* atomC;
1211 char* atomD;
1212 CubicTorsion* cTors;
1213
1214 SRI **the_sris;
1215 Atom** the_atoms;
1216 int nTorsions;
1217 the_sris = entry_plug->sr_interactions;
1218 the_atoms = entry_plug->atoms;
1219 nTorsions = entry_plug->n_torsions;
1220
1221 #ifdef IS_MPI
1222 if( worldRank == 0 ){
1223 #endif
1224
1225 // read in the torsion types.
1226
1227 headTorsionType = new LinkedType;
1228
1229 fastForward( "TorsionTypes", "initializeTorsions" );
1230
1231 // we are now at the torsionTypes section
1232
1233 eof_test = fgets( readLine, sizeof(readLine), frcFile );
1234 lineNum++;
1235
1236
1237 // read a line, and start parseing out the atom types
1238
1239 if( eof_test == NULL ){
1240 sprintf( painCave.errMsg,
1241 "Error in reading torsions from force file at line %d.\n",
1242 lineNum );
1243 painCave.isFatal = 1;
1244 simError();
1245 }
1246
1247 // stop reading at end of file, or at next section
1248 while( readLine[0] != '#' && eof_test != NULL ){
1249
1250 // toss comment lines
1251 if( readLine[0] != '!' ){
1252
1253 // the parser returns 0 if the line was blank
1254 if( parseTorsion( readLine, lineNum, info ) ){
1255 headTorsionType->add( info );
1256 }
1257 }
1258 eof_test = fgets( readLine, sizeof(readLine), frcFile );
1259 lineNum++;
1260 }
1261
1262 #ifdef IS_MPI
1263
1264 // send out the linked list to all the other processes
1265
1266 sprintf( checkPointMsg,
1267 "TraPPE_Ex torsion structures read successfully." );
1268 MPIcheckPoint();
1269
1270 currentTorsionType = headTorsionType;
1271 while( currentTorsionType != NULL ){
1272 currentTorsionType->duplicate( info );
1273 sendFrcStruct( &info, mpiTorsionStructType );
1274 currentTorsionType = currentTorsionType->next;
1275 }
1276 info.last = 1;
1277 sendFrcStruct( &info, mpiTorsionStructType );
1278
1279 }
1280
1281 else{
1282
1283 // listen for node 0 to send out the force params
1284
1285 MPIcheckPoint();
1286
1287 headTorsionType = new LinkedType;
1288 recieveFrcStruct( &info, mpiTorsionStructType );
1289 while( !info.last ){
1290
1291 headTorsionType->add( info );
1292 recieveFrcStruct( &info, mpiTorsionStructType );
1293 }
1294 }
1295 #endif // is_mpi
1296
1297 // initialize the Torsions
1298
1299 for( i=0; i<nTorsions; i++ ){
1300
1301 a = the_torsions[i].a;
1302 b = the_torsions[i].b;
1303 c = the_torsions[i].c;
1304 d = the_torsions[i].d;
1305
1306 atomA = the_atoms[a]->getType();
1307 atomB = the_atoms[b]->getType();
1308 atomC = the_atoms[c]->getType();
1309 atomD = the_atoms[d]->getType();
1310 currentTorsionType = headTorsionType->find( atomA, atomB, atomC, atomD );
1311 if( currentTorsionType == NULL ){
1312 sprintf( painCave.errMsg,
1313 "TorsionType error, %s - %s - %s - %s not found"
1314 " in force file.\n",
1315 atomA, atomB, atomC, atomD );
1316 painCave.isFatal = 1;
1317 simError();
1318 }
1319
1320 if( !strcmp( currentTorsionType->type, "cubic" ) ){
1321 index = i + entry_plug->n_bonds + entry_plug->n_bends;
1322
1323 cTors = new CubicTorsion( *the_atoms[a], *the_atoms[b],
1324 *the_atoms[c], *the_atoms[d] );
1325 cTors->setConstants( currentTorsionType->k1, currentTorsionType->k2,
1326 currentTorsionType->k3, currentTorsionType->k4 );
1327 the_sris[index] = cTors;
1328 }
1329 }
1330
1331
1332 // clean up the memory
1333
1334 delete headTorsionType;
1335
1336 #ifdef IS_MPI
1337 sprintf( checkPointMsg, "TraPPE_Ex torsions initialized succesfully" );
1338 MPIcheckPoint();
1339 #endif // is_mpi
1340
1341 }
1342
1343 void TraPPE_ExFF::fastForward( char* stopText, char* searchOwner ){
1344
1345 int foundText = 0;
1346 char* the_token;
1347
1348 rewind( frcFile );
1349 lineNum = 0;
1350
1351 eof_test = fgets( readLine, sizeof(readLine), frcFile );
1352 lineNum++;
1353 if( eof_test == NULL ){
1354 sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
1355 " file is empty.\n",
1356 searchOwner );
1357 painCave.isFatal = 1;
1358 simError();
1359 }
1360
1361
1362 while( !foundText ){
1363 while( eof_test != NULL && readLine[0] != '#' ){
1364 eof_test = fgets( readLine, sizeof(readLine), frcFile );
1365 lineNum++;
1366 }
1367 if( eof_test == NULL ){
1368 sprintf( painCave.errMsg,
1369 "Error fast forwarding force file for %s at "
1370 "line %d: file ended unexpectedly.\n",
1371 searchOwner,
1372 lineNum );
1373 painCave.isFatal = 1;
1374 simError();
1375 }
1376
1377 the_token = strtok( readLine, " ,;\t#\n" );
1378 foundText = !strcmp( stopText, the_token );
1379
1380 if( !foundText ){
1381 eof_test = fgets( readLine, sizeof(readLine), frcFile );
1382 lineNum++;
1383
1384 if( eof_test == NULL ){
1385 sprintf( painCave.errMsg,
1386 "Error fast forwarding force file for %s at "
1387 "line %d: file ended unexpectedly.\n",
1388 searchOwner,
1389 lineNum );
1390 painCave.isFatal = 1;
1391 simError();
1392 }
1393 }
1394 }
1395 }
1396
1397
1398 int TPE::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
1399
1400 char* the_token;
1401
1402 the_token = strtok( lineBuffer, " \n\t,;" );
1403 if( the_token != NULL ){
1404
1405 strcpy( info.name, the_token );
1406
1407 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1408 sprintf( painCave.errMsg,
1409 "Error parseing AtomTypes: line %d\n", lineNum );
1410 painCave.isFatal = 1;
1411 simError();
1412 }
1413
1414 info.isDipole = atoi( the_token );
1415
1416 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1417 sprintf( painCave.errMsg,
1418 "Error parseing AtomTypes: line %d\n", lineNum );
1419 painCave.isFatal = 1;
1420 simError();
1421 }
1422
1423 info.isSSD = atoi( the_token );
1424
1425 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1426 sprintf( painCave.errMsg,
1427 "Error parseing AtomTypes: line %d\n", lineNum );
1428 painCave.isFatal = 1;
1429 simError();
1430 }
1431
1432 info.mass = atof( the_token );
1433
1434 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1435 sprintf( painCave.errMsg,
1436 "Error parseing AtomTypes: line %d\n", lineNum );
1437 painCave.isFatal = 1;
1438 simError();
1439 }
1440
1441 info.epslon = atof( the_token );
1442
1443 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1444 sprintf( painCave.errMsg,
1445 "Error parseing AtomTypes: line %d\n", lineNum );
1446 painCave.isFatal = 1;
1447 simError();
1448 }
1449
1450 info.sigma = atof( the_token );
1451
1452 if( info.isDipole ){
1453
1454 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1455 sprintf( painCave.errMsg,
1456 "Error parseing AtomTypes: line %d\n", lineNum );
1457 painCave.isFatal = 1;
1458 simError();
1459 }
1460
1461 info.dipole = atof( the_token );
1462 }
1463 else info.dipole = 0.0;
1464
1465 if( info.isSSD ){
1466
1467 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1468 sprintf( painCave.errMsg,
1469 "Error parseing AtomTypes: line %d\n", lineNum );
1470 painCave.isFatal = 1;
1471 simError();
1472 }
1473
1474 info.w0 = atof( the_token );
1475
1476 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1477 sprintf( painCave.errMsg,
1478 "Error parseing AtomTypes: line %d\n", lineNum );
1479 painCave.isFatal = 1;
1480 simError();
1481 }
1482
1483 info.v0 = atof( the_token );
1484 }
1485 else info.v0 = info.w0 = 0.0;
1486
1487 return 1;
1488 }
1489 else return 0;
1490 }
1491
1492 int TPE::parseBond( char *lineBuffer, int lineNum, bondStruct &info ){
1493
1494 char* the_token;
1495
1496 the_token = strtok( lineBuffer, " \n\t,;" );
1497 if( the_token != NULL ){
1498
1499 strcpy( info.nameA, the_token );
1500
1501 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1502 sprintf( painCave.errMsg,
1503 "Error parseing BondTypes: line %d\n", lineNum );
1504 painCave.isFatal = 1;
1505 simError();
1506 }
1507
1508 strcpy( info.nameB, the_token );
1509
1510 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1511 sprintf( painCave.errMsg,
1512 "Error parseing BondTypes: line %d\n", lineNum );
1513 painCave.isFatal = 1;
1514 simError();
1515 }
1516
1517 strcpy( info.type, the_token );
1518
1519 if( !strcmp( info.type, "fixed" ) ){
1520 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1521 sprintf( painCave.errMsg,
1522 "Error parseing BondTypes: line %d\n", lineNum );
1523 painCave.isFatal = 1;
1524 simError();
1525 }
1526
1527 info.d0 = atof( the_token );
1528 }
1529 else{
1530 sprintf( painCave.errMsg,
1531 "Unknown TraPPE_Ex bond type \"%s\" at line %d\n",
1532 info.type,
1533 lineNum );
1534 painCave.isFatal = 1;
1535 simError();
1536 }
1537
1538 return 1;
1539 }
1540 else return 0;
1541 }
1542
1543
1544 int TPE::parseBend( char *lineBuffer, int lineNum, bendStruct &info ){
1545
1546 char* the_token;
1547
1548 the_token = strtok( lineBuffer, " \n\t,;" );
1549 if( the_token != NULL ){
1550
1551 strcpy( info.nameA, the_token );
1552
1553 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1554 sprintf( painCave.errMsg,
1555 "Error parseing BondTypes: line %d\n", lineNum );
1556 painCave.isFatal = 1;
1557 simError();
1558 }
1559
1560 strcpy( info.nameB, the_token );
1561
1562 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1563 sprintf( painCave.errMsg,
1564 "Error parseing BondTypes: line %d\n", lineNum );
1565 painCave.isFatal = 1;
1566 simError();
1567 }
1568
1569 strcpy( info.nameC, the_token );
1570
1571 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1572 sprintf( painCave.errMsg,
1573 "Error parseing BondTypes: line %d\n", lineNum );
1574 painCave.isFatal = 1;
1575 simError();
1576 }
1577
1578 strcpy( info.type, the_token );
1579
1580 if( !strcmp( info.type, "quadratic" ) ){
1581 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1582 sprintf( painCave.errMsg,
1583 "Error parseing BendTypes: line %d\n", lineNum );
1584 painCave.isFatal = 1;
1585 simError();
1586 }
1587
1588 info.k1 = atof( the_token );
1589
1590 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1591 sprintf( painCave.errMsg,
1592 "Error parseing BendTypes: line %d\n", lineNum );
1593 painCave.isFatal = 1;
1594 simError();
1595 }
1596
1597 info.k2 = atof( the_token );
1598
1599 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1600 sprintf( painCave.errMsg,
1601 "Error parseing BendTypes: line %d\n", lineNum );
1602 painCave.isFatal = 1;
1603 simError();
1604 }
1605
1606 info.k3 = atof( the_token );
1607
1608 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1609 sprintf( painCave.errMsg,
1610 "Error parseing BendTypes: line %d\n", lineNum );
1611 painCave.isFatal = 1;
1612 simError();
1613 }
1614
1615 info.t0 = atof( the_token );
1616 }
1617
1618 else{
1619 sprintf( painCave.errMsg,
1620 "Unknown TraPPE_Ex bend type \"%s\" at line %d\n",
1621 info.type,
1622 lineNum );
1623 painCave.isFatal = 1;
1624 simError();
1625 }
1626
1627 return 1;
1628 }
1629 else return 0;
1630 }
1631
1632 int TPE::parseTorsion( char *lineBuffer, int lineNum, torsionStruct &info ){
1633
1634 char* the_token;
1635
1636 the_token = strtok( lineBuffer, " \n\t,;" );
1637 if( the_token != NULL ){
1638
1639 strcpy( info.nameA, the_token );
1640
1641 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1642 sprintf( painCave.errMsg,
1643 "Error parseing TorsionTypes: line %d\n", lineNum );
1644 painCave.isFatal = 1;
1645 simError();
1646 }
1647
1648 strcpy( info.nameB, the_token );
1649
1650 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1651 sprintf( painCave.errMsg,
1652 "Error parseing TorsionTypes: line %d\n", lineNum );
1653 painCave.isFatal = 1;
1654 simError();
1655 }
1656
1657 strcpy( info.nameC, the_token );
1658
1659 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1660 sprintf( painCave.errMsg,
1661 "Error parseing TorsionTypes: line %d\n", lineNum );
1662 painCave.isFatal = 1;
1663 simError();
1664 }
1665
1666 strcpy( info.nameD, the_token );
1667
1668 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1669 sprintf( painCave.errMsg,
1670 "Error parseing TorsionTypes: line %d\n", lineNum );
1671 painCave.isFatal = 1;
1672 simError();
1673 }
1674
1675 strcpy( info.type, the_token );
1676
1677 if( !strcmp( info.type, "cubic" ) ){
1678 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1679 sprintf( painCave.errMsg,
1680 "Error parseing TorsionTypes: line %d\n", lineNum );
1681 painCave.isFatal = 1;
1682 simError();
1683 }
1684
1685 info.k1 = atof( the_token );
1686
1687 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1688 sprintf( painCave.errMsg,
1689 "Error parseing TorsionTypes: line %d\n", lineNum );
1690 painCave.isFatal = 1;
1691 simError();
1692 }
1693
1694 info.k2 = atof( the_token );
1695
1696 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1697 sprintf( painCave.errMsg,
1698 "Error parseing TorsionTypes: line %d\n", lineNum );
1699 painCave.isFatal = 1;
1700 simError();
1701 }
1702
1703 info.k3 = atof( the_token );
1704
1705 if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1706 sprintf( painCave.errMsg,
1707 "Error parseing TorsionTypes: line %d\n", lineNum );
1708 painCave.isFatal = 1;
1709 simError();
1710 }
1711
1712 info.k4 = atof( the_token );
1713
1714 }
1715
1716 else{
1717 sprintf( painCave.errMsg,
1718 "Unknown TraPPE_Ex torsion type \"%s\" at line %d\n",
1719 info.type,
1720 lineNum );
1721 painCave.isFatal = 1;
1722 simError();
1723 }
1724
1725 return 1;
1726 }
1727
1728 else return 0;
1729 }