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 1646 by chrisfen, Tue Oct 26 18:02:46 2004 UTC vs.
Revision 1930 by gezelter, Wed Jan 12 22:41:40 2005 UTC

# Line 1 | Line 1
1 < /*
2 < *  Shapes_FF.cpp
3 < *  oopse
1 > /*
2 > * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4 < *  Created by Chris Fennell on 10/20/04.
5 < *  Copyright 2004 __MyCompanyName__. All rights reserved.
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 <
41 >
42   #include <stdlib.h>
43   #include <stdio.h>
44   #include <string.h>
45   #include <map>
46 + #include <cmath>
47   #include <iostream>
48  
16 – using namespace std;
17 – using namespace oopse;
18 –
49   #ifdef IS_MPI
50   #include <mpi.h>
51   #endif //is_mpi
# Line 23 | Line 53 | using namespace oopse;
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 <
60 > #include "types/ShapeAtomType.hpp"
61   #include "UseTheForce/DarkSide/atype_interface.h"
62   #include "UseTheForce/DarkSide/shapes_interface.h"
63  
# Line 34 | Line 65 | Shapes_FF::Shapes_FF() {
65   #include "UseTheForce/mpiForceField.h"
66   #endif // is_mpi
67  
37 – Shapes_FF::Shapes_FF() {
38 –  Shapes_FF("");
39 – }
68  
69 < Shapes_FF::Shapes_FF(char* the_variant){
42 <  ffPath_env = "FORCE_PARAM_PATH";
69 > using namespace oopse;
70  
44 – }
45 –
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 60 | 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   }
# Line 73 | Line 100 | void Shapes_FF::readForceFile( void ){
100   }
101  
102  
103 < void Shapes_FF::readForceFile( void ){
103 > void Shapes_FF::readParams( void ){
104  
105    char readLine[1024];
106 <  char fileName[200];
107 <  char temp[200];
108 <  char* ffPath;
109 <  char *shapeFileName;
106 >
107 >   std::string fileName;
108 >   std::string shapeFileName;
109 >   std::string tempString;
110 >
111    char *nameToken;
112    char *delim = " ,;\t\n";
113    int nTokens, i;
# Line 87 | Line 115 | void Shapes_FF::readForceFile( void ){
115    int nRange = 0;
116    int nStrength = 0;
117    int myATID;
118 <  string nameString;
119 <  ShapeType* st;
120 <  map<string, AtomType*> atomTypeMap;
121 <  map<string, AtomType*>::iterator iter;
118 >  int isError;
119 >   std::string nameString;
120 >  AtomType* at;
121 >  DirectionalAtomType* dat;
122 >  ShapeAtomType* st;
123  
124 <  // 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;
124 >   std::map<string, AtomType*>::iterator iter;
125  
126 <  // build a basic file reader
127 <  ifstrsteam frcReader;
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 <  strcpy( fileName, "Shapes.frc" );
142 >  fileName = "Shapes.frc";
143    
144    // attempt to open the file in the current directory first.
145 <  frcReader.open( fileName );
146 <  
147 <  if( frcReader == NULL ){
145 >  forceFile.open( fileName.c_str() );
146 >  
147 >  if( forceFile == NULL ){
148      
149 <    // next see if the force path enviorment variable is set
149 >    tempString = ffPath;
150 >    tempString += "/";
151 >    tempString += fileName;
152 >    fileName = tempString;
153      
154 <    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 );
154 >    forceFile.open( fileName.c_str() );
155      
156 <    frcReader.open( fileName );
134 <    
135 <    if( frcFile == NULL ){
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 );
163 >               fileName.c_str() );
164        painCave.severity = OOPSE_ERROR;
165        painCave.isFatal = 1;
166        simError();
# Line 148 | Line 169 | void Shapes_FF::readForceFile( void ){
169    
170    // read in the shape types.
171    
172 <  findBegin( "ShapeTypes" );
172 >  findBegin( forceFile, "ShapeTypes" );
173    
174 <  while( !frcReader.eof() ){
175 <    frcReader.getline( readLine, sizeof(readLine) );
174 >  while( !forceFile.eof() ){
175 >    forceFile.getline( readLine, sizeof(readLine) );
176  
177      // toss comment lines
178      if( readLine[0] != '!' && readLine[0] != '#' ){
179        
180        if (isEndLine(readLine)) break;
181        
182 <      nTokens = count_tokens(readLine, " ,;\t");
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 );
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();
# Line 180 | Line 201 | void Shapes_FF::readForceFile( void ){
201          if (iter == atomTypeMap.end()) {
202            // no, it doesn't, so we may proceed:
203            
204 <          st = new ShapeType();
204 >          st = new ShapeAtomType();
205  
206            st->setName(nameString);
207 <          myATID = atomTypeMap.size();
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 string (i.e. it was
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 = (ShapeType*)(iter->second);
218 >          st = (ShapeAtomType*)(iter->second);
219            parseShapeFile(shapeFileName, st);
220          }
221        }
# Line 213 | Line 234 | void Shapes_FF::readForceFile( void ){
234    // pack up and send the necessary info to fortran
235    for (iter = atomTypeMap.begin(); iter != atomTypeMap.end(); ++iter){
236  
237 <    at = (AtomType*)(iter.second);
237 >    at = (AtomType*)(iter->second);
238  
239      if (at->isDirectional()) {
240  
# Line 232 | Line 253 | void Shapes_FF::readForceFile( void ){
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());
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();
# Line 247 | Line 268 | void Shapes_FF::readForceFile( void ){
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());
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();
# Line 262 | Line 283 | void Shapes_FF::readForceFile( void ){
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());
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 <        makeShape( &nContact, &contactL, &contactM, &contactFunc,
295 <                   &contactCoeff,
296 <                   &nRange, &rangeL, &rangeM, &rangeFunc, &rangeCoeff,
297 <                   &nStrength, &strengthL, &strengthM,
298 <                   &strengthFunc, &strengthCoeff,
294 >
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          if( isError ){
305            sprintf( painCave.errMsg,
306                     "Error initializing the \"%s\" shape in fortran\n",
307 <                   marker->first );
307 >                   (iter->first).c_str() );
308            painCave.isFatal = 1;
309            simError();
310          }
# Line 288 | Line 312 | void Shapes_FF::readForceFile( void ){
312      }
313    }
314    
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 +  
324   #ifdef IS_MPI
325    sprintf( checkPointMsg,
326             "Shapes_FF atom structures successfully sent to fortran\n" );
# Line 296 | Line 329 | void SHAPES_FF::initializeAtoms( int nAtoms, Atom** th
329  
330   }
331  
332 < void SHAPES_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
332 > void Shapes_FF::cleanMe( void ){
333 >
334 > }
335 >
336 > void Shapes_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
337    
338    int i,j,k;
339 +   std::map<string, AtomType*>::iterator iter;
340  
341    // initialize the atoms
342    DirectionalAtom* dAtom;
343    AtomType* at;
344    DirectionalAtomType* dat;
345 <  double mySigma;
345 >  ShapeAtomType* sat;
346 >  double longCutoff;
347    double ji[3];
348    double inertialMat[3][3];
349    Mat3x3d momInt;
350 <  string myTypeString;
350 >   std::string myTypeString;
351  
352 +  bigContact = 0.0;
353 +
354    for( i=0; i<nAtoms; i++ ){
355      
356 <    myTypeString = the_atoms[i]->getType();
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      } else {
# Line 327 | Line 368 | void SHAPES_FF::initializeAtoms( int nAtoms, Atom** th
368  
369        the_atoms[i]->setMass( at->getMass() );
370        the_atoms[i]->setIdent( at->getIdent() );
371 <
372 <      if( at->isLennardJones() ) {
373 <        mySigma = at->properties->getPropertyByName("sigma");
374 <        if( bigSigma < mySigma ) bigSigma = mySigma;
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        the_atoms[i]->setHasCharge(at->isCharge());
# Line 339 | Line 385 | void SHAPES_FF::initializeAtoms( int nAtoms, Atom** th
385  
386          dat = (DirectionalAtomType*)at;
387          dAtom = (DirectionalAtom *) the_atoms[i];
388 <
388 >        
389          momInt = dat->getI();
390  
391          // zero out the moments of inertia matrix
# Line 395 | Line 441 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
441    }
442   }
443  
444 < int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAtomType* st){
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 tempHarmonic;
452 <  vector<RealSphericalHarmonic> functionVector;
406 <
451 >  RealSphericalHarmonic* rsh;
452 >   std::vector<RealSphericalHarmonic*> functionVector;
453    ifstrstream shapeFile;
454 +   std::string tempString;
455 +
456 +  shapeFile.open( shapeFileName.c_str() );
457    
458 +  if( shapeFile == NULL ){
459 +    
460 +    tempString = ffPath;
461 +    tempString += "/";
462 +    tempString += shapeFileName;
463 +    shapeFileName = tempString;
464 +        
465 +    shapeFile.open( shapeFileName.c_str() );
466 +    
467 +    if( shapeFile == NULL ){
468 +      
469 +      sprintf( painCave.errMsg,
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 +               shapeFileName.c_str() );
475 +      painCave.severity = OOPSE_ERROR;
476 +      painCave.isFatal = 1;
477 +      simError();
478 +    }
479 +  }
480 +  
481 +  // read in the shape types.
482 +  
483    // first grab the values in the ShapeInfo section
484 <  findBegin(shapeFile, "ShapeInfo");
484 >  findBegin( shapeFile, "ShapeInfo");
485    
486    shapeFile.getline(inLine, MAXLEN);
487    while( !shapeFile.eof() ) {
# Line 416 | Line 490 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
490        // end marks section completion
491        if (isEndLine(inLine)) break;
492        
493 <      nTokens = count_tokens(inLine, delim);
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 <                   fileName );
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 <
502 >        } else {
503            token = strtok(inLine, delim);
504            token = strtok(NULL, delim);
505            st->setMass(atof(token));
# Line 439 | Line 513 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
513          }
514        }
515      }
516 +    shapeFile.getline(inLine, MAXLEN);
517    }
518  
519    // now grab the contact functions
# Line 451 | Line 526 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
526      if( inLine[0] != '!' && inLine[0] != '#' ){
527        // end marks section completion
528        if (isEndLine(inLine)) break;
529 <      
455 <      nTokens = count_tokens(inLine, delim);
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 <                   fileName );
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 <          
538 >        } else {
539            // read in a spherical harmonic function
540            token = strtok(inLine, delim);
541 <          tempHarmonic.setL(atoi(token));
541 >          rsh = new RealSphericalHarmonic();
542 >          rsh->setL(atoi(token));
543            token = strtok(NULL, delim);
544 <          tempHarmonic.setM(atoi(token));
544 >          rsh->setM(atoi(token));
545            token = strtok(NULL, delim);
546            if (!strcasecmp("sin",token))
547 <            tempHarmonic.makeSinFunction();
547 >            rsh->makeSinFunction();
548            else
549 <            tempHarmonic.makeCosFunction();            
549 >            rsh->makeCosFunction();            
550            token = strtok(NULL, delim);
551 <          tempHarmonic.setCoefficient(atof(token));
551 >          rsh->setCoefficient(atof(token));
552            
553 <          functionVector.push_back(tempHarmonic);
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    // now grab the range functions
# Line 495 | Line 572 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
572        // end marks section completion
573        if (isEndLine(inLine)) break;
574        
575 <      nTokens = count_tokens(inLine, delim);
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 <                   fileName );
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 <          tempHarmonic.setL(atoi(token));
588 >
589 >          rsh = new RealSphericalHarmonic();
590 >          rsh->setL(atoi(token));
591            token = strtok(NULL, delim);
592 <          tempHarmonic.setM(atoi(token));
592 >          rsh->setM(atoi(token));
593            token = strtok(NULL, delim);
594            if (!strcasecmp("sin",token))
595 <            tempHarmonic.makeSinFunction();
595 >            rsh->makeSinFunction();
596            else
597 <            tempHarmonic.makeCosFunction();            
597 >            rsh->makeCosFunction();            
598            token = strtok(NULL, delim);
599 <          tempHarmonic.setCoefficient(atof(token));
599 >          rsh->setCoefficient(atof(token));
600            
601 <          functionVector.push_back(tempHarmonic);
601 >          functionVector.push_back(rsh);
602          }
603        }
604      }
605 +    shapeFile.getline(inLine, MAXLEN);
606    }
607  
608    // pass range functions to ShapeType
# Line 538 | Line 619 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
619        // end marks section completion
620        if (isEndLine(inLine)) break;
621        
622 <      nTokens = count_tokens(inLine, delim);
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 <                   fileName );
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 <          tempHarmonic.setL(atoi(token));
635 >          rsh = new RealSphericalHarmonic();
636 >          rsh->setL(atoi(token));
637            token = strtok(NULL, delim);
638 <          tempHarmonic.setM(atoi(token));
638 >          rsh->setM(atoi(token));
639            token = strtok(NULL, delim);
640            if (!strcasecmp("sin",token))
641 <            tempHarmonic.makeSinFunction();
641 >            rsh->makeSinFunction();
642            else
643 <            tempHarmonic.makeCosFunction();            
643 >            rsh->makeCosFunction();            
644            token = strtok(NULL, delim);
645 <          tempHarmonic.setCoefficient(atof(token));
645 >          rsh->setCoefficient(atof(token));
646            
647 <          functionVector.push_back(tempHarmonic);
647 >          functionVector.push_back(rsh);
648          }
649        }
650      }
651 +    shapeFile.getline(inLine, MAXLEN);
652    }
653  
654    // pass strength functions to ShapeType
# Line 572 | Line 656 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
656  
657    // we're done reading from this file
658    shapeFile.close();
575 –  return 0;
659   }
660  
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 +  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 +    theta = thetaMin + i * thetaStep;
677 +    costheta = cos(theta);
678 +
679 +    for (j = 0; j < nSteps; j++) {
680 +
681 +      phi = j*phiStep;
682 +
683 +      sigma = st->getContactValueAt(costheta, phi);
684 +      
685 +      if (sigma > bs) bs = sigma;
686 +    }
687 +  }
688 +
689 +  return bs;  
690 + }
691 +
692 +
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 +  nSteps = 16;
701 +
702 +  thetaStep = M_PI / nSteps;
703 +  thetaMin = thetaStep / 2.0;
704 +  phiStep = thetaStep * 2.0;
705 +  bigCut = 0.0;
706 +  
707 +  for (i = 0; i < nSteps; i++) {
708 +    
709 +    theta = thetaMin + i * thetaStep;
710 +    costheta = cos(theta);
711 +
712 +    for (j = 0; j < nSteps; j++) {
713 +
714 +      phi = j*phiStep;
715 +
716 +      sigma = st->getContactValueAt(costheta, phi);
717 +      range = st->getRangeValueAt(costheta, phi);
718 +
719 +       // cutoff for a shape is taken to be (2.5*rangeVal + contactVal)
720 +      tempCut = 1.5*range + sigma;
721 +
722 +      if (tempCut > bigCut) bigCut = tempCut;
723 +    }
724 +  }
725 +
726 +  return bigCut;  
727 + }

Diff Legend

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