--- trunk/OOPSE/libmdtools/SimInfo.cpp	2003/07/15 21:34:56	618
+++ trunk/OOPSE/libmdtools/SimInfo.cpp	2004/02/06 18:58:06	1031
@@ -1,6 +1,6 @@
-#include <cstdlib>
-#include <cstring>
-#include <cmath>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
 
 #include <iostream>
 using namespace std;
@@ -20,33 +20,67 @@ inline double roundMe( double x ){
   return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
 }
 	  
+inline double min( double a, double b ){
+  return (a < b ) ? a : b;
+}
 
 SimInfo* currentInfo;
 
 SimInfo::SimInfo(){
   excludes = NULL;
   n_constraints = 0;
+  nZconstraints = 0;
   n_oriented = 0;
   n_dipoles = 0;
   ndf = 0;
   ndfRaw = 0;
+  nZconstraints = 0;
   the_integrator = NULL;
   setTemp = 0;
   thermalTime = 0.0;
+  currentTime = 0.0;
   rCut = 0.0;
   ecr = 0.0;
+  est = 0.0;
 
+  haveRcut = 0;
+  haveEcr = 0;
+  boxIsInit = 0;
+  
+  resetTime = 1e99;
+
+  orthoTolerance = 1E-6;
+  useInitXSstate = true;
+
   usePBC = 0;
   useLJ = 0; 
   useSticky = 0;
-  useDipole = 0;
+  useCharges = 0;
+  useDipoles = 0;
   useReactionField = 0;
   useGB = 0;
   useEAM = 0;
 
+  myConfiguration = new SimState();
+
+  has_minimizer = false;
+  the_minimizer =NULL;
+
   wrapMeSimInfo( this );
 }
 
