--- trunk/OOPSE/libmdtools/SimInfo.cpp 2003/04/04 19:16:11 457 +++ trunk/OOPSE/libmdtools/SimInfo.cpp 2003/04/08 22:38:43 482 @@ -16,6 +16,8 @@ SimInfo::SimInfo(){ n_constraints = 0; n_oriented = 0; n_dipoles = 0; + ndf = 0; + ndfRaw = 0; the_integrator = NULL; setTemp = 0; thermalTime = 0.0; @@ -33,10 +35,56 @@ void SimInfo::setBox(double newBox[3]) { } void SimInfo::setBox(double newBox[3]) { + double smallestBox, maxCutoff; + int status; box_x = newBox[0]; box_y = newBox[1]; box_z = newBox[2]; setFortranBoxSize(newBox); + + smallestBox = box_x; + if (box_y < smallestBox) smallestBox = box_y; + if (box_z < smallestBox) smallestBox = box_z; + + maxCutoff = smallestBox / 2.0; + + if (rList > maxCutoff) { + sprintf( painCave.errMsg, + "New Box size is forcing neighborlist radius down to %lf\n", + maxCutoff ); + painCave.isFatal = 0; + simError(); + + rList = maxCutoff; + + sprintf( painCave.errMsg, + "New Box size is forcing cutoff radius down to %lf\n", + maxCutoff - 1.0 ); + painCave.isFatal = 0; + simError(); + + rCut = rList - 1.0; + + // list radius changed so we have to refresh the simulation structure. + refreshSim(); + } + + if (rCut > maxCutoff) { + sprintf( painCave.errMsg, + "New Box size is forcing cutoff radius down to %lf\n", + maxCutoff ); + painCave.isFatal = 0; + simError(); + + status = 0; + LJ_new_rcut(&rCut, &status); + if (status != 0) { + sprintf( painCave.errMsg, + "Error in recomputing LJ shifts based on new rcut\n"); + painCave.isFatal = 1; + simError(); + } + } } void SimInfo::getBox(double theBox[3]) { @@ -44,12 +92,47 @@ void SimInfo::getBox(double theBox[3]) { theBox[1] = box_y; theBox[2] = box_z; } + +int SimInfo::getNDF(){ + int ndf_local, ndf; + ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints; + +#ifdef IS_MPI + MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); +#else + ndf = ndf_local; +#endif + + ndf = ndf - 3; + + return ndf; +} + +int SimInfo::getNDFraw() { + int ndfRaw_local, ndfRaw; + + // Raw degrees of freedom that we have to set + ndfRaw_local = 3 * n_atoms + 3 * n_oriented; + +#ifdef IS_MPI + MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); +#else + ndfRaw = ndfRaw_local; +#endif + + return ndfRaw; +} + void SimInfo::refreshSim(){ simtype fInfo; int isError; int* excl; + + fInfo.rrf = 0.0; + fInfo.rt = 0.0; + fInfo.dielect = 0.0; fInfo.box[0] = box_x; fInfo.box[1] = box_y; @@ -57,10 +140,13 @@ void SimInfo::refreshSim(){ fInfo.rlist = rList; fInfo.rcut = rCut; - fInfo.rrf = ecr; - fInfo.rt = ecr - est; - fInfo.dielect = dielectric; + if( useDipole ){ + fInfo.rrf = ecr; + fInfo.rt = ecr - est; + if( useReactionField )fInfo.dielect = dielectric; + } + fInfo.SIM_uses_PBC = usePBC; //fInfo.SIM_uses_LJ = 0; fInfo.SIM_uses_LJ = useLJ; @@ -102,5 +188,9 @@ void SimInfo::refreshSim(){ "succesfully sent the simulation information to fortran.\n"); MPIcheckPoint(); #endif // is_mpi + + this->ndf = this->getNDF(); + this->ndfRaw = this->getNDFraw(); + }