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 1261 by gezelter, Fri Jun 11 14:14:10 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 ) ){
82          sprintf( painCave.errMsg,
83 <                 "Duplicate LJ atom type \"%s\" found in "
84 <                 "the LJ_FF param file./n",
83 >                 "Duplicate EAM atom type \"%s\" found in "
84 >                 "the EAM_FF param file./n",
85                   name );
86          painCave.isFatal = 1;
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.  
156   //****************************************************************
157  
158 + EAM_FF::EAM_FF() {
159 +  EAM_FF("");
160 + }
161  
162 < EAM_FF::EAM_FF(){
162 > EAM_FF::EAM_FF(char* the_variant){
163  
164    char fileName[200];
165    char* ffPath_env = "FORCE_PARAM_PATH";
166    char* ffPath;
167    char temp[200];
126 –  char errMsg[1000];
168  
169    headAtomType = NULL;
170    currentAtomType = NULL;
171  
172 +  // Set eamRcut to 0.0
173 +  eamRcut = 0.0;
174 +
175    // do the funtion wrapping
176    wrapMeFF( this );
177  
# Line 138 | Line 182 | EAM_FF::EAM_FF(){
182    // Init the atomStruct mpi type
183  
184    atomStruct atomProto; // mpiPrototype
185 <  int atomBC[3] = {15,3,2};  // block counts
185 >  int atomBC[3] = {15,5,5};  // block counts
186    MPI_Aint atomDspls[3];           // displacements
187    MPI_Datatype atomMbrTypes[3];    // member mpi types
188  
189    MPI_Address(&atomProto.name, &atomDspls[0]);
190    MPI_Address(&atomProto.mass, &atomDspls[1]);
191 <  MPI_Address(&atomProto.ident, &atomDspls[2]);
191 >  MPI_Address(&atomProto.eam_nrho, &atomDspls[2]);
192    
193    atomMbrTypes[0] = MPI_CHAR;
194    atomMbrTypes[1] = MPI_DOUBLE;
# Line 160 | Line 204 | EAM_FF::EAM_FF(){
204    if( worldRank == 0 ){
205   #endif
206      
207 <    // generate the force file name
207 >    // generate the force file name  
208 >
209 >    strcpy( fileName, "EAM" );
210 >
211 >    if (strlen(the_variant) > 0) {
212 >      has_variant = 1;
213 >      strcpy( variant, the_variant);
214 >      strcat( fileName, ".");
215 >      strcat( fileName, variant );
216 >
217 >      sprintf( painCave.errMsg,
218 >               "Using %s variant of EAM force field.\n",
219 >               variant );
220 >      painCave.severity = OOPSE_INFO;
221 >      painCave.isFatal = 0;
222 >      simError();
223 >    }
224 >    strcat( fileName, ".frc");
225 >
226 >    //fprintf( stderr,"Trying to open %s\n", fileName );
227      
165 –    strcpy( fileName, "EAM_FF.frc" );
166 –    //    fprintf( stderr,"Trying to open %s\n", fileName );
167 –    
228      // attempt to open the file in the current directory first.
229      
230      frcFile = fopen( fileName, "r" );
# Line 187 | Line 247 | EAM_FF::EAM_FF(){
247        frcFile = fopen( fileName, "r" );
248        
249        if( frcFile == NULL ){
250 <        
250 >        
251          sprintf( painCave.errMsg,
252 <                 "Error opening the force field parameter file: %s\n"
253 <                 "Have you tried setting the FORCE_PARAM_PATH environment "
254 <                 "vairable?\n",
252 >                 "Error opening the force field parameter file:\n"
253 >                 "\t%s\n"
254 >                 "\tHave you tried setting the FORCE_PARAM_PATH environment "
255 >                 "variable?\n",
256                   fileName );
257 +        painCave.severity = OOPSE_ERROR;
258          painCave.isFatal = 1;
259          simError();
260        }
261      }
262 +
263      
264   #ifdef IS_MPI
265    }
266    
267 <  sprintf( checkPointMsg, "LJ_FF file opened sucessfully." );
267 >  sprintf( checkPointMsg, "EAM_FF file opened sucessfully." );
268    MPIcheckPoint();
269    
270   #endif // is_mpi
# Line 223 | Line 286 | EAM_FF::~EAM_FF(){
286   #endif // is_mpi
287   }
288  
289 < void EAM_FF::initForceField( int ljMixRule ){
289 >
290 > void EAM_FF::calcRcut( void ){
291    
292 + #ifdef IS_MPI
293 +  double tempEamRcut = eamRcut;
294 +  MPI_Allreduce( &tempEamRcut, &eamRcut, 1, MPI_DOUBLE, MPI_MAX,
295 +                 MPI_COMM_WORLD);
296 + #endif  //is_mpi
297 +  entry_plug->setDefaultRcut(eamRcut);
298 + }
299 +
300 +
301 + void EAM_FF::initForceField( int ljMixRule ){
302    initFortran( ljMixRule, 0 );
303   }
304  
# Line 243 | Line 317 | void EAM_FF::cleanMe( void ){
317   #endif // is_mpi
318   }
319  
320 +
321   void EAM_FF::readParams( void ){
322  
323    atomStruct info;
324    info.last = 1; // initialize last to have the last set.
325                   // if things go well, last will be set to 0
326  
252 –  int i;
327    int identNum;
328 +  double *eam_rvals;    // Z of r values
329 +  double *eam_rhovals;  // rho of r values
330 +  double *eam_Frhovals; // F of rho values
331 +  char eamPotFile[1000];
332    
333  
334    bigSigma = 0.0;
# Line 261 | Line 339 | void EAM_FF::readParams( void ){
339      // read in the atom types.
340  
341      headAtomType = new LinkedAtomType;
264 –    
265 –    fastForward( "AtomTypes", "initializeAtoms" );
342  
343 +    fastForward( "AtomTypes", "eam atom readParams" );
344 +
345      // we are now at the AtomTypes section.
346      
347      eof_test = fgets( readLine, sizeof(readLine), frcFile );
# Line 282 | Line 360 | void EAM_FF::readParams( void ){
360      
361      identNum = 1;
362      // stop reading at end of file, or at next section
363 +
364      while( readLine[0] != '#' && eof_test != NULL ){
365  
366        // toss comment lines
367        if( readLine[0] != '!' ){
368          
369          // the parser returns 0 if the line was blank
370 <        if( parseAtom( readLine, lineNum, info ) ){
370 >        if( parseAtom( readLine, lineNum, info, eamPotFile ) ){
371 >          parseEAM(info,eamPotFile, &eam_rvals,
372 >                   &eam_rhovals, &eam_Frhovals);
373            info.ident = identNum;
374 <          headAtomType->add( info );;
374 >          headAtomType->add( info, eam_rvals,
375 >                             eam_rhovals,eam_Frhovals );
376            identNum++;
377          }
378        }
379        eof_test = fgets( readLine, sizeof(readLine), frcFile );
380        lineNum++;
381      }
382 +    
383 +    
384  
385   #ifdef IS_MPI
386 +  
387      
388      // send out the linked list to all the other processes
389  
390      sprintf( checkPointMsg,
391 <             "LJ_FF atom structures read successfully." );
391 >             "EAM_FF atom structures read successfully." );
392      MPIcheckPoint();
393  
394      currentAtomType = headAtomType->next; //skip the first element who is a place holder.
# Line 314 | Line 399 | void EAM_FF::readParams( void ){
399  
400        sendFrcStruct( &info, mpiAtomStructType );
401  
402 +      // We have to now broadcast the Arrays
403 +      MPI_Bcast(currentAtomType->eam_rvals,
404 +                currentAtomType->eam_nr,
405 +                MPI_DOUBLE,0,MPI_COMM_WORLD);
406 +      MPI_Bcast(currentAtomType->eam_rhovals,
407 +                currentAtomType->eam_nr,
408 +                MPI_DOUBLE,0,MPI_COMM_WORLD);
409 +      MPI_Bcast(currentAtomType->eam_Frhovals,
410 +                currentAtomType->eam_nrho,
411 +                MPI_DOUBLE,0,MPI_COMM_WORLD);
412 +
413        sprintf( checkPointMsg,
414 <               "successfully sent lJ force type: \"%s\"\n",
414 >               "successfully sent EAM force type: \"%s\"\n",
415                 info.name );
416        MPIcheckPoint();
417  
# Line 333 | Line 429 | void EAM_FF::readParams( void ){
429      MPIcheckPoint();
430  
431      headAtomType = new LinkedAtomType;
432 <    recieveFrcStruct( &info, mpiAtomStructType );
433 <    
432 >    receiveFrcStruct( &info, mpiAtomStructType );
433 >
434      while( !info.last ){
435 +      
436 +      // allocate the arrays
437  
438 +      eam_rvals    = new double[info.eam_nr];
439 +      eam_rhovals  = new double[info.eam_nr];
440 +      eam_Frhovals = new double[info.eam_nrho];
441  
442 +      // We have to now broadcast the Arrays
443 +      MPI_Bcast(eam_rvals,
444 +                info.eam_nr,
445 +                MPI_DOUBLE,0,MPI_COMM_WORLD);
446 +      MPI_Bcast(eam_rhovals,
447 +                info.eam_nr,
448 +                MPI_DOUBLE,0,MPI_COMM_WORLD);
449 +      MPI_Bcast(eam_Frhovals,
450 +                info.eam_nrho,
451 +                MPI_DOUBLE,0,MPI_COMM_WORLD);
452 +      
453  
454 <      headAtomType->add( info );
454 >      headAtomType->add( info, eam_rvals, eam_rhovals, eam_Frhovals );
455        
456        MPIcheckPoint();
457  
458 <      recieveFrcStruct( &info, mpiAtomStructType );
458 >      receiveFrcStruct( &info, mpiAtomStructType );
459 >
460 >
461      }
462    }
463   #endif // is_mpi
# Line 353 | Line 467 | void EAM_FF::readParams( void ){
467    int isError;
468  
469    // dummy variables
470 <  int isLJ = 1;
470 >  int isLJ = 0;
471    int isDipole = 0;
472    int isSSD = 0;
473    int isGB = 0;
474 +  int isEAM = 1;
475 +  int isCharge = 0;
476    double dipole = 0.0;
477 +  double charge = 0.0;
478 +  double eamSigma = 0.0;
479 +  double eamEpslon = 0.0;
480    
481 <  currentAtomType = headAtomType;
481 >  currentAtomType = headAtomType->next;
482    while( currentAtomType != NULL ){
483      
484      if( currentAtomType->name[0] != '\0' ){
# Line 369 | Line 488 | void EAM_FF::readParams( void ){
488                   &isSSD,
489                   &isDipole,
490                   &isGB,
491 <                 &(currentAtomType->epslon),
492 <                 &(currentAtomType->sigma),
491 >                 &isEAM,
492 >                 &isCharge,
493 >                 &eamEpslon,
494 >                 &eamSigma,
495 >                 &charge,
496                   &dipole,
497                   &isError );
498        if( isError ){
# Line 384 | Line 506 | void EAM_FF::readParams( void ){
506      currentAtomType = currentAtomType->next;
507    }
508        
509 <  entry_plug->useLJ = 1;
509 >  entry_plug->useLJ = 0;
510 >  entry_plug->useEAM = 1;
511 >  // Walk down again and send out EAM type
512 >  currentAtomType = headAtomType->next;
513 >  while( currentAtomType != NULL ){
514 >    
515 >    if( currentAtomType->name[0] != '\0' ){
516 >      isError = 0;
517  
518 +      newEAMtype( &(currentAtomType->lattice_constant),
519 +                  &(currentAtomType->eam_nrho),
520 +                  &(currentAtomType->eam_drho),
521 +                  &(currentAtomType->eam_nr),
522 +                  &(currentAtomType->eam_dr),
523 +                  &(currentAtomType->eam_rcut),
524 +                  currentAtomType->eam_rvals,
525 +                  currentAtomType->eam_rhovals,
526 +                  currentAtomType->eam_Frhovals,
527 +                  &(currentAtomType->eam_ident),
528 +                  &isError);
529 +
530 +      if( isError ){
531 +        sprintf( painCave.errMsg,
532 +                 "Error initializing the \"%s\" atom type in fortran EAM\n",
533 +                 currentAtomType->name );
534 +        painCave.isFatal = 1;
535 +        simError();
536 +      }
537 +    }
538 +    currentAtomType = currentAtomType->next;
539 +  }
540 +
541 +
542 +
543   #ifdef IS_MPI
544    sprintf( checkPointMsg,
545 <           "LJ_FF atom structures successfully sent to fortran\n" );
545 >           "EAM_FF atom structures successfully sent to fortran\n" );
546    MPIcheckPoint();
547   #endif // is_mpi
548  
549 +
550 +
551   }
552  
553  
554   void EAM_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
555    
556    int i;
557 <
557 >  
558    // initialize the atoms
559    
404 –
405 –  Atom* thisAtom;
406 –
560    for( i=0; i<nAtoms; i++ ){
561      
562      currentAtomType = headAtomType->find( the_atoms[i]->getType() );
# Line 416 | Line 569 | void EAM_FF::initializeAtoms( int nAtoms, Atom** the_a
569      }
570      
571      the_atoms[i]->setMass( currentAtomType->mass );
419 –    the_atoms[i]->setEpslon( currentAtomType->epslon );
420 –    the_atoms[i]->setSigma( currentAtomType->sigma );
572      the_atoms[i]->setIdent( currentAtomType->ident );
422 –    the_atoms[i]->setLJ();
573  
574 <    if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
574 >    if (eamRcut < currentAtomType->eam_rcut) eamRcut = currentAtomType->eam_rcut;
575 >    
576    }
577   }
578  
579 < void LJ_FF::initializeBonds( int nBonds, Bond** BondArray,
579 > void EAM_FF::initializeBonds( int nBonds, Bond** BondArray,
580                               bond_pair* the_bonds ){
581    
582      if( nBonds ){
# Line 514 | Line 665 | void EAM_FF::fastForward( char* stopText, char* search
665  
666  
667  
668 < int EAM_NS::parseAtom( char *lineBuffer, int lineNum,  atomStruct &info ){
668 > int EAM_NS::parseAtom( char *lineBuffer, int lineNum,   atomStruct &info, char *eamPotFile ){
669  
670    char* the_token;
671    
# Line 539 | Line 690 | int EAM_NS::parseAtom( char *lineBuffer, int lineNum,
690        simError();
691      }
692          
693 <    info.epslon = atof( the_token );
694 <          
695 <    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
693 >    strcpy( eamPotFile, the_token );
694 >    return 1;
695 >  }
696 >  else return 0;
697 > }
698 >
699 > int EAM_NS::parseEAM(atomStruct &info, char *eamPotFile,
700 >                     double **eam_rvals,
701 >                     double **eam_rhovals,
702 >                     double **eam_Frhovals){
703 >  double* myEam_rvals;
704 >  double* myEam_rhovals;
705 >  double* myEam_Frhovals;
706 >
707 >  char* ffPath_env = "FORCE_PARAM_PATH";
708 >  char* ffPath;
709 >  char* the_token;
710 >  char* eam_eof_test;
711 >  FILE *eamFile;
712 >  const int BUFFERSIZE = 3000;
713 >
714 >  char temp[200];
715 >  int linenumber;
716 >  int nReadLines;
717 >  char eam_read_buffer[BUFFERSIZE];
718 >
719 >
720 >  int i,j;
721 >
722 >  linenumber = 0;
723 >
724 >  // Open eam file
725 >  eamFile = fopen( eamPotFile, "r" );
726 >  
727 >  
728 >  if( eamFile == NULL ){
729 >    
730 >      // next see if the force path enviorment variable is set
731 >    
732 >    ffPath = getenv( ffPath_env );
733 >    if( ffPath == NULL ) {
734 >      STR_DEFINE(ffPath, FRC_PATH );
735 >    }
736 >    
737 >    
738 >    strcpy( temp, ffPath );
739 >    strcat( temp, "/" );
740 >    strcat( temp, eamPotFile );
741 >    strcpy( eamPotFile, temp );
742 >    
743 >    eamFile = fopen( eamPotFile, "r" );
744 >
745 >    
746 >    
747 >    if( eamFile == NULL ){
748 >      
749 >      sprintf( painCave.errMsg,
750 >               "Error opening the EAM force parameter file: %s\n"
751 >               "Have you tried setting the FORCE_PARAM_PATH environment "
752 >               "variable?\n",
753 >               eamPotFile );
754 >      painCave.isFatal = 1;
755 >      simError();
756 >    }
757 >  }
758 >
759 >  // First line is a comment line, read and toss it....
760 >  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
761 >  linenumber++;
762 >  if(eam_eof_test == NULL){
763 >    sprintf( painCave.errMsg,
764 >             "error in reading commment in %s\n", eamPotFile);
765 >    painCave.isFatal = 1;
766 >    simError();
767 >  }
768 >
769 >
770 >
771 >  // The Second line contains atomic number, atomic mass and a lattice constant
772 >  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
773 >  linenumber++;
774 >  if(eam_eof_test == NULL){
775 >    sprintf( painCave.errMsg,
776 >             "error in reading Identifier line in %s\n", eamPotFile);
777 >    painCave.isFatal = 1;
778 >    simError();
779 >  }
780 >
781 >
782 >
783 >    
784 >  if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
785 >    sprintf( painCave.errMsg,
786 >             "Error parseing EAM ident  line in %s\n", eamPotFile );
787 >    painCave.isFatal = 1;
788 >    simError();
789 >  }
790 >  
791 >  info.eam_ident = atoi( the_token );
792 >
793 >  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
794 >    sprintf( painCave.errMsg,
795 >             "Error parseing EAM mass in %s\n", eamPotFile );
796 >    painCave.isFatal = 1;
797 >    simError();
798 >  }
799 >  info.mass = atof( the_token);
800 >
801 >  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
802 >    sprintf( painCave.errMsg,
803 >             "Error parseing EAM Lattice Constant %s\n", eamPotFile );
804 >    painCave.isFatal = 1;
805 >    simError();
806 >  }
807 >  info.lattice_constant = atof( the_token);
808 >
809 >  // Next line is nrho, drho, nr, dr and rcut
810 >  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
811 >  if(eam_eof_test == NULL){
812 >    sprintf( painCave.errMsg,
813 >             "error in reading number of points line in %s\n", eamPotFile);
814 >    painCave.isFatal = 1;
815 >    simError();
816 >  }
817 >
818 >  if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
819 >    sprintf( painCave.errMsg,
820 >             "Error parseing EAM nrho: line in %s\n", eamPotFile );
821 >    painCave.isFatal = 1;
822 >    simError();
823 >  }
824 >  
825 >  info.eam_nrho = atoi( the_token );
826 >  
827 >  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
828 >    sprintf( painCave.errMsg,
829 >             "Error parseing EAM drho in %s\n", eamPotFile );
830 >    painCave.isFatal = 1;
831 >    simError();
832 >  }
833 >  info.eam_drho = atof( the_token);
834 >
835 >  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
836 >    sprintf( painCave.errMsg,
837 >             "Error parseing EAM # r in %s\n", eamPotFile );
838 >    painCave.isFatal = 1;
839 >    simError();
840 >  }
841 >  info.eam_nr = atoi( the_token);
842 >  
843 >  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
844 >    sprintf( painCave.errMsg,
845 >             "Error parseing EAM dr in %s\n", eamPotFile );
846 >    painCave.isFatal = 1;
847 >    simError();
848 >  }
849 >  info.eam_dr = atof( the_token);
850 >
851 >  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
852 >    sprintf( painCave.errMsg,
853 >             "Error parseing EAM rcut in %s\n", eamPotFile );
854 >    painCave.isFatal = 1;
855 >    simError();
856 >  }
857 >  info.eam_rcut = atof( the_token);
858 >
859 >
860 >
861 >
862 >
863 >  // Ok now we have to allocate point arrays and read in number of points
864 >  // Index the arrays for fortran, starting at 1
865 >  myEam_Frhovals = new double[info.eam_nrho];
866 >  myEam_rvals    = new double[info.eam_nr];
867 >  myEam_rhovals  = new double[info.eam_nr];
868 >
869 >  // Parse F of rho vals.
870 >
871 >  // Assume for now that we have a complete number of lines
872 >  nReadLines = int(info.eam_nrho/5);
873 >  
874 >
875 >
876 >  for (i=0;i<nReadLines;i++){
877 >    j = i*5;
878 >
879 >    // Read next line
880 >    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
881 >    linenumber++;
882 >    if(eam_eof_test == NULL){
883 >      sprintf( painCave.errMsg,
884 >               "error in reading EAM file %s at line %d\n",
885 >               eamPotFile,linenumber);
886 >      painCave.isFatal = 1;
887 >      simError();
888 >    }
889 >    
890 >    // Parse 5 values on each line into array
891 >    // Value 1
892 >    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
893        sprintf( painCave.errMsg,
894 <               "Error parseing AtomTypes: line %d\n", lineNum );
894 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
895        painCave.isFatal = 1;
896        simError();
897      }
898 <        
899 <    info.sigma = atof( the_token );
898 >
899 >    myEam_Frhovals[j+0] = atof( the_token );
900      
901 <    return 1;
901 >    // Value 2
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+1] = atof( the_token );
910 >    
911 >    // Value 3
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+2] = atof( the_token );
920 >    
921 >    // Value 4
922 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
923 >      sprintf( painCave.errMsg,
924 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
925 >      painCave.isFatal = 1;
926 >      simError();
927 >    }
928 >
929 >    myEam_Frhovals[j+3] = atof( the_token );
930 >
931 >    // Value 5
932 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
933 >      sprintf( painCave.errMsg,
934 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
935 >      painCave.isFatal = 1;
936 >      simError();
937 >    }
938 >
939 >    myEam_Frhovals[j+4] = atof( the_token );
940 >    
941    }
942 <  else return 0;
942 >  // Parse Z of r vals
943 >  
944 >  // Assume for now that we have a complete number of lines
945 >  nReadLines = int(info.eam_nr/5);
946 >
947 >  for (i=0;i<nReadLines;i++){
948 >    j = i*5;
949 >
950 >    // Read next line
951 >    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
952 >    linenumber++;
953 >    if(eam_eof_test == NULL){
954 >      sprintf( painCave.errMsg,
955 >               "error in reading EAM file %s at line %d\n",
956 >               eamPotFile,linenumber);
957 >      painCave.isFatal = 1;
958 >      simError();
959 >    }
960 >    
961 >    // Parse 5 values on each line into array
962 >    // Value 1
963 >    if ( (the_token = strtok( eam_read_buffer, " \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+0] = atof( the_token );
971 >
972 >    // Value 2
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+1] = atof( the_token );
981 >
982 >    // Value 3
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+2] = atof( the_token );
991 >
992 >    // Value 4
993 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
994 >      sprintf( painCave.errMsg,
995 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
996 >      painCave.isFatal = 1;
997 >      simError();
998 >    }
999 >  
1000 >    myEam_rvals[j+3] = atof( the_token );
1001 >
1002 >    // Value 5
1003 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
1004 >      sprintf( painCave.errMsg,
1005 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
1006 >      painCave.isFatal = 1;
1007 >      simError();
1008 >    }
1009 >  
1010 >    myEam_rvals[j+4] = atof( the_token );
1011 >
1012 >  }
1013 >  // Parse rho of r vals
1014 >
1015 >  // Assume for now that we have a complete number of lines
1016 >
1017 >  for (i=0;i<nReadLines;i++){
1018 >    j = i*5;
1019 >
1020 >    // Read next line
1021 >    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
1022 >    linenumber++;
1023 >    if(eam_eof_test == NULL){
1024 >      sprintf( painCave.errMsg,
1025 >               "error in reading EAM file %s at line %d\n",
1026 >               eamPotFile,linenumber);
1027 >      painCave.isFatal = 1;
1028 >      simError();
1029 >    }
1030 >  
1031 >    // Parse 5 values on each line into array
1032 >    // Value 1
1033 >    if ( (the_token = strtok( eam_read_buffer, " \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+0] = atof( the_token );
1041 >
1042 >    // Value 2
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+1] = atof( the_token );
1051 >
1052 >    // Value 3
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+2] = atof( the_token );
1061 >
1062 >    // Value 4
1063 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
1064 >      sprintf( painCave.errMsg,
1065 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
1066 >      painCave.isFatal = 1;
1067 >      simError();
1068 >    }
1069 >  
1070 >    myEam_rhovals[j+3] = atof( the_token );
1071 >
1072 >    // Value 5
1073 >    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
1074 >      sprintf( painCave.errMsg,
1075 >               "Error parseing EAM nrho: line in %s\n", eamPotFile );
1076 >      painCave.isFatal = 1;
1077 >      simError();
1078 >    }
1079 >  
1080 >    myEam_rhovals[j+4] = atof( the_token );
1081 >
1082 >  }
1083 >  *eam_rvals = myEam_rvals;
1084 >  *eam_rhovals = myEam_rhovals;
1085 >  *eam_Frhovals = myEam_Frhovals;
1086 >
1087 >  fclose(eamFile);
1088 >  return 0;
1089   }

Diff Legend

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