--- trunk/OOPSE/libmdtools/SimInfo.cpp 2003/03/24 21:55:34 394 +++ trunk/OOPSE/libmdtools/SimInfo.cpp 2003/06/30 22:04:01 568 @@ -1,5 +1,6 @@ #include #include +#include #include "SimInfo.hpp" @@ -9,6 +10,10 @@ SimInfo* currentInfo; #include "fortranWrappers.hpp" +#ifdef IS_MPI +#include "mpiSimulation.hpp" +#endif + SimInfo* currentInfo; SimInfo::SimInfo(){ @@ -16,9 +21,12 @@ SimInfo::SimInfo(){ n_constraints = 0; n_oriented = 0; n_dipoles = 0; + ndf = 0; + ndfRaw = 0; the_integrator = NULL; setTemp = 0; thermalTime = 0.0; + rCut = 0.0; usePBC = 0; useLJ = 0; @@ -28,15 +36,266 @@ SimInfo::SimInfo(){ useGB = 0; useEAM = 0; + wrapMeSimInfo( this ); +} +void SimInfo::setBox(double newBox[3]) { - wrapMeSimInfo( this ); + double smallestBoxL, maxCutoff; + int status; + int i; + + for(i=0; i<9; i++) Hmat[i] = 0.0;; + + Hmat[0] = newBox[0]; + Hmat[4] = newBox[1]; + Hmat[8] = newBox[2]; + + calcHmatI(); + calcBoxL(); + + setFortranBoxSize(Hmat); + + smallestBoxL = boxLx; + if (boxLy < smallestBoxL) smallestBoxL = boxLy; + if (boxLz < smallestBoxL) smallestBoxL = boxLz; + + maxCutoff = smallestBoxL / 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::setBoxM( double theBox[9] ){ + + int i, status; + double smallestBoxL, maxCutoff; + + for(i=0; i<9; i++) Hmat[i] = theBox[i]; + calcHmatI(); + calcBoxL(); + + setFortranBoxSize(Hmat); + + smallestBoxL = boxLx; + if (boxLy < smallestBoxL) smallestBoxL = boxLy; + if (boxLz < smallestBoxL) smallestBoxL = boxLz; + + maxCutoff = smallestBoxL / 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[9]) { + + int i; + for(i=0; i<9; i++) theBox[i] = Hmat[i]; +} + + +void SimInfo::calcHmatI( void ) { + + double C[3][3]; + double detHmat; + int i, j, k; + + // calculate the adjunct of Hmat; + + C[0][0] = ( Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]); + C[1][0] = -( Hmat[1]*Hmat[8]) + (Hmat[7]*Hmat[2]); + C[2][0] = ( Hmat[1]*Hmat[5]) - (Hmat[4]*Hmat[2]); + + C[0][1] = -( Hmat[3]*Hmat[8]) + (Hmat[6]*Hmat[5]); + C[1][1] = ( Hmat[0]*Hmat[8]) - (Hmat[6]*Hmat[2]); + C[2][1] = -( Hmat[0]*Hmat[5]) + (Hmat[3]*Hmat[2]); + + C[0][2] = ( Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]); + C[1][2] = -( Hmat[0]*Hmat[7]) + (Hmat[6]*Hmat[1]); + C[2][2] = ( Hmat[0]*Hmat[4]) - (Hmat[3]*Hmat[1]); + + // calcutlate the determinant of Hmat + + detHmat = 0.0; + for(i=0; i<3; i++) detHmat += Hmat[i] * C[i][0]; + + + // H^-1 = C^T / det(H) + + i=0; + for(j=0; j<3; j++){ + for(k=0; k<3; k++){ + + HmatI[i] = C[j][k] / detHmat; + i++; + } + } +} + +void SimInfo::calcBoxL( void ){ + + double dx, dy, dz, dsq; + int i; + + // boxVol = h1 (dot) h2 (cross) h3 + + boxVol = Hmat[0] * ( (Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]) ) + + Hmat[1] * ( (Hmat[5]*Hmat[6]) - (Hmat[8]*Hmat[3]) ) + + Hmat[2] * ( (Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]) ); + + + // boxLx + + dx = Hmat[0]; dy = Hmat[1]; dz = Hmat[2]; + dsq = dx*dx + dy*dy + dz*dz; + boxLx = sqrt( dsq ); + + // boxLy + + dx = Hmat[3]; dy = Hmat[4]; dz = Hmat[5]; + dsq = dx*dx + dy*dy + dz*dz; + boxLy = sqrt( dsq ); + + // boxLz + + dx = Hmat[6]; dy = Hmat[7]; dz = Hmat[8]; + dsq = dx*dx + dy*dy + dz*dz; + boxLz = sqrt( dsq ); + +} + + +void SimInfo::wrapVector( double thePos[3] ){ + + int i, j, k; + double scaled[3]; + + // calc the scaled coordinates. + + for(i=0; i<3; i++) + scaled[i] = thePos[0]*Hmat[i] + thePos[1]*Hat[i+3] + thePos[3]*Hmat[i+6]; + + // wrap the scaled coordinates + + for(i=0; i<3; i++) + scaled[i] -= (copysign(1,scaled[i]) * (int)(fabs(scaled[i]) + 0.5)); + + +} + + +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 n_global; + int* excl; + + fInfo.rrf = 0.0; + fInfo.rt = 0.0; + fInfo.dielect = 0.0; fInfo.box[0] = box_x; fInfo.box[1] = box_y; @@ -44,33 +303,39 @@ 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; - //fInfo.SIM_uses_sticky = useSticky; - fInfo.SIM_uses_sticky = 0; + fInfo.SIM_uses_sticky = useSticky; + //fInfo.SIM_uses_sticky = 0; fInfo.SIM_uses_dipoles = useDipole; - fInfo.SIM_uses_RF = useReactionField; + //fInfo.SIM_uses_dipoles = 0; + //fInfo.SIM_uses_RF = useReactionField; + fInfo.SIM_uses_RF = 0; fInfo.SIM_uses_GB = useGB; fInfo.SIM_uses_EAM = useEAM; + excl = Exclude::getArray(); +#ifdef IS_MPI + n_global = mpiSim->getTotAtoms(); +#else + n_global = n_atoms; +#endif + isError = 0; - fInfo; - n_atoms; - identArray; - n_exclude; - excludes; - nGlobalExcludes; - globalExcludes; - isError; + setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl, + &nGlobalExcludes, globalExcludes, molMembershipArray, + &isError ); - setFsimulation( &fInfo, &n_atoms, identArray, &n_exclude, excludes, &nGlobalExcludes, globalExcludes, &isError ); - if( isError ){ sprintf( painCave.errMsg, @@ -84,5 +349,9 @@ void SimInfo::refreshSim(){ "succesfully sent the simulation information to fortran.\n"); MPIcheckPoint(); #endif // is_mpi + + this->ndf = this->getNDF(); + this->ndfRaw = this->getNDFraw(); + }