--- trunk/mdtools/interface_implementation/LJ_FF.cpp 2002/12/19 21:59:51 215 +++ trunk/mdtools/interface_implementation/LJ_FF.cpp 2003/01/30 15:20:21 253 @@ -10,27 +10,69 @@ using namespace std; #include "simError.h" -#ifdef IS_MPI +// Declare the structures that will be passed by the parser and MPI -#include "mpiForceField.h" - - -// Declare the structures that will be passed by MPI - typedef struct{ char name[15]; double mass; double epslon; double sigma; + int ident; int last; // 0 -> default - // 1 -> tells nodes to stop listening + // 1 -> in MPI: tells nodes to stop listening } atomStruct; + +int parseAtomLJ( char *lineBuffer, int lineNum, atomStruct &info ); + +#ifdef IS_MPI + +#include "mpiForceField.h" + MPI_Datatype mpiAtomStructType; #endif +// declaration of functions needed to wrap the fortran module +extern "C" { + void forcefactory_( char* forceName, + int* status, + void (*wrapFunction)( void (*p1)( int* ident, + double* mass, + double* epslon, + double* sigma, + int* status ), + void (*p2)( int *nLocal, + int *identArray, + int *isError ), + void (*p3)( double* positionArray, + double* forceArray, + double* potentialEnergy, + short int* doPotentialCalc )), + int forceNameLength ); +} + + +void LJfunctionWrapper( void (*p1)( int* ident, double* mass, double* epslon, + double* sigma, int* status ), + void (*p2)( int *nLocal, int *identArray, int *isError ), + void (*p3)( double* positionArray,double* forceArray, + double* potentialEnergy, + short int* doPotentialCalc ) ); + +void (*newLJtype)( int* ident, double* mass, double* epslon, double* sigma, + int* status ); + +void (*initLJfortran) ( int *nLocal, int *identArray, int *isError ); + +LJ_FF* currentLJwrap; + + +//**************************************************************** +// begins the actual forcefield stuff. +//**************************************************************** + LJ_FF::LJ_FF(){ char fileName[200]; @@ -39,6 +81,10 @@ LJ_FF::LJ_FF(){ char temp[200]; char errMsg[1000]; + // do the funtion wrapping + currentLJwrap = this; + wrapMe(); + #ifdef IS_MPI int i; @@ -46,13 +92,13 @@ LJ_FF::LJ_FF(){ // Init the atomStruct mpi type atomStruct atomProto; // mpiPrototype - int atomBC[3] = {15,3,1}; // block counts + int atomBC[3] = {15,3,2}; // 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.last, &atomDspls[2]); + MPI_Address(&atomProto.ident, &atomDspls[2]); atomMbrTypes[0] = MPI_CHAR; atomMbrTypes[1] = MPI_DOUBLE; @@ -135,6 +181,45 @@ void LJ_FF::initializeAtoms( void ){ #endif // is_mpi } + +void LJ_FF::wrapMe( void ){ + + char* currentFF = "LJ"; + int isError = 0; + + forcefactory_( currentFF, &isError, LJfunctionWrapper, strlen(currentFF) ); + + if( isError ){ + + sprintf( painCave.errMsg, + "LJ_FF error: an error was returned from fortran when the " + "the functions were being wrapped.\n" ); + painCave.isFatal = 1; + simError(); + } + +#ifdef IS_MPI + sprintf( checkPointMsg, "LJ_FF functions succesfully wrapped." ); + MPIcheckPoint(); +#endif // is_mpi +} + + +void LJfunctionWrapper( void (*p1)( int* ident, double* mass, double* epslon, + double* sigma, int* status ), + void (*p2)( int*, int*, int* ), + void (*p3)( double* positionArray,double* forceArray, + double* potentialEnergy, + short int* doPotentialCalc ) ){ + + + newLJtype = p1; + initLJfortran = p2; + currentLJwrap->setLJfortran( p3 ); +} + + + void LJ_FF::initializeAtoms( void ){ class LinkedType { @@ -151,8 +236,20 @@ void LJ_FF::initializeAtoms( void ){ return NULL; } -#ifdef IS_MPI + void add( atomStruct &info ){ + + // check for duplicates + + if( !strcmp( info.name, name ) ){ + sprintf( painCave.errMsg, + "Duplicate LJ atom type \"%s\" found in " + "the LJ_FF param file./n", + name ); + painCave.isFatal = 1; + simError(); + } + if( next != NULL ) next->add(info); else{ next = new LinkedType(); @@ -160,14 +257,19 @@ void LJ_FF::initializeAtoms( void ){ next->mass = info.mass; next->epslon = info.epslon; next->sigma = info.sigma; + next->ident = info.ident; } } + +#ifdef IS_MPI + void duplicate( atomStruct &info ){ strcpy(info.name, name); info.mass = mass; info.epslon = epslon; info.sigma = sigma; + info.ident = ident; info.last = 0; } @@ -178,62 +280,19 @@ void LJ_FF::initializeAtoms( void ){ double mass; double epslon; double sigma; + int ident; LinkedType* next; }; LinkedType* headAtomType; LinkedType* currentAtomType; - LinkedType* tempAtomType; - -#ifdef IS_MPI atomStruct info; info.last = 1; // initialize last to have the last set. // if things go well, last will be set to 0 -#endif - - char readLine[500]; - char* the_token; - char* eof_test; - int foundAtom = 0; - int lineNum = 0; int i; - + int identNum; - ////////////////////////////////////////////////// - // a quick water fix - - double waterI[3][3]; - waterI[0][0] = 1.76958347772500; - waterI[0][1] = 0.0; - waterI[0][2] = 0.0; - - waterI[1][0] = 0.0; - waterI[1][1] = 0.614537057924513; - waterI[1][2] = 0.0; - - waterI[2][0] = 0.0; - waterI[2][1] = 0.0; - waterI[2][2] = 1.15504641980049; - - - double headI[3][3]; - headI[0][0] = 1125; - headI[0][1] = 0.0; - headI[0][2] = 0.0; - - headI[1][0] = 0.0; - headI[1][1] = 1125; - headI[1][2] = 0.0; - - headI[2][0] = 0.0; - headI[2][1] = 0.0; - headI[2][2] = 250; - - - - ////////////////////////////////////////////////// - Atom** the_atoms; int nAtoms; the_atoms = entry_plug->atoms; @@ -245,55 +304,19 @@ void LJ_FF::initializeAtoms( void ){ #endif // read in the atom types. - - rewind( frcFile ); + headAtomType = new LinkedType; - currentAtomType = headAtomType; - eof_test = fgets( readLine, sizeof(readLine), frcFile ); - lineNum++; - if( eof_test == NULL ){ - sprintf( painCave.errMsg, "Error in reading Atoms from force file.\n" ); - painCave.isFatal = 1; - simError(); - } - - - while( !foundAtom ){ - while( eof_test != NULL && readLine[0] != '#' ){ - eof_test = fgets( readLine, sizeof(readLine), frcFile ); - lineNum++; - } - if( eof_test == NULL ){ - sprintf( painCave.errMsg, - "Error in reading Atoms from force file at line %d.\n", - lineNum ); - painCave.isFatal = 1; - simError(); - } - - the_token = strtok( readLine, " ,;\t#\n" ); - foundAtom = !strcmp( "AtomTypes", the_token ); - - if( !foundAtom ){ - eof_test = fgets( readLine, sizeof(readLine), frcFile ); - lineNum++; - - if( eof_test == NULL ){ - sprintf( painCave.errMsg, - "Error in reading Atoms from force file at line %d.\n", - lineNum ); - painCave.isFatal = 1; - simError(); - } - } - } - + fastForward( "AtomTypes", "initializeAtoms" ); + // we are now at the AtomTypes section. eof_test = fgets( readLine, sizeof(readLine), frcFile ); lineNum++; + + // 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", @@ -302,48 +325,20 @@ void LJ_FF::initializeAtoms( void ){ simError(); } + 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_token = strtok( readLine, " \n\t,;" ); - if( the_token != NULL ){ - - strcpy( currentAtomType->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(); - } - - sscanf( the_token, "%lf", ¤tAtomType->mass ); - - if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ - sprintf( painCave.errMsg, - "Error parseing AtomTypes: line %d\n", lineNum ); - painCave.isFatal = 1; - simError(); - } - - sscanf( the_token, "%lf", ¤tAtomType->epslon ); - - if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){ - sprintf( painCave.errMsg, - "Error parseing AtomTypes: line %d\n", lineNum ); - painCave.isFatal = 1; - simError(); - } - - sscanf( the_token, "%lf", ¤tAtomType->sigma ); + // the parser returns 0 if the line was blank + if( parseAtomLJ( readLine, lineNum, info ) ){ + info.ident = identNum; + headAtomType->add( info );; + identNum++; } } - - tempAtomType = new LinkedType; - currentAtomType->next = tempAtomType; - currentAtomType = tempAtomType; - eof_test = fgets( readLine, sizeof(readLine), frcFile ); lineNum++; } @@ -383,9 +378,41 @@ void LJ_FF::initializeAtoms( void ){ } #endif // is_mpi + // call new A_types in fortran + int isError; + currentAtomType = headAtomType; + while( currentAtomType != NULL ){ + + if( currentAtomType->name[0] != '\0' ){ + isError = 0; + newLJtype( &(currentAtomType->ident), + &(currentAtomType->mass), + &(currentAtomType->epslon), + &(currentAtomType->sigma), + &isError ); + if( isError ){ + sprintf( painCave.errMsg, + "Error initializing the \"%s\" atom type in fortran\n", + currentAtomType->name ); + painCave.isFatal = 1; + simError(); + } + } + currentAtomType = currentAtomType->next; + } + +#ifdef IS_MPI + sprintf( checkPointMsg, + "LJ_FF atom structures successfully sent to fortran\n" ); + MPIcheckPoint(); +#endif // is_mpi + + + // initialize the atoms + double bigSigma = 0.0; Atom* thisAtom; for( i=0; isetMass( currentAtomType->mass ); the_atoms[i]->setEpslon( currentAtomType->epslon ); the_atoms[i]->setSigma( currentAtomType->sigma ); + the_atoms[i]->setIdent( currentAtomType->ident ); the_atoms[i]->setLJ(); + + if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma; } + +#ifdef IS_MPI + double tempBig = bigSigma; + MPI::COMM_WORLD::Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX ); +#endif //is_mpi + //calc rCut and rList + + entry_plug->rCut = 2.5 * bigSigma; + if(entry_plug->rCut > (entry_plug->box_x / 2.0)) entry_plug->rCut = entry_plug->box_x / 2.0; + if(entry_plug->rCut > (entry_plug->box_y / 2.0)) entry_plug->rCut = entry_plug->box_y / 2.0; + if(entry_plug->rCut > (entry_plug->box_z / 2.0)) entry_plug->rCut = entry_plug->box_z / 2.0; + + entry_plug->rList = entry_plug->rCut + 1.0; + // clean up the memory delete headAtomType; @@ -415,6 +459,9 @@ void LJ_FF::initializeAtoms( void ){ MPIcheckPoint(); #endif // is_mpi + initFortran(); + entry_plug->refreshSim(); + } void LJ_FF::initializeBonds( bond_pair* the_bonds ){ @@ -458,3 +505,162 @@ void LJ_FF::initializeTorsions( torsion_set* the_torsi #endif // is_mpi } + + +void LJ_FF::fastForward( char* stopText, char* searchOwner ){ + + int foundText = 0; + char* the_token; + + rewind( frcFile ); + lineNum = 0; + + 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 ); + painCave.isFatal = 1; + simError(); + } + + + 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(); + } + } + } +} + + + +int parseAtomLJ( char *lineBuffer, int lineNum, atomStruct &info ){ + + char* the_token; + + 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(); + } + + info.epslon = 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.sigma = atof( the_token ); + + return 1; + } + else return 0; +} + + +void LJ_FF::doForces( int calcPot ){ + + int i; + double* frc; + double* pos; + short int passedCalcPot = (short int)calcPot; + + // forces are zeroed here, before any are acumulated. + // NOTE: do not rezero the forces in Fortran. + + for(i=0; in_atoms; i++){ + entry_plug->atoms[i]->zeroForces(); + } + + frc = Atom::getFrcArray(); + pos = Atom::getPosArray(); + +// entry_plug->lrPot = -1; + doLJfortran( pos, frc, &(entry_plug->lrPot), &passedCalcPot ); + + + // fprintf( stderr, +// "lrPot = %lf\n", entry_plug->lrPot ); + +} + +void LJ_FF::initFortran( void ){ + + int nLocal = entry_plug->n_atoms; + int *ident; + int isError; + int i; + + ident = new int[nLocal]; + + for(i=0; iatoms[i]->getIdent(); + } + + isError = 0; + initLJfortran( &nLocal, ident, &isError ); + + if(isError){ + sprintf( painCave.errMsg, + "LJ_FF error: There was an error initializing the component list in fortran.\n" ); + painCave.isFatal = 1; + simError(); + } + + +#ifdef IS_MPI + sprintf( checkPointMsg, "LJ_FF successfully initialized the fortran component list.\n" ); + MPIcheckPoint(); +#endif // is_mpi + + delete[] ident; + +} +