ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/InitializeFromFile.cpp
Revision: 690
Committed: Tue Aug 12 21:44:06 2003 UTC (20 years, 10 months ago) by mmeineke
File size: 16530 byte(s)
Log Message:
fixed a really annoying bug in Directional Atom, where mu was getting written to pseudorandom memory location.

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