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_module.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:	1599
Committed:	Tue Oct 19 20:40:17 2004 UTC (19 years, 8 months ago) by gezelter
File size:	21113 byte(s)
Log Message:	working on shapes
#	User	Rev	Content
1	gezelter	1520	#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	gezelter	1599	#include "UseTheForce/DarkSide/shapes_module.h"
18	gezelter	1520
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			}