--- trunk/OOPSE/libmdtools/SimInfo.cpp	2003/07/10 17:10:56	588
+++ trunk/OOPSE/libmdtools/SimInfo.cpp	2003/07/16 02:11:02	621
@@ -34,6 +34,8 @@ SimInfo::SimInfo(){
   setTemp = 0;
   thermalTime = 0.0;
   rCut = 0.0;
+  ecr = 0.0;
+  est = 0.0;
 
   usePBC = 0;
   useLJ = 0; 
@@ -77,11 +79,11 @@ void SimInfo::setBoxM( double theBox[3][3] ){
   for(i=0; i < 3; i++) 
     for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
   
-  cerr 
-    << "setting Hmat ->\n"
-    << "[ " << Hmat[0][0] << ", " << Hmat[0][1] << ", " << Hmat[0][2] << " ]\n"
-    << "[ " << Hmat[1][0] << ", " << Hmat[1][1] << ", " << Hmat[1][2] << " ]\n"
-    << "[ " << Hmat[2][0] << ", " << Hmat[2][1] << ", " << Hmat[2][2] << " ]\n";
+  //  cerr 
+  // << "setting Hmat ->\n"
+  // << "[ " << Hmat[0][0] << ", " << Hmat[0][1] << ", " << Hmat[0][2] << " ]\n"
+  // << "[ " << Hmat[1][0] << ", " << Hmat[1][1] << ", " << Hmat[1][2] << " ]\n"
+  // << "[ " << Hmat[2][0] << ", " << Hmat[2][1] << ", " << Hmat[2][2] << " ]\n";
 
   calcBoxL();
   calcHmatInv();
@@ -93,11 +95,11 @@ void SimInfo::setBoxM( double theBox[3][3] ){
     }
   }
 
-  setFortranBoxSize(FortranHmat, FortranHmatI, &orthoRhombic);
+  setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic);
  
-  smallestBoxL = boxLx;
-  if (boxLy < smallestBoxL) smallestBoxL = boxLy;
-  if (boxLz < smallestBoxL) smallestBoxL = boxLz;
+  smallestBoxL = boxL[0];
+  if (boxL[1] < smallestBoxL) smallestBoxL = boxL[1];
+  if (boxL[2] > smallestBoxL) smallestBoxL = boxL[2];
 
   maxCutoff = smallestBoxL / 2.0;
 
@@ -107,37 +109,32 @@ void SimInfo::setBoxM( double theBox[3][3] ){
              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) {
+    if (rCut > (rList - 1.0)) {
       sprintf( painCave.errMsg,
-               "Error in recomputing LJ shifts based on new rcut\n");
-      painCave.isFatal = 1;
+               "New Box size is forcing LJ cutoff radius down to %lf\n",
+               rList - 1.0 );
+      painCave.isFatal = 0;
       simError();
+      rCut = rList - 1.0;
     }
-  }
+
+    if( ecr > (rList - 1.0) ){
+      sprintf( painCave.errMsg,
+               "New Box size is forcing electrostaticCutoffRadius "
+               "down to %lf\n"
+               "electrostaticSkinThickness is now %lf\n",
+               rList - 1.0, 0.05*(rList-1.0) );
+      painCave.isFatal = 0;
+      simError();      
+      ecr = maxCutoff;
+      est = 0.05 * ecr;
+    }
+
+    // At least one of the radii changed, so we need a refresh:
+    refreshSim();
+  }    
 }
  
 
@@ -153,7 +150,7 @@ void SimInfo::scaleBox(double scale) {
   double theBox[3][3];
   int i, j;
 
-  cerr << "Scaling box by " << scale << "\n";
+  // cerr << "Scaling box by " << scale << "\n";
 
   for(i=0; i<3; i++) 
     for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;
@@ -163,7 +160,8 @@ void SimInfo::calcHmatInv( void ) {
 }
 
 void SimInfo::calcHmatInv( void ) {
-
+  
+  int i,j;
   double smallDiag;
   double tol;
   double sanity[3][3];
@@ -173,12 +171,6 @@ void SimInfo::calcHmatInv( void ) {
   // Check the inverse to make sure it is sane:
 
   matMul3( Hmat, HmatInv, sanity );
-
-  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
   
@@ -271,7 +263,39 @@ void SimInfo::matVecMul3(double m[3][3], double inVec[
   outVec[1] = m[1][0]*a0 + m[1][1]*a1 + m[1][2]*a2;
   outVec[2] = m[2][0]*a0 + m[2][1]*a1 + m[2][2]*a2;
 }
+
+void SimInfo::transposeMat3(double in[3][3], double out[3][3]) {
+  double temp[3][3];
+  int i, j;
+
+  for (i = 0; i < 3; i++) {
+    for (j = 0; j < 3; j++) {
+      temp[j][i] = in[i][j];
+    }
+  }
+  for (i = 0; i < 3; i++) {
+    for (j = 0; j < 3; j++) {
+      out[i][j] = temp[i][j];
+    }
+  }
+}
   
+void SimInfo::printMat3(double A[3][3] ){
+
+  std::cerr 
+	    << "[ " << A[0][0] << ", " << A[0][1] << ", " << A[0][2] << " ]\n"
+	    << "[ " << A[1][0] << ", " << A[1][1] << ", " << A[1][2] << " ]\n"
+	    << "[ " << A[2][0] << ", " << A[2][1] << ", " << A[2][2] << " ]\n";
+}
+
+void SimInfo::printMat9(double A[9] ){
+
+  std::cerr 
+	    << "[ " << A[0] << ", " << A[1] << ", " << A[2] << " ]\n"
+	    << "[ " << A[3] << ", " << A[4] << ", " << A[5] << " ]\n"
+	    << "[ " << A[6] << ", " << A[7] << ", " << A[8] << " ]\n";
+}
+
 void SimInfo::calcBoxL( void ){
 
   double dx, dy, dz, dsq;
@@ -285,19 +309,19 @@ void SimInfo::calcBoxL( void ){
   
   dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
   dsq = dx*dx + dy*dy + dz*dz;
-  boxLx = sqrt( dsq );
+  boxL[0] = sqrt( dsq );
 
   // boxLy
   
   dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
   dsq = dx*dx + dy*dy + dz*dz;
-  boxLy = sqrt( dsq );
+  boxL[1] = sqrt( dsq );
 
   // boxLz
   
   dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
   dsq = dx*dx + dy*dy + dz*dz;
-  boxLz = sqrt( dsq );
+  boxL[2] = sqrt( dsq );
   
 }