--- trunk/OOPSE/libmdtools/ForceFields.cpp 2003/04/02 22:19:03 443 +++ trunk/OOPSE/libmdtools/ForceFields.cpp 2003/10/29 17:55:28 841 @@ -1,8 +1,12 @@ -#include +#include +using namespace std; + + +#include + #ifdef IS_MPI #include -#include #endif // is_mpi @@ -16,21 +20,14 @@ void ForceFields::calcRcut( void ){ #ifdef IS_MPI double tempBig = bigSigma; - MPI::COMM_WORLD.Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX ); + MPI_Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX, + MPI_COMM_WORLD); #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; - + entry_plug->setDefaultRcut( 2.5 * bigSigma ); + } void ForceFields::doForces( int calcPot, int calcStress ){ @@ -39,68 +36,68 @@ 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(); } 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; - + 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; + } - + std::cerr << "before\n" + << " force[0] = " << frc[0] << "\n" + << " pos[0] = " << pos[0] << "\n" + << " trq[0] = " << trq[0] << "\n" + << " A[0] = " << A[0] << "\n" + << " tau = " << entry_plug->tau[0] << ", " << entry_plug->tau[1] << ", " << entry_plug->tau[2] << "\n" + << " " << entry_plug->tau[3] << ", " << entry_plug->tau[4] << ", " << entry_plug->tau[5] << "\n" + << " " << entry_plug->tau[6] << ", " << entry_plug->tau[7] << ", " << entry_plug->tau[8] << "\n\n"; + fortranForceLoop( pos, A, u_l, frc, trq, - tau, + entry_plug->tau, &(entry_plug->lrPot), &passedCalcPot, &passedCalcStress, &isError ); - // delete[] u_l; + std::cerr << "after\n" + << " force[0] = " << frc[0] << "\n" + << " pos[0] = " << pos[0] << "\n" + << " trq[0] = " << trq[0] << "\n" + << " A[0] = " << A[0] << "\n" + << " tau = " << entry_plug->tau[0] << ", " << entry_plug->tau[1] << ", " << entry_plug->tau[2] << "\n" + << " " << entry_plug->tau[3] << ", " << entry_plug->tau[4] << ", " << entry_plug->tau[5] << "\n" + << " " << entry_plug->tau[6] << ", " << entry_plug->tau[7] << ", " << entry_plug->tau[8] << "\n\n"; if( isError ){ sprintf( painCave.errMsg, @@ -114,6 +111,7 @@ void ForceFields::doForces( int calcPot, int calcStres "returned from the force calculation.\n" ); MPIcheckPoint(); #endif // is_mpi + }