ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/InitializeFromFile.cpp
Revision: 1108
Committed: Wed Apr 14 15:37:41 2004 UTC (20 years, 2 months ago) by tim
File size: 15806 byte(s)
Log Message:
Change DumpWriter and InitFromFile

File Contents

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