src/UseTheForce/Shapes_FF.cpp

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include <iostream>
using namespace std;

#ifdef IS_MPI
#include <mpi.h>
#endif //is_mpi

#include "UseTheForce/ForceFields.hpp"
#include "primitives/SRI.hpp"
#include "utils/simError.h"

#include "UseTheForce/fortranWrappers.hpp"
#include "UseTheForce/DarkSide/shapes_interface.h"

#ifdef IS_MPI
#include "UseTheForce/mpiForceField.h"
#endif // is_mpi

Shapes_FF::Shapes_FF() {
  Shapes_FF("");
}

Shapes_FF::Shapes_FF(char* the_variant){

  char fileName[200];
  char* ffPath_env = "FORCE_PARAM_PATH";
  char* ffPath;
  char temp[200];

  // do the funtion wrapping
  wrapMeFF( this );

#ifdef IS_MPI    
  if( worldRank == 0 ){
#endif
    
    // generate the force file name   

    strcpy( fileName, "Shapes" );

    if (strlen(the_variant) > 0) {
      has_variant = 1;
      strcpy( variant, the_variant);
      strcat( fileName, ".");
      strcat( fileName, variant );

      sprintf( painCave.errMsg,
               "Using %s variant of Shapes force field.\n",
               variant );
      painCave.severity = OOPSE_INFO;
      painCave.isFatal = 0;
      simError();
    }
    strcat( fileName, ".frc");
    
    // attempt to open the file in the current directory first.
    
    frcFile = fopen( fileName, "r" );
    
    if( frcFile == NULL ){
      
      // next see if the force path enviorment variable is set
      
      ffPath = getenv( ffPath_env );
      if( ffPath == NULL ) {
        STR_DEFINE(ffPath, FRC_PATH );
      }
      
      
      strcpy( temp, ffPath );
      strcat( temp, "/" );
      strcat( temp, fileName );
      strcpy( fileName, temp );
      
      frcFile = fopen( fileName, "r" );
      
      if( frcFile == NULL ){
        
        sprintf( painCave.errMsg,
                 "Error opening the force field parameter file:\n"
                 "\t%s\n"
                 "\tHave you tried setting the FORCE_PARAM_PATH environment "
                 "variable?\n",
                 fileName );
        painCave.severity = OOPSE_ERROR;
        painCave.isFatal = 1;
        simError();
      }
    }
    
    
#ifdef IS_MPI
  }
  
  sprintf( checkPointMsg, "Shapes_FF file opened sucessfully." );
  MPIcheckPoint();
  
#endif // is_mpi
}


Shapes_FF::~Shapes_FF(){

#ifdef IS_MPI
  if( worldRank == 0 ){
#endif // is_mpi
    
    fclose( frcFile );
    
#ifdef IS_MPI
  }
#endif // is_mpi
}


void Shapes_FF::calcRcut( void ){
  
#ifdef IS_MPI
  double tempShapesRcut = shapesRcut;
  MPI_Allreduce( &tempShapesRcut, &shapesRcut, 1, MPI_DOUBLE, MPI_MAX,
                 MPI_COMM_WORLD);
#endif  //is_mpi
  entry_plug->setDefaultRcut(shapesRcut);
}


void Shapes_FF::initForceField(int ljMixRule){
  initFortran(ljMixRule, 0);
}


