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 <iostream>

using namespace std;
using namespace oopse;

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

#include "UseTheForce/ForceFields.hpp"
#include "primitives/SRI.hpp"
#include "utils/simError.h"
#include "io/basic_ifstrstream.hpp"
#include "math/RealSphericalHarmonic.hpp"
#include "math/SquareMatrix3.hpp"

#include "UseTheForce/DarkSide/atype_interface.h"
#include "UseTheForce/DarkSide/shapes_interface.h"

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

Shapes_FF::Shapes_FF() {
  Shapes_FF("");
}

Shapes_FF::Shapes_FF(char* the_variant){
  ffPath_env = "FORCE_PARAM_PATH";

}

Shapes_FF::~Shapes_FF(){

#ifdef IS_MPI
  if( worldRank == 0 ){
#endif // is_mpi
    
    fclose( frcFile );
    
#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::readForceFile( void ){

  char readLine[1024];
  char fileName[200];
  char temp[200];
  char* ffPath;
  char *shapeFileName;
  char *nameToken;
  char *delim = " ,;\t\n";
  int nTokens, i;
  int nContact = 0;
  int nRange = 0;
  int nStrength = 0;
  int myATID;
  string nameString;
  ShapeType* st;
  map<string, AtomType*> atomTypeMap;
  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;

  // build a basic file reader
  ifstrsteam frcReader;
   
  // generate the force file name   
  strcpy( fileName, "Shapes.frc" );
  
  // attempt to open the file in the current directory first.
  frcReader.open( fileName );
  
  if( frcReader == NULL ){
    
    // next see if the force path enviorment variable is set
    
    ffPath = getenv( ffPath_env );
    if( ffPath == NULL ) {
      STR_DEFINE(ffPath, FRC_PATH );
    }
        
    strcpy( temp, ffPath );
    strcat( temp, "/" );
    strcat( temp, fileName );
    strcpy( fileName, temp );
    
    frcReader.open( fileName );
    
    if( frcFile == 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 );
      painCave.severity = OOPSE_ERROR;
      painCave.isFatal = 1;
      simError();
    }
  }
  
  // read in the shape types.
  
  findBegin( "ShapeTypes" );
  
  while( !frcReader.eof() ){
    frcReader.getline( readLine, sizeof(readLine) );

    // toss comment lines
    if( readLine[0] != '!' && readLine[0] != '#' ){
      
      if (isEndLine(readLine)) break;
      
      nTokens = count_tokens(readLine, " ,;\t");
      if (nTokens != 0) {
        if (nTokens < 2) {
          sprintf( painCave.errMsg,
                   "Not enough data on a ShapeTypes line in file: %s\n"
                   fileName );
          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 ShapeType();

          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 = (ShapeType*)(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, &contactM, &contactFunc, 
                   &contactCoeff,
                   &nRange, &rangeL, &rangeM, &rangeFunc, &rangeCoeff,
                   &nStrength, &strengthL, &strengthM,
                   &strengthFunc, &strengthCoeff,
                   &myATID, 
                   &isError);
        if( isError ){
          sprintf( painCave.errMsg,
                   "Error initializing the \"%s\" shape in fortran\n",
                   marker->first );
          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::initializeAtoms( int nAtoms, Atom** the_atoms ){
  
  int i,j,k;

  // initialize the atoms
  DirectionalAtom* dAtom;
  AtomType* at;
  DirectionalAtomType* dat;
  double mySigma;
  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->isLennardJones() ) {
        mySigma = at->properties->getPropertyByName("sigma");
        if( bigSigma < mySigma ) bigSigma = mySigma;
      }

      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();
  }
}

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

  ifstrstream shapeFile;
  
  // 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 = count_tokens(inLine, delim);
      if (nTokens != 0) {
        if (nTokens < 5) {
          sprintf( painCave.errMsg,
                   "Not enough data on a ShapeInfo line in file: %s\n"
                   fileName );
          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);
        }
      }
    }
  }

  // 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 = count_tokens(inLine, delim);
      if (nTokens != 0) {
        if (nTokens < 4) {
          sprintf( painCave.errMsg,
                   "Not enough data on a ContactFunctions line in file: %s\n"
                   fileName );
          painCave.severity = OOPSE_ERROR;
          painCave.isFatal = 1;
          simError();
          
          // read in a spherical harmonic function
          token = strtok(inLine, delim);
          tempHarmonic.setL(atoi(token));
          token = strtok(NULL, delim);
          tempHarmonic.setM(atoi(token));
          token = strtok(NULL, delim);
          if (!strcasecmp("sin",token))
            tempHarmonic.makeSinFunction();
          else
            tempHarmonic.makeCosFunction();             
          token = strtok(NULL, delim);
          tempHarmonic.setCoefficient(atof(token));
          
          functionVector.push_back(tempHarmonic);
        }
      }
    }
  }

  // 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 = count_tokens(inLine, delim);
      if (nTokens != 0) {
        if (nTokens < 4) {
          sprintf( painCave.errMsg,
                   "Not enough data on a RangeFunctions line in file: %s\n"
                   fileName );
          painCave.severity = OOPSE_ERROR;
          painCave.isFatal = 1;
          simError();
          
          // read in a spherical harmonic function
          token = strtok(inLine, delim);
          tempHarmonic.setL(atoi(token));
          token = strtok(NULL, delim);
          tempHarmonic.setM(atoi(token));
          token = strtok(NULL, delim);
          if (!strcasecmp("sin",token))
            tempHarmonic.makeSinFunction();
          else
            tempHarmonic.makeCosFunction();             
          token = strtok(NULL, delim);
          tempHarmonic.setCoefficient(atof(token));
          
          functionVector.push_back(tempHarmonic);
        }
      }
    }
  }

  // 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 = count_tokens(inLine, delim);
      if (nTokens != 0) {
        if (nTokens < 4) {
          sprintf( painCave.errMsg,
                   "Not enough data on a StrengthFunctions line in file: %s\n"
                   fileName );
          painCave.severity = OOPSE_ERROR;
          painCave.isFatal = 1;
          simError();
          
          // read in a spherical harmonic function
          token = strtok(inLine, delim);
          tempHarmonic.setL(atoi(token));
          token = strtok(NULL, delim);
          tempHarmonic.setM(atoi(token));
          token = strtok(NULL, delim);
          if (!strcasecmp("sin",token))
            tempHarmonic.makeSinFunction();
          else
            tempHarmonic.makeCosFunction();             
          token = strtok(NULL, delim);
          tempHarmonic.setCoefficient(atof(token));
          
          functionVector.push_back(tempHarmonic);
        }
      }
    }
  }

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

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

