ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/InitializeFromFile.cpp
Revision: 586
Committed: Wed Jul 9 22:14:06 2003 UTC (21 years ago) by mmeineke
File size: 16075 byte(s)
Log Message:
Bug fixing NPTi and NPTf. there is some error in the caclulation of HmatInverse.

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