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 = bigContact;
  MPI_Allreduce( &tempShapesRcut, &shapesRcut, 1, MPI_DOUBLE, MPI_MAX,
                 MPI_COMM_WORLD);
#else
  shapesRcut = bigContact;
#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() + 1;
          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();
        }
      }
    }
  }
  
  isError = 0;
  completeShapeFF(&isError);
  if( isError ){
    sprintf( painCave.errMsg,
             "Error completing Shape FF in fortran\n");
    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 longCutoff;
  double ji[3];
  double inertialMat[3][3];
  Mat3x3d momInt;
  string myTypeString;

  bigContact = 0.0;

  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;
        longCutoff = findCutoffDistance(sat);
        if (longCutoff > bigContact) bigContact = longCutoff;  
        cout << bigContact << " is the cutoff value\n";
        
        entry_plug->useShapes = 1;
      }

      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();
        } else {
          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();
        } else {
          // 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();
        } else {
          
          // 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();
        } else {
          
          // 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, bs;
  
  nSteps = 16;

  thetaStep = M_PI / nSteps;
  thetaMin = thetaStep / 2.0;
  phiStep = thetaStep * 2.0;
  bs = 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 > bs) bs = sigma;
    }
  }

  return bs;  
}


double Shapes_FF::findCutoffDistance(ShapeAtomType* st) {
  int i, j,  nSteps;
  double theta, thetaStep, thetaMin, costheta;
  double phi, phiStep;
  double sigma, range;
  double bigCut, tempCut;

  nSteps = 16;

  thetaStep = M_PI / nSteps;
  thetaMin = thetaStep / 2.0;
  phiStep = thetaStep * 2.0;
  bigCut = 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);
      range = st->getRangeValueAt(costheta, phi);

       // cutoff for a shape is taken to be (2.5*rangeVal + contactVal)
      tempCut = 1.5*range + sigma;

      if (tempCut > bigCut) bigCut = tempCut; 
    }
  }
 
  return bigCut;  
}
Revision:	1729
Committed:	Thu Nov 11 21:46:29 2004 UTC (19 years, 7 months ago) by chrisfen
File size:	16295 byte(s)
Log Message:	Got rid of some write statements
#	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	chrisfen	1688	double tempShapesRcut = bigContact;
57	gezelter	1520	MPI_Allreduce( &tempShapesRcut, &shapesRcut, 1, MPI_DOUBLE, MPI_MAX,
58			MPI_COMM_WORLD);
59	chrisfen	1688	#else
60			shapesRcut = bigContact;
61	gezelter	1520	#endif //is_mpi
62			entry_plug->setDefaultRcut(shapesRcut);
63			}
64
65
66	gezelter	1628	void Shapes_FF::initForceField(){
67			initFortran(0);
68	gezelter	1520	}
69
70
71	gezelter	1650	void Shapes_FF::readParams( void ){
72	gezelter	1520
73	chrisfen	1646	char readLine[1024];
74	gezelter	1650
75			string fileName;
76			string shapeFileName;
77			string tempString;
78
79	chrisfen	1646	char *nameToken;
80			char *delim = " ,;\t\n";
81			int nTokens, i;
82			int nContact = 0;
83			int nRange = 0;
84			int nStrength = 0;
85			int myATID;
86	gezelter	1650	int isError;
87	chrisfen	1646	string nameString;
88	gezelter	1650	AtomType* at;
89			DirectionalAtomType* dat;
90			ShapeAtomType* st;
91
92	chrisfen	1646	map<string, AtomType*>::iterator iter;
93	gezelter	1520
94	chrisfen	1646	// vectors for shape transfer to fortran
95	gezelter	1650	vector<RealSphericalHarmonic*> tempSHVector;
96	chrisfen	1646	vector<int> contactL;
97			vector<int> contactM;
98			vector<int> contactFunc;
99			vector<double> contactCoeff;
100			vector<int> rangeL;
101			vector<int> rangeM;
102			vector<int> rangeFunc;
103			vector<double> rangeCoeff;
104			vector<int> strengthL;
105			vector<int> strengthM;
106			vector<int> strengthFunc;
107			vector<double> strengthCoeff;
108
109			// generate the force file name
110	gezelter	1650	fileName = "Shapes.frc";
111	chrisfen	1646
112			// attempt to open the file in the current directory first.
113	gezelter	1650	forceFile.open( fileName.c_str() );
114	chrisfen	1646
115	gezelter	1650	if( forceFile == NULL ){
116	gezelter	1520
117	gezelter	1650	tempString = ffPath;
118			tempString += "/";
119			tempString += fileName;
120			fileName = tempString;
121	gezelter	1520
122	gezelter	1650	forceFile.open( fileName.c_str() );
123	gezelter	1520
124	gezelter	1650	if( forceFile == NULL ){
125	chrisfen	1646
126			sprintf( painCave.errMsg,
127			"Error opening the force field parameter file:\n"
128			"\t%s\n"
129			"\tHave you tried setting the FORCE_PARAM_PATH environment "
130			"variable?\n",
131	gezelter	1650	fileName.c_str() );
132	chrisfen	1646	painCave.severity = OOPSE_ERROR;
133	gezelter	1520	painCave.isFatal = 1;
134			simError();
135			}
136	chrisfen	1646	}
137
138			// read in the shape types.
139
140	gezelter	1650	findBegin( forceFile, "ShapeTypes" );
141	chrisfen	1646
142	gezelter	1650	while( !forceFile.eof() ){
143			forceFile.getline( readLine, sizeof(readLine) );
144	gezelter	1520
145	chrisfen	1646	// toss comment lines
146			if( readLine[0] != '!' && readLine[0] != '#' ){
147
148			if (isEndLine(readLine)) break;
149
150	gezelter	1650	nTokens = countTokens(readLine, " ,;\t");
151	chrisfen	1646	if (nTokens != 0) {
152			if (nTokens < 2) {
153			sprintf( painCave.errMsg,
154	gezelter	1650	"Not enough data on a ShapeTypes line in file: %s\n",
155			fileName.c_str() );
156	chrisfen	1646	painCave.severity = OOPSE_ERROR;
157			painCave.isFatal = 1;
158			simError();
159			}
160	gezelter	1520
161	chrisfen	1646	nameToken = strtok( readLine, delim );
162			shapeFileName = strtok( NULL, delim );
163	gezelter	1520
164	chrisfen	1646	// strings are not char arrays!
165			nameString = nameToken;
166	gezelter	1520
167	chrisfen	1646	// does this AtomType name already exist in the map?
168			iter = atomTypeMap.find(nameString);
169			if (iter == atomTypeMap.end()) {
170			// no, it doesn't, so we may proceed:
171
172	gezelter	1650	st = new ShapeAtomType();
173	gezelter	1520
174	chrisfen	1646	st->setName(nameString);
175	chrisfen	1688	myATID = atomTypeMap.size() + 1;
176	chrisfen	1646	st->setIdent(myATID);
177			parseShapeFile(shapeFileName, st);
178			st->complete();
179			atomTypeMap.insert(make_pair(nameString, st));
180
181			} else {
182			// atomType map already contained this string (i.e. it was
183			// declared in a previous block, and we just need to add
184			// the shape-specific information for this AtomType:
185	gezelter	1520
186	gezelter	1650	st = (ShapeAtomType*)(iter->second);
187	chrisfen	1646	parseShapeFile(shapeFileName, st);
188			}
189			}
190	gezelter	1520	}
191			}
192
193	chrisfen	1646	#ifdef IS_MPI
194	gezelter	1520
195	chrisfen	1646	// looks like all the processors have their ShapeType vectors ready...
196			sprintf( checkPointMsg,
197			"Shapes_FF shape objects read successfully." );
198			MPIcheckPoint();
199	gezelter	1520
200	chrisfen	1646	#endif // is_mpi
201	gezelter	1520
202	chrisfen	1646	// pack up and send the necessary info to fortran
203			for (iter = atomTypeMap.begin(); iter != atomTypeMap.end(); ++iter){
204	gezelter	1520
205	gezelter	1650	at = (AtomType*)(iter->second);
206	gezelter	1520
207	chrisfen	1646	if (at->isDirectional()) {
208	gezelter	1520
209	chrisfen	1646	dat = (DirectionalAtomType*)at;
210	gezelter	1520
211	chrisfen	1646	if (dat->isShape()) {
212	gezelter	1520
213	chrisfen	1646	st = (ShapeAtomType*)at;
214
215			contactL.clear();
216			contactM.clear();
217			contactFunc.clear();
218			contactCoeff.clear();
219
220			tempSHVector = st->getContactFuncs();
221
222			nContact = tempSHVector.size();
223			for (i=0; i<nContact; i++){
224	gezelter	1650	contactL.push_back(tempSHVector[i]->getL());
225			contactM.push_back(tempSHVector[i]->getM());
226			contactFunc.push_back(tempSHVector[i]->getFunctionType());
227			contactCoeff.push_back(tempSHVector[i]->getCoefficient());
228	chrisfen	1646	}
229
230			rangeL.clear();
231			rangeM.clear();
232			rangeFunc.clear();
233			rangeCoeff.clear();
234
235			tempSHVector = st->getRangeFuncs();
236
237			nRange = tempSHVector.size();
238			for (i=0; i<nRange; i++){
239	gezelter	1650	rangeL.push_back(tempSHVector[i]->getL());
240			rangeM.push_back(tempSHVector[i]->getM());
241			rangeFunc.push_back(tempSHVector[i]->getFunctionType());
242			rangeCoeff.push_back(tempSHVector[i]->getCoefficient());
243	chrisfen	1646	}
244
245			strengthL.clear();
246			strengthM.clear();
247			strengthFunc.clear();
248			strengthCoeff.clear();
249
250			tempSHVector = st->getStrengthFuncs();
251
252			nStrength = tempSHVector.size();
253			for (i=0; i<nStrength; i++){
254	gezelter	1650	strengthL.push_back(tempSHVector[i]->getL());
255			strengthM.push_back(tempSHVector[i]->getM());
256			strengthFunc.push_back(tempSHVector[i]->getFunctionType());
257			strengthCoeff.push_back(tempSHVector[i]->getCoefficient());
258	chrisfen	1646	}
259
260			isError = 0;
261			myATID = at->getIdent();
262	chrisfen	1688
263	gezelter	1650	makeShape( &nContact, &contactL[0], &contactM[0], &contactFunc[0],
264			&contactCoeff[0],
265			&nRange, &rangeL[0], &rangeM[0], &rangeFunc[0],
266			&rangeCoeff[0],
267			&nStrength, &strengthL[0], &strengthM[0],
268			&strengthFunc[0], &strengthCoeff[0],
269	chrisfen	1646	&myATID,
270			&isError);
271	chrisfen	1688
272	chrisfen	1646	if( isError ){
273			sprintf( painCave.errMsg,
274			"Error initializing the \"%s\" shape in fortran\n",
275	gezelter	1650	(iter->first).c_str() );
276	chrisfen	1646	painCave.isFatal = 1;
277			simError();
278			}
279			}
280	gezelter	1520	}
281			}
282
283	chrisfen	1688	isError = 0;
284			completeShapeFF(&isError);
285			if( isError ){
286			sprintf( painCave.errMsg,
287			"Error completing Shape FF in fortran\n");
288			painCave.isFatal = 1;
289			simError();
290			}
291
292	gezelter	1520	#ifdef IS_MPI
293			sprintf( checkPointMsg,
294			"Shapes_FF atom structures successfully sent to fortran\n" );
295			MPIcheckPoint();
296			#endif // is_mpi
297
298			}
299
300	gezelter	1650	void Shapes_FF::cleanMe( void ){
301
302			}
303
304			void Shapes_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
305	chrisfen	1646
306			int i,j,k;
307	gezelter	1650	map<string, AtomType*>::iterator iter;
308	gezelter	1520
309			// initialize the atoms
310	chrisfen	1646	DirectionalAtom* dAtom;
311			AtomType* at;
312			DirectionalAtomType* dat;
313	gezelter	1650	ShapeAtomType* sat;
314	chrisfen	1688	double longCutoff;
315	chrisfen	1646	double ji[3];
316			double inertialMat[3][3];
317			Mat3x3d momInt;
318			string myTypeString;
319
320	chrisfen	1688	bigContact = 0.0;
321
322	gezelter	1520	for( i=0; i<nAtoms; i++ ){
323
324	chrisfen	1646	myTypeString = the_atoms[i]->getType();
325			iter = atomTypeMap.find(myTypeString);
326
327			if (iter == atomTypeMap.end()) {
328	gezelter	1520	sprintf( painCave.errMsg,
329			"AtomType error, %s not found in force file.\n",
330			the_atoms[i]->getType() );
331			painCave.isFatal = 1;
332			simError();
333	chrisfen	1646	} else {
334	gezelter	1520
335	chrisfen	1646	at = (AtomType*)(iter->second);
336	gezelter	1520
337	chrisfen	1646	the_atoms[i]->setMass( at->getMass() );
338			the_atoms[i]->setIdent( at->getIdent() );
339	gezelter	1650
340			if ( at->isShape() ) {
341
342			sat = (ShapeAtomType*)at;
343	chrisfen	1688	longCutoff = findCutoffDistance(sat);
344			if (longCutoff > bigContact) bigContact = longCutoff;
345			cout << bigContact << " is the cutoff value\n";
346
347			entry_plug->useShapes = 1;
348	chrisfen	1646	}
349	gezelter	1520
350	chrisfen	1646	the_atoms[i]->setHasCharge(at->isCharge());
351	gezelter	1520
352	chrisfen	1646	if( at->isDirectional() ){
353	gezelter	1520
354	chrisfen	1646	dat = (DirectionalAtomType*)at;
355			dAtom = (DirectionalAtom *) the_atoms[i];
356	gezelter	1672
357	chrisfen	1646	momInt = dat->getI();
358	gezelter	1520
359	chrisfen	1646	// zero out the moments of inertia matrix
360			for( j=0; j<3; j++ )
361			for( k=0; k<3; k++ )
362			inertialMat[j][k] = momInt(j,k);
363			dAtom->setI( inertialMat );
364	gezelter	1520
365	chrisfen	1646	ji[0] = 0.0;
366			ji[1] = 0.0;
367			ji[2] = 0.0;
368			dAtom->setJ( ji );
369	gezelter	1520
370	chrisfen	1646	if (dat->isDipole()) {
371			dAtom->setHasDipole( dat->isDipole() );
372			entry_plug->n_dipoles++;
373			}
374			}
375	gezelter	1520	}
376			}
377			}
378
379	chrisfen	1646	void Shapes_FF::initializeBonds( int nBonds, Bond** BondArray,
380			bond_pair* the_bonds ){
381	gezelter	1520
382	chrisfen	1646	if( nBonds ){
383	gezelter	1520	sprintf( painCave.errMsg,
384	chrisfen	1646	"Shapes_FF does not support bonds.\n" );
385	gezelter	1520	painCave.isFatal = 1;
386			simError();
387			}
388	chrisfen	1646	}
389	gezelter	1520
390	chrisfen	1646	void Shapes_FF::initializeBends( int nBends, Bend** bendArray,
391			bend_set* the_bends ){
392
393			if( nBends ){
394	gezelter	1520	sprintf( painCave.errMsg,
395	chrisfen	1646	"Shapes_FF does not support bends.\n" );
396	gezelter	1520	painCave.isFatal = 1;
397			simError();
398			}
399	chrisfen	1646	}
400	gezelter	1520
401	chrisfen	1646	void Shapes_FF::initializeTorsions( int nTorsions, Torsion** torsionArray,
402			torsion_set* the_torsions ){
403	gezelter	1520
404	chrisfen	1646	if( nTorsions ){
405	gezelter	1520	sprintf( painCave.errMsg,
406	chrisfen	1646	"Shapes_FF does not support torsions.\n" );
407	gezelter	1520	painCave.isFatal = 1;
408			simError();
409			}
410	chrisfen	1646	}
411	gezelter	1520
412	gezelter	1650	void Shapes_FF::parseShapeFile(string shapeFileName, ShapeAtomType* st){
413	chrisfen	1646	const int MAXLEN = 1024;
414			char inLine[MAXLEN];
415			char *token;
416			char *delim = " ,;\t\n";
417	gezelter	1650	int nTokens;
418	chrisfen	1646	Mat3x3d momInert;
419	gezelter	1650	RealSphericalHarmonic* rsh;
420			vector<RealSphericalHarmonic*> functionVector;
421			ifstrstream shapeFile;
422			string tempString;
423	gezelter	1520
424	gezelter	1650	shapeFile.open( shapeFileName.c_str() );
425	gezelter	1520
426	gezelter	1650	if( shapeFile == NULL ){
427
428			tempString = ffPath;
429			tempString += "/";
430			tempString += shapeFileName;
431			shapeFileName = tempString;
432
433			shapeFile.open( shapeFileName.c_str() );
434
435			if( shapeFile == NULL ){
436
437			sprintf( painCave.errMsg,
438			"Error opening the shape file:\n"
439			"\t%s\n"
440			"\tHave you tried setting the FORCE_PARAM_PATH environment "
441			"variable?\n",
442			shapeFileName.c_str() );
443			painCave.severity = OOPSE_ERROR;
444			painCave.isFatal = 1;
445			simError();
446			}
447			}
448
449			// read in the shape types.
450
451	chrisfen	1646	// first grab the values in the ShapeInfo section
452	gezelter	1650	findBegin( shapeFile, "ShapeInfo");
453	gezelter	1520
454	chrisfen	1646	shapeFile.getline(inLine, MAXLEN);
455			while( !shapeFile.eof() ) {
456			// toss comment lines
457			if( inLine[0] != '!' && inLine[0] != '#' ){
458			// end marks section completion
459			if (isEndLine(inLine)) break;
460
461	gezelter	1650	nTokens = countTokens(inLine, delim);
462	chrisfen	1646	if (nTokens != 0) {
463			if (nTokens < 5) {
464			sprintf( painCave.errMsg,
465	gezelter	1650	"Not enough data on a ShapeInfo line in file: %s\n",
466			shapeFileName.c_str() );
467	chrisfen	1646	painCave.severity = OOPSE_ERROR;
468			painCave.isFatal = 1;
469			simError();
470	chrisfen	1688	} else {
471	chrisfen	1646	token = strtok(inLine, delim);
472			token = strtok(NULL, delim);
473			st->setMass(atof(token));
474			token = strtok(NULL, delim);
475			momInert(0,0) = atof(token);
476			token = strtok(NULL, delim);
477			momInert(1,1) = atof(token);
478			token = strtok(NULL, delim);
479			momInert(2,2) = atof(token);
480			st->setI(momInert);
481			}
482			}
483	gezelter	1520	}
484	gezelter	1672	shapeFile.getline(inLine, MAXLEN);
485	gezelter	1520	}
486
487	chrisfen	1646	// now grab the contact functions
488			findBegin(shapeFile, "ContactFunctions");
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	gezelter	1650	nTokens = countTokens(inLine, delim);
498	chrisfen	1646	if (nTokens != 0) {
499			if (nTokens < 4) {
500			sprintf( painCave.errMsg,
501	gezelter	1650	"Not enough data on a ContactFunctions line in file: %s\n",
502			shapeFileName.c_str() );
503	chrisfen	1646	painCave.severity = OOPSE_ERROR;
504			painCave.isFatal = 1;
505			simError();
506	chrisfen	1688	} else {
507	chrisfen	1646	// read in a spherical harmonic function
508			token = strtok(inLine, delim);
509	gezelter	1650	rsh = new RealSphericalHarmonic();
510			rsh->setL(atoi(token));
511	chrisfen	1646	token = strtok(NULL, delim);
512	gezelter	1650	rsh->setM(atoi(token));
513	chrisfen	1646	token = strtok(NULL, delim);
514			if (!strcasecmp("sin",token))
515	gezelter	1650	rsh->makeSinFunction();
516	chrisfen	1646	else
517	gezelter	1650	rsh->makeCosFunction();
518	chrisfen	1646	token = strtok(NULL, delim);
519	gezelter	1650	rsh->setCoefficient(atof(token));
520	chrisfen	1646
521	gezelter	1650	functionVector.push_back(rsh);
522	chrisfen	1646	}
523			}
524	gezelter	1520	}
525	gezelter	1672	shapeFile.getline(inLine, MAXLEN);
526	gezelter	1520	}
527
528	chrisfen	1646	// pass contact functions to ShapeType
529	chrisfen	1688
530	chrisfen	1646	st->setContactFuncs(functionVector);
531	gezelter	1520
532	chrisfen	1646	// now grab the range functions
533			findBegin(shapeFile, "RangeFunctions");
534			functionVector.clear();
535	gezelter	1520
536	chrisfen	1646	shapeFile.getline(inLine, MAXLEN);
537			while( !shapeFile.eof() ) {
538			// toss comment lines
539			if( inLine[0] != '!' && inLine[0] != '#' ){
540			// end marks section completion
541			if (isEndLine(inLine)) break;
542
543	gezelter	1650	nTokens = countTokens(inLine, delim);
544	chrisfen	1646	if (nTokens != 0) {
545			if (nTokens < 4) {
546			sprintf( painCave.errMsg,
547	gezelter	1650	"Not enough data on a RangeFunctions line in file: %s\n",
548			shapeFileName.c_str() );
549	chrisfen	1646	painCave.severity = OOPSE_ERROR;
550			painCave.isFatal = 1;
551			simError();
552	chrisfen	1688	} else {
553	chrisfen	1646
554			// read in a spherical harmonic function
555			token = strtok(inLine, delim);
556	gezelter	1650
557			rsh = new RealSphericalHarmonic();
558			rsh->setL(atoi(token));
559	chrisfen	1646	token = strtok(NULL, delim);
560	gezelter	1650	rsh->setM(atoi(token));
561	chrisfen	1646	token = strtok(NULL, delim);
562			if (!strcasecmp("sin",token))
563	gezelter	1650	rsh->makeSinFunction();
564	chrisfen	1646	else
565	gezelter	1650	rsh->makeCosFunction();
566	chrisfen	1646	token = strtok(NULL, delim);
567	gezelter	1650	rsh->setCoefficient(atof(token));
568	chrisfen	1646
569	gezelter	1650	functionVector.push_back(rsh);
570	chrisfen	1646	}
571			}
572	gezelter	1520	}
573	gezelter	1672	shapeFile.getline(inLine, MAXLEN);
574	chrisfen	1646	}
575	gezelter	1520
576	chrisfen	1646	// pass range functions to ShapeType
577			st->setRangeFuncs(functionVector);
578	gezelter	1520
579	chrisfen	1646	// finally grab the strength functions
580			findBegin(shapeFile, "StrengthFunctions");
581			functionVector.clear();
582	gezelter	1520
583	chrisfen	1646	shapeFile.getline(inLine, MAXLEN);
584			while( !shapeFile.eof() ) {
585			// toss comment lines
586			if( inLine[0] != '!' && inLine[0] != '#' ){
587			// end marks section completion
588			if (isEndLine(inLine)) break;
589
590	gezelter	1650	nTokens = countTokens(inLine, delim);
591	chrisfen	1646	if (nTokens != 0) {
592			if (nTokens < 4) {
593			sprintf( painCave.errMsg,
594	gezelter	1650	"Not enough data on a StrengthFunctions line in file: %s\n",
595			shapeFileName.c_str() );
596	chrisfen	1646	painCave.severity = OOPSE_ERROR;
597			painCave.isFatal = 1;
598			simError();
599	chrisfen	1688	} else {
600	chrisfen	1646
601			// read in a spherical harmonic function
602			token = strtok(inLine, delim);
603	gezelter	1650	rsh = new RealSphericalHarmonic();
604			rsh->setL(atoi(token));
605	chrisfen	1646	token = strtok(NULL, delim);
606	gezelter	1650	rsh->setM(atoi(token));
607	chrisfen	1646	token = strtok(NULL, delim);
608			if (!strcasecmp("sin",token))
609	gezelter	1650	rsh->makeSinFunction();
610	chrisfen	1646	else
611	gezelter	1650	rsh->makeCosFunction();
612	chrisfen	1646	token = strtok(NULL, delim);
613	gezelter	1650	rsh->setCoefficient(atof(token));
614	chrisfen	1646
615	gezelter	1650	functionVector.push_back(rsh);
616	chrisfen	1646	}
617			}
618	gezelter	1520	}
619	gezelter	1672	shapeFile.getline(inLine, MAXLEN);
620	gezelter	1520	}
621
622	chrisfen	1646	// pass strength functions to ShapeType
623			st->setStrengthFuncs(functionVector);
624
625			// we're done reading from this file
626			shapeFile.close();
627	gezelter	1520	}
628	chrisfen	1646
629	gezelter	1650	double Shapes_FF::findLargestContactDistance(ShapeAtomType* st) {
630			int i, j, nSteps;
631			double theta, thetaStep, thetaMin, costheta;
632			double phi, phiStep;
633	chrisfen	1688	double sigma, bs;
634
635	gezelter	1650	nSteps = 16;
636
637			thetaStep = M_PI / nSteps;
638			thetaMin = thetaStep / 2.0;
639			phiStep = thetaStep * 2.0;
640	chrisfen	1688	bs = 0.0;
641	gezelter	1650
642			for (i = 0; i < nSteps; i++) {
643
644			theta = thetaMin + i * thetaStep;
645			costheta = cos(theta);
646
647			for (j = 0; j < nSteps; j++) {
648
649			phi = j*phiStep;
650
651			sigma = st->getContactValueAt(costheta, phi);
652
653	chrisfen	1688	if (sigma > bs) bs = sigma;
654	gezelter	1650	}
655			}
656
657	chrisfen	1688	return bs;
658	gezelter	1650	}
659	chrisfen	1688
660
661			double Shapes_FF::findCutoffDistance(ShapeAtomType* st) {
662			int i, j, nSteps;
663			double theta, thetaStep, thetaMin, costheta;
664			double phi, phiStep;
665			double sigma, range;
666			double bigCut, tempCut;
667
668			nSteps = 16;
669
670			thetaStep = M_PI / nSteps;
671			thetaMin = thetaStep / 2.0;
672			phiStep = thetaStep * 2.0;
673			bigCut = 0.0;
674
675			for (i = 0; i < nSteps; i++) {
676
677			theta = thetaMin + i * thetaStep;
678			costheta = cos(theta);
679
680			for (j = 0; j < nSteps; j++) {
681
682			phi = j*phiStep;
683
684			sigma = st->getContactValueAt(costheta, phi);
685			range = st->getRangeValueAt(costheta, phi);
686
687			// cutoff for a shape is taken to be (2.5*rangeVal + contactVal)
688	chrisfen	1729	tempCut = 1.5*range + sigma;
689	chrisfen	1688
690			if (tempCut > bigCut) bigCut = tempCut;
691			}
692			}
693
694			return bigCut;
695			}