src/UseTheForce/Shapes_FF.cpp

 /*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
 
#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 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];

   std::string fileName;
   std::string shapeFileName;
   std::string tempString;

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

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

  // vectors for shape transfer to fortran
   std::vector<RealSphericalHarmonic*> tempSHVector;
   std::vector<int> contactL;
   std::vector<int> contactM;
   std::vector<int> contactFunc;
   std::vector<double> contactCoeff;
   std::vector<int> rangeL;
   std::vector<int> rangeM;
   std::vector<int> rangeFunc;
   std::vector<double> rangeCoeff;
   std::vector<int> strengthL;
   std::vector<int> strengthM;
   std::vector<int> strengthFunc;
   std::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  std::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;
   std::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;
   std::string myTypeString;

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

    if (iter == atomTypeMap.end()) {
      sprintf( painCave.errMsg, 
               "AtomType error, %s not found in force file.\n",
               the_atoms[i]->getType().c_str() );
      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;
   std::vector<RealSphericalHarmonic*> functionVector;
  ifstrstream shapeFile;
   std::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:	1927
Committed:	Wed Jan 12 17:14:35 2005 UTC (21 years ago) by tim
File size:	17238 byte(s)
Log Message:	change license
#	User	Rev	Content
1	tim	1927	/*
2			* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	chrisfen	1646	*
4	tim	1927	* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8	chrisfen	1646	*
9	tim	1927	* 1. Acknowledgement of the program authors must be made in any
10			* publication of scientific results based in part on use of the
11			* program. An acceptable form of acknowledgement is citation of
12			* the article in which the program was described (Matthew
13			* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			* Parallel Simulation Engine for Molecular Dynamics,"
16			* J. Comput. Chem. 26, pp. 252-271 (2005))
17			*
18			* 2. Redistributions of source code must retain the above copyright
19			* notice, this list of conditions and the following disclaimer.
20			*
21			* 3. Redistributions in binary form must reproduce the above copyright
22			* notice, this list of conditions and the following disclaimer in the
23			* documentation and/or other materials provided with the
24			* distribution.
25			*
26			* This software is provided "AS IS," without a warranty of any
27			* kind. All express or implied conditions, representations and
28			* warranties, including any implied warranty of merchantability,
29			* fitness for a particular purpose or non-infringement, are hereby
30			* excluded. The University of Notre Dame and its licensors shall not
31			* be liable for any damages suffered by licensee as a result of
32			* using, modifying or distributing the software or its
33			* derivatives. In no event will the University of Notre Dame or its
34			* licensors be liable for any lost revenue, profit or data, or for
35			* direct, indirect, special, consequential, incidental or punitive
36			* damages, however caused and regardless of the theory of liability,
37			* arising out of the use of or inability to use software, even if the
38			* University of Notre Dame has been advised of the possibility of
39			* such damages.
40	chrisfen	1646	*/
41	tim	1927
42	gezelter	1520	#include <stdlib.h>
43			#include <stdio.h>
44			#include <string.h>
45	chrisfen	1646	#include <map>
46	gezelter	1650	#include <cmath>
47	chrisfen	1646	#include <iostream>
48	gezelter	1520
49			#ifdef IS_MPI
50			#include <mpi.h>
51			#endif //is_mpi
52
53			#include "UseTheForce/ForceFields.hpp"
54			#include "primitives/SRI.hpp"
55			#include "utils/simError.h"
56	gezelter	1650	#include "utils/StringUtils.hpp"
57	chrisfen	1646	#include "io/basic_ifstrstream.hpp"
58			#include "math/RealSphericalHarmonic.hpp"
59			#include "math/SquareMatrix3.hpp"
60	gezelter	1650	#include "types/ShapeAtomType.hpp"
61	chrisfen	1646	#include "UseTheForce/DarkSide/atype_interface.h"
62	gezelter	1613	#include "UseTheForce/DarkSide/shapes_interface.h"
63	gezelter	1520
64			#ifdef IS_MPI
65			#include "UseTheForce/mpiForceField.h"
66			#endif // is_mpi
67
68	tim	1827
69	gezelter	1650	using namespace oopse;
70	gezelter	1520
71			Shapes_FF::~Shapes_FF(){
72
73			#ifdef IS_MPI
74			if( worldRank == 0 ){
75			#endif // is_mpi
76
77	gezelter	1650	forceFile.close();
78	gezelter	1520
79			#ifdef IS_MPI
80			}
81			#endif // is_mpi
82			}
83
84
85			void Shapes_FF::calcRcut( void ){
86
87			#ifdef IS_MPI
88			double tempShapesRcut = shapesRcut;
89			MPI_Allreduce( &tempShapesRcut, &shapesRcut, 1, MPI_DOUBLE, MPI_MAX,
90			MPI_COMM_WORLD);
91			#endif //is_mpi
92			entry_plug->setDefaultRcut(shapesRcut);
93			}
94
95
96	gezelter	1628	void Shapes_FF::initForceField(){
97			initFortran(0);
98	gezelter	1520	}
99
100
101	gezelter	1650	void Shapes_FF::readParams( void ){
102	gezelter	1520
103	chrisfen	1646	char readLine[1024];
104	gezelter	1650
105	tim	1829	std::string fileName;
106			std::string shapeFileName;
107			std::string tempString;
108	gezelter	1650
109	chrisfen	1646	char *nameToken;
110			char *delim = " ,;\t\n";
111			int nTokens, i;
112			int nContact = 0;
113			int nRange = 0;
114			int nStrength = 0;
115			int myATID;
116	gezelter	1650	int isError;
117	tim	1829	std::string nameString;
118	gezelter	1650	AtomType* at;
119			DirectionalAtomType* dat;
120			ShapeAtomType* st;
121
122	tim	1826	std::map<string, AtomType*>::iterator iter;
123	gezelter	1520
124	chrisfen	1646	// vectors for shape transfer to fortran
125	tim	1826	std::vector<RealSphericalHarmonic*> tempSHVector;
126			std::vector<int> contactL;
127			std::vector<int> contactM;
128			std::vector<int> contactFunc;
129			std::vector<double> contactCoeff;
130			std::vector<int> rangeL;
131			std::vector<int> rangeM;
132			std::vector<int> rangeFunc;
133			std::vector<double> rangeCoeff;
134			std::vector<int> strengthL;
135			std::vector<int> strengthM;
136			std::vector<int> strengthFunc;
137			std::vector<double> strengthCoeff;
138	chrisfen	1646
139			// generate the force file name
140	gezelter	1650	fileName = "Shapes.frc";
141	chrisfen	1646
142			// attempt to open the file in the current directory first.
143	gezelter	1650	forceFile.open( fileName.c_str() );
144	chrisfen	1646
145	gezelter	1650	if( forceFile == NULL ){
146	gezelter	1520
147	gezelter	1650	tempString = ffPath;
148			tempString += "/";
149			tempString += fileName;
150			fileName = tempString;
151	gezelter	1520
152	gezelter	1650	forceFile.open( fileName.c_str() );
153	gezelter	1520
154	gezelter	1650	if( forceFile == NULL ){
155	chrisfen	1646
156			sprintf( painCave.errMsg,
157			"Error opening the force field parameter file:\n"
158			"\t%s\n"
159			"\tHave you tried setting the FORCE_PARAM_PATH environment "
160			"variable?\n",
161	gezelter	1650	fileName.c_str() );
162	chrisfen	1646	painCave.severity = OOPSE_ERROR;
163	gezelter	1520	painCave.isFatal = 1;
164			simError();
165			}
166	chrisfen	1646	}
167
168			// read in the shape types.
169
170	gezelter	1650	findBegin( forceFile, "ShapeTypes" );
171	chrisfen	1646
172	gezelter	1650	while( !forceFile.eof() ){
173			forceFile.getline( readLine, sizeof(readLine) );
174	gezelter	1520
175	chrisfen	1646	// toss comment lines
176			if( readLine[0] != '!' && readLine[0] != '#' ){
177
178			if (isEndLine(readLine)) break;
179
180	gezelter	1650	nTokens = countTokens(readLine, " ,;\t");
181	chrisfen	1646	if (nTokens != 0) {
182			if (nTokens < 2) {
183			sprintf( painCave.errMsg,
184	gezelter	1650	"Not enough data on a ShapeTypes line in file: %s\n",
185			fileName.c_str() );
186	chrisfen	1646	painCave.severity = OOPSE_ERROR;
187			painCave.isFatal = 1;
188			simError();
189			}
190	gezelter	1520
191	chrisfen	1646	nameToken = strtok( readLine, delim );
192			shapeFileName = strtok( NULL, delim );
193	gezelter	1520
194	chrisfen	1646	// strings are not char arrays!
195			nameString = nameToken;
196	gezelter	1520
197	chrisfen	1646	// does this AtomType name already exist in the map?
198			iter = atomTypeMap.find(nameString);
199			if (iter == atomTypeMap.end()) {
200			// no, it doesn't, so we may proceed:
201
202	gezelter	1650	st = new ShapeAtomType();
203	gezelter	1520
204	chrisfen	1646	st->setName(nameString);
205			myATID = atomTypeMap.size();
206			st->setIdent(myATID);
207			parseShapeFile(shapeFileName, st);
208			st->complete();
209			atomTypeMap.insert(make_pair(nameString, st));
210
211			} else {
212	tim	1829	// atomType map already contained this std::string (i.e. it was
213	chrisfen	1646	// declared in a previous block, and we just need to add
214			// the shape-specific information for this AtomType:
215	gezelter	1520
216	gezelter	1650	st = (ShapeAtomType*)(iter->second);
217	chrisfen	1646	parseShapeFile(shapeFileName, st);
218			}
219			}
220	gezelter	1520	}
221			}
222
223	chrisfen	1646	#ifdef IS_MPI
224	gezelter	1520
225	chrisfen	1646	// looks like all the processors have their ShapeType vectors ready...
226			sprintf( checkPointMsg,
227			"Shapes_FF shape objects read successfully." );
228			MPIcheckPoint();
229	gezelter	1520
230	chrisfen	1646	#endif // is_mpi
231	gezelter	1520
232	chrisfen	1646	// pack up and send the necessary info to fortran
233			for (iter = atomTypeMap.begin(); iter != atomTypeMap.end(); ++iter){
234	gezelter	1520
235	gezelter	1650	at = (AtomType*)(iter->second);
236	gezelter	1520
237	chrisfen	1646	if (at->isDirectional()) {
238	gezelter	1520
239	chrisfen	1646	dat = (DirectionalAtomType*)at;
240	gezelter	1520
241	chrisfen	1646	if (dat->isShape()) {
242	gezelter	1520
243	chrisfen	1646	st = (ShapeAtomType*)at;
244
245			contactL.clear();
246			contactM.clear();
247			contactFunc.clear();
248			contactCoeff.clear();
249
250			tempSHVector = st->getContactFuncs();
251
252			nContact = tempSHVector.size();
253			for (i=0; i<nContact; i++){
254	gezelter	1650	contactL.push_back(tempSHVector[i]->getL());
255			contactM.push_back(tempSHVector[i]->getM());
256			contactFunc.push_back(tempSHVector[i]->getFunctionType());
257			contactCoeff.push_back(tempSHVector[i]->getCoefficient());
258	chrisfen	1646	}
259
260			rangeL.clear();
261			rangeM.clear();
262			rangeFunc.clear();
263			rangeCoeff.clear();
264
265			tempSHVector = st->getRangeFuncs();
266
267			nRange = tempSHVector.size();
268			for (i=0; i<nRange; i++){
269	gezelter	1650	rangeL.push_back(tempSHVector[i]->getL());
270			rangeM.push_back(tempSHVector[i]->getM());
271			rangeFunc.push_back(tempSHVector[i]->getFunctionType());
272			rangeCoeff.push_back(tempSHVector[i]->getCoefficient());
273	chrisfen	1646	}
274
275			strengthL.clear();
276			strengthM.clear();
277			strengthFunc.clear();
278			strengthCoeff.clear();
279
280			tempSHVector = st->getStrengthFuncs();
281
282			nStrength = tempSHVector.size();
283			for (i=0; i<nStrength; i++){
284	gezelter	1650	strengthL.push_back(tempSHVector[i]->getL());
285			strengthM.push_back(tempSHVector[i]->getM());
286			strengthFunc.push_back(tempSHVector[i]->getFunctionType());
287			strengthCoeff.push_back(tempSHVector[i]->getCoefficient());
288	chrisfen	1646	}
289
290			isError = 0;
291			myATID = at->getIdent();
292	gezelter	1650	makeShape( &nContact, &contactL[0], &contactM[0], &contactFunc[0],
293			&contactCoeff[0],
294			&nRange, &rangeL[0], &rangeM[0], &rangeFunc[0],
295			&rangeCoeff[0],
296			&nStrength, &strengthL[0], &strengthM[0],
297			&strengthFunc[0], &strengthCoeff[0],
298	chrisfen	1646	&myATID,
299			&isError);
300			if( isError ){
301			sprintf( painCave.errMsg,
302			"Error initializing the \"%s\" shape in fortran\n",
303	gezelter	1650	(iter->first).c_str() );
304	chrisfen	1646	painCave.isFatal = 1;
305			simError();
306			}
307			}
308	gezelter	1520	}
309			}
310
311			#ifdef IS_MPI
312			sprintf( checkPointMsg,
313			"Shapes_FF atom structures successfully sent to fortran\n" );
314			MPIcheckPoint();
315			#endif // is_mpi
316
317			}
318
319	gezelter	1650	void Shapes_FF::cleanMe( void ){
320
321			}
322
323			void Shapes_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
324	chrisfen	1646
325			int i,j,k;
326	tim	1826	std::map<string, AtomType*>::iterator iter;
327	gezelter	1520
328			// initialize the atoms
329	chrisfen	1646	DirectionalAtom* dAtom;
330			AtomType* at;
331			DirectionalAtomType* dat;
332	gezelter	1650	ShapeAtomType* sat;
333			double sigma;
334	chrisfen	1646	double ji[3];
335			double inertialMat[3][3];
336			Mat3x3d momInt;
337	tim	1829	std::string myTypeString;
338	chrisfen	1646
339	gezelter	1520	for( i=0; i<nAtoms; i++ ){
340
341	tim	1702	myTypeString = the_atoms[i]->getType().c_str();
342	chrisfen	1646	iter = atomTypeMap.find(myTypeString);
343
344			if (iter == atomTypeMap.end()) {
345	gezelter	1520	sprintf( painCave.errMsg,
346			"AtomType error, %s not found in force file.\n",
347	tim	1702	the_atoms[i]->getType().c_str() );
348	gezelter	1520	painCave.isFatal = 1;
349			simError();
350	chrisfen	1646	} else {
351	gezelter	1520
352	chrisfen	1646	at = (AtomType*)(iter->second);
353	gezelter	1520
354	chrisfen	1646	the_atoms[i]->setMass( at->getMass() );
355			the_atoms[i]->setIdent( at->getIdent() );
356	gezelter	1650
357			if ( at->isShape() ) {
358
359			sat = (ShapeAtomType*)at;
360			sigma = findLargestContactDistance(sat);
361			if (sigma > bigSigma) bigSigma = sigma;
362	gezelter	1520
363	chrisfen	1646	}
364	gezelter	1520
365	chrisfen	1646	the_atoms[i]->setHasCharge(at->isCharge());
366	gezelter	1520
367	chrisfen	1646	if( at->isDirectional() ){
368	gezelter	1520
369	chrisfen	1646	dat = (DirectionalAtomType*)at;
370			dAtom = (DirectionalAtom *) the_atoms[i];
371	gezelter	1672
372	chrisfen	1646	momInt = dat->getI();
373	gezelter	1520
374	chrisfen	1646	// zero out the moments of inertia matrix
375			for( j=0; j<3; j++ )
376			for( k=0; k<3; k++ )
377			inertialMat[j][k] = momInt(j,k);
378			dAtom->setI( inertialMat );
379	gezelter	1520
380	chrisfen	1646	ji[0] = 0.0;
381			ji[1] = 0.0;
382			ji[2] = 0.0;
383			dAtom->setJ( ji );
384	gezelter	1520
385	chrisfen	1646	if (dat->isDipole()) {
386			dAtom->setHasDipole( dat->isDipole() );
387			entry_plug->n_dipoles++;
388			}
389			}
390	gezelter	1520	}
391			}
392			}
393
394	chrisfen	1646	void Shapes_FF::initializeBonds( int nBonds, Bond** BondArray,
395			bond_pair* the_bonds ){
396	gezelter	1520
397	chrisfen	1646	if( nBonds ){
398	gezelter	1520	sprintf( painCave.errMsg,
399	chrisfen	1646	"Shapes_FF does not support bonds.\n" );
400	gezelter	1520	painCave.isFatal = 1;
401			simError();
402			}
403	chrisfen	1646	}
404	gezelter	1520
405	chrisfen	1646	void Shapes_FF::initializeBends( int nBends, Bend** bendArray,
406			bend_set* the_bends ){
407
408			if( nBends ){
409	gezelter	1520	sprintf( painCave.errMsg,
410	chrisfen	1646	"Shapes_FF does not support bends.\n" );
411	gezelter	1520	painCave.isFatal = 1;
412			simError();
413			}
414	chrisfen	1646	}
415	gezelter	1520
416	chrisfen	1646	void Shapes_FF::initializeTorsions( int nTorsions, Torsion** torsionArray,
417			torsion_set* the_torsions ){
418	gezelter	1520
419	chrisfen	1646	if( nTorsions ){
420	gezelter	1520	sprintf( painCave.errMsg,
421	chrisfen	1646	"Shapes_FF does not support torsions.\n" );
422	gezelter	1520	painCave.isFatal = 1;
423			simError();
424			}
425	chrisfen	1646	}
426	gezelter	1520
427	gezelter	1650	void Shapes_FF::parseShapeFile(string shapeFileName, ShapeAtomType* st){
428	chrisfen	1646	const int MAXLEN = 1024;
429			char inLine[MAXLEN];
430			char *token;
431			char *delim = " ,;\t\n";
432	gezelter	1650	int nTokens;
433	chrisfen	1646	Mat3x3d momInert;
434	gezelter	1650	RealSphericalHarmonic* rsh;
435	tim	1826	std::vector<RealSphericalHarmonic*> functionVector;
436	gezelter	1650	ifstrstream shapeFile;
437	tim	1829	std::string tempString;
438	gezelter	1520
439	gezelter	1650	shapeFile.open( shapeFileName.c_str() );
440	gezelter	1520
441	gezelter	1650	if( shapeFile == NULL ){
442
443			tempString = ffPath;
444			tempString += "/";
445			tempString += shapeFileName;
446			shapeFileName = tempString;
447
448			shapeFile.open( shapeFileName.c_str() );
449
450			if( shapeFile == NULL ){
451
452			sprintf( painCave.errMsg,
453			"Error opening the shape file:\n"
454			"\t%s\n"
455			"\tHave you tried setting the FORCE_PARAM_PATH environment "
456			"variable?\n",
457			shapeFileName.c_str() );
458			painCave.severity = OOPSE_ERROR;
459			painCave.isFatal = 1;
460			simError();
461			}
462			}
463
464			// read in the shape types.
465
466	chrisfen	1646	// first grab the values in the ShapeInfo section
467	gezelter	1650	findBegin( shapeFile, "ShapeInfo");
468	gezelter	1520
469	chrisfen	1646	shapeFile.getline(inLine, MAXLEN);
470			while( !shapeFile.eof() ) {
471			// toss comment lines
472			if( inLine[0] != '!' && inLine[0] != '#' ){
473			// end marks section completion
474			if (isEndLine(inLine)) break;
475
476	gezelter	1650	nTokens = countTokens(inLine, delim);
477	chrisfen	1646	if (nTokens != 0) {
478			if (nTokens < 5) {
479			sprintf( painCave.errMsg,
480	gezelter	1650	"Not enough data on a ShapeInfo line in file: %s\n",
481			shapeFileName.c_str() );
482	chrisfen	1646	painCave.severity = OOPSE_ERROR;
483			painCave.isFatal = 1;
484			simError();
485	gezelter	1520
486	chrisfen	1646	token = strtok(inLine, delim);
487			token = strtok(NULL, delim);
488			st->setMass(atof(token));
489			token = strtok(NULL, delim);
490			momInert(0,0) = atof(token);
491			token = strtok(NULL, delim);
492			momInert(1,1) = atof(token);
493			token = strtok(NULL, delim);
494			momInert(2,2) = atof(token);
495			st->setI(momInert);
496			}
497			}
498	gezelter	1520	}
499	gezelter	1672	shapeFile.getline(inLine, MAXLEN);
500	gezelter	1520	}
501
502	chrisfen	1646	// now grab the contact functions
503			findBegin(shapeFile, "ContactFunctions");
504			functionVector.clear();
505	gezelter	1520
506	chrisfen	1646	shapeFile.getline(inLine, MAXLEN);
507			while( !shapeFile.eof() ) {
508			// toss comment lines
509			if( inLine[0] != '!' && inLine[0] != '#' ){
510			// end marks section completion
511			if (isEndLine(inLine)) break;
512
513	gezelter	1650	nTokens = countTokens(inLine, delim);
514	chrisfen	1646	if (nTokens != 0) {
515			if (nTokens < 4) {
516			sprintf( painCave.errMsg,
517	gezelter	1650	"Not enough data on a ContactFunctions line in file: %s\n",
518			shapeFileName.c_str() );
519	chrisfen	1646	painCave.severity = OOPSE_ERROR;
520			painCave.isFatal = 1;
521			simError();
522
523			// read in a spherical harmonic function
524			token = strtok(inLine, delim);
525	gezelter	1650	rsh = new RealSphericalHarmonic();
526			rsh->setL(atoi(token));
527	chrisfen	1646	token = strtok(NULL, delim);
528	gezelter	1650	rsh->setM(atoi(token));
529	chrisfen	1646	token = strtok(NULL, delim);
530			if (!strcasecmp("sin",token))
531	gezelter	1650	rsh->makeSinFunction();
532	chrisfen	1646	else
533	gezelter	1650	rsh->makeCosFunction();
534	chrisfen	1646	token = strtok(NULL, delim);
535	gezelter	1650	rsh->setCoefficient(atof(token));
536	chrisfen	1646
537	gezelter	1650	functionVector.push_back(rsh);
538	chrisfen	1646	}
539			}
540	gezelter	1520	}
541	gezelter	1672	shapeFile.getline(inLine, MAXLEN);
542	gezelter	1520	}
543
544	chrisfen	1646	// pass contact functions to ShapeType
545			st->setContactFuncs(functionVector);
546	gezelter	1520
547	chrisfen	1646	// now grab the range functions
548			findBegin(shapeFile, "RangeFunctions");
549			functionVector.clear();
550	gezelter	1520
551	chrisfen	1646	shapeFile.getline(inLine, MAXLEN);
552			while( !shapeFile.eof() ) {
553			// toss comment lines
554			if( inLine[0] != '!' && inLine[0] != '#' ){
555			// end marks section completion
556			if (isEndLine(inLine)) break;
557
558	gezelter	1650	nTokens = countTokens(inLine, delim);
559	chrisfen	1646	if (nTokens != 0) {
560			if (nTokens < 4) {
561			sprintf( painCave.errMsg,
562	gezelter	1650	"Not enough data on a RangeFunctions line in file: %s\n",
563			shapeFileName.c_str() );
564	chrisfen	1646	painCave.severity = OOPSE_ERROR;
565			painCave.isFatal = 1;
566			simError();
567
568			// read in a spherical harmonic function
569			token = strtok(inLine, delim);
570	gezelter	1650
571			rsh = new RealSphericalHarmonic();
572			rsh->setL(atoi(token));
573	chrisfen	1646	token = strtok(NULL, delim);
574	gezelter	1650	rsh->setM(atoi(token));
575	chrisfen	1646	token = strtok(NULL, delim);
576			if (!strcasecmp("sin",token))
577	gezelter	1650	rsh->makeSinFunction();
578	chrisfen	1646	else
579	gezelter	1650	rsh->makeCosFunction();
580	chrisfen	1646	token = strtok(NULL, delim);
581	gezelter	1650	rsh->setCoefficient(atof(token));
582	chrisfen	1646
583	gezelter	1650	functionVector.push_back(rsh);
584	chrisfen	1646	}
585			}
586	gezelter	1520	}
587	gezelter	1672	shapeFile.getline(inLine, MAXLEN);
588	chrisfen	1646	}
589	gezelter	1520
590	chrisfen	1646	// pass range functions to ShapeType
591			st->setRangeFuncs(functionVector);
592	gezelter	1520
593	chrisfen	1646	// finally grab the strength functions
594			findBegin(shapeFile, "StrengthFunctions");
595			functionVector.clear();
596	gezelter	1520
597	chrisfen	1646	shapeFile.getline(inLine, MAXLEN);
598			while( !shapeFile.eof() ) {
599			// toss comment lines
600			if( inLine[0] != '!' && inLine[0] != '#' ){
601			// end marks section completion
602			if (isEndLine(inLine)) break;
603
604	gezelter	1650	nTokens = countTokens(inLine, delim);
605	chrisfen	1646	if (nTokens != 0) {
606			if (nTokens < 4) {
607			sprintf( painCave.errMsg,
608	gezelter	1650	"Not enough data on a StrengthFunctions line in file: %s\n",
609			shapeFileName.c_str() );
610	chrisfen	1646	painCave.severity = OOPSE_ERROR;
611			painCave.isFatal = 1;
612			simError();
613
614			// read in a spherical harmonic function
615			token = strtok(inLine, delim);
616	gezelter	1650	rsh = new RealSphericalHarmonic();
617			rsh->setL(atoi(token));
618	chrisfen	1646	token = strtok(NULL, delim);
619	gezelter	1650	rsh->setM(atoi(token));
620	chrisfen	1646	token = strtok(NULL, delim);
621			if (!strcasecmp("sin",token))
622	gezelter	1650	rsh->makeSinFunction();
623	chrisfen	1646	else
624	gezelter	1650	rsh->makeCosFunction();
625	chrisfen	1646	token = strtok(NULL, delim);
626	gezelter	1650	rsh->setCoefficient(atof(token));
627	chrisfen	1646
628	gezelter	1650	functionVector.push_back(rsh);
629	chrisfen	1646	}
630			}
631	gezelter	1520	}
632	gezelter	1672	shapeFile.getline(inLine, MAXLEN);
633	gezelter	1520	}
634
635	chrisfen	1646	// pass strength functions to ShapeType
636			st->setStrengthFuncs(functionVector);
637
638			// we're done reading from this file
639			shapeFile.close();
640	gezelter	1520	}
641	chrisfen	1646
642	gezelter	1650	double Shapes_FF::findLargestContactDistance(ShapeAtomType* st) {
643			int i, j, nSteps;
644			double theta, thetaStep, thetaMin, costheta;
645			double phi, phiStep;
646			double sigma, bigSigma;
647
648			nSteps = 16;
649
650			thetaStep = M_PI / nSteps;
651			thetaMin = thetaStep / 2.0;
652			phiStep = thetaStep * 2.0;
653			bigSigma = 0.0;
654
655			for (i = 0; i < nSteps; i++) {
656
657			theta = thetaMin + i * thetaStep;
658			costheta = cos(theta);
659
660			for (j = 0; j < nSteps; j++) {
661
662			phi = j*phiStep;
663
664			sigma = st->getContactValueAt(costheta, phi);
665
666			if (sigma > bigSigma) bigSigma = sigma;
667			}
668			}
669
670			return bigSigma;
671			}