ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/InitializeFromFile.cpp
Revision: 787
Committed: Thu Sep 25 19:27:15 2003 UTC (20 years, 9 months ago) by mmeineke
File size: 16375 byte(s)
Log Message:
cleaned things with gcc -Wall and g++ -Wall

File Contents

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