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 1650 by gezelter, Tue Oct 26 22:24:52 2004 UTC

# Line 11 | Line 11
11   #include <stdio.h>
12   #include <string.h>
13   #include <map>
14 + #include <cmath>
15   #include <iostream>
16  
16 – using namespace std;
17 – using namespace oopse;
18 –
17   #ifdef IS_MPI
18   #include <mpi.h>
19   #endif //is_mpi
# Line 23 | Line 21 | using namespace oopse;
21   #include "UseTheForce/ForceFields.hpp"
22   #include "primitives/SRI.hpp"
23   #include "utils/simError.h"
24 + #include "utils/StringUtils.hpp"
25   #include "io/basic_ifstrstream.hpp"
26   #include "math/RealSphericalHarmonic.hpp"
27   #include "math/SquareMatrix3.hpp"
28 <
28 > #include "types/ShapeAtomType.hpp"
29   #include "UseTheForce/DarkSide/atype_interface.h"
30   #include "UseTheForce/DarkSide/shapes_interface.h"
31  
# Line 34 | Line 33 | Shapes_FF::Shapes_FF() {
33   #include "UseTheForce/mpiForceField.h"
34   #endif // is_mpi
35  
36 < Shapes_FF::Shapes_FF() {
37 <  Shapes_FF("");
39 < }
36 > using namespace std;
37 > using namespace oopse;
38  
41 – Shapes_FF::Shapes_FF(char* the_variant){
42 –  ffPath_env = "FORCE_PARAM_PATH";
43 –
44 – }
45 –
39   Shapes_FF::~Shapes_FF(){
40  
41   #ifdef IS_MPI
42    if( worldRank == 0 ){
43   #endif // is_mpi
44      
45 <    fclose( frcFile );
45 >    forceFile.close();
46      
47   #ifdef IS_MPI
48    }
# Line 73 | Line 66 | void Shapes_FF::readForceFile( void ){
66   }
67  
68  
69 < void Shapes_FF::readForceFile( void ){
69 > void Shapes_FF::readParams( void ){
70  
71    char readLine[1024];
72 <  char fileName[200];
73 <  char temp[200];
74 <  char* ffPath;
75 <  char *shapeFileName;
72 >
73 >  string fileName;
74 >  string shapeFileName;
75 >  string tempString;
76 >
77    char *nameToken;
78    char *delim = " ,;\t\n";
79    int nTokens, i;
# Line 87 | Line 81 | void Shapes_FF::readForceFile( void ){
81    int nRange = 0;
82    int nStrength = 0;
83    int myATID;
84 +  int isError;
85    string nameString;
86 <  ShapeType* st;
87 <  map<string, AtomType*> atomTypeMap;
86 >  AtomType* at;
87 >  DirectionalAtomType* dat;
88 >  ShapeAtomType* st;
89 >
90    map<string, AtomType*>::iterator iter;
91  
92    // vectors for shape transfer to fortran
93 <  vector<RealSphericalHarmonic> tempSHVector;
93 >  vector<RealSphericalHarmonic*> tempSHVector;
94    vector<int> contactL;
95    vector<int> contactM;
96    vector<int> contactFunc;
# Line 106 | Line 103 | void Shapes_FF::readForceFile( void ){
103    vector<int> strengthM;
104    vector<int> strengthFunc;
105    vector<double> strengthCoeff;
109 –
110 –  // build a basic file reader
111 –  ifstrsteam frcReader;
106    
107    // generate the force file name  
108 <  strcpy( fileName, "Shapes.frc" );
108 >  fileName = "Shapes.frc";
109    
110    // attempt to open the file in the current directory first.
111 <  frcReader.open( fileName );
111 >  forceFile.open( fileName.c_str() );
112    
113 <  if( frcReader == NULL ){
113 >  if( forceFile == NULL ){
114      
115 <    // next see if the force path enviorment variable is set
115 >    tempString = ffPath;
116 >    tempString += "/";
117 >    tempString += fileName;
118 >    fileName = tempString;
119      
120 <    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 );
120 >    forceFile.open( fileName.c_str() );
121      
122 <    frcReader.open( fileName );
134 <    
135 <    if( frcFile == NULL ){
122 >    if( forceFile == NULL ){
123        
124        sprintf( painCave.errMsg,
125                 "Error opening the force field parameter file:\n"
126                 "\t%s\n"
127                 "\tHave you tried setting the FORCE_PARAM_PATH environment "
128                 "variable?\n",
129 <               fileName );
129 >               fileName.c_str() );
130        painCave.severity = OOPSE_ERROR;
131        painCave.isFatal = 1;
132        simError();
# Line 148 | Line 135 | void Shapes_FF::readForceFile( void ){
135    
136    // read in the shape types.
137    
138 <  findBegin( "ShapeTypes" );
138 >  findBegin( forceFile, "ShapeTypes" );
139    
140 <  while( !frcReader.eof() ){
141 <    frcReader.getline( readLine, sizeof(readLine) );
140 >  while( !forceFile.eof() ){
141 >    forceFile.getline( readLine, sizeof(readLine) );
142  
143      // toss comment lines
144      if( readLine[0] != '!' && readLine[0] != '#' ){
145        
146        if (isEndLine(readLine)) break;
147        
148 <      nTokens = count_tokens(readLine, " ,;\t");
148 >      nTokens = countTokens(readLine, " ,;\t");
149        if (nTokens != 0) {
150          if (nTokens < 2) {
151            sprintf( painCave.errMsg,
152 <                   "Not enough data on a ShapeTypes line in file: %s\n"
153 <                   fileName );
152 >                   "Not enough data on a ShapeTypes line in file: %s\n",
153 >                   fileName.c_str() );
154            painCave.severity = OOPSE_ERROR;
155            painCave.isFatal = 1;
156            simError();
# Line 180 | Line 167 | void Shapes_FF::readForceFile( void ){
167          if (iter == atomTypeMap.end()) {
168            // no, it doesn't, so we may proceed:
169            
170 <          st = new ShapeType();
170 >          st = new ShapeAtomType();
171  
172            st->setName(nameString);
173            myATID = atomTypeMap.size();
# Line 194 | Line 181 | void Shapes_FF::readForceFile( void ){
181            // declared in a previous block, and we just need to add
182            // the shape-specific information for this AtomType:
183  
184 <          st = (ShapeType*)(iter->second);
184 >          st = (ShapeAtomType*)(iter->second);
185            parseShapeFile(shapeFileName, st);
186          }
187        }
# Line 213 | Line 200 | void Shapes_FF::readForceFile( void ){
200    // pack up and send the necessary info to fortran
201    for (iter = atomTypeMap.begin(); iter != atomTypeMap.end(); ++iter){
202  
203 <    at = (AtomType*)(iter.second);
203 >    at = (AtomType*)(iter->second);
204  
205      if (at->isDirectional()) {
206  
# Line 232 | Line 219 | void Shapes_FF::readForceFile( void ){
219          
220          nContact = tempSHVector.size();
221          for (i=0; i<nContact; i++){
222 <          contactL.push_back(tempSHVector[i].getL());
223 <          contactM.push_back(tempSHVector[i].getM());
224 <          contactFunc.push_back(tempSHVector[i].getFunctionType());
225 <          contactCoeff.push_back(tempSHVector[i].getCoefficient());
222 >          contactL.push_back(tempSHVector[i]->getL());
223 >          contactM.push_back(tempSHVector[i]->getM());
224 >          contactFunc.push_back(tempSHVector[i]->getFunctionType());
225 >          contactCoeff.push_back(tempSHVector[i]->getCoefficient());
226          }
227          
228          rangeL.clear();
# Line 247 | Line 234 | void Shapes_FF::readForceFile( void ){
234          
235          nRange = tempSHVector.size();
236          for (i=0; i<nRange; i++){
237 <          rangeL.push_back(tempSHVector[i].getL());
238 <          rangeM.push_back(tempSHVector[i].getM());
239 <          rangeFunc.push_back(tempSHVector[i].getFunctionType());
240 <          rangeCoeff.push_back(tempSHVector[i].getCoefficient());
237 >          rangeL.push_back(tempSHVector[i]->getL());
238 >          rangeM.push_back(tempSHVector[i]->getM());
239 >          rangeFunc.push_back(tempSHVector[i]->getFunctionType());
240 >          rangeCoeff.push_back(tempSHVector[i]->getCoefficient());
241          }
242          
243          strengthL.clear();
# Line 262 | Line 249 | void Shapes_FF::readForceFile( void ){
249          
250          nStrength = tempSHVector.size();
251          for (i=0; i<nStrength; i++){
252 <          strengthL.push_back(tempSHVector[i].getL());
253 <          strengthM.push_back(tempSHVector[i].getM());
254 <          strengthFunc.push_back(tempSHVector[i].getFunctionType());
255 <          strengthCoeff.push_back(tempSHVector[i].getCoefficient());
252 >          strengthL.push_back(tempSHVector[i]->getL());
253 >          strengthM.push_back(tempSHVector[i]->getM());
254 >          strengthFunc.push_back(tempSHVector[i]->getFunctionType());
255 >          strengthCoeff.push_back(tempSHVector[i]->getCoefficient());
256          }
257          
258          isError = 0;
259          myATID = at->getIdent();
260 <        makeShape( &nContact, &contactL, &contactM, &contactFunc,
261 <                   &contactCoeff,
262 <                   &nRange, &rangeL, &rangeM, &rangeFunc, &rangeCoeff,
263 <                   &nStrength, &strengthL, &strengthM,
264 <                   &strengthFunc, &strengthCoeff,
260 >        makeShape( &nContact, &contactL[0], &contactM[0], &contactFunc[0],
261 >                   &contactCoeff[0],
262 >                   &nRange, &rangeL[0], &rangeM[0], &rangeFunc[0],
263 >                   &rangeCoeff[0],
264 >                   &nStrength, &strengthL[0], &strengthM[0],
265 >                   &strengthFunc[0], &strengthCoeff[0],
266                     &myATID,
267                     &isError);
268          if( isError ){
269            sprintf( painCave.errMsg,
270                     "Error initializing the \"%s\" shape in fortran\n",
271 <                   marker->first );
271 >                   (iter->first).c_str() );
272            painCave.isFatal = 1;
273            simError();
274          }
# Line 296 | Line 284 | void SHAPES_FF::initializeAtoms( int nAtoms, Atom** th
284  
285   }
286  
287 < void SHAPES_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
287 > void Shapes_FF::cleanMe( void ){
288 >
289 > }
290 >
291 > void Shapes_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
292    
293    int i,j,k;
294 +  map<string, AtomType*>::iterator iter;
295  
296    // initialize the atoms
297    DirectionalAtom* dAtom;
298    AtomType* at;
299    DirectionalAtomType* dat;
300 <  double mySigma;
300 >  ShapeAtomType* sat;
301 >  double sigma;
302    double ji[3];
303    double inertialMat[3][3];
304    Mat3x3d momInt;
# Line 327 | Line 321 | void SHAPES_FF::initializeAtoms( int nAtoms, Atom** th
321  
322        the_atoms[i]->setMass( at->getMass() );
323        the_atoms[i]->setIdent( at->getIdent() );
324 +      
325 +      if ( at->isShape() ) {
326 +        
327 +        sat = (ShapeAtomType*)at;
328 +        sigma = findLargestContactDistance(sat);
329 +        if (sigma > bigSigma) bigSigma = sigma;
330  
331 –      if( at->isLennardJones() ) {
332 –        mySigma = at->properties->getPropertyByName("sigma");
333 –        if( bigSigma < mySigma ) bigSigma = mySigma;
331        }
332  
333        the_atoms[i]->setHasCharge(at->isCharge());
# Line 395 | Line 392 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
392    }
393   }
394  
395 < int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAtomType* st){
395 > void Shapes_FF::parseShapeFile(string shapeFileName, ShapeAtomType* st){
396    const int MAXLEN = 1024;
397    char inLine[MAXLEN];  
398    char *token;
399    char *delim = " ,;\t\n";
400 +  int nTokens;
401    Mat3x3d momInert;
402 <  RealSphericalHarmonic tempHarmonic;
403 <  vector<RealSphericalHarmonic> functionVector;
406 <
402 >  RealSphericalHarmonic* rsh;
403 >  vector<RealSphericalHarmonic*> functionVector;
404    ifstrstream shapeFile;
405 +  string tempString;
406 +
407 +  shapeFile.open( shapeFileName.c_str() );
408 +  
409 +  if( shapeFile == NULL ){
410 +    
411 +    tempString = ffPath;
412 +    tempString += "/";
413 +    tempString += shapeFileName;
414 +    shapeFileName = tempString;
415 +        
416 +    shapeFile.open( shapeFileName.c_str() );
417 +    
418 +    if( shapeFile == NULL ){
419 +      
420 +      sprintf( painCave.errMsg,
421 +               "Error opening the shape file:\n"
422 +               "\t%s\n"
423 +               "\tHave you tried setting the FORCE_PARAM_PATH environment "
424 +               "variable?\n",
425 +               shapeFileName.c_str() );
426 +      painCave.severity = OOPSE_ERROR;
427 +      painCave.isFatal = 1;
428 +      simError();
429 +    }
430 +  }
431    
432 +  // read in the shape types.
433 +  
434    // first grab the values in the ShapeInfo section
435 <  findBegin(shapeFile, "ShapeInfo");
435 >  findBegin( shapeFile, "ShapeInfo");
436    
437    shapeFile.getline(inLine, MAXLEN);
438    while( !shapeFile.eof() ) {
# Line 416 | Line 441 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
441        // end marks section completion
442        if (isEndLine(inLine)) break;
443        
444 <      nTokens = count_tokens(inLine, delim);
444 >      nTokens = countTokens(inLine, delim);
445        if (nTokens != 0) {
446          if (nTokens < 5) {
447            sprintf( painCave.errMsg,
448 <                   "Not enough data on a ShapeInfo line in file: %s\n"
449 <                   fileName );
448 >                   "Not enough data on a ShapeInfo line in file: %s\n",
449 >                   shapeFileName.c_str() );
450            painCave.severity = OOPSE_ERROR;
451            painCave.isFatal = 1;
452            simError();
# Line 452 | Line 477 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
477        // end marks section completion
478        if (isEndLine(inLine)) break;
479        
480 <      nTokens = count_tokens(inLine, delim);
480 >      nTokens = countTokens(inLine, delim);
481        if (nTokens != 0) {
482          if (nTokens < 4) {
483            sprintf( painCave.errMsg,
484 <                   "Not enough data on a ContactFunctions line in file: %s\n"
485 <                   fileName );
484 >                   "Not enough data on a ContactFunctions line in file: %s\n",
485 >                   shapeFileName.c_str() );
486            painCave.severity = OOPSE_ERROR;
487            painCave.isFatal = 1;
488            simError();
489            
490            // read in a spherical harmonic function
491            token = strtok(inLine, delim);
492 <          tempHarmonic.setL(atoi(token));
492 >          rsh = new RealSphericalHarmonic();
493 >          rsh->setL(atoi(token));
494            token = strtok(NULL, delim);
495 <          tempHarmonic.setM(atoi(token));
495 >          rsh->setM(atoi(token));
496            token = strtok(NULL, delim);
497            if (!strcasecmp("sin",token))
498 <            tempHarmonic.makeSinFunction();
498 >            rsh->makeSinFunction();
499            else
500 <            tempHarmonic.makeCosFunction();            
500 >            rsh->makeCosFunction();            
501            token = strtok(NULL, delim);
502 <          tempHarmonic.setCoefficient(atof(token));
502 >          rsh->setCoefficient(atof(token));
503            
504 <          functionVector.push_back(tempHarmonic);
504 >          functionVector.push_back(rsh);
505          }
506        }
507      }
# Line 495 | Line 521 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
521        // end marks section completion
522        if (isEndLine(inLine)) break;
523        
524 <      nTokens = count_tokens(inLine, delim);
524 >      nTokens = countTokens(inLine, delim);
525        if (nTokens != 0) {
526          if (nTokens < 4) {
527            sprintf( painCave.errMsg,
528 <                   "Not enough data on a RangeFunctions line in file: %s\n"
529 <                   fileName );
528 >                   "Not enough data on a RangeFunctions line in file: %s\n",
529 >                   shapeFileName.c_str() );
530            painCave.severity = OOPSE_ERROR;
531            painCave.isFatal = 1;
532            simError();
533            
534            // read in a spherical harmonic function
535            token = strtok(inLine, delim);
536 <          tempHarmonic.setL(atoi(token));
536 >
537 >          rsh = new RealSphericalHarmonic();
538 >          rsh->setL(atoi(token));
539            token = strtok(NULL, delim);
540 <          tempHarmonic.setM(atoi(token));
540 >          rsh->setM(atoi(token));
541            token = strtok(NULL, delim);
542            if (!strcasecmp("sin",token))
543 <            tempHarmonic.makeSinFunction();
543 >            rsh->makeSinFunction();
544            else
545 <            tempHarmonic.makeCosFunction();            
545 >            rsh->makeCosFunction();            
546            token = strtok(NULL, delim);
547 <          tempHarmonic.setCoefficient(atof(token));
547 >          rsh->setCoefficient(atof(token));
548            
549 <          functionVector.push_back(tempHarmonic);
549 >          functionVector.push_back(rsh);
550          }
551        }
552      }
# Line 538 | Line 566 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
566        // end marks section completion
567        if (isEndLine(inLine)) break;
568        
569 <      nTokens = count_tokens(inLine, delim);
569 >      nTokens = countTokens(inLine, delim);
570        if (nTokens != 0) {
571          if (nTokens < 4) {
572            sprintf( painCave.errMsg,
573 <                   "Not enough data on a StrengthFunctions line in file: %s\n"
574 <                   fileName );
573 >                   "Not enough data on a StrengthFunctions line in file: %s\n",
574 >                   shapeFileName.c_str() );
575            painCave.severity = OOPSE_ERROR;
576            painCave.isFatal = 1;
577            simError();
578            
579            // read in a spherical harmonic function
580            token = strtok(inLine, delim);
581 <          tempHarmonic.setL(atoi(token));
581 >          rsh = new RealSphericalHarmonic();
582 >          rsh->setL(atoi(token));
583            token = strtok(NULL, delim);
584 <          tempHarmonic.setM(atoi(token));
584 >          rsh->setM(atoi(token));
585            token = strtok(NULL, delim);
586            if (!strcasecmp("sin",token))
587 <            tempHarmonic.makeSinFunction();
587 >            rsh->makeSinFunction();
588            else
589 <            tempHarmonic.makeCosFunction();            
589 >            rsh->makeCosFunction();            
590            token = strtok(NULL, delim);
591 <          tempHarmonic.setCoefficient(atof(token));
591 >          rsh->setCoefficient(atof(token));
592            
593 <          functionVector.push_back(tempHarmonic);
593 >          functionVector.push_back(rsh);
594          }
595        }
596      }
# Line 572 | Line 601 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
601  
602    // we're done reading from this file
603    shapeFile.close();
575 –  return 0;
604   }
605  
606 + double Shapes_FF::findLargestContactDistance(ShapeAtomType* st) {
607 +  int i, j,  nSteps;
608 +  double theta, thetaStep, thetaMin, costheta;
609 +  double phi, phiStep;
610 +  double sigma, bigSigma;
611 +
612 +  nSteps = 16;
613 +
614 +  thetaStep = M_PI / nSteps;
615 +  thetaMin = thetaStep / 2.0;
616 +  phiStep = thetaStep * 2.0;
617 +  bigSigma = 0.0;
618 +  
619 +  for (i = 0; i < nSteps; i++) {
620 +    
621 +    theta = thetaMin + i * thetaStep;
622 +    costheta = cos(theta);
623 +
624 +    for (j = 0; j < nSteps; j++) {
625 +
626 +      phi = j*phiStep;
627 +
628 +      sigma = st->getContactValueAt(costheta, phi);
629 +      
630 +      if (sigma > bigSigma) bigSigma = sigma;
631 +    }
632 +  }
633 +
634 +  return bigSigma;  
635 + }

Diff Legend

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