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

Comparing trunk/OOPSE-2.0/src/UseTheForce/Shapes_FF.cpp (file contents):
Revision 1520 by gezelter, Mon Oct 4 15:27:35 2004 UTC vs.
Revision 1646 by chrisfen, Tue Oct 26 18:02:46 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 <iostream>
15 +
16   using namespace std;
17 + using namespace oopse;
18  
19   #ifdef IS_MPI
20   #include <mpi.h>
# Line 12 | Line 23 | using namespace std;
23   #include "UseTheForce/ForceFields.hpp"
24   #include "primitives/SRI.hpp"
25   #include "utils/simError.h"
26 + #include "io/basic_ifstrstream.hpp"
27 + #include "math/RealSphericalHarmonic.hpp"
28 + #include "math/SquareMatrix3.hpp"
29  
30 < #include "UseTheForce/fortranWrappers.hpp"
30 > #include "UseTheForce/DarkSide/atype_interface.h"
31 > #include "UseTheForce/DarkSide/shapes_interface.h"
32  
33   #ifdef IS_MPI
34   #include "UseTheForce/mpiForceField.h"
# Line 24 | Line 39 | Shapes_FF::Shapes_FF(char* the_variant){
39   }
40  
41   Shapes_FF::Shapes_FF(char* the_variant){
42 +  ffPath_env = "FORCE_PARAM_PATH";
43  
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
44   }
45  
104 –
46   Shapes_FF::~Shapes_FF(){
47  
48   #ifdef IS_MPI
# Line 127 | Line 68 | void Shapes_FF::initForceField(int ljMixRule){
68   }
69  
70  
71 < void Shapes_FF::initForceField(int ljMixRule){
72 <  initFortran(ljMixRule, 0);
71 > void Shapes_FF::initForceField(){
72 >  initFortran(0);
73   }
74  
75  
76 < void Shapes_FF::readParams( void ){
76 > void Shapes_FF::readForceFile( void ){
77  
78 <  char shapePotFile[1000];
79 <  
80 < #ifdef IS_MPI
81 <  if( worldRank == 0 ){
82 < #endif
83 <    
84 <    // read in the atom types.
78 >  char readLine[1024];
79 >  char fileName[200];
80 >  char temp[200];
81 >  char* ffPath;
82 >  char *shapeFileName;
83 >  char *nameToken;
84 >  char *delim = " ,;\t\n";
85 >  int nTokens, i;
86 >  int nContact = 0;
87 >  int nRange = 0;
88 >  int nStrength = 0;
89 >  int myATID;
90 >  string nameString;
91 >  ShapeType* st;
92 >  map<string, AtomType*> atomTypeMap;
93 >  map<string, AtomType*>::iterator iter;
94  
95 <    fastForward( "AtomTypes", "eam atom readParams" );
95 >  // vectors for shape transfer to fortran
96 >  vector<RealSphericalHarmonic> tempSHVector;
97 >  vector<int> contactL;
98 >  vector<int> contactM;
99 >  vector<int> contactFunc;
100 >  vector<double> contactCoeff;
101 >  vector<int> rangeL;
102 >  vector<int> rangeM;
103 >  vector<int> rangeFunc;
104 >  vector<double> rangeCoeff;
105 >  vector<int> strengthL;
106 >  vector<int> strengthM;
107 >  vector<int> strengthFunc;
108 >  vector<double> strengthCoeff;
109  
110 <    // we are now at the AtomTypes section.
110 >  // build a basic file reader
111 >  ifstrsteam frcReader;
112 >  
113 >  // generate the force file name  
114 >  strcpy( fileName, "Shapes.frc" );
115 >  
116 >  // attempt to open the file in the current directory first.
117 >  frcReader.open( fileName );
118 >  
119 >  if( frcReader == NULL ){
120      
121 <    eof_test = fgets( readLine, sizeof(readLine), frcFile );
150 <    lineNum++;
121 >    // next see if the force path enviorment variable is set
122      
123 +    ffPath = getenv( ffPath_env );
124 +    if( ffPath == NULL ) {
125 +      STR_DEFINE(ffPath, FRC_PATH );
126 +    }
127 +        
128 +    strcpy( temp, ffPath );
129 +    strcat( temp, "/" );
130 +    strcat( temp, fileName );
131 +    strcpy( fileName, temp );
132      
133 <    // read a line, and start parseing out the atom types
134 <
135 <    if( eof_test == NULL ){
136 <      sprintf( painCave.errMsg,
137 <               "Error in reading Atoms from force file at line %d.\n",
138 <               lineNum );
133 >    frcReader.open( fileName );
134 >    
135 >    if( frcFile == NULL ){
136 >      
137 >      sprintf( painCave.errMsg,
138 >               "Error opening the force field parameter file:\n"
139 >               "\t%s\n"
140 >               "\tHave you tried setting the FORCE_PARAM_PATH environment "
141 >               "variable?\n",
142 >               fileName );
143 >      painCave.severity = OOPSE_ERROR;
144        painCave.isFatal = 1;
145        simError();
146      }
147 <    
148 <    identNum = 1;
149 <    // stop reading at end of file, or at next section
147 >  }
148 >  
149 >  // read in the shape types.
150 >  
151 >  findBegin( "ShapeTypes" );
152 >  
153 >  while( !frcReader.eof() ){
154 >    frcReader.getline( readLine, sizeof(readLine) );
155  
156 <    while( readLine[0] != '#' && eof_test != NULL ){
157 <
158 <      // toss comment lines
159 <      if( readLine[0] != '!' ){
156 >    // toss comment lines
157 >    if( readLine[0] != '!' && readLine[0] != '#' ){
158 >      
159 >      if (isEndLine(readLine)) break;
160 >      
161 >      nTokens = count_tokens(readLine, " ,;\t");
162 >      if (nTokens != 0) {
163 >        if (nTokens < 2) {
164 >          sprintf( painCave.errMsg,
165 >                   "Not enough data on a ShapeTypes line in file: %s\n"
166 >                   fileName );
167 >          painCave.severity = OOPSE_ERROR;
168 >          painCave.isFatal = 1;
169 >          simError();
170 >        }
171          
172 <        // the parser returns 0 if the line was blank
173 <        if( parseAtom( readLine, lineNum, info, eamPotFile ) ){
174 <          parseEAM(info,eamPotFile, &eam_rvals,
175 <                   &eam_rhovals, &eam_Frhovals);
176 <          info.ident = identNum;
177 <          headAtomType->add( info, eam_rvals,
178 <                             eam_rhovals,eam_Frhovals );
179 <          identNum++;
172 >        nameToken = strtok( readLine, delim );
173 >        shapeFileName = strtok( NULL, delim );
174 >
175 >        // strings are not char arrays!
176 >        nameString = nameToken;
177 >
178 >        // does this AtomType name already exist in the map?
179 >        iter = atomTypeMap.find(nameString);    
180 >        if (iter == atomTypeMap.end()) {
181 >          // no, it doesn't, so we may proceed:
182 >          
183 >          st = new ShapeType();
184 >
185 >          st->setName(nameString);
186 >          myATID = atomTypeMap.size();
187 >          st->setIdent(myATID);
188 >          parseShapeFile(shapeFileName, st);
189 >          st->complete();
190 >          atomTypeMap.insert(make_pair(nameString, st));
191 >          
192 >        } else {
193 >          // atomType map already contained this string (i.e. it was
194 >          // declared in a previous block, and we just need to add
195 >          // the shape-specific information for this AtomType:
196 >
197 >          st = (ShapeType*)(iter->second);
198 >          parseShapeFile(shapeFileName, st);
199          }
200        }
181 –      eof_test = fgets( readLine, sizeof(readLine), frcFile );
182 –      lineNum++;
201      }
202 <    
185 <    
202 >  }
203  
204   #ifdef IS_MPI
188 –  
189 –    
190 –    // send out the linked list to all the other processes
205  
206 <    sprintf( checkPointMsg,
207 <             "Shapes_FF atom structures read successfully." );
208 <    MPIcheckPoint();
206 >  // looks like all the processors have their ShapeType vectors ready...
207 >  sprintf( checkPointMsg,
208 >           "Shapes_FF shape objects read successfully." );
209 >  MPIcheckPoint();
210  
211 <    currentAtomType = headAtomType->next; //skip the first element who is a place holder.
197 <    while( currentAtomType != NULL ){
198 <      currentAtomType->duplicate( info );
211 > #endif // is_mpi
212  
213 <
213 >  // pack up and send the necessary info to fortran
214 >  for (iter = atomTypeMap.begin(); iter != atomTypeMap.end(); ++iter){
215  
216 <      sendFrcStruct( &info, mpiAtomStructType );
216 >    at = (AtomType*)(iter.second);
217  
218 <      // 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);
218 >    if (at->isDirectional()) {
219  
220 <      sprintf( checkPointMsg,
216 <               "successfully sent EAM force type: \"%s\"\n",
217 <               info.name );
218 <      MPIcheckPoint();
220 >      dat = (DirectionalAtomType*)at;
221  
222 <      currentAtomType = currentAtomType->next;
221 <    }
222 <    info.last = 1;
223 <    sendFrcStruct( &info, mpiAtomStructType );
224 <    
225 <  }
222 >      if (dat->isShape()) {
223  
224 <  else{
225 <    
226 <    // listen for node 0 to send out the force params
227 <    
228 <    MPIcheckPoint();
229 <
230 <    headAtomType = new LinkedAtomType;
231 <    receiveFrcStruct( &info, mpiAtomStructType );
232 <
233 <    while( !info.last ){
234 <      
235 <      // allocate the arrays
236 <
237 <      eam_rvals    = new double[info.eam_nr];
238 <      eam_rhovals  = new double[info.eam_nr];
239 <      eam_Frhovals = new double[info.eam_nrho];
240 <
241 <      // We have to now broadcast the Arrays
242 <      MPI_Bcast(eam_rvals,
243 <                info.eam_nr,
244 <                MPI_DOUBLE,0,MPI_COMM_WORLD);
245 <      MPI_Bcast(eam_rhovals,
246 <                info.eam_nr,
247 <                MPI_DOUBLE,0,MPI_COMM_WORLD);
248 <      MPI_Bcast(eam_Frhovals,
249 <                info.eam_nrho,
250 <                MPI_DOUBLE,0,MPI_COMM_WORLD);
251 <      
252 <
253 <      headAtomType->add( info, eam_rvals, eam_rhovals, eam_Frhovals );
254 <      
255 <      MPIcheckPoint();
256 <
257 <      receiveFrcStruct( &info, mpiAtomStructType );
258 <
259 <
260 <    }
261 <  }
262 < #endif // is_mpi
263 <
264 <  // call new A_types in fortran
265 <  
266 <  int isError;
267 <
268 <  // dummy variables
269 <  int isLJ = 0;
270 <  int isDipole = 0;
271 <  int isSSD = 0;
272 <  int isGB = 0;
273 <  int isEAM = 1;
274 <  int isCharge = 0;
275 <  double dipole = 0.0;
276 <  double charge = 0.0;
277 <  double eamSigma = 0.0;
278 <  double eamEpslon = 0.0;
279 <  
280 <  currentAtomType = headAtomType->next;
281 <  while( currentAtomType != NULL ){
282 <    
283 <    if( currentAtomType->name[0] != '\0' ){
284 <      isError = 0;
285 <      makeAtype( &(currentAtomType->ident),
286 <                 &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();
224 >        st = (ShapeAtomType*)at;
225 >        
226 >        contactL.clear();
227 >        contactM.clear();
228 >        contactFunc.clear();
229 >        contactCoeff.clear();
230 >        
231 >        tempSHVector = st->getContactFuncs();
232 >        
233 >        nContact = tempSHVector.size();
234 >        for (i=0; i<nContact; i++){
235 >          contactL.push_back(tempSHVector[i].getL());
236 >          contactM.push_back(tempSHVector[i].getM());
237 >          contactFunc.push_back(tempSHVector[i].getFunctionType());
238 >          contactCoeff.push_back(tempSHVector[i].getCoefficient());
239 >        }
240 >        
241 >        rangeL.clear();
242 >        rangeM.clear();
243 >        rangeFunc.clear();
244 >        rangeCoeff.clear();
245 >        
246 >        tempSHVector = st->getRangeFuncs();
247 >        
248 >        nRange = tempSHVector.size();
249 >        for (i=0; i<nRange; i++){
250 >          rangeL.push_back(tempSHVector[i].getL());
251 >          rangeM.push_back(tempSHVector[i].getM());
252 >          rangeFunc.push_back(tempSHVector[i].getFunctionType());
253 >          rangeCoeff.push_back(tempSHVector[i].getCoefficient());
254 >        }
255 >        
256 >        strengthL.clear();
257 >        strengthM.clear();
258 >        strengthFunc.clear();
259 >        strengthCoeff.clear();
260 >        
261 >        tempSHVector = st->getStrengthFuncs();
262 >        
263 >        nStrength = tempSHVector.size();
264 >        for (i=0; i<nStrength; i++){
265 >          strengthL.push_back(tempSHVector[i].getL());
266 >          strengthM.push_back(tempSHVector[i].getM());
267 >          strengthFunc.push_back(tempSHVector[i].getFunctionType());
268 >          strengthCoeff.push_back(tempSHVector[i].getCoefficient());
269 >        }
270 >        
271 >        isError = 0;
272 >        myATID = at->getIdent();
273 >        makeShape( &nContact, &contactL, &contactM, &contactFunc,
274 >                   &contactCoeff,
275 >                   &nRange, &rangeL, &rangeM, &rangeFunc, &rangeCoeff,
276 >                   &nStrength, &strengthL, &strengthM,
277 >                   &strengthFunc, &strengthCoeff,
278 >                   &myATID,
279 >                   &isError);
280 >        if( isError ){
281 >          sprintf( painCave.errMsg,
282 >                   "Error initializing the \"%s\" shape in fortran\n",
283 >                   marker->first );
284 >          painCave.isFatal = 1;
285 >          simError();
286 >        }
287        }
288      }
308 –    currentAtomType = currentAtomType->next;
289    }
290 <      
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 <
290 >  
291   #ifdef IS_MPI
292    sprintf( checkPointMsg,
293             "Shapes_FF atom structures successfully sent to fortran\n" );
294    MPIcheckPoint();
295   #endif // is_mpi
296  
351 –
352 –
297   }
298  
299 <
356 < void Shapes_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
299 > void SHAPES_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
300    
301 <  int i;
302 <  
301 >  int i,j,k;
302 >
303    // initialize the atoms
304 <  
304 >  DirectionalAtom* dAtom;
305 >  AtomType* at;
306 >  DirectionalAtomType* dat;
307 >  double mySigma;
308 >  double ji[3];
309 >  double inertialMat[3][3];
310 >  Mat3x3d momInt;
311 >  string myTypeString;
312 >
313    for( i=0; i<nAtoms; i++ ){
314      
315 <    currentAtomType = headAtomType->find( the_atoms[i]->getType() );
316 <    if( currentAtomType == NULL ){
315 >    myTypeString = the_atoms[i]->getType();
316 >    iter = atomTypeMap.find(myTypeString);
317 >
318 >    if (iter == atomTypeMap.end()) {
319        sprintf( painCave.errMsg,
320                 "AtomType error, %s not found in force file.\n",
321                 the_atoms[i]->getType() );
322        painCave.isFatal = 1;
323        simError();
324 <    }
372 <    
373 <    the_atoms[i]->setMass( currentAtomType->mass );
374 <    the_atoms[i]->setIdent( currentAtomType->ident );
324 >    } else {
325  
326 <    if (eamRcut < currentAtomType->eam_rcut) eamRcut = currentAtomType->eam_rcut;
377 <    
378 <  }
379 < }
326 >      at = (AtomType*)(iter->second);
327  
328 < void Shapes_FF::initializeBonds( int nBonds, Bond** BondArray,
329 <                             bond_pair* the_bonds ){
383 <  
384 <    if( nBonds ){
385 <      sprintf( painCave.errMsg,
386 <               "LJ_FF does not support bonds.\n" );
387 <      painCave.isFatal = 1;
388 <      simError();
389 <    }
390 < }
328 >      the_atoms[i]->setMass( at->getMass() );
329 >      the_atoms[i]->setIdent( at->getIdent() );
330  
331 < void Shapes_FF::initializeBends( int nBends, Bend** bendArray,
332 <                             bend_set* the_bends ){
331 >      if( at->isLennardJones() ) {
332 >        mySigma = at->properties->getPropertyByName("sigma");
333 >        if( bigSigma < mySigma ) bigSigma = mySigma;
334 >      }
335  
336 <    if( nBends ){
396 <      sprintf( painCave.errMsg,
397 <               "LJ_FF does not support bends.\n" );
398 <      painCave.isFatal = 1;
399 <      simError();
400 <    }
401 < }
336 >      the_atoms[i]->setHasCharge(at->isCharge());
337  
338 < void Shapes_FF::initializeTorsions( int nTorsions, Torsion** torsionArray,
404 <                                torsion_set* the_torsions ){
338 >      if( at->isDirectional() ){
339  
340 <    if( nTorsions ){
341 <      sprintf( painCave.errMsg,
408 <               "LJ_FF does not support torsions.\n" );
409 <      painCave.isFatal = 1;
410 <      simError();
411 <    }
412 < }
340 >        dat = (DirectionalAtomType*)at;
341 >        dAtom = (DirectionalAtom *) the_atoms[i];
342  
343 < void Shapes_FF::fastForward( char* stopText, char* searchOwner ){
343 >        momInt = dat->getI();
344  
345 <  int foundText = 0;
346 <  char* the_token;
347 <
348 <  rewind( frcFile );
349 <  lineNum = 0;
421 <
422 <  eof_test = fgets( readLine, sizeof(readLine), frcFile );
423 <  lineNum++;
424 <  if( eof_test == NULL ){
425 <    sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
426 <             " file is empty.\n",
427 <             searchOwner );
428 <    painCave.isFatal = 1;
429 <    simError();
430 <  }
431 <  
432 <  
433 <  while( !foundText ){
434 <    while( eof_test != NULL && readLine[0] != '#' ){
435 <      eof_test = fgets( readLine, sizeof(readLine), frcFile );
436 <      lineNum++;
437 <    }
438 <    if( eof_test == NULL ){
439 <      sprintf( painCave.errMsg,
440 <               "Error fast forwarding force file for %s at "
441 <               "line %d: file ended unexpectedly.\n",
442 <               searchOwner,
443 <               lineNum );
444 <      painCave.isFatal = 1;
445 <      simError();
446 <    }
447 <    
448 <    the_token = strtok( readLine, " ,;\t#\n" );
449 <    foundText = !strcmp( stopText, the_token );
450 <    
451 <    if( !foundText ){
452 <      eof_test = fgets( readLine, sizeof(readLine), frcFile );
453 <      lineNum++;
454 <      
455 <      if( eof_test == NULL ){
456 <        sprintf( painCave.errMsg,
457 <                 "Error fast forwarding force file for %s at "
458 <                 "line %d: file ended unexpectedly.\n",
459 <                 searchOwner,
460 <                 lineNum );
461 <        painCave.isFatal = 1;
462 <        simError();
463 <      }
464 <    }
465 <  }  
466 < }
467 <
345 >        // zero out the moments of inertia matrix
346 >        for( j=0; j<3; j++ )
347 >          for( k=0; k<3; k++ )
348 >            inertialMat[j][k] = momInt(j,k);
349 >        dAtom->setI( inertialMat );            
350  
351 +        ji[0] = 0.0;
352 +        ji[1] = 0.0;
353 +        ji[2] = 0.0;
354 +        dAtom->setJ( ji );
355  
356 < int EAM_NS::parseAtom( char *lineBuffer, int lineNum,   atomStruct &info, char *eamPotFile ){
357 <
358 <  char* the_token;
359 <  
360 <  the_token = strtok( lineBuffer, " \n\t,;" );
475 <  if( the_token != NULL ){
476 <    
477 <    strcpy( info.name, the_token );
478 <          
479 <    if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
480 <      sprintf( painCave.errMsg,
481 <               "Error parseing AtomTypes: line %d\n", lineNum );
482 <      painCave.isFatal = 1;
483 <      simError();
356 >        if (dat->isDipole()) {
357 >          dAtom->setHasDipole( dat->isDipole() );
358 >          entry_plug->n_dipoles++;
359 >        }              
360 >      }
361      }
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;
362    }
498 –  else return 0;
363   }
364  
365 < int EAM_NS::parseEAM(atomStruct &info, char *eamPotFile,
366 <                     double **eam_rvals,
503 <                     double **eam_rhovals,
504 <                     double **eam_Frhovals){
505 <  double* myEam_rvals;
506 <  double* myEam_rhovals;
507 <  double* myEam_Frhovals;
508 <
509 <  char* ffPath_env = "FORCE_PARAM_PATH";
510 <  char* ffPath;
511 <  char* the_token;
512 <  char* eam_eof_test;
513 <  FILE *eamFile;
514 <  const int BUFFERSIZE = 3000;
515 <
516 <  char temp[200];
517 <  int linenumber;
518 <  int nReadLines;
519 <  char eam_read_buffer[BUFFERSIZE];
520 <
521 <
522 <  int i,j;
523 <
524 <  linenumber = 0;
525 <
526 <  // Open eam file
527 <  eamFile = fopen( eamPotFile, "r" );
365 > void Shapes_FF::initializeBonds( int nBonds, Bond** BondArray,
366 >                                 bond_pair* the_bonds ){
367    
368 <  
530 <  if( eamFile == NULL ){
531 <    
532 <      // next see if the force path enviorment variable is set
533 <    
534 <    ffPath = getenv( ffPath_env );
535 <    if( ffPath == NULL ) {
536 <      STR_DEFINE(ffPath, FRC_PATH );
537 <    }
538 <    
539 <    
540 <    strcpy( temp, ffPath );
541 <    strcat( temp, "/" );
542 <    strcat( temp, eamPotFile );
543 <    strcpy( eamPotFile, temp );
544 <    
545 <    eamFile = fopen( eamPotFile, "r" );
546 <
547 <    
548 <    
549 <    if( eamFile == NULL ){
550 <      
551 <      sprintf( painCave.errMsg,
552 <               "Error opening the EAM force parameter file: %s\n"
553 <               "Have you tried setting the FORCE_PARAM_PATH environment "
554 <               "variable?\n",
555 <               eamPotFile );
556 <      painCave.isFatal = 1;
557 <      simError();
558 <    }
559 <  }
560 <
561 <  // First line is a comment line, read and toss it....
562 <  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
563 <  linenumber++;
564 <  if(eam_eof_test == NULL){
368 >  if( nBonds ){
369      sprintf( painCave.errMsg,
370 <             "error in reading commment in %s\n", eamPotFile);
370 >             "Shapes_FF does not support bonds.\n" );
371      painCave.isFatal = 1;
372      simError();
373    }
374 + }
375  
376 <
377 <
378 <  // The Second line contains atomic number, atomic mass and a lattice constant
379 <  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
575 <  linenumber++;
576 <  if(eam_eof_test == NULL){
376 > void Shapes_FF::initializeBends( int nBends, Bend** bendArray,
377 >                                 bend_set* the_bends ){
378 >  
379 >  if( nBends ){
380      sprintf( painCave.errMsg,
381 <             "error in reading Identifier line in %s\n", eamPotFile);
381 >             "Shapes_FF does not support bends.\n" );
382      painCave.isFatal = 1;
383      simError();
384    }
385 + }
386  
387 <
388 <
585 <    
586 <  if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
587 <    sprintf( painCave.errMsg,
588 <             "Error parsing EAM ident  line in %s\n", eamPotFile );
589 <    painCave.isFatal = 1;
590 <    simError();
591 <  }
387 > void Shapes_FF::initializeTorsions( int nTorsions, Torsion** torsionArray,
388 >                                    torsion_set* the_torsions ){
389    
390 <  info.eam_ident = atoi( the_token );
594 <
595 <  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
596 <    sprintf( painCave.errMsg,
597 <             "Error parsing EAM mass in %s\n", eamPotFile );
598 <    painCave.isFatal = 1;
599 <    simError();
600 <  }
601 <  info.mass = atof( the_token);
602 <
603 <  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
604 <    sprintf( painCave.errMsg,
605 <             "Error parsing EAM Lattice Constant %s\n", eamPotFile );
606 <    painCave.isFatal = 1;
607 <    simError();
608 <  }
609 <  info.lattice_constant = atof( the_token);
610 <
611 <  // Next line is nrho, drho, nr, dr and rcut
612 <  eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
613 <  if(eam_eof_test == NULL){
390 >  if( nTorsions ){
391      sprintf( painCave.errMsg,
392 <             "error in reading number of points line in %s\n", eamPotFile);
392 >             "Shapes_FF does not support torsions.\n" );
393      painCave.isFatal = 1;
394      simError();
395    }
396 + }
397  
398 <  if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
399 <    sprintf( painCave.errMsg,
400 <             "Error parseing EAM nrho: line in %s\n", eamPotFile );
401 <    painCave.isFatal = 1;
402 <    simError();
403 <  }
398 > int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAtomType* st){
399 >  const int MAXLEN = 1024;
400 >  char inLine[MAXLEN];  
401 >  char *token;
402 >  char *delim = " ,;\t\n";
403 >  Mat3x3d momInert;
404 >  RealSphericalHarmonic tempHarmonic;
405 >  vector<RealSphericalHarmonic> functionVector;
406 >
407 >  ifstrstream shapeFile;
408    
409 <  info.eam_nrho = atoi( the_token );
409 >  // first grab the values in the ShapeInfo section
410 >  findBegin(shapeFile, "ShapeInfo");
411    
412 <  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
413 <    sprintf( painCave.errMsg,
414 <             "Error parsing EAM drho in %s\n", eamPotFile );
415 <    painCave.isFatal = 1;
416 <    simError();
417 <  }
418 <  info.eam_drho = atof( the_token);
412 >  shapeFile.getline(inLine, MAXLEN);
413 >  while( !shapeFile.eof() ) {
414 >    // toss comment lines
415 >    if( inLine[0] != '!' && inLine[0] != '#' ){
416 >      // end marks section completion
417 >      if (isEndLine(inLine)) break;
418 >      
419 >      nTokens = count_tokens(inLine, delim);
420 >      if (nTokens != 0) {
421 >        if (nTokens < 5) {
422 >          sprintf( painCave.errMsg,
423 >                   "Not enough data on a ShapeInfo line in file: %s\n"
424 >                   fileName );
425 >          painCave.severity = OOPSE_ERROR;
426 >          painCave.isFatal = 1;
427 >          simError();
428  
429 <  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
430 <    sprintf( painCave.errMsg,
431 <             "Error parsing EAM # r in %s\n", eamPotFile );
432 <    painCave.isFatal = 1;
433 <    simError();
429 >          token = strtok(inLine, delim);
430 >          token = strtok(NULL, delim);
431 >          st->setMass(atof(token));
432 >          token = strtok(NULL, delim);
433 >          momInert(0,0) = atof(token);
434 >          token = strtok(NULL, delim);
435 >          momInert(1,1) = atof(token);
436 >          token = strtok(NULL, delim);
437 >          momInert(2,2) = atof(token);
438 >          st->setI(momInert);
439 >        }
440 >      }
441 >    }
442    }
443 <  info.eam_nr = atoi( the_token);
443 >
444 >  // now grab the contact functions
445 >  findBegin(shapeFile, "ContactFunctions");
446 >  functionVector.clear();
447    
448 <  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
449 <    sprintf( painCave.errMsg,
450 <             "Error parsing EAM dr in %s\n", eamPotFile );
451 <    painCave.isFatal = 1;
452 <    simError();
448 >  shapeFile.getline(inLine, MAXLEN);
449 >  while( !shapeFile.eof() ) {
450 >    // toss comment lines
451 >    if( inLine[0] != '!' && inLine[0] != '#' ){
452 >      // end marks section completion
453 >      if (isEndLine(inLine)) break;
454 >      
455 >      nTokens = count_tokens(inLine, delim);
456 >      if (nTokens != 0) {
457 >        if (nTokens < 4) {
458 >          sprintf( painCave.errMsg,
459 >                   "Not enough data on a ContactFunctions line in file: %s\n"
460 >                   fileName );
461 >          painCave.severity = OOPSE_ERROR;
462 >          painCave.isFatal = 1;
463 >          simError();
464 >          
465 >          // read in a spherical harmonic function
466 >          token = strtok(inLine, delim);
467 >          tempHarmonic.setL(atoi(token));
468 >          token = strtok(NULL, delim);
469 >          tempHarmonic.setM(atoi(token));
470 >          token = strtok(NULL, delim);
471 >          if (!strcasecmp("sin",token))
472 >            tempHarmonic.makeSinFunction();
473 >          else
474 >            tempHarmonic.makeCosFunction();            
475 >          token = strtok(NULL, delim);
476 >          tempHarmonic.setCoefficient(atof(token));
477 >          
478 >          functionVector.push_back(tempHarmonic);
479 >        }
480 >      }
481 >    }
482    }
651 –  info.eam_dr = atof( the_token);
483  
484 <  if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
485 <    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);
484 >  // pass contact functions to ShapeType
485 >  st->setContactFuncs(functionVector);
486  
487 <
488 <
489 <
664 <
665 <  // Ok now we have to allocate point arrays and read in number of points
666 <  // Index the arrays for fortran, starting at 1
667 <  myEam_Frhovals = new double[info.eam_nrho];
668 <  myEam_rvals    = new double[info.eam_nr];
669 <  myEam_rhovals  = new double[info.eam_nr];
670 <
671 <  // Parse F of rho vals.
672 <
673 <  // Assume for now that we have a complete number of lines
674 <  nReadLines = int(info.eam_nrho/5);
487 >  // now grab the range functions
488 >  findBegin(shapeFile, "RangeFunctions");
489 >  functionVector.clear();
490    
491 <
492 <
493 <  for (i=0;i<nReadLines;i++){
494 <    j = i*5;
495 <
496 <    // Read next line
497 <    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
498 <    linenumber++;
499 <    if(eam_eof_test == NULL){
500 <      sprintf( painCave.errMsg,
501 <               "error in reading EAM file %s at line %d\n",
502 <               eamPotFile,linenumber);
503 <      painCave.isFatal = 1;
504 <      simError();
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 >      
498 >      nTokens = count_tokens(inLine, delim);
499 >      if (nTokens != 0) {
500 >        if (nTokens < 4) {
501 >          sprintf( painCave.errMsg,
502 >                   "Not enough data on a RangeFunctions line in file: %s\n"
503 >                   fileName );
504 >          painCave.severity = OOPSE_ERROR;
505 >          painCave.isFatal = 1;
506 >          simError();
507 >          
508 >          // read in a spherical harmonic function
509 >          token = strtok(inLine, delim);
510 >          tempHarmonic.setL(atoi(token));
511 >          token = strtok(NULL, delim);
512 >          tempHarmonic.setM(atoi(token));
513 >          token = strtok(NULL, delim);
514 >          if (!strcasecmp("sin",token))
515 >            tempHarmonic.makeSinFunction();
516 >          else
517 >            tempHarmonic.makeCosFunction();            
518 >          token = strtok(NULL, delim);
519 >          tempHarmonic.setCoefficient(atof(token));
520 >          
521 >          functionVector.push_back(tempHarmonic);
522 >        }
523 >      }
524      }
691 –    
692 –    // Parse 5 values on each line into array
693 –    // Value 1
694 –    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
695 –      sprintf( painCave.errMsg,
696 –               "Error parsing EAM nrho: line in %s\n", eamPotFile );
697 –      painCave.isFatal = 1;
698 –      simError();
699 –    }
700 –
701 –    myEam_Frhovals[j+0] = atof( the_token );
702 –    
703 –    // Value 2
704 –    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
705 –      sprintf( painCave.errMsg,
706 –               "Error parsing EAM nrho: line in %s\n", eamPotFile );
707 –      painCave.isFatal = 1;
708 –      simError();
709 –    }
710 –
711 –    myEam_Frhovals[j+1] = atof( the_token );
712 –    
713 –    // Value 3
714 –    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
715 –      sprintf( painCave.errMsg,
716 –               "Error parsing EAM nrho: line in %s\n", eamPotFile );
717 –      painCave.isFatal = 1;
718 –      simError();
719 –    }
720 –
721 –    myEam_Frhovals[j+2] = atof( the_token );
722 –    
723 –    // Value 4
724 –    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
725 –      sprintf( painCave.errMsg,
726 –               "Error parsing EAM nrho: line in %s\n", eamPotFile );
727 –      painCave.isFatal = 1;
728 –      simError();
729 –    }
730 –
731 –    myEam_Frhovals[j+3] = atof( the_token );
732 –
733 –    // Value 5
734 –    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
735 –      sprintf( painCave.errMsg,
736 –               "Error parsing EAM nrho: line in %s\n", eamPotFile );
737 –      painCave.isFatal = 1;
738 –      simError();
739 –    }
740 –
741 –    myEam_Frhovals[j+4] = atof( the_token );
742 –    
525    }
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);
526  
527 <  for (i=0;i<nReadLines;i++){
528 <    j = i*5;
527 >  // pass range functions to ShapeType
528 >  st->setRangeFuncs(functionVector);
529  
530 <    // Read next line
531 <    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
532 <    linenumber++;
755 <    if(eam_eof_test == NULL){
756 <      sprintf( painCave.errMsg,
757 <               "error in reading EAM file %s at line %d\n",
758 <               eamPotFile,linenumber);
759 <      painCave.isFatal = 1;
760 <      simError();
761 <    }
762 <    
763 <    // Parse 5 values on each line into array
764 <    // Value 1
765 <    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
766 <      sprintf( painCave.errMsg,
767 <               "Error parsing EAM nrho: line in %s\n", eamPotFile );
768 <      painCave.isFatal = 1;
769 <      simError();
770 <    }
771 <    
772 <    myEam_rvals[j+0] = atof( the_token );
773 <
774 <    // Value 2
775 <    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
776 <      sprintf( painCave.errMsg,
777 <               "Error parsing EAM nrho: line in %s\n", eamPotFile );
778 <      painCave.isFatal = 1;
779 <      simError();
780 <    }
530 >  // finally grab the strength functions
531 >  findBegin(shapeFile, "StrengthFunctions");
532 >  functionVector.clear();
533    
534 <    myEam_rvals[j+1] = atof( the_token );
535 <
536 <    // Value 3
537 <    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
538 <      sprintf( painCave.errMsg,
539 <               "Error parsing EAM nrho: line in %s\n", eamPotFile );
540 <      painCave.isFatal = 1;
541 <      simError();
534 >  shapeFile.getline(inLine, MAXLEN);
535 >  while( !shapeFile.eof() ) {
536 >    // toss comment lines
537 >    if( inLine[0] != '!' && inLine[0] != '#' ){
538 >      // end marks section completion
539 >      if (isEndLine(inLine)) break;
540 >      
541 >      nTokens = count_tokens(inLine, delim);
542 >      if (nTokens != 0) {
543 >        if (nTokens < 4) {
544 >          sprintf( painCave.errMsg,
545 >                   "Not enough data on a StrengthFunctions line in file: %s\n"
546 >                   fileName );
547 >          painCave.severity = OOPSE_ERROR;
548 >          painCave.isFatal = 1;
549 >          simError();
550 >          
551 >          // read in a spherical harmonic function
552 >          token = strtok(inLine, delim);
553 >          tempHarmonic.setL(atoi(token));
554 >          token = strtok(NULL, delim);
555 >          tempHarmonic.setM(atoi(token));
556 >          token = strtok(NULL, delim);
557 >          if (!strcasecmp("sin",token))
558 >            tempHarmonic.makeSinFunction();
559 >          else
560 >            tempHarmonic.makeCosFunction();            
561 >          token = strtok(NULL, delim);
562 >          tempHarmonic.setCoefficient(atof(token));
563 >          
564 >          functionVector.push_back(tempHarmonic);
565 >        }
566 >      }
567      }
791 –  
792 –    myEam_rvals[j+2] = atof( the_token );
793 –
794 –    // Value 4
795 –    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
796 –      sprintf( painCave.errMsg,
797 –               "Error parsing EAM nrho: line in %s\n", eamPotFile );
798 –      painCave.isFatal = 1;
799 –      simError();
800 –    }
801 –  
802 –    myEam_rvals[j+3] = atof( the_token );
803 –
804 –    // Value 5
805 –    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
806 –      sprintf( painCave.errMsg,
807 –               "Error parsing EAM nrho: line in %s\n", eamPotFile );
808 –      painCave.isFatal = 1;
809 –      simError();
810 –    }
811 –  
812 –    myEam_rvals[j+4] = atof( the_token );
813 –
568    }
815 –  // Parse rho of r vals
569  
570 <  // Assume for now that we have a complete number of lines
570 >  // pass strength functions to ShapeType
571 >  st->setStrengthFuncs(functionVector);
572  
573 <  for (i=0;i<nReadLines;i++){
574 <    j = i*5;
821 <
822 <    // Read next line
823 <    eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile);
824 <    linenumber++;
825 <    if(eam_eof_test == NULL){
826 <      sprintf( painCave.errMsg,
827 <               "error in reading EAM file %s at line %d\n",
828 <               eamPotFile,linenumber);
829 <      painCave.isFatal = 1;
830 <      simError();
831 <    }
832 <  
833 <    // Parse 5 values on each line into array
834 <    // Value 1
835 <    if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){
836 <      sprintf( painCave.errMsg,
837 <               "Error parsing EAM nrho: line in %s\n", eamPotFile );
838 <      painCave.isFatal = 1;
839 <      simError();
840 <    }
841 <  
842 <    myEam_rhovals[j+0] = atof( the_token );
843 <
844 <    // Value 2
845 <    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
846 <      sprintf( painCave.errMsg,
847 <               "Error parsing EAM nrho: line in %s\n", eamPotFile );
848 <      painCave.isFatal = 1;
849 <      simError();
850 <    }
851 <  
852 <    myEam_rhovals[j+1] = atof( the_token );
853 <
854 <    // Value 3
855 <    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
856 <      sprintf( painCave.errMsg,
857 <               "Error parsing EAM nrho: line in %s\n", eamPotFile );
858 <      painCave.isFatal = 1;
859 <      simError();
860 <    }
861 <  
862 <    myEam_rhovals[j+2] = atof( the_token );
863 <
864 <    // Value 4
865 <    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
866 <      sprintf( painCave.errMsg,
867 <               "Error parsing EAM nrho: line in %s\n", eamPotFile );
868 <      painCave.isFatal = 1;
869 <      simError();
870 <    }
871 <  
872 <    myEam_rhovals[j+3] = atof( the_token );
873 <
874 <    // Value 5
875 <    if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){
876 <      sprintf( painCave.errMsg,
877 <               "Error parsing EAM nrho: line in %s\n", eamPotFile );
878 <      painCave.isFatal = 1;
879 <      simError();
880 <    }
881 <  
882 <    myEam_rhovals[j+4] = atof( the_token );
883 <
884 <  }
885 <  *eam_rvals = myEam_rvals;
886 <  *eam_rhovals = myEam_rhovals;
887 <  *eam_Frhovals = myEam_Frhovals;
888 <
889 <  fclose(eamFile);
573 >  // we're done reading from this file
574 >  shapeFile.close();
575    return 0;
576   }
577 +

Diff Legend

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