--- trunk/OOPSE/libmdtools/SimInfo.cpp 2003/04/09 04:06:43 483 +++ trunk/OOPSE/libmdtools/SimInfo.cpp 2003/07/09 22:14:06 586 @@ -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(){ @@ -35,19 +47,46 @@ 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); + + int i; + double tempMat[9]; - smallestBox = box_x; - if (box_y < smallestBox) smallestBox = box_y; - if (box_z < smallestBox) smallestBox = box_z; + for(i=0; i<9; i++) tempMat[i] = 0.0;; - maxCutoff = smallestBox / 2.0; + tempMat[0] = newBox[0]; + tempMat[4] = newBox[1]; + tempMat[8] = newBox[2]; + setBoxM( tempMat ); + +} + +void SimInfo::setBoxM( double theBox[9] ){ + + int i, status; + double smallestBoxL, maxCutoff; + + for(i=0; i<9; i++) Hmat[i] = theBox[i]; + + cerr + << "setting Hmat ->\n" + << "[ " << Hmat[0] << ", " << Hmat[3] << ", " << Hmat[6] << " ]\n" + << "[ " << Hmat[1] << ", " << Hmat[4] << ", " << Hmat[7] << " ]\n" + << "[ " << Hmat[2] << ", " << Hmat[5] << ", " << Hmat[8] << " ]\n"; + + 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", @@ -86,13 +125,180 @@ void SimInfo::setBox(double newBox[3]) { } } } + -void SimInfo::getBox(double theBox[3]) { - theBox[0] = box_x; - theBox[1] = box_y; - theBox[2] = box_z; +void SimInfo::getBoxM (double theBox[9]) { + + int i; + for(i=0; i<9; i++) theBox[i] = Hmat[i]; } - + + +void SimInfo::scaleBox(double scale) { + double theBox[9]; + int i; + + cerr << "Scaling box by " << scale << "\n"; + + for(i=0; i<9; i++) theBox[i] = Hmat[i]*scale; + + setBoxM(theBox); + +} + +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; @@ -128,16 +334,13 @@ 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; @@ -161,9 +364,15 @@ void SimInfo::refreshSim(){ excl = Exclude::getArray(); +#ifdef IS_MPI + n_global = mpiSim->getTotAtoms(); +#else + n_global = n_atoms; +#endif + isError = 0; - setFsimulation( &fInfo, &n_atoms, identArray, &n_exclude, excl, + setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl, &nGlobalExcludes, globalExcludes, molMembershipArray, &isError );