--- trunk/OOPSE/libmdtools/EAM_FF.cpp 2003/07/16 21:49:59 627 +++ trunk/OOPSE/libmdtools/EAM_FF.cpp 2003/07/30 21:17:01 657 @@ -33,7 +33,7 @@ namespace EAM_NS{ double eam_dr; // The distance between each of the rho points. int eam_nrho; // Number of points indexed by rho int eam_nr; // The number of points based on r (Both Phi(r) and Rho(r)). - int eam_rcut; // The cutoff radius for eam. + double eam_rcut; // The cutoff radius for eam. int eam_ident; // Atomic number int ident; int last; // 0 -> default @@ -41,7 +41,8 @@ namespace EAM_NS{ } atomStruct; int parseAtom( char *lineBuffer, int lineNum, atomStruct &info, char *eamPotFile ); - int parseEAM( atomStruct &info, char *eamPotFile ); + int parseEAM( atomStruct &info, char *eamPotFile, double **eam_rvals, + double **eam_rhovals, double **eam_Frhovals); #ifdef IS_MPI MPI_Datatype mpiAtomStructType; @@ -135,7 +136,7 @@ namespace EAM_NS{ double eam_drho; // The distance between each of the points indexed by rho. int eam_nr; // The number of points based on r (Both Phi(r) and Rho(r)). double eam_dr; // The distance between each of the rho points. - int eam_rcut; // The cutoff radius for eam. + double eam_rcut; // The cutoff radius for eam. double *eam_rvals; // Z of r values double *eam_rhovals; // rho of r values @@ -150,7 +151,7 @@ namespace EAM_NS{ } -using namespace LJ_NS; +using namespace EAM_NS; //**************************************************************** // begins the actual forcefield stuff. @@ -261,10 +262,17 @@ EAM_FF::~EAM_FF(){ #ifdef IS_MPI } #endif // is_mpi +} + + +void EAM_FF::calcRcut( void ){ + double tempEamRcut; + + entry_plug->setRcut(eamRcut); } + void EAM_FF::initForceField( int ljMixRule ){ - initFortran( ljMixRule, 0 ); } @@ -283,6 +291,7 @@ void EAM_FF::cleanMe( void ){ #endif // is_mpi } + void EAM_FF::readParams( void ){ atomStruct info; @@ -335,7 +344,7 @@ void EAM_FF::readParams( void ){ // the parser returns 0 if the line was blank if( parseAtom( readLine, lineNum, info, eamPotFile ) ){ parseEAM(info,eamPotFile, &eam_rvals, - &eam_rhovals, &eam_Frhovals)){ + &eam_rhovals, &eam_Frhovals); info.ident = identNum; headAtomType->add( info, eam_rvals, eam_rhovals,eam_Frhovals ); @@ -437,8 +446,10 @@ void EAM_FF::readParams( void ){ int isDipole = 0; int isSSD = 0; int isGB = 0; - int isEam = 1; + int isEAM= 1; double dipole = 0.0; + double eamSigma = 0.0; + double eamEpslon = 0.0; currentAtomType = headAtomType->next; while( currentAtomType != NULL ){ @@ -451,8 +462,8 @@ void EAM_FF::readParams( void ){ &isDipole, &isGB, &isEAM, - &(currentAtomType->epslon), - &(currentAtomType->sigma), + &eamEpslon, + &eamSigma, &dipole, &isError ); if( isError ){ @@ -474,16 +485,19 @@ void EAM_FF::readParams( void ){ if( currentAtomType->name[0] != '\0' ){ isError = 0; + cerr << "Calling newEAMtype for type "<eam_ident <<"\n"; 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); + cerr << "Returned from newEAMtype\n"; if( isError ){ sprintf( painCave.errMsg, "Error initializing the \"%s\" atom type in fortran EAM\n", @@ -503,6 +517,8 @@ void EAM_FF::readParams( void ){ MPIcheckPoint(); #endif // is_mpi + cerr << "Done sending eamtypes to fortran\n"; + } @@ -533,7 +549,7 @@ void EAM_FF::initializeAtoms( int nAtoms, Atom** the_a } } -void LJ_FF::initializeBonds( int nBonds, Bond** BondArray, +void EAM_FF::initializeBonds( int nBonds, Bond** BondArray, bond_pair* the_bonds ){ if( nBonds ){ @@ -657,15 +673,23 @@ int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil double **eam_rvals, double **eam_rhovals, double **eam_Frhovals){ - + double* myEam_rvals; + double* myEam_rhovals; + double* myEam_Frhovals; + 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; @@ -674,7 +698,7 @@ int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil eamFile = fopen( eamPotFile, "r" ); - if( frcFile == NULL ){ + if( eamFile == NULL ){ // next see if the force path enviorment variable is set @@ -693,7 +717,7 @@ int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil - if( frcFile == NULL ){ + if( eamFile == NULL ){ sprintf( painCave.errMsg, "Error opening the EAM force parameter file: %s\n" @@ -706,9 +730,9 @@ int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil } // First line is a comment line, read and toss it.... - eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile); + eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile); linenumber++; - if(eof_test == NULL){ + if(eam_eof_test == NULL){ sprintf( painCave.errMsg, "error in reading commment in %s\n", eamPotFile); painCave.isFatal = 1; @@ -718,9 +742,9 @@ int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil // The Second line contains atomic number, atomic mass and a lattice constant - eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile); + eam_eof_test = fgets(eam_read_buffer, sizeof(eam_read_buffer),eamFile); linenumber++; - if(eof_test == NULL){ + if(eam_eof_test == NULL){ sprintf( painCave.errMsg, "error in reading Identifier line in %s\n", eamPotFile); painCave.isFatal = 1; @@ -730,7 +754,7 @@ int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil - if ( (the_token = strtok( read_buffer, " \n\t,;")) == NULL){ + if ( (the_token = strtok( eam_read_buffer, " \n\t,;")) == NULL){ sprintf( painCave.errMsg, "Error parseing EAM ident line in %s\n", eamPotFile ); painCave.isFatal = 1; @@ -756,15 +780,15 @@ int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil info.lattice_constant = atof( the_token); // Next line is nrho, drho, nr, dr and rcut - eof_test = fgets(read_buffer, sizeof(read_buffer),eamFile) - if(eof_test == NULL){ + 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( read_buffer, " \n\t,;")) == NULL){ + 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; @@ -806,24 +830,29 @@ int EAM_NS::parseEAM(atomStruct &info, char *eamPotFil 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 - *eam_Frhovals = new double[info.nrho]; - *eam_rvals = new double[info.nr]; - *eam_rhovals = new double[info.nr]; + 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.nrho/5); + nReadLines = int(info.eam_nrho/5); + + for (i=0;i