ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/TraPPE_ExFF.cpp
Revision: 359
Committed: Mon Mar 17 21:38:57 2003 UTC (21 years, 3 months ago) by mmeineke
File size: 41455 byte(s)
Log Message:
adding more to the interface

File Contents

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