--- trunk/OOPSE/libmdtools/SimInfo.cpp 2003/03/27 20:36:16 424 +++ trunk/OOPSE/libmdtools/SimInfo.cpp 2003/07/02 21:26:55 572 @@ -1,6 +1,9 @@ #include #include +#include +#include +using namespace std; #include "SimInfo.hpp" #define __C @@ -9,6 +12,15 @@ #include "fortranWrappers.hpp" +#ifdef IS_MPI +#include "mpiSimulation.hpp" +#endif + +inline double roundMe( double x ){ + return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 ); +} + + SimInfo* currentInfo; SimInfo::SimInfo(){ @@ -16,6 +28,8 @@ SimInfo::SimInfo(){ n_constraints = 0; n_oriented = 0; n_dipoles = 0; + ndf = 0; + ndfRaw = 0; the_integrator = NULL; setTemp = 0; thermalTime = 0.0; @@ -29,53 +43,365 @@ 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, HmatI, &orthoRhombic); + + 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, HmatI, &orthoRhombic); + + 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::getBoxM (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; + double smallDiag; + double tol; + double sanity[3][3]; + + // 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++; + } + } + + // sanity check + + for(i=0; i<3; i++){ + for(j=0; j<3; j++){ + + sanity[i][j] = 0.0; + for(k=0; k<3; k++){ + sanity[i][j] += Hmat[3*k+i] * HmatI[3*j+k]; + } + } + } + + cerr << "sanity => \n" + << sanity[0][0] << "\t" << sanity[0][1] << "\t" << sanity [0][2] << "\n" + << sanity[1][0] << "\t" << sanity[1][1] << "\t" << sanity [1][2] << "\n" + << sanity[2][0] << "\t" << sanity[2][1] << "\t" << sanity [2][2] + << "\n"; + + + // check to see if Hmat is orthorhombic + + smallDiag = Hmat[0]; + if(smallDiag > Hmat[4]) smallDiag = Hmat[4]; + if(smallDiag > Hmat[8]) smallDiag = Hmat[8]; + tol = smallDiag * 1E-6; + + orthoRhombic = 1; + for(i=0; (i<9) && orthoRhombic; i++){ + + if( (i%4) ){ // ignore the diagonals (0, 4, and 8) + orthoRhombic = (Hmat[i] <= tol); + } + } + } +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]; + + if( !orthoRhombic ){ + // calc the scaled coordinates. + + for(i=0; i<3; i++) + scaled[i] = + thePos[0]*HmatI[i] + thePos[1]*HmatI[i+3] + thePos[3]*HmatI[i+6]; + + // wrap the scaled coordinates + + for(i=0; i<3; i++) + scaled[i] -= roundMe(scaled[i]); + + // calc the wrapped real coordinates from the wrapped scaled coordinates + + for(i=0; i<3; i++) + thePos[i] = + scaled[0]*Hmat[i] + scaled[1]*Hmat[i+3] + scaled[2]*Hmat[i+6]; + } + else{ + // calc the scaled coordinates. + + for(i=0; i<3; i++) + scaled[i] = thePos[i]*HmatI[i*4]; + + // wrap the scaled coordinates + + for(i=0; i<3; i++) + scaled[i] -= roundMe(scaled[i]); + + // calc the wrapped real coordinates from the wrapped scaled coordinates + + for(i=0; i<3; i++) + thePos[i] = scaled[i]*Hmat[i*4]; + } + + +} + + +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; - fInfo.box[2] = box_z; - 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_dipoles = 0; - fInfo.SIM_uses_RF = useReactionField; + //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, excl, - &nGlobalExcludes, globalExcludes, &isError ); - if( isError ){ sprintf( painCave.errMsg, @@ -89,5 +415,9 @@ void SimInfo::refreshSim(){ "succesfully sent the simulation information to fortran.\n"); MPIcheckPoint(); #endif // is_mpi + + this->ndf = this->getNDF(); + this->ndfRaw = this->getNDFraw(); + }