--- trunk/OOPSE/libmdtools/EAM_FF.cpp 2003/07/16 21:49:59 627 +++ trunk/OOPSE/libmdtools/EAM_FF.cpp 2003/10/23 19:57:25 814 @@ -30,10 +30,10 @@ namespace EAM_NS{ double mass; double lattice_constant; double eam_drho; // The distance between each of the points indexed by rho. + double eam_rcut; // The cutoff radius for eam. 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. 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; @@ -75,8 +76,6 @@ namespace EAM_NS{ void add( atomStruct &info, double *the_eam_rvals, double *the_eam_rhovals,double *the_eam_Frhovals ){ - int i; - // check for duplicates if( !strcmp( info.name, name ) ){ @@ -135,7 +134,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 +149,7 @@ namespace EAM_NS{ } -using namespace LJ_NS; +using namespace EAM_NS; //**************************************************************** // begins the actual forcefield stuff. @@ -163,11 +162,13 @@ EAM_FF::EAM_FF(){ char* ffPath_env = "FORCE_PARAM_PATH"; char* ffPath; char temp[200]; - char errMsg[1000]; headAtomType = NULL; currentAtomType = NULL; + // Set eamRcut to 0.0 + eamRcut = 0.0; + // do the funtion wrapping wrapMeFF( this ); @@ -178,7 +179,7 @@ EAM_FF::EAM_FF(){ // Init the atomStruct mpi type atomStruct atomProto; // mpiPrototype - int atomBC[3] = {15,4,6}; // block counts + int atomBC[3] = {15,5,5}; // block counts MPI_Aint atomDspls[3]; // displacements MPI_Datatype atomMbrTypes[3]; // member mpi types @@ -263,8 +264,19 @@ EAM_FF::~EAM_FF(){ #endif // is_mpi } -void EAM_FF::initForceField( int ljMixRule ){ + +void EAM_FF::calcRcut( void ){ + #ifdef IS_MPI + double tempEamRcut = eamRcut; + MPI_Allreduce( &tempEamRcut, &eamRcut, 1, MPI_DOUBLE, MPI_MAX, + MPI_COMM_WORLD); +#endif //is_mpi + entry_plug->setRcut(eamRcut); +} + + +void EAM_FF::initForceField( int ljMixRule ){ initFortran( ljMixRule, 0 ); } @@ -283,13 +295,13 @@ void EAM_FF::cleanMe( void ){ #endif // is_mpi } + void EAM_FF::readParams( void ){ atomStruct info; info.last = 1; // initialize last to have the last set. // if things go well, last will be set to 0 - int i; int identNum; double *eam_rvals; // Z of r values double *eam_rhovals; // rho of r values @@ -335,7 +347,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 +449,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 +465,8 @@ void EAM_FF::readParams( void ){ &isDipole, &isGB, &isEAM, - &(currentAtomType->epslon), - &(currentAtomType->sigma), + &eamEpslon, + &eamSigma, &dipole, &isError ); if( isError ){ @@ -467,23 +481,26 @@ void EAM_FF::readParams( void ){ } entry_plug->useLJ = 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", @@ -503,18 +520,17 @@ void EAM_FF::readParams( void ){ MPIcheckPoint(); #endif // is_mpi + + } void EAM_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){ int i; - + // initialize the atoms - - Atom* thisAtom; - for( i=0; ifind( the_atoms[i]->getType() ); @@ -529,11 +545,14 @@ void EAM_FF::initializeAtoms( int nAtoms, Atom** the_a the_atoms[i]->setMass( currentAtomType->mass ); the_atoms[i]->setIdent( currentAtomType->ident ); the_atoms[i]->setEAM(); + the_atoms[i]->setEamRcut( currentAtomType->eam_rcut); + if (eamRcut < currentAtomType->eam_rcut) eamRcut = currentAtomType->eam_rcut; + } } -void LJ_FF::initializeBonds( int nBonds, Bond** BondArray, +void EAM_FF::initializeBonds( int nBonds, Bond** BondArray, bond_pair* the_bonds ){ if( nBonds ){ @@ -657,15 +676,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 +701,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 +720,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 +733,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 +745,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 +757,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 +783,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 +833,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