src/UseTheForce/Shapes_FF.cpp

/*
 *  Shapes_FF.cpp
 *  oopse
 *
 *  Created by Chris Fennell on 10/20/04.
 *  Copyright 2004 __MyCompanyName__. All rights reserved.
 *
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <map>
#include <cmath>
#include <iostream>

#ifdef IS_MPI
#include <mpi.h>
#endif //is_mpi

#include "UseTheForce/ForceFields.hpp"
#include "primitives/SRI.hpp"
#include "utils/simError.h"
#include "utils/StringUtils.hpp"
#include "io/basic_ifstrstream.hpp"
#include "math/RealSphericalHarmonic.hpp"
#include "math/SquareMatrix3.hpp"
#include "types/ShapeAtomType.hpp"
#include "UseTheForce/DarkSide/atype_interface.h"
#include "UseTheForce/DarkSide/shapes_interface.h"

#ifdef IS_MPI
#include "UseTheForce/mpiForceField.h"
#endif // is_mpi

using namespace std;
using namespace oopse;

Shapes_FF::~Shapes_FF(){

#ifdef IS_MPI
  if( worldRank == 0 ){
#endif // is_mpi
    
    forceFile.close();
    
#ifdef IS_MPI
  }
#endif // is_mpi
}


void Shapes_FF::calcRcut( void ){
  
#ifdef IS_MPI
  double tempShapesRcut = shapesRcut;
  MPI_Allreduce( &tempShapesRcut, &shapesRcut, 1, MPI_DOUBLE, MPI_MAX,
                 MPI_COMM_WORLD);
#endif  //is_mpi
  entry_plug->setDefaultRcut(shapesRcut);
}


void Shapes_FF::initForceField(){
  initFortran(0);
}


void Shapes_FF::readParams( void ){

  char readLine[1024];

  string fileName;
  string shapeFileName;
  string tempString;

  char *nameToken;
  char *delim = " ,;\t\n";
  int nTokens, i;
  int nContact = 0;
  int nRange = 0;
  int nStrength = 0;
  int myATID;
  int isError;
  string nameString;
  AtomType* at;
  DirectionalAtomType* dat;
  ShapeAtomType* st;

  map<string, AtomType*>::iterator iter;

  // vectors for shape transfer to fortran
  vector<RealSphericalHarmonic*> tempSHVector;
  vector<int> contactL;
  vector<int> contactM;
  vector<int> contactFunc;
  vector<double> contactCoeff;
  vector<int> rangeL;
  vector<int> rangeM;
  vector<int> rangeFunc;
  vector<double> rangeCoeff;
  vector<int> strengthL;
  vector<int> strengthM;
  vector<int> strengthFunc;
  vector<double> strengthCoeff;
   
  // generate the force file name   
  fileName = "Shapes.frc";
  
  // attempt to open the file in the current directory first.
  forceFile.open( fileName.c_str() );
  
  if( forceFile == NULL ){
    
    tempString = ffPath;
    tempString += "/";
    tempString += fileName;
    fileName = tempString;
    
    forceFile.open( fileName.c_str() );
    
    if( forceFile == NULL ){
      
      sprintf( painCave.errMsg,
               "Error opening the force field parameter file:\n"
               "\t%s\n"
               "\tHave you tried setting the FORCE_PARAM_PATH environment "
               "variable?\n",
               fileName.c_str() );
      painCave.severity = OOPSE_ERROR;
      painCave.isFatal = 1;
      simError();
    }
  }
  
  // read in the shape types.
  
  findBegin( forceFile, "ShapeTypes" );
  
  while( !forceFile.eof() ){
    forceFile.getline( readLine, sizeof(readLine) );

    // toss comment lines
    if( readLine[0] != '!' && readLine[0] != '#' ){
      
      if (isEndLine(readLine)) break;
      
      nTokens = countTokens(readLine, " ,;\t");
      if (nTokens != 0) {
        if (nTokens < 2) {
          sprintf( painCave.errMsg,
                   "Not enough data on a ShapeTypes line in file: %s\n",
                   fileName.c_str() );
          painCave.severity = OOPSE_ERROR;
          painCave.isFatal = 1;
          simError();
        }
        
        nameToken = strtok( readLine, delim );
        shapeFileName = strtok( NULL, delim );

        // strings are not char arrays!
        nameString = nameToken;

        // does this AtomType name already exist in the map?
        iter = atomTypeMap.find(nameString);    
        if (iter == atomTypeMap.end()) {
          // no, it doesn't, so we may proceed:
          
          st = new ShapeAtomType();

          st->setName(nameString);
          myATID = atomTypeMap.size();
          st->setIdent(myATID);
          parseShapeFile(shapeFileName, st);
          st->complete();
          atomTypeMap.insert(make_pair(nameString, st));
          
        } else {
          // atomType map already contained this string (i.e. it was
          // declared in a previous block, and we just need to add
          // the shape-specific information for this AtomType:

          st = (ShapeAtomType*)(iter->second);
          parseShapeFile(shapeFileName, st);
        }
      }
    }
  }

#ifdef IS_MPI

  // looks like all the processors have their ShapeType vectors ready...
  sprintf( checkPointMsg,
           "Shapes_FF shape objects read successfully." );
  MPIcheckPoint();

#endif // is_mpi

  // pack up and send the necessary info to fortran
  for (iter = atomTypeMap.begin(); iter != atomTypeMap.end(); ++iter){

    at = (AtomType*)(iter->second);

    if (at->isDirectional()) {

      dat = (DirectionalAtomType*)at;

      if (dat->isShape()) {

        st = (ShapeAtomType*)at;
        
        contactL.clear();
        contactM.clear();
        contactFunc.clear();
        contactCoeff.clear();
        
        tempSHVector = st->getContactFuncs();
        
        nContact = tempSHVector.size();
        for (i=0; i<nContact; i++){
          contactL.push_back(tempSHVector[i]->getL());
          contactM.push_back(tempSHVector[i]->getM());
          contactFunc.push_back(tempSHVector[i]->getFunctionType());
          contactCoeff.push_back(tempSHVector[i]->getCoefficient());
        }
        
        rangeL.clear();
        rangeM.clear();
        rangeFunc.clear();
        rangeCoeff.clear();
        
        tempSHVector = st->getRangeFuncs();
        
        nRange = tempSHVector.size();
        for (i=0; i<nRange; i++){
          rangeL.push_back(tempSHVector[i]->getL());
          rangeM.push_back(tempSHVector[i]->getM());
          rangeFunc.push_back(tempSHVector[i]->getFunctionType());
          rangeCoeff.push_back(tempSHVector[i]->getCoefficient());
        }
        
        strengthL.clear();
        strengthM.clear();
        strengthFunc.clear();
        strengthCoeff.clear();
        
        tempSHVector = st->getStrengthFuncs();
        
        nStrength = tempSHVector.size();
        for (i=0; i<nStrength; i++){
          strengthL.push_back(tempSHVector[i]->getL());
          strengthM.push_back(tempSHVector[i]->getM());
          strengthFunc.push_back(tempSHVector[i]->getFunctionType());
          strengthCoeff.push_back(tempSHVector[i]->getCoefficient());
        }
        
        isError = 0;
        myATID = at->getIdent();
        makeShape( &nContact, &contactL[0], &contactM[0], &contactFunc[0], 
                   &contactCoeff[0],
                   &nRange, &rangeL[0], &rangeM[0], &rangeFunc[0], 
                   &rangeCoeff[0],
                   &nStrength, &strengthL[0], &strengthM[0],
                   &strengthFunc[0], &strengthCoeff[0],
                   &myATID, 
                   &isError);
        if( isError ){
          sprintf( painCave.errMsg,
                   "Error initializing the \"%s\" shape in fortran\n",
                   (iter->first).c_str() );
          painCave.isFatal = 1;
          simError();
        }
      }
    }
  }
  
#ifdef IS_MPI
  sprintf( checkPointMsg,
           "Shapes_FF atom structures successfully sent to fortran\n" );
  MPIcheckPoint();
#endif // is_mpi

}

void Shapes_FF::cleanMe( void ){

}

void Shapes_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
  
  int i,j,k;
  map<string, AtomType*>::iterator iter;

  // initialize the atoms
  DirectionalAtom* dAtom;
  AtomType* at;
  DirectionalAtomType* dat;
  ShapeAtomType* sat;
  double sigma;
  double ji[3];
  double inertialMat[3][3];
  Mat3x3d momInt;
  string myTypeString;

  for( i=0; i<nAtoms; i++ ){
    
    myTypeString = the_atoms[i]->getType();
    iter = atomTypeMap.find(myTypeString);

    if (iter == atomTypeMap.end()) {
      sprintf( painCave.errMsg, 
               "AtomType error, %s not found in force file.\n",
               the_atoms[i]->getType() );
      painCave.isFatal = 1;
      simError();
    } else {

      at = (AtomType*)(iter->second);

      the_atoms[i]->setMass( at->getMass() );
      the_atoms[i]->setIdent( at->getIdent() );
      
      if ( at->isShape() ) {
        
        sat = (ShapeAtomType*)at;
        sigma = findLargestContactDistance(sat);
        if (sigma > bigSigma) bigSigma = sigma;

      }

      the_atoms[i]->setHasCharge(at->isCharge());

      if( at->isDirectional() ){

        dat = (DirectionalAtomType*)at;
        dAtom = (DirectionalAtom *) the_atoms[i];
        
        momInt = dat->getI();

        // zero out the moments of inertia matrix
        for( j=0; j<3; j++ )
          for( k=0; k<3; k++ )
            inertialMat[j][k] = momInt(j,k);
        dAtom->setI( inertialMat );             

        ji[0] = 0.0;
        ji[1] = 0.0;
        ji[2] = 0.0;
        dAtom->setJ( ji );

        if (dat->isDipole()) {
          dAtom->setHasDipole( dat->isDipole() );
          entry_plug->n_dipoles++;
        }               
      }
    }
  }
}

void Shapes_FF::initializeBonds( int nBonds, Bond** BondArray,
                                 bond_pair* the_bonds ){
  
  if( nBonds ){
    sprintf( painCave.errMsg,
             "Shapes_FF does not support bonds.\n" );
    painCave.isFatal = 1;
    simError();
  }
}

void Shapes_FF::initializeBends( int nBends, Bend** bendArray,
                                 bend_set* the_bends ){
  
  if( nBends ){
    sprintf( painCave.errMsg,
             "Shapes_FF does not support bends.\n" );
    painCave.isFatal = 1;
    simError();
  }
}

void Shapes_FF::initializeTorsions( int nTorsions, Torsion** torsionArray,
                                    torsion_set* the_torsions ){
  
  if( nTorsions ){
    sprintf( painCave.errMsg,
             "Shapes_FF does not support torsions.\n" );
    painCave.isFatal = 1;
    simError();
  }
}

void Shapes_FF::parseShapeFile(string shapeFileName, ShapeAtomType* st){
  const int MAXLEN = 1024;
  char inLine[MAXLEN];   
  char *token;
  char *delim = " ,;\t\n";
  int nTokens;
  Mat3x3d momInert;
  RealSphericalHarmonic* rsh;
  vector<RealSphericalHarmonic*> functionVector;
  ifstrstream shapeFile;
  string tempString;

  shapeFile.open( shapeFileName.c_str() );
  
  if( shapeFile == NULL ){
    
    tempString = ffPath;
    tempString += "/";
    tempString += shapeFileName;
    shapeFileName = tempString;
        
    shapeFile.open( shapeFileName.c_str() );
    
    if( shapeFile == NULL ){
      
      sprintf( painCave.errMsg,
               "Error opening the shape file:\n"
               "\t%s\n"
               "\tHave you tried setting the FORCE_PARAM_PATH environment "
               "variable?\n",
               shapeFileName.c_str() );
      painCave.severity = OOPSE_ERROR;
      painCave.isFatal = 1;
      simError();
    }
  }
  
  // read in the shape types. 
  
  // first grab the values in the ShapeInfo section
  findBegin( shapeFile, "ShapeInfo");
  
  shapeFile.getline(inLine, MAXLEN);
  while( !shapeFile.eof() ) {
    // toss comment lines
    if( inLine[0] != '!' && inLine[0] != '#' ){
      // end marks section completion
      if (isEndLine(inLine)) break;
      
      nTokens = countTokens(inLine, delim);
      if (nTokens != 0) {
        if (nTokens < 5) {
          sprintf( painCave.errMsg,
                   "Not enough data on a ShapeInfo line in file: %s\n",
                   shapeFileName.c_str() );
          painCave.severity = OOPSE_ERROR;
          painCave.isFatal = 1;
          simError();

          token = strtok(inLine, delim);
          token = strtok(NULL, delim);
          st->setMass(atof(token));
          token = strtok(NULL, delim);
          momInert(0,0) = atof(token);
          token = strtok(NULL, delim);
          momInert(1,1) = atof(token);
          token = strtok(NULL, delim);
          momInert(2,2) = atof(token);
          st->setI(momInert);
        }
      }
    }
    shapeFile.getline(inLine, MAXLEN); 
  }

  // now grab the contact functions
  findBegin(shapeFile, "ContactFunctions");
  functionVector.clear();
  
  shapeFile.getline(inLine, MAXLEN);
  while( !shapeFile.eof() ) {
    // toss comment lines
    if( inLine[0] != '!' && inLine[0] != '#' ){
      // end marks section completion
      if (isEndLine(inLine)) break;
      
      nTokens = countTokens(inLine, delim);
      if (nTokens != 0) {
        if (nTokens < 4) {
          sprintf( painCave.errMsg,
                   "Not enough data on a ContactFunctions line in file: %s\n",
                   shapeFileName.c_str() );
          painCave.severity = OOPSE_ERROR;
          painCave.isFatal = 1;
          simError();
          
          // read in a spherical harmonic function
          token = strtok(inLine, delim);
          rsh = new RealSphericalHarmonic();
          rsh->setL(atoi(token));
          token = strtok(NULL, delim);
          rsh->setM(atoi(token));
          token = strtok(NULL, delim);
          if (!strcasecmp("sin",token))
            rsh->makeSinFunction();
          else
            rsh->makeCosFunction();             
          token = strtok(NULL, delim);
          rsh->setCoefficient(atof(token));
          
          functionVector.push_back(rsh);
        }
      }
    }
    shapeFile.getline(inLine, MAXLEN);
  }

  // pass contact functions to ShapeType
  st->setContactFuncs(functionVector);

  // now grab the range functions
  findBegin(shapeFile, "RangeFunctions");
  functionVector.clear();
  
  shapeFile.getline(inLine, MAXLEN);
  while( !shapeFile.eof() ) {
    // toss comment lines
    if( inLine[0] != '!' && inLine[0] != '#' ){
      // end marks section completion
      if (isEndLine(inLine)) break;
      
      nTokens = countTokens(inLine, delim);
      if (nTokens != 0) {
        if (nTokens < 4) {
          sprintf( painCave.errMsg,
                   "Not enough data on a RangeFunctions line in file: %s\n",
                   shapeFileName.c_str() );
          painCave.severity = OOPSE_ERROR;
          painCave.isFatal = 1;
          simError();
          
          // read in a spherical harmonic function
          token = strtok(inLine, delim);

          rsh = new RealSphericalHarmonic();
          rsh->setL(atoi(token));
          token = strtok(NULL, delim);
          rsh->setM(atoi(token));
          token = strtok(NULL, delim);
          if (!strcasecmp("sin",token))
            rsh->makeSinFunction();
          else
            rsh->makeCosFunction();             
          token = strtok(NULL, delim);
          rsh->setCoefficient(atof(token));
          
          functionVector.push_back(rsh);
        }
      }
    }
    shapeFile.getline(inLine, MAXLEN);
  }

  // pass range functions to ShapeType
  st->setRangeFuncs(functionVector);

  // finally grab the strength functions
  findBegin(shapeFile, "StrengthFunctions");
  functionVector.clear();
  
  shapeFile.getline(inLine, MAXLEN);
  while( !shapeFile.eof() ) {
    // toss comment lines
    if( inLine[0] != '!' && inLine[0] != '#' ){
      // end marks section completion
      if (isEndLine(inLine)) break;
      
      nTokens = countTokens(inLine, delim);
      if (nTokens != 0) {
        if (nTokens < 4) {
          sprintf( painCave.errMsg,
                   "Not enough data on a StrengthFunctions line in file: %s\n",
                   shapeFileName.c_str() );
          painCave.severity = OOPSE_ERROR;
          painCave.isFatal = 1;
          simError();
          
          // read in a spherical harmonic function
          token = strtok(inLine, delim);
          rsh = new RealSphericalHarmonic();
          rsh->setL(atoi(token));
          token = strtok(NULL, delim);
          rsh->setM(atoi(token));
          token = strtok(NULL, delim);
          if (!strcasecmp("sin",token))
            rsh->makeSinFunction();
          else
            rsh->makeCosFunction();             
          token = strtok(NULL, delim);
          rsh->setCoefficient(atof(token));
          
          functionVector.push_back(rsh);
        }
      }
    }
    shapeFile.getline(inLine, MAXLEN);
  }

  // pass strength functions to ShapeType
  st->setStrengthFuncs(functionVector);

  // we're done reading from this file
  shapeFile.close();
}

double Shapes_FF::findLargestContactDistance(ShapeAtomType* st) {
  int i, j,  nSteps;
  double theta, thetaStep, thetaMin, costheta;
  double phi, phiStep;
  double sigma, bigSigma;

  nSteps = 16;

  thetaStep = M_PI / nSteps;
  thetaMin = thetaStep / 2.0;
  phiStep = thetaStep * 2.0;
  bigSigma = 0.0;
  
  for (i = 0; i < nSteps; i++) {
    
    theta = thetaMin + i * thetaStep;
    costheta = cos(theta);

    for (j = 0; j < nSteps; j++) {

      phi = j*phiStep;

      sigma = st->getContactValueAt(costheta, phi);
      
      if (sigma > bigSigma) bigSigma = sigma;
    }
  }

  return bigSigma;  
}
Revision:	1672
Committed:	Thu Oct 28 17:20:22 2004 UTC (21 years ago) by gezelter
File size:	15170 byte(s)
Log Message:	forgot to advance lines in the shapeFile
#	User	Rev	Content
1	chrisfen	1646	/*
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	gezelter	1520	#include <stdlib.h>
11			#include <stdio.h>
12			#include <string.h>
13	chrisfen	1646	#include <map>
14	gezelter	1650	#include <cmath>
15	chrisfen	1646	#include <iostream>
16	gezelter	1520
17			#ifdef IS_MPI
18			#include <mpi.h>
19			#endif //is_mpi
20
21			#include "UseTheForce/ForceFields.hpp"
22			#include "primitives/SRI.hpp"
23			#include "utils/simError.h"
24	gezelter	1650	#include "utils/StringUtils.hpp"
25	chrisfen	1646	#include "io/basic_ifstrstream.hpp"
26			#include "math/RealSphericalHarmonic.hpp"
27			#include "math/SquareMatrix3.hpp"
28	gezelter	1650	#include "types/ShapeAtomType.hpp"
29	chrisfen	1646	#include "UseTheForce/DarkSide/atype_interface.h"
30	gezelter	1613	#include "UseTheForce/DarkSide/shapes_interface.h"
31	gezelter	1520
32			#ifdef IS_MPI
33			#include "UseTheForce/mpiForceField.h"
34			#endif // is_mpi
35
36	gezelter	1650	using namespace std;
37			using namespace oopse;
38	gezelter	1520
39			Shapes_FF::~Shapes_FF(){
40
41			#ifdef IS_MPI
42			if( worldRank == 0 ){
43			#endif // is_mpi
44
45	gezelter	1650	forceFile.close();
46	gezelter	1520
47			#ifdef IS_MPI
48			}
49			#endif // is_mpi
50			}
51
52
53			void Shapes_FF::calcRcut( void ){
54
55			#ifdef IS_MPI
56			double tempShapesRcut = shapesRcut;
57			MPI_Allreduce( &tempShapesRcut, &shapesRcut, 1, MPI_DOUBLE, MPI_MAX,
58			MPI_COMM_WORLD);
59			#endif //is_mpi
60			entry_plug->setDefaultRcut(shapesRcut);
61			}
62
63
64	gezelter	1628	void Shapes_FF::initForceField(){
65			initFortran(0);
66	gezelter	1520	}
67
68
69	gezelter	1650	void Shapes_FF::readParams( void ){
70	gezelter	1520
71	chrisfen	1646	char readLine[1024];
72	gezelter	1650
73			string fileName;
74			string shapeFileName;
75			string tempString;
76
77	chrisfen	1646	char *nameToken;
78			char *delim = " ,;\t\n";
79			int nTokens, i;
80			int nContact = 0;
81			int nRange = 0;
82			int nStrength = 0;
83			int myATID;
84	gezelter	1650	int isError;
85	chrisfen	1646	string nameString;
86	gezelter	1650	AtomType* at;
87			DirectionalAtomType* dat;
88			ShapeAtomType* st;
89
90	chrisfen	1646	map<string, AtomType*>::iterator iter;
91	gezelter	1520
92	chrisfen	1646	// vectors for shape transfer to fortran
93	gezelter	1650	vector<RealSphericalHarmonic*> tempSHVector;
94	chrisfen	1646	vector<int> contactL;
95			vector<int> contactM;
96			vector<int> contactFunc;
97			vector<double> contactCoeff;
98			vector<int> rangeL;
99			vector<int> rangeM;
100			vector<int> rangeFunc;
101			vector<double> rangeCoeff;
102			vector<int> strengthL;
103			vector<int> strengthM;
104			vector<int> strengthFunc;
105			vector<double> strengthCoeff;
106
107			// generate the force file name
108	gezelter	1650	fileName = "Shapes.frc";
109	chrisfen	1646
110			// attempt to open the file in the current directory first.
111	gezelter	1650	forceFile.open( fileName.c_str() );
112	chrisfen	1646
113	gezelter	1650	if( forceFile == NULL ){
114	gezelter	1520
115	gezelter	1650	tempString = ffPath;
116			tempString += "/";
117			tempString += fileName;
118			fileName = tempString;
119	gezelter	1520
120	gezelter	1650	forceFile.open( fileName.c_str() );
121	gezelter	1520
122	gezelter	1650	if( forceFile == NULL ){
123	chrisfen	1646
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	gezelter	1650	fileName.c_str() );
130	chrisfen	1646	painCave.severity = OOPSE_ERROR;
131	gezelter	1520	painCave.isFatal = 1;
132			simError();
133			}
134	chrisfen	1646	}
135
136			// read in the shape types.
137
138	gezelter	1650	findBegin( forceFile, "ShapeTypes" );
139	chrisfen	1646
140	gezelter	1650	while( !forceFile.eof() ){
141			forceFile.getline( readLine, sizeof(readLine) );
142	gezelter	1520
143	chrisfen	1646	// toss comment lines
144			if( readLine[0] != '!' && readLine[0] != '#' ){
145
146			if (isEndLine(readLine)) break;
147
148	gezelter	1650	nTokens = countTokens(readLine, " ,;\t");
149	chrisfen	1646	if (nTokens != 0) {
150			if (nTokens < 2) {
151			sprintf( painCave.errMsg,
152	gezelter	1650	"Not enough data on a ShapeTypes line in file: %s\n",
153			fileName.c_str() );
154	chrisfen	1646	painCave.severity = OOPSE_ERROR;
155			painCave.isFatal = 1;
156			simError();
157			}
158	gezelter	1520
159	chrisfen	1646	nameToken = strtok( readLine, delim );
160			shapeFileName = strtok( NULL, delim );
161	gezelter	1520
162	chrisfen	1646	// strings are not char arrays!
163			nameString = nameToken;
164	gezelter	1520
165	chrisfen	1646	// does this AtomType name already exist in the map?
166			iter = atomTypeMap.find(nameString);
167			if (iter == atomTypeMap.end()) {
168			// no, it doesn't, so we may proceed:
169
170	gezelter	1650	st = new ShapeAtomType();
171	gezelter	1520
172	chrisfen	1646	st->setName(nameString);
173			myATID = atomTypeMap.size();
174			st->setIdent(myATID);
175			parseShapeFile(shapeFileName, st);
176			st->complete();
177			atomTypeMap.insert(make_pair(nameString, st));
178
179			} else {
180			// atomType map already contained this string (i.e. it was
181			// declared in a previous block, and we just need to add
182			// the shape-specific information for this AtomType:
183	gezelter	1520
184	gezelter	1650	st = (ShapeAtomType*)(iter->second);
185	chrisfen	1646	parseShapeFile(shapeFileName, st);
186			}
187			}
188	gezelter	1520	}
189			}
190
191	chrisfen	1646	#ifdef IS_MPI
192	gezelter	1520
193	chrisfen	1646	// looks like all the processors have their ShapeType vectors ready...
194			sprintf( checkPointMsg,
195			"Shapes_FF shape objects read successfully." );
196			MPIcheckPoint();
197	gezelter	1520
198	chrisfen	1646	#endif // is_mpi
199	gezelter	1520
200	chrisfen	1646	// pack up and send the necessary info to fortran
201			for (iter = atomTypeMap.begin(); iter != atomTypeMap.end(); ++iter){
202	gezelter	1520
203	gezelter	1650	at = (AtomType*)(iter->second);
204	gezelter	1520
205	chrisfen	1646	if (at->isDirectional()) {
206	gezelter	1520
207	chrisfen	1646	dat = (DirectionalAtomType*)at;
208	gezelter	1520
209	chrisfen	1646	if (dat->isShape()) {
210	gezelter	1520
211	chrisfen	1646	st = (ShapeAtomType*)at;
212
213			contactL.clear();
214			contactM.clear();
215			contactFunc.clear();
216			contactCoeff.clear();
217
218			tempSHVector = st->getContactFuncs();
219
220			nContact = tempSHVector.size();
221			for (i=0; i<nContact; i++){
222	gezelter	1650	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	chrisfen	1646	}
227
228			rangeL.clear();
229			rangeM.clear();
230			rangeFunc.clear();
231			rangeCoeff.clear();
232
233			tempSHVector = st->getRangeFuncs();
234
235			nRange = tempSHVector.size();
236			for (i=0; i<nRange; i++){
237	gezelter	1650	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	chrisfen	1646	}
242
243			strengthL.clear();
244			strengthM.clear();
245			strengthFunc.clear();
246			strengthCoeff.clear();
247
248			tempSHVector = st->getStrengthFuncs();
249
250			nStrength = tempSHVector.size();
251			for (i=0; i<nStrength; i++){
252	gezelter	1650	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	chrisfen	1646	}
257
258			isError = 0;
259			myATID = at->getIdent();
260	gezelter	1650	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	chrisfen	1646	&myATID,
267			&isError);
268			if( isError ){
269			sprintf( painCave.errMsg,
270			"Error initializing the \"%s\" shape in fortran\n",
271	gezelter	1650	(iter->first).c_str() );
272	chrisfen	1646	painCave.isFatal = 1;
273			simError();
274			}
275			}
276	gezelter	1520	}
277			}
278
279			#ifdef IS_MPI
280			sprintf( checkPointMsg,
281			"Shapes_FF atom structures successfully sent to fortran\n" );
282			MPIcheckPoint();
283			#endif // is_mpi
284
285			}
286
287	gezelter	1650	void Shapes_FF::cleanMe( void ){
288
289			}
290
291			void Shapes_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
292	chrisfen	1646
293			int i,j,k;
294	gezelter	1650	map<string, AtomType*>::iterator iter;
295	gezelter	1520
296			// initialize the atoms
297	chrisfen	1646	DirectionalAtom* dAtom;
298			AtomType* at;
299			DirectionalAtomType* dat;
300	gezelter	1650	ShapeAtomType* sat;
301			double sigma;
302	chrisfen	1646	double ji[3];
303			double inertialMat[3][3];
304			Mat3x3d momInt;
305			string myTypeString;
306
307	gezelter	1520	for( i=0; i<nAtoms; i++ ){
308
309	chrisfen	1646	myTypeString = the_atoms[i]->getType();
310			iter = atomTypeMap.find(myTypeString);
311
312			if (iter == atomTypeMap.end()) {
313	gezelter	1520	sprintf( painCave.errMsg,
314			"AtomType error, %s not found in force file.\n",
315			the_atoms[i]->getType() );
316			painCave.isFatal = 1;
317			simError();
318	chrisfen	1646	} else {
319	gezelter	1520
320	chrisfen	1646	at = (AtomType*)(iter->second);
321	gezelter	1520
322	chrisfen	1646	the_atoms[i]->setMass( at->getMass() );
323			the_atoms[i]->setIdent( at->getIdent() );
324	gezelter	1650
325			if ( at->isShape() ) {
326
327			sat = (ShapeAtomType*)at;
328			sigma = findLargestContactDistance(sat);
329			if (sigma > bigSigma) bigSigma = sigma;
330	gezelter	1520
331	chrisfen	1646	}
332	gezelter	1520
333	chrisfen	1646	the_atoms[i]->setHasCharge(at->isCharge());
334	gezelter	1520
335	chrisfen	1646	if( at->isDirectional() ){
336	gezelter	1520
337	chrisfen	1646	dat = (DirectionalAtomType*)at;
338			dAtom = (DirectionalAtom *) the_atoms[i];
339	gezelter	1672
340	chrisfen	1646	momInt = dat->getI();
341	gezelter	1520
342	chrisfen	1646	// zero out the moments of inertia matrix
343			for( j=0; j<3; j++ )
344			for( k=0; k<3; k++ )
345			inertialMat[j][k] = momInt(j,k);
346			dAtom->setI( inertialMat );
347	gezelter	1520
348	chrisfen	1646	ji[0] = 0.0;
349			ji[1] = 0.0;
350			ji[2] = 0.0;
351			dAtom->setJ( ji );
352	gezelter	1520
353	chrisfen	1646	if (dat->isDipole()) {
354			dAtom->setHasDipole( dat->isDipole() );
355			entry_plug->n_dipoles++;
356			}
357			}
358	gezelter	1520	}
359			}
360			}
361
362	chrisfen	1646	void Shapes_FF::initializeBonds( int nBonds, Bond** BondArray,
363			bond_pair* the_bonds ){
364	gezelter	1520
365	chrisfen	1646	if( nBonds ){
366	gezelter	1520	sprintf( painCave.errMsg,
367	chrisfen	1646	"Shapes_FF does not support bonds.\n" );
368	gezelter	1520	painCave.isFatal = 1;
369			simError();
370			}
371	chrisfen	1646	}
372	gezelter	1520
373	chrisfen	1646	void Shapes_FF::initializeBends( int nBends, Bend** bendArray,
374			bend_set* the_bends ){
375
376			if( nBends ){
377	gezelter	1520	sprintf( painCave.errMsg,
378	chrisfen	1646	"Shapes_FF does not support bends.\n" );
379	gezelter	1520	painCave.isFatal = 1;
380			simError();
381			}
382	chrisfen	1646	}
383	gezelter	1520
384	chrisfen	1646	void Shapes_FF::initializeTorsions( int nTorsions, Torsion** torsionArray,
385			torsion_set* the_torsions ){
386	gezelter	1520
387	chrisfen	1646	if( nTorsions ){
388	gezelter	1520	sprintf( painCave.errMsg,
389	chrisfen	1646	"Shapes_FF does not support torsions.\n" );
390	gezelter	1520	painCave.isFatal = 1;
391			simError();
392			}
393	chrisfen	1646	}
394	gezelter	1520
395	gezelter	1650	void Shapes_FF::parseShapeFile(string shapeFileName, ShapeAtomType* st){
396	chrisfen	1646	const int MAXLEN = 1024;
397			char inLine[MAXLEN];
398			char *token;
399			char *delim = " ,;\t\n";
400	gezelter	1650	int nTokens;
401	chrisfen	1646	Mat3x3d momInert;
402	gezelter	1650	RealSphericalHarmonic* rsh;
403			vector<RealSphericalHarmonic*> functionVector;
404			ifstrstream shapeFile;
405			string tempString;
406	gezelter	1520
407	gezelter	1650	shapeFile.open( shapeFileName.c_str() );
408	gezelter	1520
409	gezelter	1650	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	chrisfen	1646	// first grab the values in the ShapeInfo section
435	gezelter	1650	findBegin( shapeFile, "ShapeInfo");
436	gezelter	1520
437	chrisfen	1646	shapeFile.getline(inLine, MAXLEN);
438			while( !shapeFile.eof() ) {
439			// toss comment lines
440			if( inLine[0] != '!' && inLine[0] != '#' ){
441			// end marks section completion
442			if (isEndLine(inLine)) break;
443
444	gezelter	1650	nTokens = countTokens(inLine, delim);
445	chrisfen	1646	if (nTokens != 0) {
446			if (nTokens < 5) {
447			sprintf( painCave.errMsg,
448	gezelter	1650	"Not enough data on a ShapeInfo line in file: %s\n",
449			shapeFileName.c_str() );
450	chrisfen	1646	painCave.severity = OOPSE_ERROR;
451			painCave.isFatal = 1;
452			simError();
453	gezelter	1520
454	chrisfen	1646	token = strtok(inLine, delim);
455			token = strtok(NULL, delim);
456			st->setMass(atof(token));
457			token = strtok(NULL, delim);
458			momInert(0,0) = atof(token);
459			token = strtok(NULL, delim);
460			momInert(1,1) = atof(token);
461			token = strtok(NULL, delim);
462			momInert(2,2) = atof(token);
463			st->setI(momInert);
464			}
465			}
466	gezelter	1520	}
467	gezelter	1672	shapeFile.getline(inLine, MAXLEN);
468	gezelter	1520	}
469
470	chrisfen	1646	// now grab the contact functions
471			findBegin(shapeFile, "ContactFunctions");
472			functionVector.clear();
473	gezelter	1520
474	chrisfen	1646	shapeFile.getline(inLine, MAXLEN);
475			while( !shapeFile.eof() ) {
476			// toss comment lines
477			if( inLine[0] != '!' && inLine[0] != '#' ){
478			// end marks section completion
479			if (isEndLine(inLine)) break;
480
481	gezelter	1650	nTokens = countTokens(inLine, delim);
482	chrisfen	1646	if (nTokens != 0) {
483			if (nTokens < 4) {
484			sprintf( painCave.errMsg,
485	gezelter	1650	"Not enough data on a ContactFunctions line in file: %s\n",
486			shapeFileName.c_str() );
487	chrisfen	1646	painCave.severity = OOPSE_ERROR;
488			painCave.isFatal = 1;
489			simError();
490
491			// read in a spherical harmonic function
492			token = strtok(inLine, delim);
493	gezelter	1650	rsh = new RealSphericalHarmonic();
494			rsh->setL(atoi(token));
495	chrisfen	1646	token = strtok(NULL, delim);
496	gezelter	1650	rsh->setM(atoi(token));
497	chrisfen	1646	token = strtok(NULL, delim);
498			if (!strcasecmp("sin",token))
499	gezelter	1650	rsh->makeSinFunction();
500	chrisfen	1646	else
501	gezelter	1650	rsh->makeCosFunction();
502	chrisfen	1646	token = strtok(NULL, delim);
503	gezelter	1650	rsh->setCoefficient(atof(token));
504	chrisfen	1646
505	gezelter	1650	functionVector.push_back(rsh);
506	chrisfen	1646	}
507			}
508	gezelter	1520	}
509	gezelter	1672	shapeFile.getline(inLine, MAXLEN);
510	gezelter	1520	}
511
512	chrisfen	1646	// pass contact functions to ShapeType
513			st->setContactFuncs(functionVector);
514	gezelter	1520
515	chrisfen	1646	// now grab the range functions
516			findBegin(shapeFile, "RangeFunctions");
517			functionVector.clear();
518	gezelter	1520
519	chrisfen	1646	shapeFile.getline(inLine, MAXLEN);
520			while( !shapeFile.eof() ) {
521			// toss comment lines
522			if( inLine[0] != '!' && inLine[0] != '#' ){
523			// end marks section completion
524			if (isEndLine(inLine)) break;
525
526	gezelter	1650	nTokens = countTokens(inLine, delim);
527	chrisfen	1646	if (nTokens != 0) {
528			if (nTokens < 4) {
529			sprintf( painCave.errMsg,
530	gezelter	1650	"Not enough data on a RangeFunctions line in file: %s\n",
531			shapeFileName.c_str() );
532	chrisfen	1646	painCave.severity = OOPSE_ERROR;
533			painCave.isFatal = 1;
534			simError();
535
536			// read in a spherical harmonic function
537			token = strtok(inLine, delim);
538	gezelter	1650
539			rsh = new RealSphericalHarmonic();
540			rsh->setL(atoi(token));
541	chrisfen	1646	token = strtok(NULL, delim);
542	gezelter	1650	rsh->setM(atoi(token));
543	chrisfen	1646	token = strtok(NULL, delim);
544			if (!strcasecmp("sin",token))
545	gezelter	1650	rsh->makeSinFunction();
546	chrisfen	1646	else
547	gezelter	1650	rsh->makeCosFunction();
548	chrisfen	1646	token = strtok(NULL, delim);
549	gezelter	1650	rsh->setCoefficient(atof(token));
550	chrisfen	1646
551	gezelter	1650	functionVector.push_back(rsh);
552	chrisfen	1646	}
553			}
554	gezelter	1520	}
555	gezelter	1672	shapeFile.getline(inLine, MAXLEN);
556	chrisfen	1646	}
557	gezelter	1520
558	chrisfen	1646	// pass range functions to ShapeType
559			st->setRangeFuncs(functionVector);
560	gezelter	1520
561	chrisfen	1646	// finally grab the strength functions
562			findBegin(shapeFile, "StrengthFunctions");
563			functionVector.clear();
564	gezelter	1520
565	chrisfen	1646	shapeFile.getline(inLine, MAXLEN);
566			while( !shapeFile.eof() ) {
567			// toss comment lines
568			if( inLine[0] != '!' && inLine[0] != '#' ){
569			// end marks section completion
570			if (isEndLine(inLine)) break;
571
572	gezelter	1650	nTokens = countTokens(inLine, delim);
573	chrisfen	1646	if (nTokens != 0) {
574			if (nTokens < 4) {
575			sprintf( painCave.errMsg,
576	gezelter	1650	"Not enough data on a StrengthFunctions line in file: %s\n",
577			shapeFileName.c_str() );
578	chrisfen	1646	painCave.severity = OOPSE_ERROR;
579			painCave.isFatal = 1;
580			simError();
581
582			// read in a spherical harmonic function
583			token = strtok(inLine, delim);
584	gezelter	1650	rsh = new RealSphericalHarmonic();
585			rsh->setL(atoi(token));
586	chrisfen	1646	token = strtok(NULL, delim);
587	gezelter	1650	rsh->setM(atoi(token));
588	chrisfen	1646	token = strtok(NULL, delim);
589			if (!strcasecmp("sin",token))
590	gezelter	1650	rsh->makeSinFunction();
591	chrisfen	1646	else
592	gezelter	1650	rsh->makeCosFunction();
593	chrisfen	1646	token = strtok(NULL, delim);
594	gezelter	1650	rsh->setCoefficient(atof(token));
595	chrisfen	1646
596	gezelter	1650	functionVector.push_back(rsh);
597	chrisfen	1646	}
598			}
599	gezelter	1520	}
600	gezelter	1672	shapeFile.getline(inLine, MAXLEN);
601	gezelter	1520	}
602
603	chrisfen	1646	// pass strength functions to ShapeType
604			st->setStrengthFuncs(functionVector);
605
606			// we're done reading from this file
607			shapeFile.close();
608	gezelter	1520	}
609	chrisfen	1646
610	gezelter	1650	double Shapes_FF::findLargestContactDistance(ShapeAtomType* st) {
611			int i, j, nSteps;
612			double theta, thetaStep, thetaMin, costheta;
613			double phi, phiStep;
614			double sigma, bigSigma;
615
616			nSteps = 16;
617
618			thetaStep = M_PI / nSteps;
619			thetaMin = thetaStep / 2.0;
620			phiStep = thetaStep * 2.0;
621			bigSigma = 0.0;
622
623			for (i = 0; i < nSteps; i++) {
624
625			theta = thetaMin + i * thetaStep;
626			costheta = cos(theta);
627
628			for (j = 0; j < nSteps; j++) {
629
630			phi = j*phiStep;
631
632			sigma = st->getContactValueAt(costheta, phi);
633
634			if (sigma > bigSigma) bigSigma = sigma;
635			}
636			}
637
638			return bigSigma;
639			}