+
+SimInfo::~SimInfo(){
+
+  delete myConfiguration;
+
+  map<string, GenericData*>::iterator i;
+  
+  for(i = properties.begin(); i != properties.end(); i++)
+    delete (*i).second;
+    
+}
+
 void SimInfo::setBox(double newBox[3]) {
   
   int i, j;
@@ -65,8 +99,7 @@ void SimInfo::setBoxM( double theBox[3][3] ){
 
 void SimInfo::setBoxM( double theBox[3][3] ){
   
-  int i, j, status;
-  double smallestBoxL, maxCutoff;
+  int i, j;
   double FortranHmat[9]; // to preserve compatibility with Fortran the
                          // ordering in the array is as follows:
                          // [ 0 3 6 ]
@@ -74,16 +107,11 @@ void SimInfo::setBoxM( double theBox[3][3] ){
                          // [ 2 5 8 ]
   double FortranHmatInv[9]; // the inverted Hmat (for Fortran);
 
+  if( !boxIsInit ) boxIsInit = 1;
 
   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";
-
   calcBoxL();
   calcHmatInv();
 
@@ -96,48 +124,6 @@ void SimInfo::setBoxM( double theBox[3][3] ){
 
   setFortranBoxSize(FortranHmat, FortranHmatInv, &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( ecr > maxCutoff ){
-
-    sprintf( painCave.errMsg,
-             "New Box size is forcing electrostatic cutoff radius "
-	     "down to %lf\n",
-             maxCutoff );
-    painCave.isFatal = 0;
-    simError();
-
-    ecr = maxCutoff;
-    est = 0.05 * ecr;
-
-    refreshSim();
-  }
-    
 }
  
 
@@ -164,6 +150,7 @@ void SimInfo::calcHmatInv( void ) {
 
 void SimInfo::calcHmatInv( void ) {
   
+  int oldOrtho;
   int i,j;
   double smallDiag;
   double tol;
@@ -171,28 +158,46 @@ void SimInfo::calcHmatInv( void ) {
 
   invertMat3( Hmat, HmatInv );
 
-  // Check the inverse to make sure it is sane:
-
-  matMul3( Hmat, HmatInv, sanity );
-    
   // check to see if Hmat is orthorhombic
   
-  smallDiag = Hmat[0][0];
-  if(smallDiag > Hmat[1][1]) smallDiag = Hmat[1][1];
-  if(smallDiag > Hmat[2][2]) smallDiag = Hmat[2][2];
-  tol = smallDiag * 1E-6;
+  oldOrtho = orthoRhombic;
 
+  smallDiag = fabs(Hmat[0][0]);
+  if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]);
+  if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]);
+  tol = smallDiag * orthoTolerance;
+
   orthoRhombic = 1;
   
   for (i = 0; i < 3; i++ ) {
     for (j = 0 ; j < 3; j++) {
       if (i != j) {
         if (orthoRhombic) {
-          if (Hmat[i][j] >= tol) orthoRhombic = 0;
+          if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0;
         }        
       }
     }
   }
+
+  if( oldOrtho != orthoRhombic ){
+    
+    if( orthoRhombic ){
+      sprintf( painCave.errMsg,
+	       "Hmat is switching from Non-Orthorhombic to Orthorhombic Box.\n"
+	       "\tIf this is a bad thing, change the orthoBoxTolerance\n"
+               "\tvariable ( currently set to %G ).\n",
+	       orthoTolerance);
+      simError();
+    }
+    else {
+      sprintf( painCave.errMsg,
+	       "Hmat is switching from Orthorhombic to Non-Orthorhombic Box.\n"
+	       "\tIf this is a bad thing, change the orthoBoxTolerance\n"
+               "\tvariable ( currently set to %G ).\n",
+	       orthoTolerance);
+      simError();
+    }
+  }
 }
 
 double SimInfo::matDet3(double a[3][3]) {
@@ -299,10 +304,26 @@ void SimInfo::printMat9(double A[9] ){
 	    << "[ " << A[6] << ", " << A[7] << ", " << A[8] << " ]\n";
 }
 
+
+void SimInfo::crossProduct3(double a[3],double b[3], double out[3]){
+
+      out[0] = a[1] * b[2] - a[2] * b[1];
+      out[1] = a[2] * b[0] - a[0] * b[2] ;
+      out[2] = a[0] * b[1] - a[1] * b[0];
+      
+}
+
+double SimInfo::dotProduct3(double a[3], double b[3]){
+  return a[0]*b[0] + a[1]*b[1]+ a[2]*b[2];
+}
+
+double SimInfo::length3(double a[3]){
+  return sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]);
+}
+
 void SimInfo::calcBoxL( void ){
 
   double dx, dy, dz, dsq;
-  int i;
 
   // boxVol = Determinant of Hmat
 
@@ -312,26 +333,67 @@ 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 );
+  //maxCutoff = 0.5 * boxL[0];
 
   // 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 );
+  //if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
 
+
   // 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 );
+  //if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
+
+  //calculate the max cutoff
+  maxCutoff =  calcMaxCutOff(); 
   
+  checkCutOffs();
+
 }
+
+
+double SimInfo::calcMaxCutOff(){
+
+  double ri[3], rj[3], rk[3];
+  double rij[3], rjk[3], rki[3];
+  double minDist;
+
+  ri[0] = Hmat[0][0];
+  ri[1] = Hmat[1][0];
+  ri[2] = Hmat[2][0];
+
+  rj[0] = Hmat[0][1];
+  rj[1] = Hmat[1][1];
+  rj[2] = Hmat[2][1];
+
+  rk[0] = Hmat[0][2];
+  rk[1] = Hmat[1][2];
+  rk[2] = Hmat[2][2];
+  
+  crossProduct3(ri,rj, rij);
+  distXY = dotProduct3(rk,rij) / length3(rij);
 
+  crossProduct3(rj,rk, rjk);
+  distYZ = dotProduct3(ri,rjk) / length3(rjk);
 
+  crossProduct3(rk,ri, rki);
+  distZX = dotProduct3(rj,rki) / length3(rki);
+
+  minDist = min(min(distXY, distYZ), distZX);
+  return minDist/2;
+  
+}
+
 void SimInfo::wrapVector( double thePos[3] ){
 
-  int i, j, k;
+  int i;
   double scaled[3];
 
   if( !orthoRhombic ){
@@ -369,7 +431,7 @@ int SimInfo::getNDF(){
 
 
 int SimInfo::getNDF(){
-  int ndf_local, ndf;
+  int ndf_local;
   
   ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints;
 
@@ -379,13 +441,13 @@ int SimInfo::getNDF(){
   ndf = ndf_local;
 #endif
 
-  ndf = ndf - 3;
+  ndf = ndf - 3 - nZconstraints;
 
   return ndf;
 }
 
 int SimInfo::getNDFraw() {
-  int ndfRaw_local, ndfRaw;
+  int ndfRaw_local;
 
   // Raw degrees of freedom that we have to set
   ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
@@ -398,24 +460,33 @@ int SimInfo::getNDFraw() {
 
   return ndfRaw;
 }
- 
+
+int SimInfo::getNDFtranslational() {
+  int ndfTrans_local;
+
+  ndfTrans_local = 3 * n_atoms - n_constraints;
+
+#ifdef IS_MPI
+  MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
+#else
+  ndfTrans = ndfTrans_local;
+#endif
+
+  ndfTrans = ndfTrans - 3 - nZconstraints;
+
+  return ndfTrans;
+}
+
 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.rlist = rList;
-  fInfo.rcut = rCut;
-
-  if( useDipole ){
-    fInfo.rrf = ecr;
-    fInfo.rt = ecr - est;
+  if( useDipoles ){
     if( useReactionField )fInfo.dielect = dielectric;
   }
 
@@ -424,10 +495,11 @@ 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_charges = useCharges;
+  fInfo.SIM_uses_dipoles = useDipoles;
   //fInfo.SIM_uses_dipoles = 0;
-  //fInfo.SIM_uses_RF = useReactionField;
-  fInfo.SIM_uses_RF = 0;
+  fInfo.SIM_uses_RF = useReactionField;
+  //fInfo.SIM_uses_RF = 0;
   fInfo.SIM_uses_GB = useGB;
   fInfo.SIM_uses_EAM = useEAM;
 
@@ -461,6 +533,134 @@ void SimInfo::refreshSim(){
 
   this->ndf = this->getNDF();
   this->ndfRaw = this->getNDFraw();
+  this->ndfTrans = this->getNDFtranslational();
+}
 
+void SimInfo::setDefaultRcut( double theRcut ){
+
+  haveRcut = 1;
+  rCut = theRcut;
+
+  ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
+
+  notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
 }
 
+void SimInfo::setDefaultEcr( double theEcr ){
+
+  haveEcr = 1;
+  ecr = theEcr;
+  
+  ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
+
+  notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
+}
+
+void SimInfo::setDefaultEcr( double theEcr, double theEst ){
+
+  est = theEst;
+  setDefaultEcr( theEcr );
+}
+
+
+void SimInfo::checkCutOffs( void ){
+  
+  if( boxIsInit ){
+    
+    //we need to check cutOffs against the box
+    
+    if( rCut > maxCutoff ){
+      sprintf( painCave.errMsg,
+	       "Box size is too small for the long range cutoff radius, "
+	       "%G, at time %G\n"
+	       "\t[ %G %G %G ]\n"
+	       "\t[ %G %G %G ]\n"
+	       "\t[ %G %G %G ]\n",
+	       rCut, currentTime,
+	       Hmat[0][0], Hmat[0][1], Hmat[0][2],
+	       Hmat[1][0], Hmat[1][1], Hmat[1][2],
+	       Hmat[2][0], Hmat[2][1], Hmat[2][2]);
+      painCave.isFatal = 1;
+      simError();
+    }
+    
+    if( haveEcr ){
+      if( ecr > maxCutoff ){
+	sprintf( painCave.errMsg,
+		 "Box size is too small for the electrostatic cutoff radius, "
+		 "%G, at time %G\n"
+		 "\t[ %G %G %G ]\n"
+		 "\t[ %G %G %G ]\n"
+		 "\t[ %G %G %G ]\n",
+		 ecr, currentTime,
+		 Hmat[0][0], Hmat[0][1], Hmat[0][2],
+		 Hmat[1][0], Hmat[1][1], Hmat[1][2],
+		 Hmat[2][0], Hmat[2][1], Hmat[2][2]);
+	painCave.isFatal = 1;
+	simError();
+      }
+    }
+  } else {
+    // initialize this stuff before using it, OK?
+    sprintf( painCave.errMsg,
+             "Trying to check cutoffs without a box.\n"
+             "\tOOPSE should have better programmers than that.\n" );
+    painCave.isFatal = 1;
+    simError();      
+  }
+  
+}
+
+void SimInfo::addProperty(GenericData* prop){
+
+  map<string, GenericData*>::iterator result;
+  result = properties.find(prop->getID());
+  
+  //we can't simply use  properties[prop->getID()] = prop,
+  //it will cause memory leak if we already contain a propery which has the same name of prop
+  
+  if(result != properties.end()){
+    
+    delete (*result).second;
+    (*result).second = prop;
+      
+  }
+  else{
+
+    properties[prop->getID()] = prop;
+
+  }
+    
+}
+
+GenericData* SimInfo::getProperty(const string& propName){
+ 
+  map<string, GenericData*>::iterator result;
+  
+  //string lowerCaseName = ();
+  
+  result = properties.find(propName);
+  
+  if(result != properties.end()) 
+    return (*result).second;  
+  else   
+    return NULL;  
+}
+
+vector<GenericData*> SimInfo::getProperties(){
+
+  vector<GenericData*> result;
+  map<string, GenericData*>::iterator i;
+  
+  for(i = properties.begin(); i != properties.end(); i++)
+    result.push_back((*i).second);
+    
+  return result;
+}
+
+double SimInfo::matTrace3(double m[3][3]){
+  double trace;
+  trace = m[0][0] + m[1][1] + m[2][2];
+
+  return trace;
+}