--- trunk/OOPSE/libmdtools/ForceFields.cpp 2003/04/03 20:21:54 447 +++ trunk/OOPSE/libmdtools/ForceFields.cpp 2003/12/22 21:27:04 892 @@ -1,10 +1,19 @@ -#include +#include +using namespace std; + + +#include + #ifdef IS_MPI #include #endif // is_mpi +#ifdef PROFILE +#include "mdProfile.hpp" +#endif + #include "simError.h" #include "ForceFields.hpp" #include "Atom.hpp" @@ -21,16 +30,8 @@ void ForceFields::calcRcut( void ){ //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; - + entry_plug->setDefaultRcut( 2.5 * bigSigma ); + } void ForceFields::doForces( int calcPot, int calcStress ){ @@ -39,68 +40,67 @@ void ForceFields::doForces( int calcPot, int calcStres double* frc; double* pos; double* trq; - double* tau; double* A; double* u_l;; - DirectionalAtom* dAtom; + SimState* config; - double ut[3]; - - //u_l = new double[entry_plug->n_atoms]; - short int passedCalcPot = (short int)calcPot; short int passedCalcStress = (short int)calcStress; - // forces are zeroed here, before any are acumulated. + // forces are zeroed here, before any are accumulated. // NOTE: do not rezero the forces in Fortran. - for(i=0; in_atoms; i++){ - entry_plug->atoms[i]->zeroForces(); - -// if( entry_plug->atoms[i]->isDirectional() ){ -// dAtom = (DirectionalAtom *)entry_plug->atoms[i]; -// dAtom->getU(ut); - - -// if(dAtom->getIndex()== 1){ -// std::cerr << "atom 2's u_l = " << ut[0] << ", " << ut[1] -// << ", " << ut[2] << "\n"; -// } -// } - + entry_plug->atoms[i]->zeroForces(); } +#ifdef PROFILE + startProfile(pro7); +#endif + for(i=0; in_mol; i++ ){ entry_plug->molecules[i].calcForces(); } - frc = Atom::getFrcArray(); - pos = Atom::getPosArray(); - trq = Atom::getTrqArray(); - A = Atom::getAmatArray(); - u_l = Atom::getUlArray(); - tau = entry_plug->tau; +#ifdef PROFILE + endProfile( pro7 ); +#endif + config = entry_plug->getConfiguration(); + frc = config->getFrcArray(); + pos = config->getPosArray(); + trq = config->getTrqArray(); + A = config->getAmatArray(); + u_l = config->getUlArray(); + isError = 0; entry_plug->lrPot = 0.0; + for (i=0; i<9; i++) { + entry_plug->tau[i] = 0.0; + } - + +#ifdef PROFILE + startProfile(pro8); +#endif + fortranForceLoop( pos, A, u_l, frc, trq, - tau, + entry_plug->tau, &(entry_plug->lrPot), &passedCalcPot, &passedCalcStress, &isError ); +#ifdef PROFILE + endProfile(pro8); +#endif - // delete[] u_l; if( isError ){ sprintf( painCave.errMsg, @@ -114,6 +114,7 @@ void ForceFields::doForces( int calcPot, int calcStres "returned from the force calculation.\n" ); MPIcheckPoint(); #endif // is_mpi + }