--- trunk/OOPSE/libmdtools/SimInfo.cpp 2003/04/07 20:06:31 469 +++ 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(){ @@ -35,19 +47,28 @@ void SimInfo::setBox(double newBox[3]) { } void SimInfo::setBox(double newBox[3]) { - double smallestBox, maxCutoff; + + double smallestBoxL, maxCutoff; int status; - box_x = newBox[0]; - box_y = newBox[1]; - box_z = newBox[2]; - setFortranBoxSize(newBox); + int i; - smallestBox = box_x; - if (box_y < smallestBox) smallestBox = box_y; - if (box_z < smallestBox) smallestBox = box_z; + for(i=0; i<9; i++) Hmat[i] = 0.0;; - maxCutoff = smallestBox / 2.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", @@ -87,12 +108,223 @@ 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::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; @@ -128,16 +360,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; @@ -152,8 +381,8 @@ void SimInfo::refreshSim(){ fInfo.SIM_uses_LJ = useLJ; fInfo.SIM_uses_sticky = useSticky; //fInfo.SIM_uses_sticky = 0; - //fInfo.SIM_uses_dipoles = useDipole; - fInfo.SIM_uses_dipoles = 0; + fInfo.SIM_uses_dipoles = useDipole; + //fInfo.SIM_uses_dipoles = 0; //fInfo.SIM_uses_RF = useReactionField; fInfo.SIM_uses_RF = 0; fInfo.SIM_uses_GB = useGB; @@ -161,20 +390,18 @@ void SimInfo::refreshSim(){ 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, @@ -189,8 +416,8 @@ void SimInfo::refreshSim(){ MPIcheckPoint(); #endif // is_mpi - ndf = this->getNDF(); - ndfRaw = this->getNDFraw(); + this->ndf = this->getNDF(); + this->ndfRaw = this->getNDFraw(); }