ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/InitializeFromFile.cpp
Revision: 830
Committed: Tue Oct 28 16:20:28 2003 UTC (20 years, 8 months ago) by gezelter
File size: 16405 byte(s)
Log Message:
fixes for compatibility

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