ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/InitializeFromFile.cpp
Revision: 837
Committed: Wed Oct 29 00:19:10 2003 UTC (20 years, 8 months ago) by tim
File size: 16072 byte(s)
Log Message:
add chi and eta to the comment line of dump file.

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 strcpy( checkPointMsg, "Infile opened for reading successfully." );
49 MPIcheckPoint();
50 #endif
51 return;
52 }
53
54 InitializeFromFile::~InitializeFromFile( ){
55 #ifdef IS_MPI
56 if (worldRank == 0) {
57 #endif
58 int error;
59 error = fclose( c_in_file );
60 if( error ){
61 sprintf( painCave.errMsg,
62 "Error closing %s\n", c_in_name );
63 simError();
64 }
65 #ifdef IS_MPI
66 }
67 strcpy( checkPointMsg, "Infile closed successfully." );
68 MPIcheckPoint();
69 #endif
70
71 return;
72 }
73
74
75 void InitializeFromFile :: readInit( SimInfo* the_simnfo ){
76
77 int i, j;
78
79 #ifdef IS_MPI
80 int done, which_node, which_atom; // loop counter
81 #endif //is_mpi
82
83 const int BUFFERSIZE = 2000; // size of the read buffer
84 int n_atoms; // the number of atoms
85 char read_buffer[BUFFERSIZE]; //the line buffer for reading
86
87 char *eof_test; // ptr to see when we reach the end of the file
88 char *parseErr;
89
90 simnfo = the_simnfo;
91
92
93 #ifndef IS_MPI
94 eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
95 if( eof_test == NULL ){
96 sprintf( painCave.errMsg,
97 "InitializeFromFile error: error reading 1st line of \"%s\"\n",
98 c_in_name );
99 painCave.isFatal = 1;
100 simError();
101 }
102
103 n_atoms = atoi( read_buffer );
104
105 if( n_atoms != simnfo->n_atoms ){
106 sprintf( painCave.errMsg,
107 "Initialize from File error. %s n_atoms, %d, "
108 "does not match the BASS file's n_atoms, %d.\n",
109 c_in_name, n_atoms, simnfo->n_atoms );
110 painCave.isFatal = 1;
111 simError();
112 }
113
114 //read the box mat from the comment line
115
116 eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
117 if(eof_test == NULL){
118 sprintf( painCave.errMsg,
119 "error in reading commment in %s\n", c_in_name);
120 painCave.isFatal = 1;
121 simError();
122 }
123
124
125
126 parseErr = parseCommentLine( read_buffer, simnfo);
127 if( parseErr != NULL ){
128 strcpy( painCave.errMsg, parseErr );
129 painCave.isFatal = 1;
130 simError();
131 }
132
133 //parse dump lines
134
135 for( i=0; i < n_atoms; i++){
136
137 eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
138 if(eof_test == NULL){
139 sprintf(painCave.errMsg,
140 "error in reading file %s\n"
141 "natoms = %d; index = %d\n"
142 "error reading the line from the file.\n",
143 c_in_name, n_atoms, i );
144 painCave.isFatal = 1;
145 simError();
146 }
147
148
149 parseErr = parseDumpLine( read_buffer, i );
150 if( parseErr != NULL ){
151 strcpy( painCave.errMsg, parseErr );
152 painCave.isFatal = 1;
153 simError();
154 }
155 }
156
157
158 // MPI Section of code..........
159 #else //IS_MPI
160
161 // first thing first, suspend fatalities.
162 painCave.isEventLoop = 1;
163
164 int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
165 int haveError;
166
167 MPI_Status istatus;
168 int *AtomToProcMap = mpiSim->getAtomToProcMap();
169
170
171 haveError = 0;
172 if (worldRank == 0) {
173
174 eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
175 if( eof_test == NULL ){
176 sprintf( painCave.errMsg,
177 "Error reading 1st line of %s \n ",c_in_name);
178 haveError = 1;
179 simError();
180 }
181
182 n_atoms = atoi( read_buffer );
183
184 // Check to see that the number of atoms in the intial configuration file is the
185 // same as declared in simBass.
186
187 if( n_atoms != mpiSim->getTotAtoms() ){
188 sprintf( painCave.errMsg,
189 "Initialize from File error. %s n_atoms, %d, "
190 "does not match the BASS file's n_atoms, %d.\n",
191 c_in_name, n_atoms, simnfo->n_atoms );
192 haveError= 1;
193 simError();
194 }
195
196 //read the boxMat from the comment line
197
198 eof_test = fgets(read_buffer, sizeof(read_buffer), c_in_file);
199 if(eof_test == NULL){
200 sprintf( painCave.errMsg,
201 "error in reading commment in %s\n", c_in_name);
202 haveError = 1;
203 simError();
204 }
205
206 //Every single processor will parse the comment line by itself
207 //By using this way, we might lose some efficiency, but if we want to add
208 //more parameters into comment line, we only need to modify function
209 //parseCommentLine
210
211 MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, 0, MPI_COMM_WORLD);
212
213 parseErr = parseCommentLine( read_buffer, simnfo);
214
215 if( parseErr != NULL ){
216 strcpy( painCave.errMsg, parseErr );
217 haveError = 1;
218 simError();
219 }
220
221 for (i=0 ; i < mpiSim->getTotAtoms(); i++) {
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 file %s\n"
227 "natoms = %d; index = %d\n"
228 "error reading the line from the file.\n",
229 c_in_name, n_atoms, i );
230 haveError= 1;
231 simError();
232 }
233
234 if(haveError) nodeZeroError();
235
236 // Get the Node number which wants this atom:
237 which_node = AtomToProcMap[i];
238 if (which_node == 0) {
239 parseErr = parseDumpLine( read_buffer, i );
240 if( parseErr != NULL ){
241 strcpy( painCave.errMsg, parseErr );
242 haveError = 1;
243 simError();
244 }
245 if(haveError) nodeZeroError();
246 }
247
248 else {
249
250 myStatus = 1;
251 MPI_Send(&myStatus, 1, MPI_INT, which_node,
252 TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
253 MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node,
254 TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
255 MPI_Send(&i, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
256 MPI_COMM_WORLD);
257 MPI_Recv(&myStatus, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
258 MPI_COMM_WORLD, &istatus);
259
260 if(!myStatus) nodeZeroError();
261 }
262 }
263 myStatus = -1;
264 for (j = 0; j < mpiSim->getNumberProcessors(); j++) {
265 MPI_Send( &myStatus, 1, MPI_INT, j,
266 TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
267 }
268
269 } else {
270
271 MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, 0, MPI_COMM_WORLD);
272
273 parseErr = parseCommentLine( read_buffer, simnfo);
274
275 if( parseErr != NULL ){
276 strcpy( painCave.errMsg, parseErr );
277 haveError = 1;
278 simError();
279 }
280
281
282 done = 0;
283 while (!done) {
284
285 MPI_Recv(&myStatus, 1, MPI_INT, 0,
286 TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
287
288 if(!myStatus) anonymousNodeDie();
289
290 if(myStatus < 0) break;
291
292 MPI_Recv(read_buffer, BUFFERSIZE, MPI_CHAR, 0,
293 TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
294 MPI_Recv(&which_atom, 1, MPI_INT, 0,
295 TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
296
297 myStatus = 1;
298 parseErr = parseDumpLine( read_buffer, which_atom );
299 if( parseErr != NULL ){
300 strcpy( painCave.errMsg, parseErr );
301 myStatus = 0;;
302 simError();
303 }
304
305 MPI_Send( &myStatus, 1, MPI_INT, 0,
306 TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
307
308 }
309 }
310
311 // last thing last, enable fatalities.
312 painCave.isEventLoop = 0;
313
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::parseCommentLine(char* readLine, SimInfo* entry_plug){
542
543 double currTime;
544 double boxMat[9];
545 double theBoxMat3[3][3];
546 double chi;
547 double integralOfChidt;
548 double eta[9];
549
550 char *foo; // the pointer to the current string token
551
552 // set the string tokenizer
553
554 foo = strtok(readLine, " ,;\t");
555 // set the timeToken.
556
557 if(foo == NULL){
558 sprintf( painCave.errMsg,
559 "error in reading Time from %s\n",
560 c_in_name );
561 return strdup( painCave.errMsg );
562 }
563
564 currTime = atof( foo );
565 entry_plug->setTime( currTime );
566
567 //get H-Matrix
568
569 for(int i = 0 ; i < 9; i++){
570 foo = strtok(NULL, " ,;\t");
571 if(foo == NULL){
572 sprintf( painCave.errMsg,
573 "error in reading H[%d] from %s\n", i, c_in_name );
574 return strdup( painCave.errMsg );
575 }
576 boxMat[i] = atof( foo );
577 }
578
579 for(int i=0;i<3;i++)
580 for(int j=0;j<3;j++) theBoxMat3[i][j] = boxMat[3*j+i];
581
582 //set H-Matrix
583 entry_plug->setBoxM( theBoxMat3 );
584
585 //get chi and integralOfChidt, they should appear by pair
586 foo = strtok(NULL, " ,;\t\n");
587 if(foo != NULL){
588 chi = atof(foo);
589
590 foo = strtok(NULL, " ,;\t\n");
591 if(foo == NULL){
592 sprintf( painCave.errMsg,
593 "chi and integralOfChidt should appear by pair in %s\n", c_in_name );
594 return strdup( painCave.errMsg );
595 }
596 integralOfChidt = atof( foo );
597
598 //push chi and integralOfChidt into SimInfo::properties which can be
599 //retrieved by integrator later
600 DoubleData* chiValue = new DoubleData();
601 chiValue->setID(CHIVALUE_ID);
602 chiValue->setData(chi);
603 entry_plug->addProperty(chiValue);
604
605 DoubleData* integralOfChidtValue = new DoubleData();
606 integralOfChidtValue->setID(INTEGRALOFCHIDT_ID);
607 integralOfChidtValue->setData(integralOfChidt);
608 entry_plug->addProperty(integralOfChidtValue);
609
610 }
611 else
612 return NULL;
613
614 //get eta
615 for(int i = 0 ; i < 9; i++){
616 foo = strtok(NULL, " ,;\t");
617 if(foo == NULL){
618 sprintf( painCave.errMsg,
619 "error in reading eta[%d] from %s\n", i, c_in_name );
620 return strdup( painCave.errMsg );
621 }
622 eta[i] = atof( foo );
623 }
624
625 //push eta into SimInfo::properties which can be
626 //retrieved by integrator later
627 //entry_plug->setBoxM( theBoxMat3 );
628 DoubleArrayData* etaValue = new DoubleArrayData();
629 etaValue->setID(ETAVALUE_ID);
630 etaValue->setData(eta, 9);
631 entry_plug->addProperty(etaValue);
632
633
634 return NULL;
635 }
636
637 #ifdef IS_MPI
638
639 // a couple of functions to let us escape the read loop
640
641 void initFile::nodeZeroError( void ){
642 int j, myStatus;
643
644 myStatus = 0;
645 for (j = 0; j < mpiSim->getNumberProcessors(); j++) {
646 MPI_Send( &myStatus, 1, MPI_INT, j,
647 TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
648 }
649
650
651 MPI_Finalize();
652 exit (0);
653
654 }
655
656 void initFile::anonymousNodeDie( void ){
657
658 MPI_Finalize();
659 exit (0);
660 }
661
662 #endif //is_mpi