--- trunk/OOPSE/libmdtools/SimInfo.cpp	2003/06/30 22:04:01	568
+++ trunk/OOPSE/libmdtools/SimInfo.cpp	2003/07/09 22:14:06	586
@@ -2,6 +2,8 @@
 #include <cstring>
 #include <cmath>
 
+#include <iostream>
+using namespace std;
 
 #include "SimInfo.hpp"
 #define __C
@@ -14,6 +16,11 @@
 #include "mpiSimulation.hpp"
 #endif
 
+inline double roundMe( double x ){
+  return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
+}
+	  
+
 SimInfo* currentInfo;
 
 SimInfo::SimInfo(){
@@ -40,65 +47,18 @@ void SimInfo::setBox(double newBox[3]) {
 }
 
 void SimInfo::setBox(double newBox[3]) {
-
-  double smallestBoxL, maxCutoff;
-  int status;
+  
   int i;
+  double tempMat[9];
 
-  for(i=0; i<9; i++) Hmat[i] = 0.0;;
+  for(i=0; i<9; i++) tempMat[i] = 0.0;;
 
-  Hmat[0] = newBox[0];
-  Hmat[4] = newBox[1];
-  Hmat[8] = newBox[2];
+  tempMat[0] = newBox[0];
+  tempMat[4] = newBox[1];
+  tempMat[8] = newBox[2];
 
-  calcHmatI();
-  calcBoxL();
+  setBoxM( tempMat );
 
-  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] ){
@@ -107,10 +67,19 @@ void SimInfo::setBoxM( double theBox[9] ){
   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);
+
+
+  setFortranBoxSize(Hmat, HmatI, &orthoRhombic);
  
   smallestBoxL = boxLx;
   if (boxLy < smallestBoxL) smallestBoxL = boxLy;
@@ -158,18 +127,33 @@ void SimInfo::setBoxM( double theBox[9] ){
 }
  
 
-void SimInfo::getBox(double theBox[9]) {
+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;
 
@@ -201,6 +185,41 @@ void SimInfo::calcHmatI( void ) {
       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 ){
@@ -241,17 +260,42 @@ 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));
-  
-
+  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];
+  }
+    
+    
 }
 
 
@@ -297,10 +341,6 @@ void SimInfo::refreshSim(){
   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;