void Shapes_FF::readParams( void ){

  char shapePotFile[1000];
  
#ifdef IS_MPI
  if( worldRank == 0 ){
#endif
    
    // read in the atom types.

    fastForward( "AtomTypes", "eam atom readParams" );

    // we are now at the AtomTypes section.
    
    eof_test = fgets( readLine, sizeof(readLine), frcFile );
    lineNum++;
    
    
    // read a line, and start parseing out the atom types 

    if( eof_test == NULL ){
      sprintf( painCave.errMsg, 
               "Error in reading Atoms from force file at line %d.\n",
               lineNum );
      painCave.isFatal = 1;
      simError();
    }
    
    identNum = 1;
    // stop reading at end of file, or at next section

    while( readLine[0] != '#' && eof_test != NULL ){

      // toss comment lines
      if( readLine[0] != '!' ){
        
        // the parser returns 0 if the line was blank
        if( parseAtom( readLine, lineNum, info, eamPotFile ) ){
          parseEAM(info,eamPotFile, &eam_rvals,
                   &eam_rhovals, &eam_Frhovals);
          info.ident = identNum;
          headAtomType->add( info, eam_rvals,
                             eam_rhovals,eam_Frhovals );
          identNum++;
        }
      }
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;
    }
    
    
#ifdef IS_MPI
  
    
    // send out the linked list to all the other processes

    sprintf( checkPointMsg,
             "Shapes_FF atom structures read successfully." );
    MPIcheckPoint();

    currentAtomType = headAtomType->next; //skip the first element who is a place holder.
    while( currentAtomType != NULL ){
      currentAtomType->duplicate( info );

 
      sendFrcStruct( &info, mpiAtomStructType );

      // We have to now broadcast the Arrays
      MPI_Bcast(currentAtomType->eam_rvals,
                currentAtomType->eam_nr,
                MPI_DOUBLE,0,MPI_COMM_WORLD);
      MPI_Bcast(currentAtomType->eam_rhovals,
                currentAtomType->eam_nr,
                MPI_DOUBLE,0,MPI_COMM_WORLD);
      MPI_Bcast(currentAtomType->eam_Frhovals,
                currentAtomType->eam_nrho,
                MPI_DOUBLE,0,MPI_COMM_WORLD);

      sprintf( checkPointMsg,
               "successfully sent EAM force type: \"%s\"\n",
               info.name );
      MPIcheckPoint();

      currentAtomType = currentAtomType->next;
    }
    info.last = 1;
    sendFrcStruct( &info, mpiAtomStructType );
    
  }

  else{
    
    // listen for node 0 to send out the force params
    
    MPIcheckPoint();

    headAtomType = new LinkedAtomType;
    receiveFrcStruct( &info, mpiAtomStructType );

    while( !info.last ){
      
      // allocate the arrays

      eam_rvals    = new double[info.eam_nr];
      eam_rhovals  = new double[info.eam_nr];
      eam_Frhovals = new double[info.eam_nrho];

      // We have to now broadcast the Arrays
      MPI_Bcast(eam_rvals,
                info.eam_nr,
                MPI_DOUBLE,0,MPI_COMM_WORLD);
      MPI_Bcast(eam_rhovals,
                info.eam_nr,
                MPI_DOUBLE,0,MPI_COMM_WORLD);
      MPI_Bcast(eam_Frhovals,
                info.eam_nrho,
                MPI_DOUBLE,0,MPI_COMM_WORLD);
      

      headAtomType->add( info, eam_rvals, eam_rhovals, eam_Frhovals );
      
      MPIcheckPoint();

      receiveFrcStruct( &info, mpiAtomStructType );


    }
  }
#endif // is_mpi

  // call new A_types in fortran
  
  int isError;

  // dummy variables
  int isLJ = 0;
  int isDipole = 0;
  int isSSD = 0;
  int isGB = 0;
  int isEAM = 1;
  int isCharge = 0;
  double dipole = 0.0;
  double charge = 0.0;
  double eamSigma = 0.0;
  double eamEpslon = 0.0;
  
  currentAtomType = headAtomType->next;
  while( currentAtomType != NULL ){
    
    if( currentAtomType->name[0] != '\0' ){
      isError = 0;
      makeAtype( &(currentAtomType->ident),
                 &isLJ,
                 &isSSD,
                 &isDipole,
                 &isGB,
                 &isEAM,
                 &isCharge,
                 &eamEpslon,
                 &eamSigma,
                 &charge,
                 &dipole,
                 &isError );
      if( isError ){
        sprintf( painCave.errMsg,
                 "Error initializing the \"%s\" atom type in fortran\n",
                 currentAtomType->name );
        painCave.isFatal = 1;
        simError();
      }
    }
    currentAtomType = currentAtomType->next;
  }
      
  entry_plug->useLJ = 0;
  entry_plug->useEAM = 1;
  // Walk down again and send out EAM type
  currentAtomType = headAtomType->next;
  while( currentAtomType != NULL ){
    
    if( currentAtomType->name[0] != '\0' ){
      isError = 0;

      newEAMtype( &(currentAtomType->lattice_constant),
                  &(currentAtomType->eam_nrho),
                  &(currentAtomType->eam_drho),
                  &(currentAtomType->eam_nr),
                  &(currentAtomType->eam_dr),
                  &(currentAtomType->eam_rcut),
                  currentAtomType->eam_rvals,
                  currentAtomType->eam_rhovals,
                  currentAtomType->eam_Frhovals,
                  &(currentAtomType->eam_ident),
                  &isError);

      if( isError ){
        sprintf( painCave.errMsg,
                 "Error initializing the \"%s\" atom type in fortran EAM\n",
                 currentAtomType->name );
        painCave.isFatal = 1;
        simError();
      }
    }
    currentAtomType = currentAtomType->next;
  }


#ifdef IS_MPI
  sprintf( checkPointMsg,
           "Shapes_FF atom structures successfully sent to fortran\n" );
  MPIcheckPoint();
#endif // is_mpi


}


void Shapes_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
  
  int i;
  
  // initialize the atoms
  
  for( i=0; i<nAtoms; i++ ){
    
    currentAtomType = headAtomType->find( the_atoms[i]->getType() );
    if( currentAtomType == NULL ){
      sprintf( painCave.errMsg, 
               "AtomType error, %s not found in force file.\n",
               the_atoms[i]->getType() );
      painCave.isFatal = 1;
      simError();
    }
    
    the_atoms[i]->setMass( currentAtomType->mass );
    the_atoms[i]->setIdent( currentAtomType->ident );

    if (eamRcut < currentAtomType->eam_rcut) eamRcut = currentAtomType->eam_rcut;
    
  }
}

void Shapes_FF::initializeBonds( int nBonds, Bond** BondArray,
                             bond_pair* the_bonds ){
  
    if( nBonds ){
      sprintf( painCave.errMsg,
               "LJ_FF does not support bonds.\n" );
      painCave.isFatal = 1;
      simError();
    }
}

void Shapes_FF::initializeBends( int nBends, Bend** bendArray,
                             bend_set* the_bends ){

    if( nBends ){
      sprintf( painCave.errMsg,
               "LJ_FF does not support bends.\n" );
      painCave.isFatal = 1;
      simError();
    }
}

void Shapes_FF::initializeTorsions( int nTorsions, Torsion** torsionArray,
                                torsion_set* the_torsions ){

    if( nTorsions ){
      sprintf( painCave.errMsg,
               "LJ_FF does not support torsions.\n" );
      painCave.isFatal = 1;
      simError();
    }
}

void Shapes_FF::fastForward( char* stopText, char* searchOwner ){

  int foundText = 0;
  char* the_token;

  rewind( frcFile );
  lineNum = 0;

  eof_test = fgets( readLine, sizeof(readLine), frcFile );
  lineNum++;
  if( eof_test == NULL ){
    sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
             " file is empty.\n",
             searchOwner );
    painCave.isFatal = 1;
    simError();
  }
  
  
  while( !foundText ){
    while( eof_test != NULL && readLine[0] != '#' ){
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;
    }
    if( eof_test == NULL ){
      sprintf( painCave.errMsg,
               "Error fast forwarding force file for %s at "
               "line %d: file ended unexpectedly.\n",
               searchOwner,
               lineNum );
      painCave.isFatal = 1;
      simError();
    }
    
    the_token = strtok( readLine, " ,;\t#\n" );
    foundText = !strcmp( stopText, the_token );
    
    if( !foundText ){
      eof_test = fgets( readLine, sizeof(readLine), frcFile );
      lineNum++;
      
      if( eof_test == NULL ){
        sprintf( painCave.errMsg, 
                 "Error fast forwarding force file for %s at "
                 "line %d: file ended unexpectedly.\n",
                 searchOwner,
                 lineNum );
        painCave.isFatal = 1;
        simError();
      } 
    }
  }  
}


