--- trunk/OOPSE-2.0/src/UseTheForce/Shapes_FF.cpp 2004/10/22 22:54:01 1636 +++ trunk/OOPSE-2.0/src/UseTheForce/Shapes_FF.cpp 2004/10/28 17:20:22 1672 @@ -1,9 +1,18 @@ +/* + * Shapes_FF.cpp + * oopse + * + * Created by Chris Fennell on 10/20/04. + * Copyright 2004 __MyCompanyName__. All rights reserved. + * + */ + #include #include #include - +#include +#include #include -using namespace std; #ifdef IS_MPI #include @@ -12,104 +21,28 @@ using namespace std; #include "UseTheForce/ForceFields.hpp" #include "primitives/SRI.hpp" #include "utils/simError.h" - -#include "UseTheForce/fortranWrappers.hpp" +#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 -Shapes_FF::Shapes_FF() { - Shapes_FF(""); -} +using namespace std; +using namespace oopse; -Shapes_FF::Shapes_FF(char* the_variant){ - - char fileName[200]; - char* ffPath_env = "FORCE_PARAM_PATH"; - char* ffPath; - char temp[200]; - - // do the funtion wrapping - wrapMeFF( this ); - -#ifdef IS_MPI - if( worldRank == 0 ){ -#endif - - // generate the force file name - - strcpy( fileName, "Shapes" ); - - if (strlen(the_variant) > 0) { - has_variant = 1; - strcpy( variant, the_variant); - strcat( fileName, "."); - strcat( fileName, variant ); - - sprintf( painCave.errMsg, - "Using %s variant of Shapes force field.\n", - variant ); - painCave.severity = OOPSE_INFO; - painCave.isFatal = 0; - simError(); - } - strcat( fileName, ".frc"); - - // attempt to open the file in the current directory first. - - frcFile = fopen( fileName, "r" ); - - if( frcFile == 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 ); - - frcFile = fopen( fileName, "r" ); - - 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(); - } - } - - -#ifdef IS_MPI - } - - sprintf( checkPointMsg, "Shapes_FF file opened sucessfully." ); - MPIcheckPoint(); - -#endif // is_mpi -} - - Shapes_FF::~Shapes_FF(){ #ifdef IS_MPI if( worldRank == 0 ){ #endif // is_mpi - fclose( frcFile ); + forceFile.close(); #ifdef IS_MPI } @@ -135,758 +68,572 @@ void Shapes_FF::readParams( void ){ void Shapes_FF::readParams( void ){ - char shapePotFile[1000]; - -#ifdef IS_MPI - if( worldRank == 0 ){ -#endif - - // read in the atom types. + char readLine[1024]; - fastForward( "AtomTypes", "eam atom readParams" ); + string fileName; + string shapeFileName; + string tempString; - // we are now at the AtomTypes section. + 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::iterator iter; + + // vectors for shape transfer to fortran + vector tempSHVector; + vector contactL; + vector contactM; + vector contactFunc; + vector contactCoeff; + vector rangeL; + vector rangeM; + vector rangeFunc; + vector rangeCoeff; + vector strengthL; + vector strengthM; + vector strengthFunc; + vector 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 ){ - eof_test = fgets( readLine, sizeof(readLine), frcFile ); - lineNum++; + tempString = ffPath; + tempString += "/"; + tempString += fileName; + fileName = tempString; + forceFile.open( fileName.c_str() ); - // read a line, and start parseing out the atom types - - if( eof_test == NULL ){ - sprintf( painCave.errMsg, - "Error in reading Atoms from force file at line %d.\n", - lineNum ); + 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(); } - - identNum = 1; - // stop reading at end of file, or at next section + } + + // read in the shape types. + + findBegin( forceFile, "ShapeTypes" ); + + while( !forceFile.eof() ){ + forceFile.getline( readLine, sizeof(readLine) ); - while( readLine[0] != '#' && eof_test != NULL ){ - - // toss comment lines - if( readLine[0] != '!' ){ - - // the parser returns 0 if the line was blank - if( parseAtom( readLine, lineNum, info, eamPotFile ) ){ - parseEAM(info,eamPotFile, &eam_rvals, - &eam_rhovals, &eam_Frhovals); - info.ident = identNum; - headAtomType->add( info, eam_rvals, - eam_rhovals,eam_Frhovals ); - identNum++; + // 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(); } - } - eof_test = fgets( readLine, sizeof(readLine), frcFile ); - lineNum++; - } - - + + nameToken = strtok( readLine, delim ); + shapeFileName = strtok( NULL, delim ); -#ifdef IS_MPI - - - // send out the linked list to all the other processes + // strings are not char arrays! + nameString = nameToken; - sprintf( checkPointMsg, - "Shapes_FF atom structures read successfully." ); - MPIcheckPoint(); + // 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(); - currentAtomType = headAtomType->next; //skip the first element who is a place holder. - while( currentAtomType != NULL ){ - currentAtomType->duplicate( info ); + 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: - - - sendFrcStruct( &info, mpiAtomStructType ); - - // We have to now broadcast the Arrays - MPI_Bcast(currentAtomType->eam_rvals, - currentAtomType->eam_nr, - MPI_DOUBLE,0,MPI_COMM_WORLD); - MPI_Bcast(currentAtomType->eam_rhovals, - currentAtomType->eam_nr, - MPI_DOUBLE,0,MPI_COMM_WORLD); - MPI_Bcast(currentAtomType->eam_Frhovals, - currentAtomType->eam_nrho, - MPI_DOUBLE,0,MPI_COMM_WORLD); - - sprintf( checkPointMsg, - "successfully sent EAM force type: \"%s\"\n", - info.name ); - MPIcheckPoint(); - - currentAtomType = currentAtomType->next; + st = (ShapeAtomType*)(iter->second); + parseShapeFile(shapeFileName, st); + } + } } - info.last = 1; - sendFrcStruct( &info, mpiAtomStructType ); - } - else{ - - // listen for node 0 to send out the force params - - MPIcheckPoint(); +#ifdef IS_MPI - headAtomType = new LinkedAtomType; - receiveFrcStruct( &info, mpiAtomStructType ); + // looks like all the processors have their ShapeType vectors ready... + sprintf( checkPointMsg, + "Shapes_FF shape objects read successfully." ); + MPIcheckPoint(); - while( !info.last ){ - - // allocate the arrays +#endif // is_mpi - eam_rvals = new double[info.eam_nr]; - eam_rhovals = new double[info.eam_nr]; - eam_Frhovals = new double[info.eam_nrho]; + // pack up and send the necessary info to fortran + for (iter = atomTypeMap.begin(); iter != atomTypeMap.end(); ++iter){ - // We have to now broadcast the Arrays - MPI_Bcast(eam_rvals, - info.eam_nr, - MPI_DOUBLE,0,MPI_COMM_WORLD); - MPI_Bcast(eam_rhovals, - info.eam_nr, - MPI_DOUBLE,0,MPI_COMM_WORLD); - MPI_Bcast(eam_Frhovals, - info.eam_nrho, - MPI_DOUBLE,0,MPI_COMM_WORLD); - + at = (AtomType*)(iter->second); - headAtomType->add( info, eam_rvals, eam_rhovals, eam_Frhovals ); - - MPIcheckPoint(); + if (at->isDirectional()) { + + dat = (DirectionalAtomType*)at; - receiveFrcStruct( &info, mpiAtomStructType ); + if (dat->isShape()) { - - } - } -#endif // is_mpi - - // call new A_types in fortran - - int isError; - - // dummy variables - int isLJ = 0; - int isDipole = 0; - int isSSD = 0; - int isGB = 0; - int isEAM = 1; - int isCharge = 0; - double dipole = 0.0; - double charge = 0.0; - double eamSigma = 0.0; - double eamEpslon = 0.0; - - currentAtomType = headAtomType->next; - while( currentAtomType != NULL ){ - - if( currentAtomType->name[0] != '\0' ){ - isError = 0; - makeAtype( &(currentAtomType->ident), - &isLJ, - &isSSD, - &isDipole, - &isGB, - &isEAM, - &isCharge, - &eamEpslon, - &eamSigma, - &charge, - &dipole, - &isError ); - if( isError ){ - sprintf( painCave.errMsg, - "Error initializing the \"%s\" atom type in fortran\n", - currentAtomType->name ); - painCave.isFatal = 1; - simError(); + st = (ShapeAtomType*)at; + + contactL.clear(); + contactM.clear(); + contactFunc.clear(); + contactCoeff.clear(); + + tempSHVector = st->getContactFuncs(); + + nContact = tempSHVector.size(); + for (i=0; igetL()); + 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; igetL()); + 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; igetL()); + 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(); + } } } - currentAtomType = currentAtomType->next; } - - entry_plug->useLennardJones = 0; - entry_plug->useEAM = 1; - // Walk down again and send out EAM type - currentAtomType = headAtomType->next; - while( currentAtomType != NULL ){ - - if( currentAtomType->name[0] != '\0' ){ - isError = 0; - - newEAMtype( &(currentAtomType->lattice_constant), - &(currentAtomType->eam_nrho), - &(currentAtomType->eam_drho), - &(currentAtomType->eam_nr), - &(currentAtomType->eam_dr), - &(currentAtomType->eam_rcut), - currentAtomType->eam_rvals, - currentAtomType->eam_rhovals, - currentAtomType->eam_Frhovals, - &(currentAtomType->eam_ident), - &isError); - - if( isError ){ - sprintf( painCave.errMsg, - "Error initializing the \"%s\" atom type in fortran EAM\n", - currentAtomType->name ); - painCave.isFatal = 1; - simError(); - } - } - currentAtomType = currentAtomType->next; - } - - - + #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; - + int i,j,k; + map::iterator iter; + // initialize the atoms - + DirectionalAtom* dAtom; + AtomType* at; + DirectionalAtomType* dat; + ShapeAtomType* sat; + double sigma; + double ji[3]; + double inertialMat[3][3]; + Mat3x3d momInt; + string myTypeString; + for( i=0; ifind( the_atoms[i]->getType() ); - if( currentAtomType == NULL ){ + 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(); - } - - the_atoms[i]->setMass( currentAtomType->mass ); - the_atoms[i]->setIdent( currentAtomType->ident ); + } else { - if (eamRcut < currentAtomType->eam_rcut) eamRcut = currentAtomType->eam_rcut; - - } -} + at = (AtomType*)(iter->second); -void Shapes_FF::initializeBonds( int nBonds, Bond** BondArray, - bond_pair* the_bonds ){ - - if( nBonds ){ - sprintf( painCave.errMsg, - "LJ_FF does not support bonds.\n" ); - painCave.isFatal = 1; - simError(); - } -} + 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; -void Shapes_FF::initializeBends( int nBends, Bend** bendArray, - bend_set* the_bends ){ + } - if( nBends ){ - sprintf( painCave.errMsg, - "LJ_FF does not support bends.\n" ); - painCave.isFatal = 1; - simError(); - } -} + the_atoms[i]->setHasCharge(at->isCharge()); -void Shapes_FF::initializeTorsions( int nTorsions, Torsion** torsionArray, - torsion_set* the_torsions ){ + if( at->isDirectional() ){ - if( nTorsions ){ - sprintf( painCave.errMsg, - "LJ_FF does not support torsions.\n" ); - painCave.isFatal = 1; - simError(); - } -} + dat = (DirectionalAtomType*)at; + dAtom = (DirectionalAtom *) the_atoms[i]; + + momInt = dat->getI(); -void Shapes_FF::fastForward( char* stopText, char* searchOwner ){ + // 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 ); - int foundText = 0; - char* the_token; + ji[0] = 0.0; + ji[1] = 0.0; + ji[2] = 0.0; + dAtom->setJ( ji ); - rewind( frcFile ); - lineNum = 0; + if (dat->isDipole()) { + dAtom->setHasDipole( dat->isDipole() ); + entry_plug->n_dipoles++; + } + } + } + } +} - eof_test = fgets( readLine, sizeof(readLine), frcFile ); - lineNum++; - if( eof_test == NULL ){ - sprintf( painCave.errMsg, "Error fast forwarding force file for %s: " - " file is empty.\n", - searchOwner ); +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 ){ - - while( !foundText ){ - while( eof_test != NULL && readLine[0] != '#' ){ - eof_test = fgets( readLine, sizeof(readLine), frcFile ); - lineNum++; - } - if( eof_test == NULL ){ - sprintf( painCave.errMsg, - "Error fast forwarding force file for %s at " - "line %d: file ended unexpectedly.\n", - searchOwner, - lineNum ); - painCave.isFatal = 1; - simError(); - } - - the_token = strtok( readLine, " ,;\t#\n" ); - foundText = !strcmp( stopText, the_token ); - - if( !foundText ){ - eof_test = fgets( readLine, sizeof(readLine), frcFile ); - lineNum++; - - if( eof_test == NULL ){ - sprintf( painCave.errMsg, - "Error fast forwarding force file for %s at " - "line %d: file ended unexpectedly.\n", - searchOwner, - lineNum ); - painCave.isFatal = 1; - simError(); - } - } - } + if( nBends ){ + sprintf( painCave.errMsg, + "Shapes_FF does not support bends.\n" ); + painCave.isFatal = 1; + simError(); + } } - - -int EAM_NS::parseAtom( char *lineBuffer, int lineNum, atomStruct &info, char *eamPotFile ){ - - char* the_token; +void Shapes_FF::initializeTorsions( int nTorsions, Torsion** torsionArray, + torsion_set* the_torsions ){ - the_token = strtok( lineBuffer, " \n\t,;" ); - if( the_token != NULL ){ - - strcpy( info.name, the_token ); - - if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ - sprintf( painCave.errMsg, - "Error parseing AtomTypes: line %d\n", lineNum ); - painCave.isFatal = 1; - simError(); - } - - info.mass = atof( the_token ); - - if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ - sprintf( painCave.errMsg, - "Error parseing AtomTypes: line %d\n", lineNum ); - painCave.isFatal = 1; - simError(); - } - - strcpy( eamPotFile, the_token ); - return 1; + if( nTorsions ){ + sprintf( painCave.errMsg, + "Shapes_FF does not support torsions.\n" ); + painCave.isFatal = 1; + simError(); } - else return 0; } -int EAM_NS::parseEAM(atomStruct &info, char *eamPotFile, - double **eam_rvals, - double **eam_rhovals, - double **eam_Frhovals){ - double* myEam_rvals; - double* myEam_rhovals; - double* myEam_Frhovals; +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 functionVector; + ifstrstream shapeFile; + string tempString; - char* ffPath_env = "FORCE_PARAM_PATH"; - char* ffPath; - char* the_token; - char* eam_eof_test; - FILE *eamFile; - const int BUFFERSIZE = 3000; - - char temp[200]; - int linenumber; - int nReadLines; - char eam_read_buffer[BUFFERSIZE]; - - - int i,j; - - linenumber = 0; - - // Open eam file - eamFile = fopen( eamPotFile, "r" ); + shapeFile.open( shapeFileName.c_str() ); - - if( eamFile == NULL ){ + if( shapeFile == NULL ){ - // next see if the force path enviorment variable is set + tempString = ffPath; + tempString += "/"; + tempString += shapeFileName; + shapeFileName = tempString; + + shapeFile.open( shapeFileName.c_str() ); - ffPath = getenv( ffPath_env ); - if( ffPath == NULL ) { - STR_DEFINE(ffPath, FRC_PATH ); - } - - - strcpy( temp, ffPath ); - strcat( temp, "/" ); - strcat( temp, eamPotFile ); - strcpy( eamPotFile, temp ); - - eamFile = fopen( eamPotFile, "r" ); - - - - if( eamFile == NULL ){ + if( shapeFile == NULL ){ sprintf( painCave.errMsg, - "Error opening the EAM force parameter file: %s\n" - "Have you tried setting the FORCE_PARAM_PATH environment " + "Error opening the shape file:\n" + "\t%s\n" + "\tHave you tried setting the FORCE_PARAM_PATH environment " "variable?\n", - eamPotFile ); + 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(); - // First line is a comment line, read and toss it.... - eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile); - linenumber++; - if(eam_eof_test == NULL){ - sprintf( painCave.errMsg, - "error in reading commment in %s\n", eamPotFile); - 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); } - - - // The Second line contains atomic number, atomic mass and a lattice constant - eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile); - linenumber++; - if(eam_eof_test == NULL){ - sprintf( painCave.errMsg, - "error in reading Identifier line in %s\n", eamPotFile); - painCave.isFatal = 1; - simError(); - } - - - - - if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){ - sprintf( painCave.errMsg, - "Error parsing EAM ident line in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } + // now grab the contact functions + findBegin(shapeFile, "ContactFunctions"); + functionVector.clear(); - info.eam_ident = atoi( the_token ); - - if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ - sprintf( painCave.errMsg, - "Error parsing EAM mass in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); + 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); } - info.mass = atof( the_token); - - if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ - sprintf( painCave.errMsg, - "Error parsing EAM Lattice Constant %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } - info.lattice_constant = atof( the_token); - - // Next line is nrho, drho, nr, dr and rcut - eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile); - if(eam_eof_test == NULL){ - sprintf( painCave.errMsg, - "error in reading number of points line in %s\n", eamPotFile); - painCave.isFatal = 1; - simError(); - } - if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){ - sprintf( painCave.errMsg, - "Error parseing EAM nrho: line in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } - - info.eam_nrho = atoi( the_token ); - - if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ - sprintf( painCave.errMsg, - "Error parsing EAM drho in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } - info.eam_drho = atof( the_token); + // pass contact functions to ShapeType + st->setContactFuncs(functionVector); - if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ - sprintf( painCave.errMsg, - "Error parsing EAM # r in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } - info.eam_nr = atoi( the_token); + // now grab the range functions + findBegin(shapeFile, "RangeFunctions"); + functionVector.clear(); - if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ - sprintf( painCave.errMsg, - "Error parsing EAM dr in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } - info.eam_dr = atof( the_token); + 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); - if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ - sprintf( painCave.errMsg, - "Error parsing EAM rcut in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } - info.eam_rcut = atof( the_token); - - - - - - // Ok now we have to allocate point arrays and read in number of points - // Index the arrays for fortran, starting at 1 - myEam_Frhovals = new double[info.eam_nrho]; - myEam_rvals = new double[info.eam_nr]; - myEam_rhovals = new double[info.eam_nr]; - - // Parse F of rho vals. - - // Assume for now that we have a complete number of lines - nReadLines = int(info.eam_nrho/5); - - - - for (i=0;isetL(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); + } + } } - - // Parse 5 values on each line into array - // Value 1 - if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){ - sprintf( painCave.errMsg, - "Error parsing EAM nrho: line in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } - - myEam_Frhovals[j+0] = atof( the_token ); - - // Value 2 - if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){ - sprintf( painCave.errMsg, - "Error parsing EAM nrho: line in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } - - myEam_Frhovals[j+1] = atof( the_token ); - - // Value 3 - if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){ - sprintf( painCave.errMsg, - "Error parsing EAM nrho: line in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } - - myEam_Frhovals[j+2] = atof( the_token ); - - // Value 4 - if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){ - sprintf( painCave.errMsg, - "Error parsing EAM nrho: line in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } - - myEam_Frhovals[j+3] = atof( the_token ); - - // Value 5 - if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){ - sprintf( painCave.errMsg, - "Error parsing EAM nrho: line in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } - - myEam_Frhovals[j+4] = atof( the_token ); - + shapeFile.getline(inLine, MAXLEN); } - // Parse Z of r vals - - // Assume for now that we have a complete number of lines - nReadLines = int(info.eam_nr/5); - for (i=0;isetRangeFuncs(functionVector); - // Read next line - eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile); - linenumber++; - if(eam_eof_test == NULL){ - sprintf( painCave.errMsg, - "error in reading EAM file %s at line %d\n", - eamPotFile,linenumber); - painCave.isFatal = 1; - simError(); - } - - // Parse 5 values on each line into array - // Value 1 - if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){ - sprintf( painCave.errMsg, - "Error parsing EAM nrho: line in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } - - myEam_rvals[j+0] = atof( the_token ); - - // Value 2 - if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){ - sprintf( painCave.errMsg, - "Error parsing EAM nrho: line in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } + // finally grab the strength functions + findBegin(shapeFile, "StrengthFunctions"); + functionVector.clear(); - myEam_rvals[j+1] = atof( the_token ); - - // Value 3 - if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){ - sprintf( painCave.errMsg, - "Error parsing EAM nrho: line in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); + 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); + } + } } - - myEam_rvals[j+2] = atof( the_token ); + shapeFile.getline(inLine, MAXLEN); + } - // Value 4 - if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){ - sprintf( painCave.errMsg, - "Error parsing EAM nrho: line in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } - - myEam_rvals[j+3] = atof( the_token ); + // pass strength functions to ShapeType + st->setStrengthFuncs(functionVector); - // Value 5 - if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){ - sprintf( painCave.errMsg, - "Error parsing EAM nrho: line in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } - - myEam_rvals[j+4] = atof( the_token ); + // we're done reading from this file + shapeFile.close(); +} - } - // Parse rho of r vals +double Shapes_FF::findLargestContactDistance(ShapeAtomType* st) { + int i, j, nSteps; + double theta, thetaStep, thetaMin, costheta; + double phi, phiStep; + double sigma, bigSigma; - // Assume for now that we have a complete number of lines + nSteps = 16; - for (i=0;igetContactValueAt(costheta, phi); + + if (sigma > bigSigma) bigSigma = sigma; } - - myEam_rhovals[j+3] = atof( the_token ); - - // Value 5 - if ( (the_token = strtok( NULL, " \n\t,;")) == NULL){ - sprintf( painCave.errMsg, - "Error parsing EAM nrho: line in %s\n", eamPotFile ); - painCave.isFatal = 1; - simError(); - } - - myEam_rhovals[j+4] = atof( the_token ); - } - *eam_rvals = myEam_rvals; - *eam_rhovals = myEam_rhovals; - *eam_Frhovals = myEam_Frhovals; - fclose(eamFile); - return 0; + return bigSigma; }