ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/TraPPE_ExFF.cpp
Revision: 373
Committed: Thu Mar 20 19:01:44 2003 UTC (21 years, 3 months ago) by mmeineke
File size: 41376 byte(s)
Log Message:
fixed a bug in the torsions.

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