--- trunk/OOPSE/libmdtools/EAM_FF.cpp 2003/04/15 16:37:59 499 +++ trunk/OOPSE/libmdtools/EAM_FF.cpp 2004/01/16 21:55:39 952 @@ -1,6 +1,6 @@ -#include -#include -#include +#include +#include +#include #include using namespace std; @@ -28,15 +28,21 @@ namespace EAM_NS{ typedef struct{ char name[15]; double mass; - double epslon; - double sigma; + 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_ident; // Atomic number int ident; int last; // 0 -> default // 1 -> in MPI: tells nodes to stop listening } atomStruct; - int parseAtom( char *lineBuffer, int lineNum, atomStruct &info ); - + int parseAtom( char *lineBuffer, int lineNum, 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; @@ -48,9 +54,18 @@ namespace EAM_NS{ LinkedAtomType(){ next = NULL; name[0] = '\0'; + eam_rvals = NULL; + eam_rhovals = NULL; + eam_Frhovals = NULL; } - ~LinkedAtomType(){ if( next != NULL ) delete next; } + ~LinkedAtomType(){ + if( next != NULL ) delete next; + if( eam_rvals != NULL ) delete[] eam_rvals; + if( eam_rhovals != NULL ) delete[] eam_rhovals; + if( eam_Frhovals != NULL ) delete[] eam_Frhovals; + } + LinkedAtomType* find(char* key){ if( !strcmp(name, key) ) return this; if( next != NULL ) return next->find(key); @@ -58,8 +73,9 @@ namespace EAM_NS{ } - void add( atomStruct &info ){ - + void add( atomStruct &info, double *the_eam_rvals, + double *the_eam_rhovals,double *the_eam_Frhovals ){ + // check for duplicates if( !strcmp( info.name, name ) ){ @@ -71,14 +87,23 @@ namespace EAM_NS{ simError(); } - if( next != NULL ) next->add(info); + if( next != NULL ) next->add(info, the_eam_rvals, the_eam_rhovals, the_eam_Frhovals); else{ next = new LinkedAtomType(); strcpy(next->name, info.name); - next->mass = info.mass; - next->epslon = info.epslon; - next->sigma = info.sigma; - next->ident = info.ident; + next->mass = info.mass; + next->lattice_constant = info.lattice_constant; + next->eam_nrho = info.eam_nrho; + next->eam_drho = info.eam_drho; + next->eam_nr = info.eam_nr; + next->eam_dr = info.eam_dr; + next->eam_rcut = info.eam_rcut; + next->eam_ident = info.eam_ident; + next->ident = info.ident; + + next->eam_rvals = the_eam_rvals; + next->eam_rhovals = the_eam_rhovals; + next->eam_Frhovals = the_eam_Frhovals; } } @@ -87,11 +112,16 @@ namespace EAM_NS{ void duplicate( atomStruct &info ){ strcpy(info.name, name); - info.mass = mass; - info.epslon = epslon; - info.sigma = sigma; - info.ident = ident; - info.last = 0; + info.mass = mass; + info.lattice_constant = lattice_constant; + info.eam_nrho = eam_nrho; + info.eam_drho = eam_drho; + info.eam_nr = eam_nr; + info.eam_dr = eam_dr; + info.eam_rcut = eam_rcut; + info.eam_ident = eam_ident; + info.ident = ident; + info.last = 0; } @@ -99,8 +129,17 @@ namespace EAM_NS{ char name[15]; double mass; - double epslon; - double sigma; + double lattice_constant; + int eam_nrho; // Number of points indexed by rho + 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. + double eam_rcut; // The cutoff radius for eam. + + double *eam_rvals; // Z of r values + double *eam_rhovals; // rho of r values + double *eam_Frhovals; // F of rho values + int eam_ident; // eam identity (atomic number) int ident; LinkedAtomType* next; }; @@ -110,7 +149,7 @@ using namespace LJ_NS; } -using namespace LJ_NS; +using namespace EAM_NS; //**************************************************************** // begins the actual forcefield stuff. @@ -123,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 ); @@ -138,13 +179,13 @@ EAM_FF::EAM_FF(){ // Init the atomStruct mpi type atomStruct atomProto; // mpiPrototype - int atomBC[3] = {15,3,2}; // block counts + int atomBC[3] = {15,5,5}; // block counts MPI_Aint atomDspls[3]; // displacements MPI_Datatype atomMbrTypes[3]; // member mpi types MPI_Address(&atomProto.name, &atomDspls[0]); MPI_Address(&atomProto.mass, &atomDspls[1]); - MPI_Address(&atomProto.ident, &atomDspls[2]); + MPI_Address(&atomProto.eam_nrho, &atomDspls[2]); atomMbrTypes[0] = MPI_CHAR; atomMbrTypes[1] = MPI_DOUBLE; @@ -201,7 +242,7 @@ EAM_FF::EAM_FF(){ #ifdef IS_MPI } - sprintf( checkPointMsg, "LJ_FF file opened sucessfully." ); + sprintf( checkPointMsg, "EAM_FF file opened sucessfully." ); MPIcheckPoint(); #endif // is_mpi @@ -223,8 +264,19 @@ void EAM_FF::initForceField( int ljMixRule ){ #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->setDefaultRcut(eamRcut); +} + + +void EAM_FF::initForceField( int ljMixRule ){ initFortran( ljMixRule, 0 ); } @@ -243,14 +295,18 @@ void EAM_FF::readParams( 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 + double *eam_Frhovals; // F of rho values + char eamPotFile[1000]; bigSigma = 0.0; @@ -262,7 +318,7 @@ void EAM_FF::readParams( void ){ headAtomType = new LinkedAtomType; - fastForward( "AtomTypes", "initializeAtoms" ); + fastForward( "AtomTypes", "eam atom readParams" ); // we are now at the AtomTypes section. @@ -282,28 +338,35 @@ void EAM_FF::readParams( void ){ identNum = 1; // stop reading at end of file, or at next section + 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 ) ){ + if( parseAtom( readLine, lineNum, info, eamPotFile ) ){ + parseEAM(info,eamPotFile, &eam_rvals, + &eam_rhovals, &eam_Frhovals); info.ident = identNum; - headAtomType->add( info );; + headAtomType->add( info, eam_rvals, + eam_rhovals,eam_Frhovals ); identNum++; } } eof_test = fgets( readLine, sizeof(readLine), frcFile ); lineNum++; } + + #ifdef IS_MPI + // send out the linked list to all the other processes sprintf( checkPointMsg, - "LJ_FF atom structures read successfully." ); + "EAM_FF atom structures read successfully." ); MPIcheckPoint(); currentAtomType = headAtomType->next; //skip the first element who is a place holder. @@ -314,8 +377,19 @@ void EAM_FF::readParams( void ){ 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 lJ force type: \"%s\"\n", + "successfully sent EAM force type: \"%s\"\n", info.name ); MPIcheckPoint(); @@ -334,16 +408,34 @@ void EAM_FF::readParams( void ){ headAtomType = new LinkedAtomType; recieveFrcStruct( &info, mpiAtomStructType ); - + while( !info.last ){ + + // allocate the arrays + eam_rvals = new double[info.eam_nr]; + eam_rhovals = new double[info.eam_nr]; + eam_Frhovals = new double[info.eam_nrho]; + // 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); + - headAtomType->add( info ); + headAtomType->add( info, eam_rvals, eam_rhovals, eam_Frhovals ); MPIcheckPoint(); recieveFrcStruct( &info, mpiAtomStructType ); + + } } #endif // is_mpi @@ -353,13 +445,18 @@ void EAM_FF::readParams( void ){ int isError; // dummy variables - int isLJ = 1; + 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; + currentAtomType = headAtomType->next; while( currentAtomType != NULL ){ if( currentAtomType->name[0] != '\0' ){ @@ -369,8 +466,11 @@ void EAM_FF::readParams( void ){ &isSSD, &isDipole, &isGB, - &(currentAtomType->epslon), - &(currentAtomType->sigma), + &isEAM, + &isCharge, + &eamEpslon, + &eamSigma, + &charge, &dipole, &isError ); if( isError ){ @@ -384,26 +484,57 @@ void EAM_FF::readParams( void ){ currentAtomType = currentAtomType->next; } - entry_plug->useLJ = 1; + 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", + currentAtomType->name ); + painCave.isFatal = 1; + simError(); + } + } + currentAtomType = currentAtomType->next; + } + + + #ifdef IS_MPI sprintf( checkPointMsg, - "LJ_FF atom structures successfully sent to fortran\n" ); + "EAM_FF atom structures successfully sent to fortran\n" ); 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() ); @@ -416,16 +547,16 @@ void EAM_FF::initializeAtoms( int nAtoms, Atom** the_a } the_atoms[i]->setMass( currentAtomType->mass ); - the_atoms[i]->setEpslon( currentAtomType->epslon ); - the_atoms[i]->setSigma( currentAtomType->sigma ); the_atoms[i]->setIdent( currentAtomType->ident ); - the_atoms[i]->setLJ(); + the_atoms[i]->setEAM(); + the_atoms[i]->setEamRcut( currentAtomType->eam_rcut); - if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma; + 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 ){ @@ -514,7 +645,7 @@ int EAM_NS::parseAtom( char *lineBuffer, int lineNum, -int EAM_NS::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){ +int EAM_NS::parseAtom( char *lineBuffer, int lineNum, atomStruct &info, char *eamPotFile ){ char* the_token; @@ -532,15 +663,6 @@ int EAM_NS::parseAtom( char *lineBuffer, int lineNum, 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(); - } - - info.epslon = atof( the_token ); - if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ sprintf( painCave.errMsg, "Error parseing AtomTypes: line %d\n", lineNum ); @@ -548,9 +670,400 @@ int EAM_NS::parseAtom( char *lineBuffer, int lineNum, simError(); } - info.sigma = atof( the_token ); - + strcpy( eamPotFile, the_token ); return 1; } 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; + + 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" ); + + + if( eamFile == 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, eamPotFile ); + strcpy( eamPotFile, temp ); + + eamFile = fopen( eamPotFile, "r" ); + + + + if( eamFile == NULL ){ + + sprintf( painCave.errMsg, + "Error opening the EAM force parameter file: %s\n" + "Have you tried setting the FORCE_PARAM_PATH environment " + "vairable?\n", + eamPotFile ); + 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(); + } + + + + // 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 parseing EAM ident line in %s\n", eamPotFile ); + painCave.isFatal = 1; + simError(); + } + + info.eam_ident = atoi( the_token ); + + if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ + sprintf( painCave.errMsg, + "Error parseing EAM mass in %s\n", eamPotFile ); + painCave.isFatal = 1; + simError(); + } + info.mass = atof( the_token); + + if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ + sprintf( painCave.errMsg, + "Error parseing 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 parseing EAM drho in %s\n", eamPotFile ); + painCave.isFatal = 1; + simError(); + } + info.eam_drho = atof( the_token); + + if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ + sprintf( painCave.errMsg, + "Error parseing EAM # r in %s\n", eamPotFile ); + painCave.isFatal = 1; + simError(); + } + info.eam_nr = atoi( the_token); + + if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ + sprintf( painCave.errMsg, + "Error parseing EAM dr in %s\n", eamPotFile ); + painCave.isFatal = 1; + simError(); + } + info.eam_dr = atof( the_token); + + if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ + sprintf( painCave.errMsg, + "Error parseing 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;i