ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/EAM_FF.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/EAM_FF.cpp (file contents):
Revision 499 by chuckv, Tue Apr 15 16:37:59 2003 UTC vs.
Revision 1224 by gezelter, Wed Jun 2 18:27:52 2004 UTC

# Line 1 | Line 1
1 < #include <cstdlib>
2 < #include <cstdio>
3 < #include <cstring>
1 > #include <stdlib.h>
2 > #include <stdio.h>
3 > #include <string.h>
4  
5   #include <iostream>
6   using namespace std;
# Line 28 | Line 28 | namespace EAM_NS{
28    typedef struct{
29      char name[15];
30      double mass;
31 <    double epslon;
32 <    double sigma;
31 >    double lattice_constant;
32 >    double eam_drho;  // The distance between each of the points indexed by rho.
33 >    double eam_rcut;  // The cutoff radius for eam.
34 >    double eam_dr;    // The distance between each of the rho points.    
35 >    int eam_nrho;  // Number of points indexed by rho
36 >    int eam_nr;    // The number of points based on r (Both Phi(r) and Rho(r)).
37 >    int eam_ident; // Atomic number
38      int ident;
39      int last;      //  0  -> default
40                     //  1  -> in MPI: tells nodes to stop listening
41    } atomStruct;
42  
43 <  int parseAtom( char *lineBuffer, int lineNum, atomStruct &info );
44 <  
43 >  int parseAtom( char *lineBuffer, int lineNum, atomStruct &info, char *eamPotFile );
44 >  int parseEAM( atomStruct &info, char *eamPotFile, double **eam_rvals,
45 >                double **eam_rhovals, double **eam_Frhovals);
46   #ifdef IS_MPI
47    
48    MPI_Datatype mpiAtomStructType;
# Line 48 | Line 54 | namespace EAM_NS{
54      LinkedAtomType(){
55        next = NULL;
56        name[0] = '\0';
57 +      eam_rvals    = NULL;
58 +      eam_rhovals  = NULL;
59 +      eam_Frhovals = NULL;
60      }
52 –    ~LinkedAtomType(){ if( next != NULL ) delete next; }
61  
62 +    ~LinkedAtomType(){
63 +      if( next != NULL ) delete next;
64 +      if( eam_rvals != NULL ) delete[] eam_rvals;
65 +      if( eam_rhovals != NULL ) delete[] eam_rhovals;
66 +      if( eam_Frhovals != NULL ) delete[] eam_Frhovals;
67 +    }
68 +
69      LinkedAtomType* find(char* key){
70        if( !strcmp(name, key) ) return this;
71        if( next != NULL ) return next->find(key);
# Line 58 | Line 73 | namespace EAM_NS{
73      }
74      
75  
76 <    void add( atomStruct &info ){
77 <    
76 >    void add( atomStruct &info, double *the_eam_rvals,
77 >              double *the_eam_rhovals,double *the_eam_Frhovals ){
78 >
79        // check for duplicates
80        
81        if( !strcmp( info.name, name ) ){
# Line 71 | Line 87 | namespace EAM_NS{
87          simError();
88        }
89        
90 <      if( next != NULL ) next->add(info);
90 >      if( next != NULL ) next->add(info, the_eam_rvals, the_eam_rhovals, the_eam_Frhovals);
91        else{
92          next = new LinkedAtomType();
93          strcpy(next->name, info.name);
94 <        next->mass     = info.mass;
95 <        next->epslon   = info.epslon;
96 <        next->sigma    = info.sigma;
97 <        next->ident    = info.ident;
94 >        next->mass             = info.mass;
95 >        next->lattice_constant = info.lattice_constant;
96 >        next->eam_nrho         = info.eam_nrho;
97 >        next->eam_drho         = info.eam_drho;
98 >        next->eam_nr           = info.eam_nr;
99 >        next->eam_dr           = info.eam_dr;
100 >        next->eam_rcut         = info.eam_rcut;
101 >        next->eam_ident        = info.eam_ident;
102 >        next->ident            = info.ident;
103 >
104 >        next->eam_rvals    = the_eam_rvals;
105 >        next->eam_rhovals  = the_eam_rhovals;
106 >        next->eam_Frhovals = the_eam_Frhovals;
107        }
108      }
109      
# Line 87 | Line 112 | namespace EAM_NS{
112      
113      void duplicate( atomStruct &info ){
114        strcpy(info.name, name);
115 <      info.mass     = mass;
116 <      info.epslon   = epslon;
117 <      info.sigma    = sigma;
118 <      info.ident    = ident;
119 <      info.last     = 0;
115 >      info.mass             = mass;
116 >      info.lattice_constant = lattice_constant;
117 >      info.eam_nrho         = eam_nrho;
118 >      info.eam_drho         = eam_drho;
119 >      info.eam_nr           = eam_nr;
120 >      info.eam_dr           = eam_dr;
121 >      info.eam_rcut         = eam_rcut;
122 >      info.eam_ident        = eam_ident;
123 >      info.ident            = ident;
124 >      info.last             = 0;
125      }
126  
127  
# Line 99 | Line 129 | namespace EAM_NS{
129  
130      char name[15];
131      double mass;
132 <    double epslon;
133 <    double sigma;
132 >    double lattice_constant;
133 >    int eam_nrho; // Number of points indexed by rho
134 >    double eam_drho; // The distance between each of the points indexed by rho.
135 >    int eam_nr;   // The number of points based on r (Both Phi(r) and Rho(r)).
136 >    double eam_dr;   // The distance between each of the rho points.
137 >    double eam_rcut; // The cutoff radius for eam.
138 >
139 >    double *eam_rvals;    // Z of r values
140 >    double *eam_rhovals;  // rho of r values
141 >    double *eam_Frhovals; // F of rho values
142 >    int eam_ident;        // eam identity (atomic number)
143      int ident;
144      LinkedAtomType* next;
145    };
# Line 110 | Line 149 | namespace EAM_NS{
149  
150   }
151  
152 < using namespace LJ_NS;
152 > using namespace EAM_NS;
153  
154   //****************************************************************
155   // begins the actual forcefield stuff.  
# Line 123 | Line 162 | EAM_FF::EAM_FF(){
162    char* ffPath_env = "FORCE_PARAM_PATH";
163    char* ffPath;
164    char temp[200];
126 –  char errMsg[1000];
165  
166    headAtomType = NULL;
167    currentAtomType = NULL;
168  
169 +  // Set eamRcut to 0.0
170 +  eamRcut = 0.0;
171 +
172    // do the funtion wrapping
173    wrapMeFF( this );
174  
# Line 138 | Line 179 | EAM_FF::EAM_FF(){
179    // Init the atomStruct mpi type
180  
181    atomStruct atomProto; // mpiPrototype
182 <  int atomBC[3] = {15,3,2};  // block counts
182 >  int atomBC[3] = {15,5,5};  // block counts
183    MPI_Aint atomDspls[3];           // displacements
184    MPI_Datatype atomMbrTypes[3];    // member mpi types
185  
186    MPI_Address(&atomProto.name, &atomDspls[0]);
187    MPI_Address(&atomProto.mass, &atomDspls[1]);
188 <  MPI_Address(&atomProto.ident, &atomDspls[2]);
188 >  MPI_Address(&atomProto.eam_nrho, &atomDspls[2]);
189    
190    atomMbrTypes[0] = MPI_CHAR;
191    atomMbrTypes[1] = MPI_DOUBLE;
# Line 187 | Line 228 | EAM_FF::EAM_FF(){
228        frcFile = fopen( fileName, "r" );
229        
230        if( frcFile == NULL ){
231 <        
231 >
232          sprintf( painCave.errMsg,
233 <                 "Error opening the force field parameter file: %s\n"
234 <                 "Have you tried setting the FORCE_PARAM_PATH environment "
235 <                 "vairable?\n",
233 >                 "Error opening the force field parameter file:\n"
234 >                 "\t%s\n"
235 >                 "\tHave you tried setting the FORCE_PARAM_PATH environment "
236 >                 "variable?\n",
237                   fileName );
238 +        painCave.severity = OOPSE_ERROR;
239          painCave.isFatal = 1;
240          simError();
241        }
# Line 201 | Line 244 | EAM_FF::EAM_FF(){
244   #ifdef IS_MPI
245    }
246    
247 <  sprintf( checkPointMsg, "LJ_FF file opened sucessfully." );
247 >  sprintf( checkPointMsg, "EAM_FF file opened sucessfully." );
248    MPIcheckPoint();
249    
250   #endif // is_mpi
# Line 223 | Line 266 | EAM_FF::~EAM_FF(){
266   #endif // is_mpi
267   }
268  
269 < void EAM_FF::initForceField( int ljMixRule ){
269 >
270 > void EAM_FF::calcRcut( void ){
271    
272 + #ifdef IS_MPI
273 +  double tempEamRcut = eamRcut;
274 +  MPI_Allreduce( &tempEamRcut, &eamRcut, 1, MPI_DOUBLE, MPI_MAX,
275 +                 MPI_COMM_WORLD);
276 + #endif  //is_mpi
277 +  entry_plug->setDefaultRcut(eamRcut);
278 + }
279 +
280 +
281 + void EAM_FF::initForceField( int ljMixRule ){
282    initFortran( ljMixRule, 0 );
283   }
284  
# Line 243 | Line 297 | void EAM_FF::cleanMe( void ){
297   #endif // is_mpi
298   }
299  
300 +
301   void EAM_FF::readParams( void ){
302  
303    atomStruct info;
304    info.last = 1; // initialize last to have the last set.
305                   // if things go well, last will be set to 0
306  
252 –  int i;
307    int identNum;
308 +  double *eam_rvals;    // Z of r values
309 +  double *eam_rhovals;  // rho of r values
310 +  double *eam_Frhovals; // F of rho values
311 +  char eamPotFile[1000];
312    
313  
314    bigSigma = 0.0;
# Line 262 | Line 320 | void EAM_FF::readParams( void ){
320  
321      headAtomType = new LinkedAtomType;
322      
323 <    fastForward( "AtomTypes", "initializeAtoms" );
323 >    fastForward( "AtomTypes", "eam atom readParams" );
324  
325      // we are now at the AtomTypes section.
326      
# Line 282 | Line 340 | void EAM_FF::readParams( void ){
340      
341      identNum = 1;
342      // stop reading at end of file, or at next section
343 +
344      while( readLine[0] != '#' && eof_test != NULL ){
345  
346        // toss comment lines
347        if( readLine[0] != '!' ){
348          
349          // the parser returns 0 if the line was blank
350 <        if( parseAtom( readLine, lineNum, info ) ){
350 >        if( parseAtom( readLine, lineNum, info, eamPotFile ) ){
351 >          parseEAM(info,eamPotFile, &eam_rvals,
352 >                   &eam_rhovals, &eam_Frhovals);
353            info.ident = identNum;
354 <          headAtomType->add( info );;
354 >          headAtomType->add( info, eam_rvals,
355 >                             eam_rhovals,eam_Frhovals );
356            identNum++;
357          }
358        }
359        eof_test = fgets( readLine, sizeof(readLine), frcFile );
360        lineNum++;
361      }
362 +    
363 +    
364  
365   #ifdef IS_MPI
366 +  
367      
368      // send out the linked list to all the other processes
369  
370      sprintf( checkPointMsg,
371 <             "LJ_FF atom structures read successfully." );
371 >             "EAM_FF atom structures read successfully." );
372      MPIcheckPoint();
373  
374      currentAtomType = headAtomType->next; //skip the first element who is a place holder.
# Line 314 | Line 379 | void EAM_FF::readParams( void ){
379  
380        sendFrcStruct( &info, mpiAtomStructType );
381  
382 +      // We have to now broadcast the Arrays
383 +      MPI_Bcast(currentAtomType->eam_rvals,
384 +                currentAtomType->eam_nr,
385 +                MPI_DOUBLE,0,MPI_COMM_WORLD);
386 +      MPI_Bcast(currentAtomType->eam_rhovals,
387 +                currentAtomType->eam_nr,
388 +                MPI_DOUBLE,0,MPI_COMM_WORLD);
389 +      MPI_Bcast(currentAtomType->eam_Frhovals,
390 +                currentAtomType->eam_nrho,
391 +                MPI_DOUBLE,0,MPI_COMM_WORLD);
392 +
393        sprintf( checkPointMsg,
394 <               "successfully sent lJ force type: \"%s\"\n",
394 >               "successfully sent EAM force type: \"%s\"\n",
395                 info.name );
396        MPIcheckPoint();
397  
# Line 333 | Line 409 | void EAM_FF::readParams( void ){
409      MPIcheckPoint();
410  
411      headAtomType = new LinkedAtomType;
412 <    recieveFrcStruct( &info, mpiAtomStructType );
413 <    
412 >    receiveFrcStruct( &info, mpiAtomStructType );
413 >
414      while( !info.last ){
415 +      
416 +      // allocate the arrays
417  
418 +      eam_rvals    = new double[info.eam_nr];
419 +      eam_rhovals  = new double[info.eam_nr];
420 +      eam_Frhovals = new double[info.eam_nrho];
421  
422 +      // We have to now broadcast the Arrays
423 +      MPI_Bcast(eam_rvals,
424 +                info.eam_nr,
425 +                MPI_DOUBLE,0,MPI_COMM_WORLD);
426 +      MPI_Bcast(eam_rhovals,
427 +                info.eam_nr,
428 +                MPI_DOUBLE,0,MPI_COMM_WORLD);
429 +      MPI_Bcast(eam_Frhovals,
430 +                info.eam_nrho,
431 +                MPI_DOUBLE,0,MPI_COMM_WORLD);
432 +      
433  
434 <      headAtomType->add( info );
434 >      headAtomType->add( info, eam_rvals, eam_rhovals, eam_Frhovals );
435        
436        MPIcheckPoint();
437  
438 <      recieveFrcStruct( &info, mpiAtomStructType );
438 >      receiveFrcStruct( &info, mpiAtomStructType );
439 >
440 >
441      }
442    }
443   #endif // is_mpi
# Line 353 | Line 447 | void EAM_FF::readParams( void ){
447    int isError;
448  
449    // dummy variables
450 <  int isLJ = 1;
450 >  int isLJ = 0;
451    int isDipole = 0;
452    int isSSD = 0;
453    int isGB = 0;
454 +  int isEAM = 1;
455 +  int isCharge = 0;
456    double dipole = 0.0;
457 +  double charge = 0.0;
458 +  double eamSigma = 0.0;
459 +  double eamEpslon = 0.0;
460    
461 <  currentAtomType = headAtomType;
461 >  currentAtomType = headAtomType->next;
462    while( currentAtomType != NULL ){
463      
464      if( currentAtomType->name[0] != '\0' ){
# Line 369 | Line 468 | void EAM_FF::readParams( void ){
468                   &isSSD,
469                   &isDipole,
470                   &isGB,
471 <                 &(currentAtomType->epslon),
472 <                 &(currentAtomType->sigma),
471 >                 &isEAM,
472 >                 &isCharge,
473 >                 &eamEpslon,
474 >                 &eamSigma,
475 >                 &charge,
476                   &dipole,
477                   &isError );
478        if( isError ){
# Line 384 | Line 486 | void EAM_FF::readParams( void ){
486      currentAtomType = currentAtomType->next;
487    }
488        
489 <  entry_plug->useLJ = 1;
489 >  entry_plug->useLJ = 0;
490 >  entry_plug->useEAM = 1;
491 >  // Walk down again and send out EAM type
492 >  currentAtomType = headAtomType->next;
493 >  while( currentAtomType != NULL ){
494 >    
495 >    if( currentAtomType->name[0] != '\0' ){
496 >      isError = 0;
497  
498 +      newEAMtype( &(currentAtomType->lattice_constant),
499 +                  &(currentAtomType->eam_nrho),
500 +                  &(currentAtomType->eam_drho),
501 +                  &(currentAtomType->eam_nr),
502 +                  &(currentAtomType->eam_dr),
503 +                  &(currentAtomType->eam_rcut),
504 +                  currentAtomType->eam_rvals,
505 +                  currentAtomType->eam_rhovals,
506 +                  currentAtomType->eam_Frhovals,
507 +                  &(currentAtomType->eam_ident),
508 +                  &isError);
509 +
510 +      if( isError ){
511 +        sprintf( painCave.errMsg,
512 +                 "Error initializing the \"%s\" atom type in fortran EAM\n",
513 +                 currentAtomType->name );
514 +        painCave.isFatal = 1;
515 +        simError();
516 +      }
517 +    }
518 +    currentAtomType = currentAtomType->next;
519 +  }
520 +
521 +
522 +
523   #ifdef IS_MPI
524    sprintf( checkPointMsg,
525 <           "LJ_FF atom structures successfully sent to fortran\n" );
525 >           "EAM_FF atom structures successfully sent to fortran\n" );
526    MPIcheckPoint();
527   #endif // is_mpi
528  
529 +
530 +
531   }
532  
533  
534   void EAM_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
535    
536    int i;
537 <
537 >  
538    // initialize the atoms
539    
404 –
405 –  Atom* thisAtom;
406 –
540    for( i=0; i<nAtoms; i++ ){
541      
542      currentAtomType = headAtomType->find( the_atoms[i]->getType() );
# Line 416 | Line 549 | void EAM_FF::initializeAtoms( int nAtoms, Atom** the_a
549      }
550      
551      the_atoms[i]->setMass( currentAtomType->mass );
419 –    the_atoms[i]->setEpslon( currentAtomType->epslon );
420 –    the_atoms[i]->setSigma( currentAtomType->sigma );
552      the_atoms[i]->setIdent( currentAtomType->ident );
422 –    the_atoms[i]->setLJ();
553  
554 <    if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
554 >    if (eamRcut < currentAtomType->eam_rcut) eamRcut = currentAtomType->eam_rcut;
555 >    
556    }
557   }
558  
559 < void LJ_FF::initializeBonds( int nBonds, Bond** BondArray,
559 > void EAM_FF::initializeBonds( int nBonds, Bond** BondArray,
560                               bond_pair* the_bonds ){
561    
562      if( nBonds ){
# Line 514 | Line 645 | void EAM_FF::fastForward( char* stopText, char* search
645  
646  
647  
648 < int EAM_NS::parseAtom( char *lineBuffer, int lineNum,  atomStruct &info ){
648 > int EAM_NS::parseAtom( char *lineBuffer, int lineNum,   atomStruct &info, char *eamPotFile ){
649  
650    char* the_token;
651    
# Line 539 | Line 670 | int EAM_NS::parseAtom( char *lineBuffer, int lineNum,
670        simError();
671      }
672          
673 <    info.epslon = atof( the_token );
674 <          
675 <    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
676 <      sprintf( painCave.errMsg,
677 <               "Error parseing AtomTypes: line %d\n", lineNum );
673 >    strcpy( eamPotFile, the_token );
674 >    return 1;
675 >  }
676 >  else return 0;
677 > }
678 >
679 > int EAM_NS::parseEAM(atomStruct &info, char *eamPotFile,
680 >                     double **eam_rvals,
681 >                     double **eam_rhovals,
682 >                     double **eam_Frhovals){
683 >  double* myEam_rvals;
684 >  double* myEam_rhovals;
685 >  double* myEam_Frhovals;
686 >
687 >  char* ffPath_env = "FORCE_PARAM_PATH";
688 >  char* ffPath;
689 >  char* the_token;
690 >  char* eam_eof_test;
691 >  FILE *eamFile;
692 >  const int BUFFERSIZE = 3000;
693 >
694 >  char temp[200];
695 >  int linenumber;
696 >  int nReadLines;
697 >  char eam_read_buffer[BUFFERSIZE];
698 >
699 >
700 >  int i,j;
701 >
702 >  linenumber = 0;
703 >
704 >  // Open eam file
705 >  eamFile = fopen( eamPotFile, "r" );
706 >  
707 >  
708 >  if( eamFile == NULL ){
709 >    
710 >      // next see if the force path enviorment variable is set
711 >    
712 >    ffPath = getenv( ffPath_env );
713 >    if( ffPath == NULL ) {
714 >      STR_DEFINE(ffPath, FRC_PATH );
715 >    }
716 >    
717 >    
718 >    strcpy( temp, ffPath );
719 >    strcat( temp, "/" );
720 >    strcat( temp, eamPotFile );
721 >    strcpy( eamPotFile, temp );
722 >    
723 >    eamFile = fopen( eamPotFile, "r" );
724 >
725 >    
726 >    
727 >    if( eamFile == NULL ){
728 >      
729 >      sprintf( painCave.errMsg,
730 >               "Error opening the EAM force parameter file: %s\n"
731 >               "Have you tried setting the FORCE_PARAM_PATH environment "
732 >               "variable?\n",
733 >               eamPotFile );
734        painCave.isFatal = 1;
735        simError();
736      }
550 –        
551 –    info.sigma = atof( the_token );
552 –    
553 –    return 1;
737    }
738 <  else return 0;
738 >
739 >  // First line is a comment line, read and toss it....
740 >  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
741 >  linenumber++;
742 >  if(eam_eof_test == NULL){
743 >    sprintf( painCave.errMsg,
744 >             "error in reading commment in %s\n", eamPotFile);
745 >    painCave.isFatal = 1;
746 >    simError();
747 >  }
748 >
749 >
750 >
751 >  // The Second line contains atomic number, atomic mass and a lattice constant
752 >  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
753 >  linenumber++;
754 >  if(eam_eof_test == NULL){
755 >    sprintf( painCave.errMsg,
756 >             "error in reading Identifier line in %s\n", eamPotFile);
757 >    painCave.isFatal = 1;
758 >    simError();
759 >  }
760 >
761 >
762 >
763 >    
764 >  if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
765 >    sprintf( painCave.errMsg,
766 >             "Error parseing EAM ident  line in %s\n", eamPotFile );
767 >    painCave.isFatal = 1;
768 >    simError();
769 >  }
770 >  
771 >  info.eam_ident = atoi( the_token );
772 >
773 >  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
774 >    sprintf( painCave.errMsg,
775 >             "Error parseing EAM mass in %s\n", eamPotFile );
776 >    painCave.isFatal = 1;
777 >    simError();
778 >  }
779 >  info.mass = atof( the_token);
780 >
781 >  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
782 >    sprintf( painCave.errMsg,
783 >             "Error parseing EAM Lattice Constant %s\n", eamPotFile );
784 >    painCave.isFatal = 1;
785 >    simError();
786 >  }
787 >  info.lattice_constant = atof( the_token);
788 >
789 >  // Next line is nrho, drho, nr, dr and rcut
790 >  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
791 >  if(eam_eof_test == NULL){
792 >    sprintf( painCave.errMsg,
793 >             "error in reading number of points line in %s\n", eamPotFile);
794 >    painCave.isFatal = 1;
795 >    simError();
796 >  }
797 >
798 >  if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
799 >    sprintf( painCave.errMsg,
800 >             "Error parseing EAM nrho: line in %s\n", eamPotFile );
801 >    painCave.isFatal = 1;
802 >    simError();
803 >  }
804 >  
805 >  info.eam_nrho = atoi( the_token );
806 >  
807 >  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
808 >    sprintf( painCave.errMsg,
809 >             "Error parseing EAM drho in %s\n", eamPotFile );
810 >    painCave.isFatal = 1;
811 >    simError();
812 >  }
813 >  info.eam_drho = atof( the_token);
814 >
815 >  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
816 >    sprintf( painCave.errMsg,
817 >             "Error parseing EAM # r in %s\n", eamPotFile );
818 >    painCave.isFatal = 1;
819 >    simError();
820 >  }
821 >  info.eam_nr = atoi( the_token);
822 >  
823 >  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
824 >    sprintf( painCave.errMsg,
825 >             "Error parseing EAM dr in %s\n", eamPotFile );
826 >    painCave.isFatal = 1;
827 >    simError();
828 >  }
829 >  info.eam_dr = atof( the_token);
830 >
831 >  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
832 >    sprintf( painCave.errMsg,
833 >             "Error parseing EAM rcut in %s\n", eamPotFile );
834 >    painCave.isFatal = 1;
835 >    simError();
836 >  }
837 >  info.eam_rcut = atof( the_token);
838 >
839 >
840 >
841 >
842 >
843 >  // Ok now we have to allocate point arrays and read in number of points
844 >  // Index the arrays for fortran, starting at 1
845 >  myEam_Frhovals = new double[info.eam_nrho];
846 >  myEam_rvals    = new double[info.eam_nr];
847 >  myEam_rhovals  = new double[info.eam_nr];
848 >
849 >  // Parse F of rho vals.
850 >
851 >  // Assume for now that we have a complete number of lines
852 >  nReadLines = int(info.eam_nrho/5);
853 >  
854 >
855 >
856 >  for (i=0;i<nReadLines;i++){
857 >    j = i*5;
858 >
859 >    // Read next line
860 >    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
861 >    linenumber++;
862 >    if(eam_eof_test == NULL){
863 >      sprintf( painCave.errMsg,
864 >               "error in reading EAM file %s at line %d\n",
865 >               eamPotFile,linenumber);
866 >      painCave.isFatal = 1;
867 >      simError();
868 >    }
869 >    
870 >    // Parse 5 values on each line into array
871 >    // Value 1
872 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
873 >      sprintf( painCave.errMsg,
874 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
875 >      painCave.isFatal = 1;
876 >      simError();
877 >    }
878 >
879 >    myEam_Frhovals[j+0] = atof( the_token );
880 >    
881 >    // Value 2
882 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
883 >      sprintf( painCave.errMsg,
884 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
885 >      painCave.isFatal = 1;
886 >      simError();
887 >    }
888 >
889 >    myEam_Frhovals[j+1] = atof( the_token );
890 >    
891 >    // Value 3
892 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
893 >      sprintf( painCave.errMsg,
894 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
895 >      painCave.isFatal = 1;
896 >      simError();
897 >    }
898 >
899 >    myEam_Frhovals[j+2] = atof( the_token );
900 >    
901 >    // Value 4
902 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
903 >      sprintf( painCave.errMsg,
904 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
905 >      painCave.isFatal = 1;
906 >      simError();
907 >    }
908 >
909 >    myEam_Frhovals[j+3] = atof( the_token );
910 >
911 >    // Value 5
912 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
913 >      sprintf( painCave.errMsg,
914 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
915 >      painCave.isFatal = 1;
916 >      simError();
917 >    }
918 >
919 >    myEam_Frhovals[j+4] = atof( the_token );
920 >    
921 >  }
922 >  // Parse Z of r vals
923 >  
924 >  // Assume for now that we have a complete number of lines
925 >  nReadLines = int(info.eam_nr/5);
926 >
927 >  for (i=0;i<nReadLines;i++){
928 >    j = i*5;
929 >
930 >    // Read next line
931 >    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
932 >    linenumber++;
933 >    if(eam_eof_test == NULL){
934 >      sprintf( painCave.errMsg,
935 >               "error in reading EAM file %s at line %d\n",
936 >               eamPotFile,linenumber);
937 >      painCave.isFatal = 1;
938 >      simError();
939 >    }
940 >    
941 >    // Parse 5 values on each line into array
942 >    // Value 1
943 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
944 >      sprintf( painCave.errMsg,
945 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
946 >      painCave.isFatal = 1;
947 >      simError();
948 >    }
949 >    
950 >    myEam_rvals[j+0] = atof( the_token );
951 >
952 >    // Value 2
953 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
954 >      sprintf( painCave.errMsg,
955 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
956 >      painCave.isFatal = 1;
957 >      simError();
958 >    }
959 >  
960 >    myEam_rvals[j+1] = atof( the_token );
961 >
962 >    // Value 3
963 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
964 >      sprintf( painCave.errMsg,
965 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
966 >      painCave.isFatal = 1;
967 >      simError();
968 >    }
969 >  
970 >    myEam_rvals[j+2] = atof( the_token );
971 >
972 >    // Value 4
973 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
974 >      sprintf( painCave.errMsg,
975 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
976 >      painCave.isFatal = 1;
977 >      simError();
978 >    }
979 >  
980 >    myEam_rvals[j+3] = atof( the_token );
981 >
982 >    // Value 5
983 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
984 >      sprintf( painCave.errMsg,
985 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
986 >      painCave.isFatal = 1;
987 >      simError();
988 >    }
989 >  
990 >    myEam_rvals[j+4] = atof( the_token );
991 >
992 >  }
993 >  // Parse rho of r vals
994 >
995 >  // Assume for now that we have a complete number of lines
996 >
997 >  for (i=0;i<nReadLines;i++){
998 >    j = i*5;
999 >
1000 >    // Read next line
1001 >    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
1002 >    linenumber++;
1003 >    if(eam_eof_test == NULL){
1004 >      sprintf( painCave.errMsg,
1005 >               "error in reading EAM file %s at line %d\n",
1006 >               eamPotFile,linenumber);
1007 >      painCave.isFatal = 1;
1008 >      simError();
1009 >    }
1010 >  
1011 >    // Parse 5 values on each line into array
1012 >    // Value 1
1013 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
1014 >      sprintf( painCave.errMsg,
1015 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
1016 >      painCave.isFatal = 1;
1017 >      simError();
1018 >    }
1019 >  
1020 >    myEam_rhovals[j+0] = atof( the_token );
1021 >
1022 >    // Value 2
1023 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
1024 >      sprintf( painCave.errMsg,
1025 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
1026 >      painCave.isFatal = 1;
1027 >      simError();
1028 >    }
1029 >  
1030 >    myEam_rhovals[j+1] = atof( the_token );
1031 >
1032 >    // Value 3
1033 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
1034 >      sprintf( painCave.errMsg,
1035 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
1036 >      painCave.isFatal = 1;
1037 >      simError();
1038 >    }
1039 >  
1040 >    myEam_rhovals[j+2] = atof( the_token );
1041 >
1042 >    // Value 4
1043 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
1044 >      sprintf( painCave.errMsg,
1045 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
1046 >      painCave.isFatal = 1;
1047 >      simError();
1048 >    }
1049 >  
1050 >    myEam_rhovals[j+3] = atof( the_token );
1051 >
1052 >    // Value 5
1053 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
1054 >      sprintf( painCave.errMsg,
1055 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
1056 >      painCave.isFatal = 1;
1057 >      simError();
1058 >    }
1059 >  
1060 >    myEam_rhovals[j+4] = atof( the_token );
1061 >
1062 >  }
1063 >  *eam_rvals = myEam_rvals;
1064 >  *eam_rhovals = myEam_rhovals;
1065 >  *eam_Frhovals = myEam_Frhovals;
1066 >
1067 >  fclose(eamFile);
1068 >  return 0;
1069   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines