--- trunk/OOPSE/libmdtools/SimInfo.cpp	2003/04/11 15:16:59	490
+++ trunk/OOPSE/libmdtools/SimInfo.cpp	2003/06/30 22:04:01	568
@@ -1,5 +1,6 @@
 #include <cstdlib>
 #include <cstring>
+#include <cmath>
 
 
 #include "SimInfo.hpp"
@@ -39,19 +40,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);
+
+  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",
@@ -91,12 +101,160 @@ void SimInfo::getBox(double theBox[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);
+ 
+  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::getBox(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;
+
+  // 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++;
+    }
+  }
+}
+
+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];
+
+  // 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));
+  
+
+}
+
+
 int SimInfo::getNDF(){
   int ndf_local, ndf;