ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DumpReader.cpp
Revision: 646
Committed: Tue Jul 22 20:05:30 2003 UTC (20 years, 11 months ago) by mmeineke
File size: 16732 byte(s)
Log Message:
making some changes to read dump files

File Contents

# Content
1 #define _FILE_OFFSET_BITS 64
2
3 #include <iostream>
4 #include <cmath>
5 #include <vector>
6
7 #include <stdio.h>
8 #include <stdlib.h>
9 #include <string.h>
10 #include <unistd.h>
11 #include <sys/types.h>
12 #include <sys/stat.h>
13
14 #include "ReadWrite.hpp"
15 #include "simError.h"
16
17 #ifdef IS_MPI
18 #include <mpi.h>
19 #include "mpiSimulation.hpp"
20 #define TAKE_THIS_TAG_CHAR 0
21 #define TAKE_THIS_TAG_INT 1
22
23 namespace dumpRead{
24 void nodeZeroError( void );
25 void anonymousNodeDie( void );
26
27
28 class FilePos{
29 public:
30
31
32 private:
33
34
35 };
36
37 vector<FilePos> frameStart;
38 }
39
40 using namespace dumpRead;
41
42 #endif // is_mpi
43
44 DumpReader :: DumpReader( char *in_name ){
45 #ifdef IS_MPI
46 if (worldRank == 0) {
47 #endif
48
49 c_in_file = fopen(in_name, "r");
50 if(c_in_file == NULL){
51 sprintf(painCave.errMsg,
52 "Cannot open file: %s\n", in_name);
53 painCave.isFatal = 1;
54 simError();
55 }
56
57 strcpy( c_in_name, in_name);
58 #ifdef IS_MPI
59 }
60 strcpy( checkPointMsg, "Dump file opened for reading successfully." );
61 MPIcheckPoint();
62 #endif
63 return;
64 }
65
66 DumpReader :: ~DumpReader( ){
67 #ifdef IS_MPI
68 if (worldRank == 0) {
69 #endif
70 int error;
71 error = fclose( c_in_file );
72 if( error ){
73 sprintf( painCave.errMsg,
74 "Error closing %s\n", c_in_name );
75 simError();
76 }
77 #ifdef IS_MPI
78 }
79 strcpy( checkPointMsg, "Dump file closed successfully." );
80 MPIcheckPoint();
81 #endif
82
83 return;
84 }
85
86
87
88 void DumpReader :: getFrame( SimInfo* the_simnfo, int whichFrame){
89
90 int i, j, done, which_node, which_atom; // loop counter
91
92 const int BUFFERSIZE = 2000; // size of the read buffer
93 int n_atoms; // the number of atoms
94 char read_buffer[BUFFERSIZE]; //the line buffer for reading
95 #ifdef IS_MPI
96 char send_buffer[BUFFERSIZE];
97 #endif
98
99 char *eof_test; // ptr to see when we reach the end of the file
100 char *parseErr;
101 int procIndex;
102 double boxMat[9];
103 double theBoxMat3[3][3];
104
105 simnfo = the_simnfo;
106
107 #ifndef IS_MPI
108 eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
109 if( eof_test == NULL ){
110 sprintf( painCave.errMsg,
111 "DumpReader error: error reading 1st line of \"%s\"\n",
112 c_in_name );
113 painCave.isFatal = 1;
114 simError();
115 }
116
117 n_atoms = atoi( read_buffer );
118
119 Atom **atoms = simnfo->atoms;
120 DirectionalAtom* dAtom;
121
122 if( n_atoms != simnfo->n_atoms ){
123 sprintf( painCave.errMsg,
124 "DumpReader error. %s n_atoms, %d, "
125 "does not match the BASS file's n_atoms, %d.\n",
126 c_in_name, n_atoms, simnfo->n_atoms );
127 painCave.isFatal = 1;
128 simError();
129 }
130
131 //read the box mat from the comment line
132
133 eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
134 if(eof_test == NULL){
135 sprintf( painCave.errMsg,
136 "error in reading commment in %s\n", c_in_name);
137 painCave.isFatal = 1;
138 simError();
139 }
140
141 parseErr = parseCommentLine( read_buffer, time, boxMat );
142 if( parseErr != NULL ){
143 strcpy( painCave.errMsg, parseErr );
144 painCave.isFatal = 1;
145 simError();
146 }
147
148 simnfo->setTime( time );
149
150 for(i=0;i<3;i++)
151 for(j=0;j<3;j++) theBoxMat3[i][j] = boxMat[3*j+i];
152
153
154 simnfo->setBoxM( theBoxMat3 );
155
156
157 for( i=0; i < n_atoms; i++){
158
159 eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
160 if(eof_test == NULL){
161 sprintf(painCave.errMsg,
162 "error in reading file %s\n"
163 "natoms = %d; index = %d\n"
164 "error reading the line from the file.\n",
165 c_in_name, n_atoms, i );
166 painCave.isFatal = 1;
167 simError();
168 }
169
170
171 parseErr = parseDumpLine( read_buffer, i );
172 if( parseErr != NULL ){
173 strcpy( painCave.errMsg, parseErr );
174 painCave.isFatal = 1;
175 simError();
176 }
177 }
178
179
180 // MPI Section of code..........
181 #else //IS_MPI
182
183 // first thing first, suspend fatalities.
184 painCave.isEventLoop = 1;
185
186 int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
187 int haveError;
188
189 MPI_Status istatus;
190 int *AtomToProcMap = mpiSim->getAtomToProcMap();
191
192
193 haveError = 0;
194 if (worldRank == 0) {
195
196 eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
197 if( eof_test == NULL ){
198 sprintf( painCave.errMsg,
199 "Error reading 1st line of %d \n ",c_in_name);
200 haveError = 1;
201 simError();
202 }
203
204 n_atoms = atoi( read_buffer );
205
206 Atom **atoms = simnfo->atoms;
207 DirectionalAtom* dAtom;
208
209 // Check to see that the number of atoms in the intial configuration file is the
210 // same as declared in simBass.
211
212 if( n_atoms != mpiSim->getTotAtoms() ){
213 sprintf( painCave.errMsg,
214 "Initialize from File error. %s n_atoms, %d, "
215 "does not match the BASS file's n_atoms, %d.\n",
216 c_in_name, n_atoms, simnfo->n_atoms );
217 haveError= 1;
218 simError();
219 }
220
221 //read the time and boxMat from the comment line
222
223 eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
224 if(eof_test == NULL){
225 sprintf( painCave.errMsg,
226 "error in reading commment in %s\n", c_in_name);
227 haveError = 1;
228 simError();
229 }
230
231 parseErr = parseCommentLine( read_buffer, time, boxMat );
232 if( parseErr != NULL ){
233 strcpy( painCave.errMsg, parseErr );
234 haveError = 1;
235 simError();
236 }
237
238 MPI_Bcast(time, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD );
239
240 MPI_Bcast(boxMat, 9, MPI_DOUBLE, 0, MPI_COMM_WORLD );
241
242 if(haveError) nodeZeroError();
243
244 for (i=0 ; i < mpiSim->getTotAtoms(); i++) {
245
246 eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
247 if(eof_test == NULL){
248 sprintf(painCave.errMsg,
249 "error in reading file %s\n"
250 "natoms = %d; index = %d\n"
251 "error reading the line from the file.\n",
252 c_in_name, n_atoms, i );
253 haveError= 1;
254 simError();
255 }
256
257 if(haveError) nodeZeroError();
258
259 // Get the Node number which wants this atom:
260 which_node = AtomToProcMap[i];
261 if (which_node == 0) {
262 parseErr = parseDumpLine( read_buffer, i );
263 if( parseErr != NULL ){
264 strcpy( painCave.errMsg, parseErr );
265 haveError = 1;
266 simError();
267 }
268 if(haveError) nodeZeroError();
269 }
270
271 else {
272
273 myStatus = 1;
274 MPI_Send(&myStatus, 1, MPI_INT, which_node,
275 TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
276 MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node,
277 TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
278 MPI_Send(&i, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
279 MPI_COMM_WORLD);
280 MPI_Recv(&myStatus, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
281 MPI_COMM_WORLD, &istatus);
282
283 if(!myStatus) nodeZeroError();
284 }
285 }
286 myStatus = -1;
287 for (j = 0; j < mpiSim->getNumberProcessors(); j++) {
288 MPI_Send( &myStatus, 1, MPI_INT, j,
289 TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
290 }
291
292 } else {
293
294 MPI_Bcast(time, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
295 MPI_Bcast(boxMat, 9, MPI_DOUBLE, 0, MPI_COMM_WORLD);
296
297 done = 0;
298 while (!done) {
299
300 MPI_Recv(&myStatus, 1, MPI_INT, 0,
301 TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
302
303 if(!myStatus) anonymousNodeDie();
304
305 if(myStatus < 0) break;
306
307 MPI_Recv(read_buffer, BUFFERSIZE, MPI_CHAR, 0,
308 TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
309 MPI_Recv(&which_atom, 1, MPI_INT, 0,
310 TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
311
312 myStatus = 1;
313 parseErr = parseDumpLine( read_buffer, which_atom );
314 if( parseErr != NULL ){
315 strcpy( painCave.errMsg, parseErr );
316 myStatus = 0;;
317 simError();
318 }
319
320 MPI_Send( &myStatus, 1, MPI_INT, 0,
321 TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
322
323 }
324 }
325
326 // last thing last, enable fatalities.
327 painCave.isEventLoop = 0;
328
329 simnfo->setTime( time );
330
331 for(i=0;i<3;i++)
332 for(j=0;j<3;j++) theBoxMat3[i][j] = boxMat[3*j+i];
333
334 simnfo->setBoxM( theBoxMat3 );
335
336
337 #endif
338 }
339
340 char* DumpReader::parseDumpLine(char* readLine, int globalIndex){
341
342 char *foo; // the pointer to the current string token
343
344 double rx, ry, rz; // position place holders
345 double vx, vy, vz; // velocity placeholders
346 double q[4]; // the quaternions
347 double jx, jy, jz; // angular velocity placeholders;
348 double qSqr, qLength; // needed to normalize the quaternion vector.
349
350 Atom **atoms = simnfo->atoms;
351 DirectionalAtom* dAtom;
352
353 int j, n_atoms, atomIndex;
354
355 #ifdef IS_MPI
356 n_atoms = mpiSim->getTotAtoms();
357 atomIndex=-1;
358 for (j=0; j < mpiSim->getMyNlocal(); j++) {
359 if (atoms[j]->getGlobalIndex() == globalIndex) atomIndex = j;
360 }
361 if (atomIndex == -1) {
362 sprintf( painCave.errMsg,
363 "Initialize from file error. Atom at index %d "
364 "in file %s does not exist on processor %d .\n",
365 globalIndex, c_in_name, mpiSim->getMyNode() );
366 return strdup( painCave.errMsg );
367 }
368 #else
369 n_atoms = simnfo->n_atoms;
370 atomIndex = globalIndex;
371 #endif // is_mpi
372
373 // set the string tokenizer
374
375 foo = strtok(readLine, " ,;\t");
376
377 // check the atom name to the current atom
378
379 if( strcmp( foo, atoms[atomIndex]->getType() ) ){
380 sprintf( painCave.errMsg,
381 "Initialize from file error. Atom %s at index %d "
382 "in file %s does not"
383 " match the BASS atom %s.\n",
384 foo, atomIndex, c_in_name, atoms[atomIndex]->getType() );
385 return strdup( painCave.errMsg );
386 }
387
388 // get the positions
389
390 foo = strtok(NULL, " ,;\t");
391 if(foo == NULL){
392 sprintf( painCave.errMsg,
393 "error in reading postition x from %s\n"
394 "natoms = %d, index = %d\n",
395 c_in_name, n_atoms, atomIndex );
396 return strdup( painCave.errMsg );
397 }
398 rx = atof( foo );
399
400 foo = strtok(NULL, " ,;\t");
401 if(foo == NULL){
402 sprintf( painCave.errMsg,
403 "error in reading postition y from %s\n"
404 "natoms = %d, index = %d\n",
405 c_in_name, n_atoms, atomIndex );
406 return strdup( painCave.errMsg );
407 }
408 ry = atof( foo );
409
410 foo = strtok(NULL, " ,;\t");
411 if(foo == NULL){
412 sprintf( painCave.errMsg,
413 "error in reading postition z from %s\n"
414 "natoms = %d, index = %d\n",
415 c_in_name, n_atoms, atomIndex );
416 return strdup( painCave.errMsg );
417 }
418 rz = atof( foo );
419
420
421 // get the velocities
422
423 foo = strtok(NULL, " ,;\t");
424 if(foo == NULL){
425 sprintf( painCave.errMsg,
426 "error in reading velocity x from %s\n"
427 "natoms = %d, index = %d\n",
428 c_in_name, n_atoms, atomIndex );
429 return strdup( painCave.errMsg );
430 }
431 vx = atof( foo );
432
433 foo = strtok(NULL, " ,;\t");
434 if(foo == NULL){
435 sprintf( painCave.errMsg,
436 "error in reading velocity y from %s\n"
437 "natoms = %d, index = %d\n",
438 c_in_name, n_atoms, atomIndex );
439 return strdup( painCave.errMsg );
440 }
441 vy = atof( foo );
442
443 foo = strtok(NULL, " ,;\t");
444 if(foo == NULL){
445 sprintf( painCave.errMsg,
446 "error in reading velocity z from %s\n"
447 "natoms = %d, index = %d\n",
448 c_in_name, n_atoms, atomIndex );
449 return strdup( painCave.errMsg );
450 }
451 vz = atof( foo );
452
453
454 // get the quaternions
455
456 if( atoms[atomIndex]->isDirectional() ){
457
458 foo = strtok(NULL, " ,;\t");
459 if(foo == NULL){
460 sprintf(painCave.errMsg,
461 "error in reading quaternion 0 from %s\n"
462 "natoms = %d, index = %d\n",
463 c_in_name, n_atoms, atomIndex );
464 return strdup( painCave.errMsg );
465 }
466 q[0] = atof( foo );
467
468 foo = strtok(NULL, " ,;\t");
469 if(foo == NULL){
470 sprintf( painCave.errMsg,
471 "error in reading quaternion 1 from %s\n"
472 "natoms = %d, index = %d\n",
473 c_in_name, n_atoms, atomIndex );
474 return strdup( painCave.errMsg );
475 }
476 q[1] = atof( foo );
477
478 foo = strtok(NULL, " ,;\t");
479 if(foo == NULL){
480 sprintf( painCave.errMsg,
481 "error in reading quaternion 2 from %s\n"
482 "natoms = %d, index = %d\n",
483 c_in_name, n_atoms, atomIndex );
484 return strdup( painCave.errMsg );
485 }
486 q[2] = atof( foo );
487
488 foo = strtok(NULL, " ,;\t");
489 if(foo == NULL){
490 sprintf( painCave.errMsg,
491 "error in reading quaternion 3 from %s\n"
492 "natoms = %d, index = %d\n",
493 c_in_name, n_atoms, atomIndex );
494 return strdup( painCave.errMsg );
495 }
496 q[3] = atof( foo );
497
498 // get the angular velocities
499
500 foo = strtok(NULL, " ,;\t");
501 if(foo == NULL){
502 sprintf( painCave.errMsg,
503 "error in reading angular momentum jx from %s\n"
504 "natoms = %d, index = %d\n",
505 c_in_name, n_atoms, atomIndex );
506 return strdup( painCave.errMsg );
507 }
508 jx = atof( foo );
509
510 foo = strtok(NULL, " ,;\t");
511 if(foo == NULL){
512 sprintf( painCave.errMsg,
513 "error in reading angular momentum jy from %s\n"
514 "natoms = %d, index = %d\n",
515 c_in_name, n_atoms, atomIndex );
516 return strdup( painCave.errMsg );
517 }
518 jy = atof(foo );
519
520 foo = strtok(NULL, " ,;\t");
521 if(foo == NULL){
522 sprintf( painCave.errMsg,
523 "error in reading angular momentum jz from %s\n"
524 "natoms = %d, index = %d\n",
525 c_in_name, n_atoms, atomIndex );
526 return strdup( painCave.errMsg );
527 }
528 jz = atof( foo );
529
530 dAtom = ( DirectionalAtom* )atoms[atomIndex];
531
532 // check that the quaternion vector is normalized
533
534 qSqr = (q[0] * q[0]) + (q[1] * q[1]) + (q[2] * q[2]) + (q[3] * q[3]);
535
536 qLength = sqrt( qSqr );
537 q[0] = q[0] / qLength;
538 q[1] = q[1] / qLength;
539 q[2] = q[2] / qLength;
540 q[3] = q[3] / qLength;
541
542 dAtom->setQ( q );
543
544 // add the angular velocities
545
546 dAtom->setJx( jx );
547 dAtom->setJy( jy );
548 dAtom->setJz( jz );
549 }
550
551 // add the positions and velocities to the atom
552
553 atoms[atomIndex]->setX( rx );
554 atoms[atomIndex]->setY( ry );
555 atoms[atomIndex]->setZ( rz );
556
557 atoms[atomIndex]->set_vx( vx );
558 atoms[atomIndex]->set_vy( vy );
559 atoms[atomIndex]->set_vz( vz );
560
561 return NULL;
562 }
563
564
565 char* DumpReader::parseCommentLine(char* readLine, double time,
566 double boxMat[9]){
567
568 char *foo; // the pointer to the current string token
569 int j;
570
571 // set the string tokenizer
572
573 foo = strtok(readLine, " ,;\t");
574 if(foo == NULL){
575 sprintf( painCave.errMsg,
576 "error in reading time from %s\n",
577 c_in_name );
578 return strdup( painCave.errMsg );
579 }
580 time = atof( foo );
581
582 // get the Hx vector
583
584 foo = strtok(NULL, " ,;\t");
585 if(foo == NULL){
586 sprintf( painCave.errMsg,
587 "error in reading Hx[0] from %s\n",
588 c_in_name );
589 return strdup( painCave.errMsg );
590 }
591 boxMat[0] = atof( foo );
592
593 foo = strtok(NULL, " ,;\t");
594 if(foo == NULL){
595 sprintf( painCave.errMsg,
596 "error in reading Hx[1] from %s\n",
597 c_in_name );
598 return strdup( painCave.errMsg );
599 }
600 boxMat[1] = atof( foo );
601
602 foo = strtok(NULL, " ,;\t");
603 if(foo == NULL){
604 sprintf( painCave.errMsg,
605 "error in reading Hx[2] from %s\n",
606 c_in_name );
607 return strdup( painCave.errMsg );
608 }
609 boxMat[2] = atof( foo );
610
611 // get the Hy vector
612
613 foo = strtok(NULL, " ,;\t");
614 if(foo == NULL){
615 sprintf( painCave.errMsg,
616 "error in reading Hy[0] from %s\n",
617 c_in_name );
618 return strdup( painCave.errMsg );
619 }
620 boxMat[3] = atof( foo );
621
622 foo = strtok(NULL, " ,;\t");
623 if(foo == NULL){
624 sprintf( painCave.errMsg,
625 "error in reading Hy[1] from %s\n",
626 c_in_name );
627 return strdup( painCave.errMsg );
628 }
629 boxMat[4] = atof( foo );
630
631 foo = strtok(NULL, " ,;\t");
632 if(foo == NULL){
633 sprintf( painCave.errMsg,
634 "error in reading Hy[2] from %s\n",
635 c_in_name );
636 return strdup( painCave.errMsg );
637 }
638 boxMat[5] = atof( foo );
639
640 // get the Hz vector
641
642 foo = strtok(NULL, " ,;\t");
643 if(foo == NULL){
644 sprintf( painCave.errMsg,
645 "error in reading Hz[0] from %s\n",
646 c_in_name );
647 return strdup( painCave.errMsg );
648 }
649 boxMat[6] = atof( foo );
650
651 foo = strtok(NULL, " ,;\t");
652 if(foo == NULL){
653 sprintf( painCave.errMsg,
654 "error in reading Hz[1] from %s\n",
655 c_in_name );
656 return strdup( painCave.errMsg );
657 }
658 boxMat[7] = atof( foo );
659
660 foo = strtok(NULL, " ,;\t");
661 if(foo == NULL){
662 sprintf( painCave.errMsg,
663 "error in reading Hz[2] from %s\n",
664 c_in_name );
665 return strdup( painCave.errMsg );
666 }
667 boxMat[8] = atof( foo );
668
669 return NULL;
670 }
671
672
673 #ifdef IS_MPI
674
675 // a couple of functions to let us escape the read loop
676
677 void initFile::nodeZeroError( void ){
678 int j, myStatus;
679
680 myStatus = 0;
681 for (j = 0; j < mpiSim->getNumberProcessors(); j++) {
682 MPI_Send( &myStatus, 1, MPI_INT, j,
683 TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
684 }
685
686
687 MPI_Finalize();
688 exit (0);
689
690 }
691
692 void initFile::anonymousNodeDie( void ){
693
694 MPI_Finalize();
695 exit (0);
696 }
697
698 #endif //is_mpi