--- trunk/mdtools/interface_implementation/LJ_FF.cpp 2003/01/09 19:40:38 230 +++ trunk/mdtools/interface_implementation/LJ_FF.cpp 2003/02/04 20:15:48 263 @@ -9,7 +9,6 @@ using namespace std; #include "SRI.hpp" #include "simError.h" - // Declare the structures that will be passed by the parser and MPI typedef struct{ @@ -25,7 +24,6 @@ int parseAtomLJ( char *lineBuffer, int lineNum, atomSt int parseAtomLJ( char *lineBuffer, int lineNum, atomStruct &info ); #ifdef IS_MPI - #include "mpiForceField.h" MPI_Datatype mpiAtomStructType; @@ -33,6 +31,47 @@ LJ_FF::LJ_FF(){ #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]; @@ -41,6 +80,10 @@ LJ_FF::LJ_FF(){ char temp[200]; char errMsg[1000]; + // do the funtion wrapping + currentLJwrap = this; + wrapMe(); + #ifdef IS_MPI int i; @@ -133,10 +176,49 @@ LJ_FF::~LJ_FF(){ fclose( frcFile ); #ifdef IS_MPI + } +#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 { @@ -159,7 +241,7 @@ void LJ_FF::initializeAtoms( void ){ // check for duplicates if( !strcmp( info.name, name ) ){ - sprintf( simError.painCave, + sprintf( painCave.errMsg, "Duplicate LJ atom type \"%s\" found in " "the LJ_FF param file./n", name ); @@ -224,7 +306,7 @@ void LJ_FF::initializeAtoms( void ){ headAtomType = new LinkedType; - fastFoward( "AtomTypes", "initializeAtoms" ); + fastForward( "AtomTypes", "initializeAtoms" ); // we are now at the AtomTypes section. @@ -268,10 +350,19 @@ void LJ_FF::initializeAtoms( void ){ "LJ_FF atom structures read successfully." ); MPIcheckPoint(); - currentAtomType = headAtomType; + currentAtomType = headAtomType->next; //skip the first element who is a place holder. while( currentAtomType != NULL ){ currentAtomType->duplicate( info ); + + + sendFrcStruct( &info, mpiAtomStructType ); + + sprintf( checkPointMsg, + "successfully sent lJ force type: \"%s\"\n", + info.name ); + MPIcheckPoint(); + currentAtomType = currentAtomType->next; } info.last = 1; @@ -287,19 +378,55 @@ void LJ_FF::initializeAtoms( void ){ headAtomType = new LinkedType; recieveFrcStruct( &info, mpiAtomStructType ); + while( !info.last ){ + + headAtomType->add( info ); + + MPIcheckPoint(); + recieveFrcStruct( &info, mpiAtomStructType ); } } #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; isetSigma( 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; @@ -330,6 +473,9 @@ void LJ_FF::initializeAtoms( void ){ MPIcheckPoint(); #endif // is_mpi + initFortran(); + entry_plug->refreshSim(); + } void LJ_FF::initializeBonds( bond_pair* the_bonds ){ @@ -471,3 +617,64 @@ int parseAtomLJ( char *lineBuffer, int lineNum, atomS } 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; + +} +