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

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

Diff Legend

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