Revision:	1646
Committed:	Tue Oct 26 18:02:46 2004 UTC (19 years, 8 months ago) by chrisfen
Original Path:	*trunk/OOPSE-4/src/UseTheForce/Shapes_FF.cpp*
File size:	13714 byte(s)
Log Message:	Changes to Shapes force field reader
#	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			#include <iostream>
15	gezelter	1520
16			using namespace std;
17	chrisfen	1646	using namespace oopse;
18	gezelter	1520
19			#ifdef IS_MPI
20			#include <mpi.h>
21			#endif //is_mpi
22
23			#include "UseTheForce/ForceFields.hpp"
24			#include "primitives/SRI.hpp"
25			#include "utils/simError.h"
26	chrisfen	1646	#include "io/basic_ifstrstream.hpp"
27			#include "math/RealSphericalHarmonic.hpp"
28			#include "math/SquareMatrix3.hpp"
29	gezelter	1520
30	chrisfen	1646	#include "UseTheForce/DarkSide/atype_interface.h"
31	gezelter	1613	#include "UseTheForce/DarkSide/shapes_interface.h"
32	gezelter	1520
33			#ifdef IS_MPI
34			#include "UseTheForce/mpiForceField.h"
35			#endif // is_mpi
36
37			Shapes_FF::Shapes_FF() {
38			Shapes_FF("");
39			}
40
41			Shapes_FF::Shapes_FF(char* the_variant){
42	chrisfen	1646	ffPath_env = "FORCE_PARAM_PATH";
43	gezelter	1520
44			}
45
46			Shapes_FF::~Shapes_FF(){
47
48			#ifdef IS_MPI
49			if( worldRank == 0 ){
50			#endif // is_mpi
51
52			fclose( frcFile );
53
54			#ifdef IS_MPI
55			}
56			#endif // is_mpi
57			}
58
59
60			void Shapes_FF::calcRcut( void ){
61
62			#ifdef IS_MPI
63			double tempShapesRcut = shapesRcut;
64			MPI_Allreduce( &tempShapesRcut, &shapesRcut, 1, MPI_DOUBLE, MPI_MAX,
65			MPI_COMM_WORLD);
66			#endif //is_mpi
67			entry_plug->setDefaultRcut(shapesRcut);
68			}
69
70
71	gezelter	1628	void Shapes_FF::initForceField(){
72			initFortran(0);
73	gezelter	1520	}
74
75
76	chrisfen	1646	void Shapes_FF::readForceFile( void ){
77	gezelter	1520
78	chrisfen	1646	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	gezelter	1520
95	chrisfen	1646	// 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	gezelter	1520
110	chrisfen	1646	// 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	gezelter	1520
121	chrisfen	1646	// next see if the force path enviorment variable is set
122	gezelter	1520
123	chrisfen	1646	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	gezelter	1520
133	chrisfen	1646	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	gezelter	1520	painCave.isFatal = 1;
145			simError();
146			}
147	chrisfen	1646	}
148
149			// read in the shape types.
150
151			findBegin( "ShapeTypes" );
152
153			while( !frcReader.eof() ){
154			frcReader.getline( readLine, sizeof(readLine) );
155	gezelter	1520
156	chrisfen	1646	// 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	gezelter	1520
172	chrisfen	1646	nameToken = strtok( readLine, delim );
173			shapeFileName = strtok( NULL, delim );
174	gezelter	1520
175	chrisfen	1646	// strings are not char arrays!
176			nameString = nameToken;
177	gezelter	1520
178	chrisfen	1646	// 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	gezelter	1520
185	chrisfen	1646	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	gezelter	1520
197	chrisfen	1646	st = (ShapeType*)(iter->second);
198			parseShapeFile(shapeFileName, st);
199			}
200			}
201	gezelter	1520	}
202			}
203
204	chrisfen	1646	#ifdef IS_MPI
205	gezelter	1520
206	chrisfen	1646	// looks like all the processors have their ShapeType vectors ready...
207			sprintf( checkPointMsg,
208			"Shapes_FF shape objects read successfully." );
209			MPIcheckPoint();
210	gezelter	1520
211	chrisfen	1646	#endif // is_mpi
212	gezelter	1520
213	chrisfen	1646	// pack up and send the necessary info to fortran
214			for (iter = atomTypeMap.begin(); iter != atomTypeMap.end(); ++iter){
215	gezelter	1520
216	chrisfen	1646	at = (AtomType*)(iter.second);
217	gezelter	1520
218	chrisfen	1646	if (at->isDirectional()) {
219	gezelter	1520
220	chrisfen	1646	dat = (DirectionalAtomType*)at;
221	gezelter	1520
222	chrisfen	1646	if (dat->isShape()) {
223	gezelter	1520
224	chrisfen	1646	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	gezelter	1520	}
289			}
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
297			}
298
299	chrisfen	1646	void SHAPES_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
300
301			int i,j,k;
302	gezelter	1520
303			// initialize the atoms
304	chrisfen	1646	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	gezelter	1520	for( i=0; i<nAtoms; i++ ){
314
315	chrisfen	1646	myTypeString = the_atoms[i]->getType();
316			iter = atomTypeMap.find(myTypeString);
317
318			if (iter == atomTypeMap.end()) {
319	gezelter	1520	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	chrisfen	1646	} else {
325	gezelter	1520
326	chrisfen	1646	at = (AtomType*)(iter->second);
327	gezelter	1520
328	chrisfen	1646	the_atoms[i]->setMass( at->getMass() );
329			the_atoms[i]->setIdent( at->getIdent() );
330	gezelter	1520
331	chrisfen	1646	if( at->isLennardJones() ) {
332			mySigma = at->properties->getPropertyByName("sigma");
333			if( bigSigma < mySigma ) bigSigma = mySigma;
334			}
335	gezelter	1520
336	chrisfen	1646	the_atoms[i]->setHasCharge(at->isCharge());
337	gezelter	1520
338	chrisfen	1646	if( at->isDirectional() ){
339	gezelter	1520
340	chrisfen	1646	dat = (DirectionalAtomType*)at;
341			dAtom = (DirectionalAtom *) the_atoms[i];
342	gezelter	1520
343	chrisfen	1646	momInt = dat->getI();
344	gezelter	1520
345	chrisfen	1646	// 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	gezelter	1520
351	chrisfen	1646	ji[0] = 0.0;
352			ji[1] = 0.0;
353			ji[2] = 0.0;
354			dAtom->setJ( ji );
355	gezelter	1520
356	chrisfen	1646	if (dat->isDipole()) {
357			dAtom->setHasDipole( dat->isDipole() );
358			entry_plug->n_dipoles++;
359			}
360			}
361	gezelter	1520	}
362			}
363			}
364
365	chrisfen	1646	void Shapes_FF::initializeBonds( int nBonds, Bond** BondArray,
366			bond_pair* the_bonds ){
367	gezelter	1520
368	chrisfen	1646	if( nBonds ){
369	gezelter	1520	sprintf( painCave.errMsg,
370	chrisfen	1646	"Shapes_FF does not support bonds.\n" );
371	gezelter	1520	painCave.isFatal = 1;
372			simError();
373			}
374	chrisfen	1646	}
375	gezelter	1520
376	chrisfen	1646	void Shapes_FF::initializeBends( int nBends, Bend** bendArray,
377			bend_set* the_bends ){
378
379			if( nBends ){
380	gezelter	1520	sprintf( painCave.errMsg,
381	chrisfen	1646	"Shapes_FF does not support bends.\n" );
382	gezelter	1520	painCave.isFatal = 1;
383			simError();
384			}
385	chrisfen	1646	}
386	gezelter	1520
387	chrisfen	1646	void Shapes_FF::initializeTorsions( int nTorsions, Torsion** torsionArray,
388			torsion_set* the_torsions ){
389	gezelter	1520
390	chrisfen	1646	if( nTorsions ){
391	gezelter	1520	sprintf( painCave.errMsg,
392	chrisfen	1646	"Shapes_FF does not support torsions.\n" );
393	gezelter	1520	painCave.isFatal = 1;
394			simError();
395			}
396	chrisfen	1646	}
397	gezelter	1520
398	chrisfen	1646	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	gezelter	1520
407	chrisfen	1646	ifstrstream shapeFile;
408	gezelter	1520
409	chrisfen	1646	// first grab the values in the ShapeInfo section
410			findBegin(shapeFile, "ShapeInfo");
411	gezelter	1520
412	chrisfen	1646	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	gezelter	1520
429	chrisfen	1646	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	gezelter	1520	}
442			}
443
444	chrisfen	1646	// now grab the contact functions
445			findBegin(shapeFile, "ContactFunctions");
446			functionVector.clear();
447	gezelter	1520
448	chrisfen	1646	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	gezelter	1520	}
482			}
483
484	chrisfen	1646	// pass contact functions to ShapeType
485			st->setContactFuncs(functionVector);
486	gezelter	1520
487	chrisfen	1646	// now grab the range functions
488			findBegin(shapeFile, "RangeFunctions");
489			functionVector.clear();
490	gezelter	1520
491	chrisfen	1646	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	gezelter	1520	}
525	chrisfen	1646	}
526	gezelter	1520
527	chrisfen	1646	// pass range functions to ShapeType
528			st->setRangeFuncs(functionVector);
529	gezelter	1520
530	chrisfen	1646	// finally grab the strength functions
531			findBegin(shapeFile, "StrengthFunctions");
532			functionVector.clear();
533	gezelter	1520
534	chrisfen	1646	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	gezelter	1520	}
568			}
569
570	chrisfen	1646	// pass strength functions to ShapeType
571			st->setStrengthFuncs(functionVector);
572
573			// we're done reading from this file
574			shapeFile.close();
575	gezelter	1520	return 0;
576			}
577	chrisfen	1646