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 2175 by chrisfen, Tue Apr 12 20:57:32 2005 UTC

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

Diff Legend

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