int EAM_NS::parseAtom( char *lineBuffer, int lineNum,   atomStruct &info, char *eamPotFile ){

  char* the_token;
  
  the_token = strtok( lineBuffer, " \n\t,;" );
  if( the_token != NULL ){
    
    strcpy( info.name, the_token );
          
    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
      sprintf( painCave.errMsg, 
               "Error parseing AtomTypes: line %d\n", lineNum );
      painCave.isFatal = 1;
      simError();
    }
    
    info.mass = atof( the_token );
    
    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
      sprintf( painCave.errMsg, 
               "Error parseing AtomTypes: line %d\n", lineNum );
      painCave.isFatal = 1;
      simError();
    }
        
    strcpy( eamPotFile, the_token );
    return 1;
  }
  else return 0;
}

int EAM_NS::parseEAM(atomStruct &info, char *eamPotFile, 
                     double **eam_rvals, 
                     double **eam_rhovals,
                     double **eam_Frhovals){
  double* myEam_rvals;
  double* myEam_rhovals;
  double* myEam_Frhovals;

  char* ffPath_env = "FORCE_PARAM_PATH";
  char* ffPath;
  char* the_token;
  char* eam_eof_test;
  FILE *eamFile;
  const int BUFFERSIZE = 3000;

  char temp[200];
  int linenumber;
  int nReadLines;
  char eam_read_buffer[BUFFERSIZE];


  int i,j;

  linenumber = 0;

  // Open eam file
  eamFile = fopen( eamPotFile, "r" );
  
  
  if( eamFile == NULL ){
    
      // next see if the force path enviorment variable is set
    
    ffPath = getenv( ffPath_env );
    if( ffPath == NULL ) {
      STR_DEFINE(ffPath, FRC_PATH );
    }
    
    
    strcpy( temp, ffPath );
    strcat( temp, "/" );
    strcat( temp, eamPotFile );
    strcpy( eamPotFile, temp );
    
    eamFile = fopen( eamPotFile, "r" );

    
    if( eamFile == NULL ){
      
      sprintf( painCave.errMsg,
               "Error opening the EAM force parameter file: %s\n"
               "Have you tried setting the FORCE_PARAM_PATH environment "
               "variable?\n",
               eamPotFile );
      painCave.isFatal = 1;
      simError();
    }
  }

  // First line is a comment line, read and toss it....
  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
  linenumber++;
  if(eam_eof_test == NULL){
    sprintf( painCave.errMsg,
             "error in reading commment in %s\n", eamPotFile);
    painCave.isFatal = 1;
    simError();
  }


  // The Second line contains atomic number, atomic mass and a lattice constant
  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
  linenumber++;
  if(eam_eof_test == NULL){
    sprintf( painCave.errMsg,
             "error in reading Identifier line in %s\n", eamPotFile);
    painCave.isFatal = 1;
    simError();
  }


  if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
    sprintf( painCave.errMsg, 
             "Error parsing EAM ident  line in %s\n", eamPotFile );
    painCave.isFatal = 1;
    simError();
  }
  
  info.eam_ident = atoi( the_token );
 
  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
    sprintf( painCave.errMsg, 
             "Error parsing EAM mass in %s\n", eamPotFile );
    painCave.isFatal = 1;
    simError();
  }
  info.mass = atof( the_token);
 
  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
    sprintf( painCave.errMsg, 
             "Error parsing EAM Lattice Constant %s\n", eamPotFile );
    painCave.isFatal = 1;
    simError();
  }
  info.lattice_constant = atof( the_token);
 
  // Next line is nrho, drho, nr, dr and rcut
  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
  if(eam_eof_test == NULL){
    sprintf( painCave.errMsg,
             "error in reading number of points line in %s\n", eamPotFile);
    painCave.isFatal = 1;
    simError();
  }

  if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
    sprintf( painCave.errMsg, 
             "Error parseing EAM nrho: line in %s\n", eamPotFile );
    painCave.isFatal = 1;
    simError();
  }
  
  info.eam_nrho = atoi( the_token );
  
  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
    sprintf( painCave.errMsg, 
             "Error parsing EAM drho in %s\n", eamPotFile );
    painCave.isFatal = 1;
    simError();
  }
  info.eam_drho = atof( the_token);

  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
    sprintf( painCave.errMsg, 
             "Error parsing EAM # r in %s\n", eamPotFile );
    painCave.isFatal = 1;
    simError();
  }
  info.eam_nr = atoi( the_token);
  
  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
    sprintf( painCave.errMsg, 
             "Error parsing EAM dr in %s\n", eamPotFile );
    painCave.isFatal = 1;
    simError();
  }
  info.eam_dr = atof( the_token);

  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
    sprintf( painCave.errMsg, 
             "Error parsing EAM rcut in %s\n", eamPotFile );
    painCave.isFatal = 1;
    simError();
  }
  info.eam_rcut = atof( the_token);


  // Ok now we have to allocate point arrays and read in number of points
  // Index the arrays for fortran, starting at 1
  myEam_Frhovals = new double[info.eam_nrho];
  myEam_rvals    = new double[info.eam_nr];
  myEam_rhovals  = new double[info.eam_nr];

  // Parse F of rho vals.

  // Assume for now that we have a complete number of lines
  nReadLines = int(info.eam_nrho/5);
  

  for (i=0;i<nReadLines;i++){
    j = i*5;

    // Read next line
    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
    linenumber++;
    if(eam_eof_test == NULL){
      sprintf( painCave.errMsg,
               "error in reading EAM file %s at line %d\n",
               eamPotFile,linenumber);
      painCave.isFatal = 1;
      simError();
    }
    
    // Parse 5 values on each line into array
    // Value 1
    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
      sprintf( painCave.errMsg, 
               "Error parsing EAM nrho: line in %s\n", eamPotFile );
      painCave.isFatal = 1;
      simError();
    }

    myEam_Frhovals[j+0] = atof( the_token );
    
    // Value 2
    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
      sprintf( painCave.errMsg, 
               "Error parsing EAM nrho: line in %s\n", eamPotFile );
      painCave.isFatal = 1;
      simError();
    }

    myEam_Frhovals[j+1] = atof( the_token );
    
    // Value 3
    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
      sprintf( painCave.errMsg, 
               "Error parsing EAM nrho: line in %s\n", eamPotFile );
      painCave.isFatal = 1;
      simError();
    }

    myEam_Frhovals[j+2] = atof( the_token );
    
    // Value 4
    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
      sprintf( painCave.errMsg, 
               "Error parsing EAM nrho: line in %s\n", eamPotFile );
      painCave.isFatal = 1;
      simError();
    }

    myEam_Frhovals[j+3] = atof( the_token );

    // Value 5
    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
      sprintf( painCave.errMsg, 
               "Error parsing EAM nrho: line in %s\n", eamPotFile );
      painCave.isFatal = 1;
      simError();
    }

    myEam_Frhovals[j+4] = atof( the_token );
    
  }
  // Parse Z of r vals
  
  // Assume for now that we have a complete number of lines
  nReadLines = int(info.eam_nr/5);

  for (i=0;i<nReadLines;i++){
    j = i*5;

    // Read next line
    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
    linenumber++;
    if(eam_eof_test == NULL){
      sprintf( painCave.errMsg,
               "error in reading EAM file %s at line %d\n",
               eamPotFile,linenumber);
      painCave.isFatal = 1;
      simError();
    }
    
    // Parse 5 values on each line into array
    // Value 1
    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
      sprintf( painCave.errMsg, 
               "Error parsing EAM nrho: line in %s\n", eamPotFile );
      painCave.isFatal = 1;
      simError();
    }
    
    myEam_rvals[j+0] = atof( the_token );

    // Value 2
    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
      sprintf( painCave.errMsg, 
               "Error parsing EAM nrho: line in %s\n", eamPotFile );
      painCave.isFatal = 1;
      simError();
    }
  
    myEam_rvals[j+1] = atof( the_token );

    // Value 3
    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
      sprintf( painCave.errMsg, 
               "Error parsing EAM nrho: line in %s\n", eamPotFile );
      painCave.isFatal = 1;
      simError();
    }
  
    myEam_rvals[j+2] = atof( the_token );

    // Value 4
    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
      sprintf( painCave.errMsg, 
               "Error parsing EAM nrho: line in %s\n", eamPotFile );
      painCave.isFatal = 1;
      simError();
    }
  
    myEam_rvals[j+3] = atof( the_token );

    // Value 5
    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
      sprintf( painCave.errMsg, 
               "Error parsing EAM nrho: line in %s\n", eamPotFile );
      painCave.isFatal = 1;
      simError();
    }
  
    myEam_rvals[j+4] = atof( the_token );

  }
  // Parse rho of r vals

  // Assume for now that we have a complete number of lines

  for (i=0;i<nReadLines;i++){
    j = i*5;

    // Read next line
    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
    linenumber++;
    if(eam_eof_test == NULL){
      sprintf( painCave.errMsg,
               "error in reading EAM file %s at line %d\n",
               eamPotFile,linenumber);
      painCave.isFatal = 1;
      simError();
    }
  
    // Parse 5 values on each line into array
    // Value 1
    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
      sprintf( painCave.errMsg, 
               "Error parsing EAM nrho: line in %s\n", eamPotFile );
      painCave.isFatal = 1;
      simError();
    }
  
    myEam_rhovals[j+0] = atof( the_token );

    // Value 2
    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
      sprintf( painCave.errMsg, 
               "Error parsing EAM nrho: line in %s\n", eamPotFile );
      painCave.isFatal = 1;
      simError();
    }
  
    myEam_rhovals[j+1] = atof( the_token );

    // Value 3
    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
      sprintf( painCave.errMsg, 
               "Error parsing EAM nrho: line in %s\n", eamPotFile );
      painCave.isFatal = 1;
      simError();
    }
  
    myEam_rhovals[j+2] = atof( the_token );

    // Value 4
    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
      sprintf( painCave.errMsg, 
               "Error parsing EAM nrho: line in %s\n", eamPotFile );
      painCave.isFatal = 1;
      simError();
    }
  
    myEam_rhovals[j+3] = atof( the_token );

    // Value 5
    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
      sprintf( painCave.errMsg, 
               "Error parsing EAM nrho: line in %s\n", eamPotFile );
      painCave.isFatal = 1;
      simError();
    }
  
    myEam_rhovals[j+4] = atof( the_token );
 
  }
  *eam_rvals = myEam_rvals;
  *eam_rhovals = myEam_rhovals;
  *eam_Frhovals = myEam_Frhovals;

  fclose(eamFile);
  return 0;
}
Revision:	1613
Committed:	Wed Oct 20 04:54:23 2004 UTC (19 years, 8 months ago) by gezelter
File size:	21116 byte(s)
Log Message:	fixed include
#	Content
1	#include <stdlib.h>
2	#include <stdio.h>
3	#include <string.h>
4
5	#include <iostream>
6	using namespace std;
7
8	#ifdef IS_MPI
9	#include <mpi.h>
10	#endif //is_mpi
11
12	#include "UseTheForce/ForceFields.hpp"
13	#include "primitives/SRI.hpp"
14	#include "utils/simError.h"
15
16	#include "UseTheForce/fortranWrappers.hpp"
17	#include "UseTheForce/DarkSide/shapes_interface.h"
18
19	#ifdef IS_MPI
20	#include "UseTheForce/mpiForceField.h"
21	#endif // is_mpi
22
23	Shapes_FF::Shapes_FF() {
24	Shapes_FF("");
25	}
26
27	Shapes_FF::Shapes_FF(char* the_variant){
28
29	char fileName[200];
30	char* ffPath_env = "FORCE_PARAM_PATH";
31	char* ffPath;
32	char temp[200];
33
34	// do the funtion wrapping
35	wrapMeFF( this );
36
37	#ifdef IS_MPI
38	if( worldRank == 0 ){
39	#endif
40
41	// generate the force file name
42
43	strcpy( fileName, "Shapes" );
44
45	if (strlen(the_variant) > 0) {
46	has_variant = 1;
47	strcpy( variant, the_variant);
48	strcat( fileName, ".");
49	strcat( fileName, variant );
50
51	sprintf( painCave.errMsg,
52	"Using %s variant of Shapes force field.\n",
53	variant );
54	painCave.severity = OOPSE_INFO;
55	painCave.isFatal = 0;
56	simError();
57	}
58	strcat( fileName, ".frc");
59
60	// attempt to open the file in the current directory first.
61
62	frcFile = fopen( fileName, "r" );
63
64	if( frcFile == NULL ){
65
66	// next see if the force path enviorment variable is set
67
68	ffPath = getenv( ffPath_env );
69	if( ffPath == NULL ) {
70	STR_DEFINE(ffPath, FRC_PATH );
71	}
72
73
74	strcpy( temp, ffPath );
75	strcat( temp, "/" );
76	strcat( temp, fileName );
77	strcpy( fileName, temp );
78
79	frcFile = fopen( fileName, "r" );
80
81	if( frcFile == NULL ){
82
83	sprintf( painCave.errMsg,
84	"Error opening the force field parameter file:\n"
85	"\t%s\n"
86	"\tHave you tried setting the FORCE_PARAM_PATH environment "
87	"variable?\n",
88	fileName );
89	painCave.severity = OOPSE_ERROR;
90	painCave.isFatal = 1;
91	simError();
92	}
93	}
94
95
96	#ifdef IS_MPI
97	}
98
99	sprintf( checkPointMsg, "Shapes_FF file opened sucessfully." );
100	MPIcheckPoint();
101
102	#endif // is_mpi
103	}
104
105
106	Shapes_FF::~Shapes_FF(){
107
108	#ifdef IS_MPI
109	if( worldRank == 0 ){
110	#endif // is_mpi
111
112	fclose( frcFile );
113
114	#ifdef IS_MPI
115	}
116	#endif // is_mpi
117	}
118
119
120	void Shapes_FF::calcRcut( void ){
121
122	#ifdef IS_MPI
123	double tempShapesRcut = shapesRcut;
124	MPI_Allreduce( &tempShapesRcut, &shapesRcut, 1, MPI_DOUBLE, MPI_MAX,
125	MPI_COMM_WORLD);
126	#endif //is_mpi
127	entry_plug->setDefaultRcut(shapesRcut);
128	}
129
130
131	void Shapes_FF::initForceField(int ljMixRule){
132	initFortran(ljMixRule, 0);
133	}
134
135
136	void Shapes_FF::readParams( void ){
137
138	char shapePotFile[1000];
139
140	#ifdef IS_MPI
141	if( worldRank == 0 ){
142	#endif
143
144	// read in the atom types.
145
146	fastForward( "AtomTypes", "eam atom readParams" );
147
148	// we are now at the AtomTypes section.
149
150	eof_test = fgets( readLine, sizeof(readLine), frcFile );
151	lineNum++;
152
153
154	// read a line, and start parseing out the atom types
155
156	if( eof_test == NULL ){
157	sprintf( painCave.errMsg,
158	"Error in reading Atoms from force file at line %d.\n",
159	lineNum );
160	painCave.isFatal = 1;
161	simError();
162	}
163
164	identNum = 1;
165	// stop reading at end of file, or at next section
166
167	while( readLine[0] != '#' && eof_test != NULL ){
168
169	// toss comment lines
170	if( readLine[0] != '!' ){
171
172	// the parser returns 0 if the line was blank
173	if( parseAtom( readLine, lineNum, info, eamPotFile ) ){
174	parseEAM(info,eamPotFile, &eam_rvals,
175	&eam_rhovals, &eam_Frhovals);
176	info.ident = identNum;
177	headAtomType->add( info, eam_rvals,
178	eam_rhovals,eam_Frhovals );
179	identNum++;
180	}
181	}
182	eof_test = fgets( readLine, sizeof(readLine), frcFile );
183	lineNum++;
184	}
185
186
187
188	#ifdef IS_MPI
189
190
191	// send out the linked list to all the other processes
192
193	sprintf( checkPointMsg,
194	"Shapes_FF atom structures read successfully." );
195	MPIcheckPoint();
196
197	currentAtomType = headAtomType->next; //skip the first element who is a place holder.
198	while( currentAtomType != NULL ){
199	currentAtomType->duplicate( info );
200
201
202
203	sendFrcStruct( &info, mpiAtomStructType );
204
205	// We have to now broadcast the Arrays
206	MPI_Bcast(currentAtomType->eam_rvals,
207	currentAtomType->eam_nr,
208	MPI_DOUBLE,0,MPI_COMM_WORLD);
209	MPI_Bcast(currentAtomType->eam_rhovals,
210	currentAtomType->eam_nr,
211	MPI_DOUBLE,0,MPI_COMM_WORLD);
212	MPI_Bcast(currentAtomType->eam_Frhovals,
213	currentAtomType->eam_nrho,
214	MPI_DOUBLE,0,MPI_COMM_WORLD);
215
216	sprintf( checkPointMsg,
217	"successfully sent EAM force type: \"%s\"\n",
218	info.name );
219	MPIcheckPoint();
220
221	currentAtomType = currentAtomType->next;
222	}
223	info.last = 1;
224	sendFrcStruct( &info, mpiAtomStructType );
225
226	}
227
228	else{
229
230	// listen for node 0 to send out the force params
231
232	MPIcheckPoint();
233
234	headAtomType = new LinkedAtomType;
235	receiveFrcStruct( &info, mpiAtomStructType );
236
237	while( !info.last ){
238
239	// allocate the arrays
240
241	eam_rvals = new double[info.eam_nr];
242	eam_rhovals = new double[info.eam_nr];
243	eam_Frhovals = new double[info.eam_nrho];
244
245	// We have to now broadcast the Arrays
246	MPI_Bcast(eam_rvals,
247	info.eam_nr,
248	MPI_DOUBLE,0,MPI_COMM_WORLD);
249	MPI_Bcast(eam_rhovals,
250	info.eam_nr,
251	MPI_DOUBLE,0,MPI_COMM_WORLD);
252	MPI_Bcast(eam_Frhovals,
253	info.eam_nrho,
254	MPI_DOUBLE,0,MPI_COMM_WORLD);
255
256
257	headAtomType->add( info, eam_rvals, eam_rhovals, eam_Frhovals );
258
259	MPIcheckPoint();
260
261	receiveFrcStruct( &info, mpiAtomStructType );
262
263
264	}
265	}
266	#endif // is_mpi
267
268	// call new A_types in fortran
269
270	int isError;
271
272	// dummy variables
273	int isLJ = 0;
274	int isDipole = 0;
275	int isSSD = 0;
276	int isGB = 0;
277	int isEAM = 1;
278	int isCharge = 0;
279	double dipole = 0.0;
280	double charge = 0.0;
281	double eamSigma = 0.0;
282	double eamEpslon = 0.0;
283
284	currentAtomType = headAtomType->next;
285	while( currentAtomType != NULL ){
286
287	if( currentAtomType->name[0] != '\0' ){
288	isError = 0;
289	makeAtype( &(currentAtomType->ident),
290	&isLJ,
291	&isSSD,
292	&isDipole,
293	&isGB,
294	&isEAM,
295	&isCharge,
296	&eamEpslon,
297	&eamSigma,
298	&charge,
299	&dipole,
300	&isError );
301	if( isError ){
302	sprintf( painCave.errMsg,
303	"Error initializing the \"%s\" atom type in fortran\n",
304	currentAtomType->name );
305	painCave.isFatal = 1;
306	simError();
307	}
308	}
309	currentAtomType = currentAtomType->next;
310	}
311
312	entry_plug->useLJ = 0;
313	entry_plug->useEAM = 1;
314	// Walk down again and send out EAM type
315	currentAtomType = headAtomType->next;
316	while( currentAtomType != NULL ){
317
318	if( currentAtomType->name[0] != '\0' ){
319	isError = 0;
320
321	newEAMtype( &(currentAtomType->lattice_constant),
322	&(currentAtomType->eam_nrho),
323	&(currentAtomType->eam_drho),
324	&(currentAtomType->eam_nr),
325	&(currentAtomType->eam_dr),
326	&(currentAtomType->eam_rcut),
327	currentAtomType->eam_rvals,
328	currentAtomType->eam_rhovals,
329	currentAtomType->eam_Frhovals,
330	&(currentAtomType->eam_ident),
331	&isError);
332
333	if( isError ){
334	sprintf( painCave.errMsg,
335	"Error initializing the \"%s\" atom type in fortran EAM\n",
336	currentAtomType->name );
337	painCave.isFatal = 1;
338	simError();
339	}
340	}
341	currentAtomType = currentAtomType->next;
342	}
343
344
345
346	#ifdef IS_MPI
347	sprintf( checkPointMsg,
348	"Shapes_FF atom structures successfully sent to fortran\n" );
349	MPIcheckPoint();
350	#endif // is_mpi
351
352
353
354	}
355
356
357	void Shapes_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
358
359	int i;
360
361	// initialize the atoms
362
363	for( i=0; i<nAtoms; i++ ){
364
365	currentAtomType = headAtomType->find( the_atoms[i]->getType() );
366	if( currentAtomType == NULL ){
367	sprintf( painCave.errMsg,
368	"AtomType error, %s not found in force file.\n",
369	the_atoms[i]->getType() );
370	painCave.isFatal = 1;
371	simError();
372	}
373
374	the_atoms[i]->setMass( currentAtomType->mass );
375	the_atoms[i]->setIdent( currentAtomType->ident );
376
377	if (eamRcut < currentAtomType->eam_rcut) eamRcut = currentAtomType->eam_rcut;
378
379	}
380	}
381
382	void Shapes_FF::initializeBonds( int nBonds, Bond** BondArray,
383	bond_pair* the_bonds ){
384
385	if( nBonds ){
386	sprintf( painCave.errMsg,
387	"LJ_FF does not support bonds.\n" );
388	painCave.isFatal = 1;
389	simError();
390	}
391	}
392
393	void Shapes_FF::initializeBends( int nBends, Bend** bendArray,
394	bend_set* the_bends ){
395
396	if( nBends ){
397	sprintf( painCave.errMsg,
398	"LJ_FF does not support bends.\n" );
399	painCave.isFatal = 1;
400	simError();
401	}
402	}
403
404	void Shapes_FF::initializeTorsions( int nTorsions, Torsion** torsionArray,
405	torsion_set* the_torsions ){
406
407	if( nTorsions ){
408	sprintf( painCave.errMsg,
409	"LJ_FF does not support torsions.\n" );
410	painCave.isFatal = 1;
411	simError();
412	}
413	}
414
415	void Shapes_FF::fastForward( char* stopText, char* searchOwner ){
416
417	int foundText = 0;
418	char* the_token;
419
420	rewind( frcFile );
421	lineNum = 0;
422
423	eof_test = fgets( readLine, sizeof(readLine), frcFile );
424	lineNum++;
425	if( eof_test == NULL ){
426	sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
427	" file is empty.\n",
428	searchOwner );
429	painCave.isFatal = 1;
430	simError();
431	}
432
433
434	while( !foundText ){
435	while( eof_test != NULL && readLine[0] != '#' ){
436	eof_test = fgets( readLine, sizeof(readLine), frcFile );
437	lineNum++;
438	}
439	if( eof_test == NULL ){
440	sprintf( painCave.errMsg,
441	"Error fast forwarding force file for %s at "
442	"line %d: file ended unexpectedly.\n",
443	searchOwner,
444	lineNum );
445	painCave.isFatal = 1;
446	simError();
447	}
448
449	the_token = strtok( readLine, " ,;\t#\n" );
450	foundText = !strcmp( stopText, the_token );
451
452	if( !foundText ){
453	eof_test = fgets( readLine, sizeof(readLine), frcFile );
454	lineNum++;
455
456	if( eof_test == NULL ){
457	sprintf( painCave.errMsg,
458	"Error fast forwarding force file for %s at "
459	"line %d: file ended unexpectedly.\n",
460	searchOwner,
461	lineNum );
462	painCave.isFatal = 1;
463	simError();
464	}
465	}
466	}
467	}
468
469
470
471	int EAM_NS::parseAtom( char lineBuffer, int lineNum, atomStruct &info, char eamPotFile ){
472
473	char* the_token;
474
475	the_token = strtok( lineBuffer, " \n\t,;" );
476	if( the_token != NULL ){
477
478	strcpy( info.name, the_token );
479
480	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
481	sprintf( painCave.errMsg,
482	"Error parseing AtomTypes: line %d\n", lineNum );
483	painCave.isFatal = 1;
484	simError();
485	}
486
487	info.mass = atof( the_token );
488
489	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
490	sprintf( painCave.errMsg,
491	"Error parseing AtomTypes: line %d\n", lineNum );
492	painCave.isFatal = 1;
493	simError();
494	}
495
496	strcpy( eamPotFile, the_token );
497	return 1;
498	}
499	else return 0;
500	}
501
502	int EAM_NS::parseEAM(atomStruct &info, char *eamPotFile,
503	double **eam_rvals,
504	double **eam_rhovals,
505	double **eam_Frhovals){
506	double* myEam_rvals;
507	double* myEam_rhovals;
508	double* myEam_Frhovals;
509
510	char* ffPath_env = "FORCE_PARAM_PATH";
511	char* ffPath;
512	char* the_token;
513	char* eam_eof_test;
514	FILE *eamFile;
515	const int BUFFERSIZE = 3000;
516
517	char temp[200];
518	int linenumber;
519	int nReadLines;
520	char eam_read_buffer[BUFFERSIZE];
521
522
523	int i,j;
524
525	linenumber = 0;
526
527	// Open eam file
528	eamFile = fopen( eamPotFile, "r" );
529
530
531	if( eamFile == NULL ){
532
533	// next see if the force path enviorment variable is set
534
535	ffPath = getenv( ffPath_env );
536	if( ffPath == NULL ) {
537	STR_DEFINE(ffPath, FRC_PATH );
538	}
539
540
541	strcpy( temp, ffPath );
542	strcat( temp, "/" );
543	strcat( temp, eamPotFile );
544	strcpy( eamPotFile, temp );
545
546	eamFile = fopen( eamPotFile, "r" );
547
548
549
550	if( eamFile == NULL ){
551
552	sprintf( painCave.errMsg,
553	"Error opening the EAM force parameter file: %s\n"
554	"Have you tried setting the FORCE_PARAM_PATH environment "
555	"variable?\n",
556	eamPotFile );
557	painCave.isFatal = 1;
558	simError();
559	}
560	}
561
562	// First line is a comment line, read and toss it....
563	eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
564	linenumber++;
565	if(eam_eof_test == NULL){
566	sprintf( painCave.errMsg,
567	"error in reading commment in %s\n", eamPotFile);
568	painCave.isFatal = 1;
569	simError();
570	}
571
572
573
574	// The Second line contains atomic number, atomic mass and a lattice constant
575	eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
576	linenumber++;
577	if(eam_eof_test == NULL){
578	sprintf( painCave.errMsg,
579	"error in reading Identifier line in %s\n", eamPotFile);
580	painCave.isFatal = 1;
581	simError();
582	}
583
584
585
586
587	if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
588	sprintf( painCave.errMsg,
589	"Error parsing EAM ident line in %s\n", eamPotFile );
590	painCave.isFatal = 1;
591	simError();
592	}
593
594	info.eam_ident = atoi( the_token );
595
596	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
597	sprintf( painCave.errMsg,
598	"Error parsing EAM mass in %s\n", eamPotFile );
599	painCave.isFatal = 1;
600	simError();
601	}
602	info.mass = atof( the_token);
603
604	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
605	sprintf( painCave.errMsg,
606	"Error parsing EAM Lattice Constant %s\n", eamPotFile );
607	painCave.isFatal = 1;
608	simError();
609	}
610	info.lattice_constant = atof( the_token);
611
612	// Next line is nrho, drho, nr, dr and rcut
613	eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
614	if(eam_eof_test == NULL){
615	sprintf( painCave.errMsg,
616	"error in reading number of points line in %s\n", eamPotFile);
617	painCave.isFatal = 1;
618	simError();
619	}
620
621	if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
622	sprintf( painCave.errMsg,
623	"Error parseing EAM nrho: line in %s\n", eamPotFile );
624	painCave.isFatal = 1;
625	simError();
626	}
627
628	info.eam_nrho = atoi( the_token );
629
630	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
631	sprintf( painCave.errMsg,
632	"Error parsing EAM drho in %s\n", eamPotFile );
633	painCave.isFatal = 1;
634	simError();
635	}
636	info.eam_drho = atof( the_token);
637
638	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
639	sprintf( painCave.errMsg,
640	"Error parsing EAM # r in %s\n", eamPotFile );
641	painCave.isFatal = 1;
642	simError();
643	}
644	info.eam_nr = atoi( the_token);
645
646	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
647	sprintf( painCave.errMsg,
648	"Error parsing EAM dr in %s\n", eamPotFile );
649	painCave.isFatal = 1;
650	simError();
651	}
652	info.eam_dr = atof( the_token);
653
654	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
655	sprintf( painCave.errMsg,
656	"Error parsing EAM rcut in %s\n", eamPotFile );
657	painCave.isFatal = 1;
658	simError();
659	}
660	info.eam_rcut = atof( the_token);
661
662
663
664
665
666	// Ok now we have to allocate point arrays and read in number of points
667	// Index the arrays for fortran, starting at 1
668	myEam_Frhovals = new double[info.eam_nrho];
669	myEam_rvals = new double[info.eam_nr];
670	myEam_rhovals = new double[info.eam_nr];
671
672	// Parse F of rho vals.
673
674	// Assume for now that we have a complete number of lines
675	nReadLines = int(info.eam_nrho/5);
676
677
678
679	for (i=0;i<nReadLines;i++){
680	j = i*5;
681
682	// Read next line
683	eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
684	linenumber++;
685	if(eam_eof_test == NULL){
686	sprintf( painCave.errMsg,
687	"error in reading EAM file %s at line %d\n",
688	eamPotFile,linenumber);
689	painCave.isFatal = 1;
690	simError();
691	}
692
693	// Parse 5 values on each line into array
694	// Value 1
695	if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
696	sprintf( painCave.errMsg,
697	"Error parsing EAM nrho: line in %s\n", eamPotFile );
698	painCave.isFatal = 1;
699	simError();
700	}
701
702	myEam_Frhovals[j+0] = atof( the_token );
703
704	// Value 2
705	if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
706	sprintf( painCave.errMsg,
707	"Error parsing EAM nrho: line in %s\n", eamPotFile );
708	painCave.isFatal = 1;
709	simError();
710	}
711
712	myEam_Frhovals[j+1] = atof( the_token );
713
714	// Value 3
715	if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
716	sprintf( painCave.errMsg,
717	"Error parsing EAM nrho: line in %s\n", eamPotFile );
718	painCave.isFatal = 1;
719	simError();
720	}
721
722	myEam_Frhovals[j+2] = atof( the_token );
723
724	// Value 4
725	if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
726	sprintf( painCave.errMsg,
727	"Error parsing EAM nrho: line in %s\n", eamPotFile );
728	painCave.isFatal = 1;
729	simError();
730	}
731
732	myEam_Frhovals[j+3] = atof( the_token );
733
734	// Value 5
735	if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
736	sprintf( painCave.errMsg,
737	"Error parsing EAM nrho: line in %s\n", eamPotFile );
738	painCave.isFatal = 1;
739	simError();
740	}
741
742	myEam_Frhovals[j+4] = atof( the_token );
743
744	}
745	// Parse Z of r vals
746
747	// Assume for now that we have a complete number of lines
748	nReadLines = int(info.eam_nr/5);
749
750	for (i=0;i<nReadLines;i++){
751	j = i*5;
752
753	// Read next line
754	eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
755	linenumber++;
756	if(eam_eof_test == NULL){
757	sprintf( painCave.errMsg,
758	"error in reading EAM file %s at line %d\n",
759	eamPotFile,linenumber);
760	painCave.isFatal = 1;
761	simError();
762	}
763
764	// Parse 5 values on each line into array
765	// Value 1
766	if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
767	sprintf( painCave.errMsg,
768	"Error parsing EAM nrho: line in %s\n", eamPotFile );
769	painCave.isFatal = 1;
770	simError();
771	}
772
773	myEam_rvals[j+0] = atof( the_token );
774
775	// Value 2
776	if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
777	sprintf( painCave.errMsg,
778	"Error parsing EAM nrho: line in %s\n", eamPotFile );
779	painCave.isFatal = 1;
780	simError();
781	}
782
783	myEam_rvals[j+1] = atof( the_token );
784
785	// Value 3
786	if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
787	sprintf( painCave.errMsg,
788	"Error parsing EAM nrho: line in %s\n", eamPotFile );
789	painCave.isFatal = 1;
790	simError();
791	}
792
793	myEam_rvals[j+2] = atof( the_token );
794
795	// Value 4
796	if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
797	sprintf( painCave.errMsg,
798	"Error parsing EAM nrho: line in %s\n", eamPotFile );
799	painCave.isFatal = 1;
800	simError();
801	}
802
803	myEam_rvals[j+3] = atof( the_token );
804
805	// Value 5
806	if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
807	sprintf( painCave.errMsg,
808	"Error parsing EAM nrho: line in %s\n", eamPotFile );
809	painCave.isFatal = 1;
810	simError();
811	}
812
813	myEam_rvals[j+4] = atof( the_token );
814
815	}
816	// Parse rho of r vals
817
818	// Assume for now that we have a complete number of lines
819
820	for (i=0;i<nReadLines;i++){
821	j = i*5;
822
823	// Read next line
824	eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
825	linenumber++;
826	if(eam_eof_test == NULL){
827	sprintf( painCave.errMsg,
828	"error in reading EAM file %s at line %d\n",
829	eamPotFile,linenumber);
830	painCave.isFatal = 1;
831	simError();
832	}
833
834	// Parse 5 values on each line into array
835	// Value 1
836	if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
837	sprintf( painCave.errMsg,
838	"Error parsing EAM nrho: line in %s\n", eamPotFile );
839	painCave.isFatal = 1;
840	simError();
841	}
842
843	myEam_rhovals[j+0] = atof( the_token );
844
845	// Value 2
846	if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
847	sprintf( painCave.errMsg,
848	"Error parsing EAM nrho: line in %s\n", eamPotFile );
849	painCave.isFatal = 1;
850	simError();
851	}
852
853	myEam_rhovals[j+1] = atof( the_token );
854
855	// Value 3
856	if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
857	sprintf( painCave.errMsg,
858	"Error parsing EAM nrho: line in %s\n", eamPotFile );
859	painCave.isFatal = 1;
860	simError();
861	}
862
863	myEam_rhovals[j+2] = atof( the_token );
864
865	// Value 4
866	if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
867	sprintf( painCave.errMsg,
868	"Error parsing EAM nrho: line in %s\n", eamPotFile );
869	painCave.isFatal = 1;
870	simError();
871	}
872
873	myEam_rhovals[j+3] = atof( the_token );
874
875	// Value 5
876	if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
877	sprintf( painCave.errMsg,
878	"Error parsing EAM nrho: line in %s\n", eamPotFile );
879	painCave.isFatal = 1;
880	simError();
881	}
882
883	myEam_rhovals[j+4] = atof( the_token );
884
885	}
886	*eam_rvals = myEam_rvals;
887	*eam_rhovals = myEam_rhovals;
888	*eam_Frhovals = myEam_Frhovals;
889
890	fclose(eamFile);
891	return 0;
892	}