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"

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