--- trunk/mdtools/interface_implementation/LJ_FF.cpp 2003/01/09 19:40:38 230 +++ trunk/mdtools/interface_implementation/LJ_FF.cpp 2003/01/28 22:16:55 252 @@ -33,6 +33,46 @@ 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 +81,10 @@ LJ_FF::LJ_FF(){ char temp[200]; char errMsg[1000]; + // do the funtion wrapping + currentLJwrap = this; + wrapMe(); + #ifdef IS_MPI int i; @@ -137,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 { @@ -159,7 +242,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 +307,7 @@ void LJ_FF::initializeAtoms( void ){ headAtomType = new LinkedType; - fastFoward( "AtomTypes", "initializeAtoms" ); + fastForward( "AtomTypes", "initializeAtoms" ); // we are now at the AtomTypes section. @@ -296,10 +379,40 @@ 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; 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 +459,9 @@ void LJ_FF::initializeAtoms( void ){ MPIcheckPoint(); #endif // is_mpi + initFortran(); + entry_plug->refreshSim(); + } void LJ_FF::initializeBonds( bond_pair* the_bonds ){ @@ -471,3 +603,58 @@ int parseAtomLJ( char *lineBuffer, int lineNum, atomS } else return 0; } + + +void LJ_FF::doForces( void ){ + + int i; + double* frc; + double* pos; + short int calcPot = 1; + + // 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(); + + doLJfortran( pos, frc, &(entry_plug->lrPot), &calcPot ); +} + +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; + +} +