ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/InitializeFromFile.cpp
Revision: 874
Committed: Fri Nov 21 20:10:02 2003 UTC (20 years, 7 months ago) by mmeineke
File size: 16144 byte(s)
Log Message:
added a more verbose error message in SimInfo. Added a more informative error message in InitializeFromFile

File Contents

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