--- trunk/OOPSE/libmdtools/ForceFields.cpp 2003/03/24 15:26:05 389 +++ trunk/OOPSE/libmdtools/ForceFields.cpp 2004/05/20 20:24:07 1180 @@ -1,62 +1,191 @@ -#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" #include "fortranWrappers.hpp" +void ForceFields::calcRcut( void ){ + +#ifdef IS_MPI + double tempBig = bigSigma; + MPI_Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX, + MPI_COMM_WORLD); +#endif //is_mpi + + //calc rCut and rList + + entry_plug->setDefaultRcut( 2.5 * bigSigma ); + +} + +void ForceFields::setRcut( double LJrcut ) { + +#ifdef IS_MPI + double tempBig = bigSigma; + MPI_Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX, + MPI_COMM_WORLD); +#endif //is_mpi + + if (LJrcut < 2.5 * bigSigma) { + sprintf( painCave.errMsg, + "Setting Lennard-Jones cutoff radius to %lf.\n" + "\tThis value is smaller than %lf, which is\n" + "\t2.5 * bigSigma, where bigSigma is the largest\n" + "\tvalue of sigma present in the simulation.\n" + "\tThis is potentially a problem since the LJ potential may\n" + "\tbe appreciable at this distance. If you don't want the\n" + "\tsmaller cutoff, change the LJrcut variable.\n", + LJrcut, 2.5*bigSigma); + painCave.isFatal = 0; + simError(); + } else { + sprintf( painCave.errMsg, + "Setting Lennard-Jones cutoff radius to %lf.\n" + "\tThis value is larger than %lf, which is\n" + "\t2.5 * bigSigma, where bigSigma is the largest\n" + "\tvalue of sigma present in the simulation. This should\n" + "\tnot be a problem, but could adversely effect performance.\n", + LJrcut, 2.5*bigSigma); + painCave.isFatal = 0; + simError(); + } + + //calc rCut and rList + + entry_plug->setDefaultRcut( LJrcut ); +} + void ForceFields::doForces( int calcPot, int calcStress ){ - int i, isError; + int i, j, isError; double* frc; double* pos; double* trq; - double* tau; double* A; double* u_l; + double* rc; + double* massRatio; + double factor; + SimState* config; + Molecule* myMols; + Atom** myAtoms; + int numAtom; + int curIndex; + double mtot; + int numMol; + int numCutoffGroups; + CutoffGroup* myCutoffGroup; + vector::iterator iterCutoff; + double com[3]; + vector rcGroup; + 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(); + entry_plug->atoms[i]->zeroForces(); } - for(i=0; in_SRI; i++ ){ - entry_plug->sr_interactions[i]->calc_forces(); +#ifdef PROFILE + startProfile(pro7); +#endif + + for(i=0; in_mol; i++ ){ + // CalcForces in molecules takes care of mapping rigid body coordinates + // into atomic coordinates + entry_plug->molecules[i].calcForces(); } - frc = Atom::getFrcArray(); - pos = Atom::getPosArray(); - trq = Atom::getTrqArray(); - A = Atom::getAmatArray(); - u_l = Atom::getUlArray(); +#ifdef PROFILE + endProfile( pro7 ); +#endif - tau = entry_plug->tau; - + config = entry_plug->getConfiguration(); + + frc = config->getFrcArray(); + pos = config->getPosArray(); + trq = config->getTrqArray(); + A = config->getAmatArray(); + u_l = config->getUlArray(); + + if(entry_plug->haveCutoffGroups){ + myMols = entry_plug->molecules; + numMol = entry_plug->n_mol; + for(int i = 0; i < numMol; i++){ + + numCutoffGroups = myMols[i].getNCutoffGroups(); + for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff); myCutoffGroup != NULL; + myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){ + //get center of mass of the cutoff group + myCutoffGroup->getCOM(com); + + rcGroup.push_back(com[0]); + rcGroup.push_back(com[1]); + rcGroup.push_back(com[2]); + + }// end for(myCutoffGroup) + + }//end for(int i = 0) + + rc = &rcGroup[0]; + } + else{ + // center of mass of the group is the same as position of the atom if cutoff group does not exist + rc = pos; + } + + + 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, + rc, A, u_l, frc, trq, - tau, + entry_plug->tau, &(entry_plug->lrPot), &passedCalcPot, &passedCalcStress, &isError ); +#ifdef PROFILE + endProfile(pro8); +#endif + + if( isError ){ sprintf( painCave.errMsg, "Error returned from the fortran force calculation.\n" ); @@ -64,11 +193,35 @@ void ForceFields::doForces( int calcPot, int calcStres simError(); } + for(i=0; in_mol; i++ ){ + entry_plug->molecules[i].atoms2rigidBodies(); + } + + + if (entry_plug->useThermInt) { + myStunts = entry_plug->integrableObjects; + + factor = pow(entry_plug->thermIntLambda, entry_plug->thermIntK); + for (i=0; i < myStunts.size(); i++) { + for (j=0; j< 3; j++) + frc[3*i + j] *= factor; + if (myStunts[i]->isDirectional()) { + for (j=0; j< 3; j++) + trq[3*i + j] *= factor; + } + } + entry_plug->vRaw = entry_plug->lrPot; + entry_plug->lrPot *= factor; + entry_plug->lrPot += entry_plug->restraint->Calc_Restraint_Forces(myStunts); + entry_plug->vHarm = entry_plug->restraint->getVharm(); + } + #ifdef IS_MPI sprintf( checkPointMsg, "returned from the force calculation.\n" ); MPIcheckPoint(); #endif // is_mpi + } @@ -94,3 +247,18 @@ void ForceFields::initFortran(int ljMixPolicy, int use #endif // is_mpi } + + +void ForceFields::initRestraints(){ + + // store the initial info. + entry_plug->restraint->Store_Init_Info(); + +} + +void ForceFields::dumpzAngle(){ + + // store the initial info. + entry_plug->restraint->Write_zAngle_File(); + +}