ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-4/src/UseTheForce/Shapes_FF.cpp
(Generate patch)

Comparing:
trunk/OOPSE-4/src/UseTheForce/Shapes_FF.cpp (file contents), Revision 1520 by gezelter, Mon Oct 4 15:27:35 2004 UTC vs.
branches/new_design/OOPSE-4/src/UseTheForce/Shapes_FF.cpp (file contents), Revision 1683, Thu Oct 28 22:34:02 2004 UTC

# Line 1 | Line 1
1 + /*
2 + *  Shapes_FF.cpp
3 + *  oopse
4 + *
5 + *  Created by Chris Fennell on 10/20/04.
6 + *  Copyright 2004 __MyCompanyName__. All rights reserved.
7 + *
8 + */
9 +
10   #include <stdlib.h>
11   #include <stdio.h>
12   #include <string.h>
13 <
13 > #include <map>
14 > #include <cmath>
15   #include <iostream>
6 – using namespace std;
16  
17   #ifdef IS_MPI
18   #include <mpi.h>
# Line 12 | Line 21 | using namespace std;
21   #include "UseTheForce/ForceFields.hpp"
22   #include "primitives/SRI.hpp"
23   #include "utils/simError.h"
24 + #include "utils/StringUtils.hpp"
25 + #include "io/basic_ifstrstream.hpp"
26 + #include "math/RealSphericalHarmonic.hpp"
27 + #include "math/SquareMatrix3.hpp"
28 + #include "types/ShapeAtomType.hpp"
29 + #include "UseTheForce/DarkSide/atype_interface.h"
30 + #include "UseTheForce/DarkSide/shapes_interface.h"
31  
16 – #include "UseTheForce/fortranWrappers.hpp"
17 –
32   #ifdef IS_MPI
33   #include "UseTheForce/mpiForceField.h"
34   #endif // is_mpi
35  
36 < Shapes_FF::Shapes_FF() {
37 <  Shapes_FF("");
24 < }
36 > using namespace std;
37 > using namespace oopse;
38  
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 –
39   Shapes_FF::~Shapes_FF(){
40  
41   #ifdef IS_MPI
42    if( worldRank == 0 ){
43   #endif // is_mpi
44      
45 <    fclose( frcFile );
45 >    forceFile.close();
46      
47   #ifdef IS_MPI
48    }
# Line 127 | Line 61 | void Shapes_FF::initForceField(int ljMixRule){
61   }
62  
63  
64 < void Shapes_FF::initForceField(int ljMixRule){
65 <  initFortran(ljMixRule, 0);
64 > void Shapes_FF::initForceField(){
65 >  initFortran(0);
66   }
67  
68  
69   void Shapes_FF::readParams( void ){
70  
71 <  char shapePotFile[1000];
138 <  
139 < #ifdef IS_MPI
140 <  if( worldRank == 0 ){
141 < #endif
142 <    
143 <    // read in the atom types.
71 >  char readLine[1024];
72  
73 <    fastForward( "AtomTypes", "eam atom readParams" );
73 >  string fileName;
74 >  string shapeFileName;
75 >  string tempString;
76  
77 <    // we are now at the AtomTypes section.
77 >  char *nameToken;
78 >  char *delim = " ,;\t\n";
79 >  int nTokens, i;
80 >  int nContact = 0;
81 >  int nRange = 0;
82 >  int nStrength = 0;
83 >  int myATID;
84 >  int isError;
85 >  string nameString;
86 >  AtomType* at;
87 >  DirectionalAtomType* dat;
88 >  ShapeAtomType* st;
89 >
90 >  map<string, AtomType*>::iterator iter;
91 >
92 >  // vectors for shape transfer to fortran
93 >  vector<RealSphericalHarmonic*> tempSHVector;
94 >  vector<int> contactL;
95 >  vector<int> contactM;
96 >  vector<int> contactFunc;
97 >  vector<double> contactCoeff;
98 >  vector<int> rangeL;
99 >  vector<int> rangeM;
100 >  vector<int> rangeFunc;
101 >  vector<double> rangeCoeff;
102 >  vector<int> strengthL;
103 >  vector<int> strengthM;
104 >  vector<int> strengthFunc;
105 >  vector<double> strengthCoeff;
106 >  
107 >  // generate the force file name  
108 >  fileName = "Shapes.frc";
109 >  
110 >  // attempt to open the file in the current directory first.
111 >  forceFile.open( fileName.c_str() );
112 >  
113 >  if( forceFile == NULL ){
114      
115 <    eof_test = fgets( readLine, sizeof(readLine), frcFile );
116 <    lineNum++;
115 >    tempString = ffPath;
116 >    tempString += "/";
117 >    tempString += fileName;
118 >    fileName = tempString;
119      
120 +    forceFile.open( fileName.c_str() );
121      
122 <    // read a line, and start parseing out the atom types
123 <
124 <    if( eof_test == NULL ){
125 <      sprintf( painCave.errMsg,
126 <               "Error in reading Atoms from force file at line %d.\n",
127 <               lineNum );
122 >    if( forceFile == NULL ){
123 >      
124 >      sprintf( painCave.errMsg,
125 >               "Error opening the force field parameter file:\n"
126 >               "\t%s\n"
127 >               "\tHave you tried setting the FORCE_PARAM_PATH environment "
128 >               "variable?\n",
129 >               fileName.c_str() );
130 >      painCave.severity = OOPSE_ERROR;
131        painCave.isFatal = 1;
132        simError();
133      }
134 <    
135 <    identNum = 1;
136 <    // stop reading at end of file, or at next section
134 >  }
135 >  
136 >  // read in the shape types.
137 >  
138 >  findBegin( forceFile, "ShapeTypes" );
139 >  
140 >  while( !forceFile.eof() ){
141 >    forceFile.getline( readLine, sizeof(readLine) );
142  
143 <    while( readLine[0] != '#' && eof_test != NULL ){
144 <
145 <      // toss comment lines
146 <      if( readLine[0] != '!' ){
147 <        
148 <        // the parser returns 0 if the line was blank
149 <        if( parseAtom( readLine, lineNum, info, eamPotFile ) ){
150 <          parseEAM(info,eamPotFile, &eam_rvals,
151 <                   &eam_rhovals, &eam_Frhovals);
152 <          info.ident = identNum;
153 <          headAtomType->add( info, eam_rvals,
154 <                             eam_rhovals,eam_Frhovals );
155 <          identNum++;
143 >    // toss comment lines
144 >    if( readLine[0] != '!' && readLine[0] != '#' ){
145 >      
146 >      if (isEndLine(readLine)) break;
147 >      
148 >      nTokens = countTokens(readLine, " ,;\t");
149 >      if (nTokens != 0) {
150 >        if (nTokens < 2) {
151 >          sprintf( painCave.errMsg,
152 >                   "Not enough data on a ShapeTypes line in file: %s\n",
153 >                   fileName.c_str() );
154 >          painCave.severity = OOPSE_ERROR;
155 >          painCave.isFatal = 1;
156 >          simError();
157          }
158 <      }
159 <      eof_test = fgets( readLine, sizeof(readLine), frcFile );
160 <      lineNum++;
183 <    }
184 <    
185 <    
158 >        
159 >        nameToken = strtok( readLine, delim );
160 >        shapeFileName = strtok( NULL, delim );
161  
162 < #ifdef IS_MPI
163 <  
189 <    
190 <    // send out the linked list to all the other processes
162 >        // strings are not char arrays!
163 >        nameString = nameToken;
164  
165 <    sprintf( checkPointMsg,
166 <             "Shapes_FF atom structures read successfully." );
167 <    MPIcheckPoint();
165 >        // does this AtomType name already exist in the map?
166 >        iter = atomTypeMap.find(nameString);    
167 >        if (iter == atomTypeMap.end()) {
168 >          // no, it doesn't, so we may proceed:
169 >          
170 >          st = new ShapeAtomType();
171  
172 <    currentAtomType = headAtomType->next; //skip the first element who is a place holder.
173 <    while( currentAtomType != NULL ){
174 <      currentAtomType->duplicate( info );
172 >          st->setName(nameString);
173 >          myATID = atomTypeMap.size();
174 >          st->setIdent(myATID);
175 >          parseShapeFile(shapeFileName, st);
176 >          st->complete();
177 >          atomTypeMap.insert(make_pair(nameString, st));
178 >          
179 >        } else {
180 >          // atomType map already contained this string (i.e. it was
181 >          // declared in a previous block, and we just need to add
182 >          // the shape-specific information for this AtomType:
183  
184 <
185 <
186 <      sendFrcStruct( &info, mpiAtomStructType );
187 <
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;
184 >          st = (ShapeAtomType*)(iter->second);
185 >          parseShapeFile(shapeFileName, st);
186 >        }
187 >      }
188      }
222 –    info.last = 1;
223 –    sendFrcStruct( &info, mpiAtomStructType );
224 –    
189    }
190  
191 <  else{
228 <    
229 <    // listen for node 0 to send out the force params
230 <    
231 <    MPIcheckPoint();
191 > #ifdef IS_MPI
192  
193 <    headAtomType = new LinkedAtomType;
194 <    receiveFrcStruct( &info, mpiAtomStructType );
193 >  // looks like all the processors have their ShapeType vectors ready...
194 >  sprintf( checkPointMsg,
195 >           "Shapes_FF shape objects read successfully." );
196 >  MPIcheckPoint();
197  
198 <    while( !info.last ){
237 <      
238 <      // allocate the arrays
198 > #endif // is_mpi
199  
200 <      eam_rvals    = new double[info.eam_nr];
201 <      eam_rhovals  = new double[info.eam_nr];
242 <      eam_Frhovals = new double[info.eam_nrho];
200 >  // pack up and send the necessary info to fortran
201 >  for (iter = atomTypeMap.begin(); iter != atomTypeMap.end(); ++iter){
202  
203 <      // 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 <      
203 >    at = (AtomType*)(iter->second);
204  
205 <      headAtomType->add( info, eam_rvals, eam_rhovals, eam_Frhovals );
257 <      
258 <      MPIcheckPoint();
205 >    if (at->isDirectional()) {
206  
207 <      receiveFrcStruct( &info, mpiAtomStructType );
207 >      dat = (DirectionalAtomType*)at;
208  
209 +      if (dat->isShape()) {
210  
211 <    }
212 <  }
213 < #endif // is_mpi
214 <
215 <  // call new A_types in fortran
216 <  
217 <  int isError;
218 <
219 <  // dummy variables
220 <  int isLJ = 0;
221 <  int isDipole = 0;
222 <  int isSSD = 0;
223 <  int isGB = 0;
224 <  int isEAM = 1;
225 <  int isCharge = 0;
226 <  double dipole = 0.0;
227 <  double charge = 0.0;
228 <  double eamSigma = 0.0;
229 <  double eamEpslon = 0.0;
230 <  
231 <  currentAtomType = headAtomType->next;
232 <  while( currentAtomType != NULL ){
233 <    
234 <    if( currentAtomType->name[0] != '\0' ){
235 <      isError = 0;
236 <      makeAtype( &(currentAtomType->ident),
237 <                 &isLJ,
238 <                 &isSSD,
239 <                 &isDipole,
240 <                 &isGB,
241 <                 &isEAM,
242 <                 &isCharge,
243 <                 &eamEpslon,
244 <                 &eamSigma,
245 <                 &charge,
246 <                 &dipole,
247 <                 &isError );
248 <      if( isError ){
249 <        sprintf( painCave.errMsg,
250 <                 "Error initializing the \"%s\" atom type in fortran\n",
251 <                 currentAtomType->name );
252 <        painCave.isFatal = 1;
253 <        simError();
211 >        st = (ShapeAtomType*)at;
212 >        
213 >        contactL.clear();
214 >        contactM.clear();
215 >        contactFunc.clear();
216 >        contactCoeff.clear();
217 >        
218 >        tempSHVector = st->getContactFuncs();
219 >        
220 >        nContact = tempSHVector.size();
221 >        for (i=0; i<nContact; i++){
222 >          contactL.push_back(tempSHVector[i]->getL());
223 >          contactM.push_back(tempSHVector[i]->getM());
224 >          contactFunc.push_back(tempSHVector[i]->getFunctionType());
225 >          contactCoeff.push_back(tempSHVector[i]->getCoefficient());
226 >        }
227 >        
228 >        rangeL.clear();
229 >        rangeM.clear();
230 >        rangeFunc.clear();
231 >        rangeCoeff.clear();
232 >        
233 >        tempSHVector = st->getRangeFuncs();
234 >        
235 >        nRange = tempSHVector.size();
236 >        for (i=0; i<nRange; i++){
237 >          rangeL.push_back(tempSHVector[i]->getL());
238 >          rangeM.push_back(tempSHVector[i]->getM());
239 >          rangeFunc.push_back(tempSHVector[i]->getFunctionType());
240 >          rangeCoeff.push_back(tempSHVector[i]->getCoefficient());
241 >        }
242 >        
243 >        strengthL.clear();
244 >        strengthM.clear();
245 >        strengthFunc.clear();
246 >        strengthCoeff.clear();
247 >        
248 >        tempSHVector = st->getStrengthFuncs();
249 >        
250 >        nStrength = tempSHVector.size();
251 >        for (i=0; i<nStrength; i++){
252 >          strengthL.push_back(tempSHVector[i]->getL());
253 >          strengthM.push_back(tempSHVector[i]->getM());
254 >          strengthFunc.push_back(tempSHVector[i]->getFunctionType());
255 >          strengthCoeff.push_back(tempSHVector[i]->getCoefficient());
256 >        }
257 >        
258 >        isError = 0;
259 >        myATID = at->getIdent();
260 >        makeShape( &nContact, &contactL[0], &contactM[0], &contactFunc[0],
261 >                   &contactCoeff[0],
262 >                   &nRange, &rangeL[0], &rangeM[0], &rangeFunc[0],
263 >                   &rangeCoeff[0],
264 >                   &nStrength, &strengthL[0], &strengthM[0],
265 >                   &strengthFunc[0], &strengthCoeff[0],
266 >                   &myATID,
267 >                   &isError);
268 >        if( isError ){
269 >          sprintf( painCave.errMsg,
270 >                   "Error initializing the \"%s\" shape in fortran\n",
271 >                   (iter->first).c_str() );
272 >          painCave.isFatal = 1;
273 >          simError();
274 >        }
275        }
276      }
308 –    currentAtomType = currentAtomType->next;
277    }
278 <      
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 <
278 >  
279   #ifdef IS_MPI
280    sprintf( checkPointMsg,
281             "Shapes_FF atom structures successfully sent to fortran\n" );
282    MPIcheckPoint();
283   #endif // is_mpi
284  
285 + }
286  
287 + void Shapes_FF::cleanMe( void ){
288  
289   }
290  
355 –
291   void Shapes_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
292    
293 <  int i;
294 <  
293 >  int i,j,k;
294 >  map<string, AtomType*>::iterator iter;
295 >
296    // initialize the atoms
297 <  
297 >  DirectionalAtom* dAtom;
298 >  AtomType* at;
299 >  DirectionalAtomType* dat;
300 >  ShapeAtomType* sat;
301 >  double sigma;
302 >  double ji[3];
303 >  double inertialMat[3][3];
304 >  Mat3x3d momInt;
305 >  string myTypeString;
306 >
307    for( i=0; i<nAtoms; i++ ){
308      
309 <    currentAtomType = headAtomType->find( the_atoms[i]->getType() );
310 <    if( currentAtomType == NULL ){
309 >    myTypeString = the_atoms[i]->getType();
310 >    iter = atomTypeMap.find(myTypeString);
311 >
312 >    if (iter == atomTypeMap.end()) {
313        sprintf( painCave.errMsg,
314                 "AtomType error, %s not found in force file.\n",
315                 the_atoms[i]->getType() );
316        painCave.isFatal = 1;
317        simError();
318 <    }
372 <    
373 <    the_atoms[i]->setMass( currentAtomType->mass );
374 <    the_atoms[i]->setIdent( currentAtomType->ident );
318 >    } else {
319  
320 <    if (eamRcut < currentAtomType->eam_rcut) eamRcut = currentAtomType->eam_rcut;
377 <    
378 <  }
379 < }
320 >      at = (AtomType*)(iter->second);
321  
322 < void Shapes_FF::initializeBonds( int nBonds, Bond** BondArray,
323 <                             bond_pair* the_bonds ){
324 <  
325 <    if( nBonds ){
326 <      sprintf( painCave.errMsg,
327 <               "LJ_FF does not support bonds.\n" );
328 <      painCave.isFatal = 1;
329 <      simError();
389 <    }
390 < }
322 >      the_atoms[i]->setMass( at->getMass() );
323 >      the_atoms[i]->setIdent( at->getIdent() );
324 >      
325 >      if ( at->isShape() ) {
326 >        
327 >        sat = (ShapeAtomType*)at;
328 >        sigma = findLargestContactDistance(sat);
329 >        if (sigma > bigSigma) bigSigma = sigma;
330  
331 < void Shapes_FF::initializeBends( int nBends, Bend** bendArray,
393 <                             bend_set* the_bends ){
331 >      }
332  
333 <    if( nBends ){
396 <      sprintf( painCave.errMsg,
397 <               "LJ_FF does not support bends.\n" );
398 <      painCave.isFatal = 1;
399 <      simError();
400 <    }
401 < }
333 >      the_atoms[i]->setHasCharge(at->isCharge());
334  
335 < void Shapes_FF::initializeTorsions( int nTorsions, Torsion** torsionArray,
404 <                                torsion_set* the_torsions ){
335 >      if( at->isDirectional() ){
336  
337 <    if( nTorsions ){
338 <      sprintf( painCave.errMsg,
339 <               "LJ_FF does not support torsions.\n" );
340 <      painCave.isFatal = 1;
410 <      simError();
411 <    }
412 < }
337 >        dat = (DirectionalAtomType*)at;
338 >        dAtom = (DirectionalAtom *) the_atoms[i];
339 >        
340 >        momInt = dat->getI();
341  
342 < void Shapes_FF::fastForward( char* stopText, char* searchOwner ){
342 >        // zero out the moments of inertia matrix
343 >        for( j=0; j<3; j++ )
344 >          for( k=0; k<3; k++ )
345 >            inertialMat[j][k] = momInt(j,k);
346 >        dAtom->setI( inertialMat );            
347  
348 <  int foundText = 0;
349 <  char* the_token;
348 >        ji[0] = 0.0;
349 >        ji[1] = 0.0;
350 >        ji[2] = 0.0;
351 >        dAtom->setJ( ji );
352  
353 <  rewind( frcFile );
354 <  lineNum = 0;
353 >        if (dat->isDipole()) {
354 >          dAtom->setHasDipole( dat->isDipole() );
355 >          entry_plug->n_dipoles++;
356 >        }              
357 >      }
358 >    }
359 >  }
360 > }
361  
362 <  eof_test = fgets( readLine, sizeof(readLine), frcFile );
363 <  lineNum++;
364 <  if( eof_test == NULL ){
365 <    sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
366 <             " file is empty.\n",
367 <             searchOwner );
362 > void Shapes_FF::initializeBonds( int nBonds, Bond** BondArray,
363 >                                 bond_pair* the_bonds ){
364 >  
365 >  if( nBonds ){
366 >    sprintf( painCave.errMsg,
367 >             "Shapes_FF does not support bonds.\n" );
368      painCave.isFatal = 1;
369      simError();
370    }
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 –  }  
371   }
372  
373 <
374 <
470 < int EAM_NS::parseAtom( char *lineBuffer, int lineNum,   atomStruct &info, char *eamPotFile ){
471 <
472 <  char* the_token;
373 > void Shapes_FF::initializeBends( int nBends, Bend** bendArray,
374 >                                 bend_set* the_bends ){
375    
376 <  the_token = strtok( lineBuffer, " \n\t,;" );
377 <  if( the_token != NULL ){
378 <    
379 <    strcpy( info.name, the_token );
380 <          
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;
376 >  if( nBends ){
377 >    sprintf( painCave.errMsg,
378 >             "Shapes_FF does not support bends.\n" );
379 >    painCave.isFatal = 1;
380 >    simError();
381    }
498 –  else return 0;
382   }
383  
384 < int EAM_NS::parseEAM(atomStruct &info, char *eamPotFile,
385 <                     double **eam_rvals,
386 <                     double **eam_rhovals,
387 <                     double **eam_Frhovals){
388 <  double* myEam_rvals;
389 <  double* myEam_rhovals;
390 <  double* myEam_Frhovals;
384 > void Shapes_FF::initializeTorsions( int nTorsions, Torsion** torsionArray,
385 >                                    torsion_set* the_torsions ){
386 >  
387 >  if( nTorsions ){
388 >    sprintf( painCave.errMsg,
389 >             "Shapes_FF does not support torsions.\n" );
390 >    painCave.isFatal = 1;
391 >    simError();
392 >  }
393 > }
394  
395 <  char* ffPath_env = "FORCE_PARAM_PATH";
396 <  char* ffPath;
397 <  char* the_token;
398 <  char* eam_eof_test;
399 <  FILE *eamFile;
400 <  const int BUFFERSIZE = 3000;
395 > void Shapes_FF::parseShapeFile(string shapeFileName, ShapeAtomType* st){
396 >  const int MAXLEN = 1024;
397 >  char inLine[MAXLEN];  
398 >  char *token;
399 >  char *delim = " ,;\t\n";
400 >  int nTokens;
401 >  Mat3x3d momInert;
402 >  RealSphericalHarmonic* rsh;
403 >  vector<RealSphericalHarmonic*> functionVector;
404 >  ifstrstream shapeFile;
405 >  string tempString;
406  
407 <  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" );
407 >  shapeFile.open( shapeFileName.c_str() );
408    
409 <  
530 <  if( eamFile == NULL ){
409 >  if( shapeFile == NULL ){
410      
411 <      // next see if the force path enviorment variable is set
411 >    tempString = ffPath;
412 >    tempString += "/";
413 >    tempString += shapeFileName;
414 >    shapeFileName = tempString;
415 >        
416 >    shapeFile.open( shapeFileName.c_str() );
417      
418 <    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 ){
418 >    if( shapeFile == NULL ){
419        
420        sprintf( painCave.errMsg,
421 <               "Error opening the EAM force parameter file: %s\n"
422 <               "Have you tried setting the FORCE_PARAM_PATH environment "
421 >               "Error opening the shape file:\n"
422 >               "\t%s\n"
423 >               "\tHave you tried setting the FORCE_PARAM_PATH environment "
424                 "variable?\n",
425 <               eamPotFile );
425 >               shapeFileName.c_str() );
426 >      painCave.severity = OOPSE_ERROR;
427        painCave.isFatal = 1;
428        simError();
429      }
430    }
431 +  
432 +  // read in the shape types.
433 +  
434 +  // first grab the values in the ShapeInfo section
435 +  findBegin( shapeFile, "ShapeInfo");
436 +  
437 +  shapeFile.getline(inLine, MAXLEN);
438 +  while( !shapeFile.eof() ) {
439 +    // toss comment lines
440 +    if( inLine[0] != '!' && inLine[0] != '#' ){
441 +      // end marks section completion
442 +      if (isEndLine(inLine)) break;
443 +      
444 +      nTokens = countTokens(inLine, delim);
445 +      if (nTokens != 0) {
446 +        if (nTokens < 5) {
447 +          sprintf( painCave.errMsg,
448 +                   "Not enough data on a ShapeInfo line in file: %s\n",
449 +                   shapeFileName.c_str() );
450 +          painCave.severity = OOPSE_ERROR;
451 +          painCave.isFatal = 1;
452 +          simError();
453  
454 <  // First line is a comment line, read and toss it....
455 <  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
456 <  linenumber++;
457 <  if(eam_eof_test == NULL){
458 <    sprintf( painCave.errMsg,
459 <             "error in reading commment in %s\n", eamPotFile);
460 <    painCave.isFatal = 1;
461 <    simError();
454 >          token = strtok(inLine, delim);
455 >          token = strtok(NULL, delim);
456 >          st->setMass(atof(token));
457 >          token = strtok(NULL, delim);
458 >          momInert(0,0) = atof(token);
459 >          token = strtok(NULL, delim);
460 >          momInert(1,1) = atof(token);
461 >          token = strtok(NULL, delim);
462 >          momInert(2,2) = atof(token);
463 >          st->setI(momInert);
464 >        }
465 >      }
466 >    }
467 >    shapeFile.getline(inLine, MAXLEN);
468    }
469  
470 <
471 <
472 <  // 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 <  }
470 >  // now grab the contact functions
471 >  findBegin(shapeFile, "ContactFunctions");
472 >  functionVector.clear();
473    
474 <  info.eam_ident = atoi( the_token );
475 <
476 <  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
477 <    sprintf( painCave.errMsg,
478 <             "Error parsing EAM mass in %s\n", eamPotFile );
479 <    painCave.isFatal = 1;
480 <    simError();
474 >  shapeFile.getline(inLine, MAXLEN);
475 >  while( !shapeFile.eof() ) {
476 >    // toss comment lines
477 >    if( inLine[0] != '!' && inLine[0] != '#' ){
478 >      // end marks section completion
479 >      if (isEndLine(inLine)) break;
480 >      
481 >      nTokens = countTokens(inLine, delim);
482 >      if (nTokens != 0) {
483 >        if (nTokens < 4) {
484 >          sprintf( painCave.errMsg,
485 >                   "Not enough data on a ContactFunctions line in file: %s\n",
486 >                   shapeFileName.c_str() );
487 >          painCave.severity = OOPSE_ERROR;
488 >          painCave.isFatal = 1;
489 >          simError();
490 >          
491 >          // read in a spherical harmonic function
492 >          token = strtok(inLine, delim);
493 >          rsh = new RealSphericalHarmonic();
494 >          rsh->setL(atoi(token));
495 >          token = strtok(NULL, delim);
496 >          rsh->setM(atoi(token));
497 >          token = strtok(NULL, delim);
498 >          if (!strcasecmp("sin",token))
499 >            rsh->makeSinFunction();
500 >          else
501 >            rsh->makeCosFunction();            
502 >          token = strtok(NULL, delim);
503 >          rsh->setCoefficient(atof(token));
504 >          
505 >          functionVector.push_back(rsh);
506 >        }
507 >      }
508 >    }
509 >    shapeFile.getline(inLine, MAXLEN);
510    }
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 –  }
511  
512 <  if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
513 <    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);
512 >  // pass contact functions to ShapeType
513 >  st->setContactFuncs(functionVector);
514  
515 <  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
516 <    sprintf( painCave.errMsg,
517 <             "Error parsing EAM # r in %s\n", eamPotFile );
640 <    painCave.isFatal = 1;
641 <    simError();
642 <  }
643 <  info.eam_nr = atoi( the_token);
515 >  // now grab the range functions
516 >  findBegin(shapeFile, "RangeFunctions");
517 >  functionVector.clear();
518    
519 <  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
520 <    sprintf( painCave.errMsg,
521 <             "Error parsing EAM dr in %s\n", eamPotFile );
522 <    painCave.isFatal = 1;
523 <    simError();
524 <  }
525 <  info.eam_dr = atof( the_token);
519 >  shapeFile.getline(inLine, MAXLEN);
520 >  while( !shapeFile.eof() ) {
521 >    // toss comment lines
522 >    if( inLine[0] != '!' && inLine[0] != '#' ){
523 >      // end marks section completion
524 >      if (isEndLine(inLine)) break;
525 >      
526 >      nTokens = countTokens(inLine, delim);
527 >      if (nTokens != 0) {
528 >        if (nTokens < 4) {
529 >          sprintf( painCave.errMsg,
530 >                   "Not enough data on a RangeFunctions line in file: %s\n",
531 >                   shapeFileName.c_str() );
532 >          painCave.severity = OOPSE_ERROR;
533 >          painCave.isFatal = 1;
534 >          simError();
535 >          
536 >          // read in a spherical harmonic function
537 >          token = strtok(inLine, delim);
538  
539 <  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
540 <    sprintf( painCave.errMsg,
541 <             "Error parsing EAM rcut in %s\n", eamPotFile );
542 <    painCave.isFatal = 1;
543 <    simError();
544 <  }
545 <  info.eam_rcut = atof( the_token);
546 <
547 <
548 <
549 <
550 <
551 <  // Ok now we have to allocate point arrays and read in number of points
552 <  // Index the arrays for fortran, starting at 1
553 <  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();
539 >          rsh = new RealSphericalHarmonic();
540 >          rsh->setL(atoi(token));
541 >          token = strtok(NULL, delim);
542 >          rsh->setM(atoi(token));
543 >          token = strtok(NULL, delim);
544 >          if (!strcasecmp("sin",token))
545 >            rsh->makeSinFunction();
546 >          else
547 >            rsh->makeCosFunction();            
548 >          token = strtok(NULL, delim);
549 >          rsh->setCoefficient(atof(token));
550 >          
551 >          functionVector.push_back(rsh);
552 >        }
553 >      }
554      }
555 <    
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 <    
555 >    shapeFile.getline(inLine, MAXLEN);
556    }
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);
557  
558 <  for (i=0;i<nReadLines;i++){
559 <    j = i*5;
558 >  // pass range functions to ShapeType
559 >  st->setRangeFuncs(functionVector);
560  
561 <    // Read next line
562 <    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
563 <    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 <    }
561 >  // finally grab the strength functions
562 >  findBegin(shapeFile, "StrengthFunctions");
563 >  functionVector.clear();
564    
565 <    myEam_rvals[j+1] = atof( the_token );
566 <
567 <    // Value 3
568 <    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
569 <      sprintf( painCave.errMsg,
570 <               "Error parsing EAM nrho: line in %s\n", eamPotFile );
571 <      painCave.isFatal = 1;
572 <      simError();
565 >  shapeFile.getline(inLine, MAXLEN);
566 >  while( !shapeFile.eof() ) {
567 >    // toss comment lines
568 >    if( inLine[0] != '!' && inLine[0] != '#' ){
569 >      // end marks section completion
570 >      if (isEndLine(inLine)) break;
571 >      
572 >      nTokens = countTokens(inLine, delim);
573 >      if (nTokens != 0) {
574 >        if (nTokens < 4) {
575 >          sprintf( painCave.errMsg,
576 >                   "Not enough data on a StrengthFunctions line in file: %s\n",
577 >                   shapeFileName.c_str() );
578 >          painCave.severity = OOPSE_ERROR;
579 >          painCave.isFatal = 1;
580 >          simError();
581 >          
582 >          // read in a spherical harmonic function
583 >          token = strtok(inLine, delim);
584 >          rsh = new RealSphericalHarmonic();
585 >          rsh->setL(atoi(token));
586 >          token = strtok(NULL, delim);
587 >          rsh->setM(atoi(token));
588 >          token = strtok(NULL, delim);
589 >          if (!strcasecmp("sin",token))
590 >            rsh->makeSinFunction();
591 >          else
592 >            rsh->makeCosFunction();            
593 >          token = strtok(NULL, delim);
594 >          rsh->setCoefficient(atof(token));
595 >          
596 >          functionVector.push_back(rsh);
597 >        }
598 >      }
599      }
600 <  
601 <    myEam_rvals[j+2] = atof( the_token );
600 >    shapeFile.getline(inLine, MAXLEN);
601 >  }
602  
603 <    // Value 4
604 <    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 );
603 >  // pass strength functions to ShapeType
604 >  st->setStrengthFuncs(functionVector);
605  
606 <    // Value 5
607 <    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
608 <      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 );
606 >  // we're done reading from this file
607 >  shapeFile.close();
608 > }
609  
610 <  }
611 <  // Parse rho of r vals
610 > double Shapes_FF::findLargestContactDistance(ShapeAtomType* st) {
611 >  int i, j,  nSteps;
612 >  double theta, thetaStep, thetaMin, costheta;
613 >  double phi, phiStep;
614 >  double sigma, bigSigma;
615  
616 <  // Assume for now that we have a complete number of lines
616 >  nSteps = 16;
617  
618 <  for (i=0;i<nReadLines;i++){
619 <    j = i*5;
620 <
621 <    // 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 <    }
618 >  thetaStep = M_PI / nSteps;
619 >  thetaMin = thetaStep / 2.0;
620 >  phiStep = thetaStep * 2.0;
621 >  bigSigma = 0.0;
622    
623 <    // Parse 5 values on each line into array
624 <    // Value 1
625 <    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
626 <      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 );
623 >  for (i = 0; i < nSteps; i++) {
624 >    
625 >    theta = thetaMin + i * thetaStep;
626 >    costheta = cos(theta);
627  
628 <    // 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 );
628 >    for (j = 0; j < nSteps; j++) {
629  
630 <    // 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 );
630 >      phi = j*phiStep;
631  
632 <    // Value 4
633 <    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
634 <      sprintf( painCave.errMsg,
867 <               "Error parsing EAM nrho: line in %s\n", eamPotFile );
868 <      painCave.isFatal = 1;
869 <      simError();
632 >      sigma = st->getContactValueAt(costheta, phi);
633 >      
634 >      if (sigma > bigSigma) bigSigma = sigma;
635      }
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 –
636    }
885 –  *eam_rvals = myEam_rvals;
886 –  *eam_rhovals = myEam_rhovals;
887 –  *eam_Frhovals = myEam_Frhovals;
637  
638 <  fclose(eamFile);
890 <  return 0;
638 >  return bigSigma;  
639   }

Diff Legend

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