--- trunk/OOPSE-2.0/src/restraints/Restraints.cpp	2004/09/24 04:16:43	1490
+++ trunk/OOPSE-2.0/src/restraints/Restraints.cpp	2004/11/23 22:48:31	1772
@@ -1,3 +1,5 @@
+// Thermodynamic integration is not multiprocessor friendly right now
+
 #include <iostream>
 #include <stdlib.h>
 #include <cstdio>
@@ -9,17 +11,23 @@ using namespace std;
 
 using namespace std;
 
-#include "Restraints.hpp"
-#include "SimInfo.hpp"
-#include "simError.h"
+#include "restraints/Restraints.hpp"
+#include "brains/SimInfo.hpp"
+#include "utils/simError.h"
+#include "io/basic_ifstrstream.hpp"
 
+#ifdef IS_MPI
+#include<mpi.h>
+#include "brains/mpiSimulation.hpp"
+#endif // is_mpi
+
 #define PI 3.14159265359
 #define TWO_PI 6.28318530718
 
 Restraints::Restraints(double lambdaVal, double lambdaExp){
   lambdaValue = lambdaVal;
   lambdaK = lambdaExp;
-
+  vector<double> resConsts;
   const char *jolt = " \t\n;,";
 
 #ifdef IS_MPI
@@ -32,10 +40,10 @@ Restraints::Restraints(double lambdaVal, double lambda
     
     if (!springs) { 
       sprintf(painCave.errMsg,
-	      "In Restraints: Unable to open HarmSpringConsts.txt for reading.\n"
-	      "\tDefault spring constants will be loaded. If you want to specify\n"
-	      "\tspring constants, include a three line HarmSpringConsts.txt file\n"
-	      "\tin the current directory.\n");
+	      "Unable to open HarmSpringConsts.txt for reading, so the\n"
+	      "\tdefault spring constants will be loaded. If you want\n"
+	      "\tto specify spring constants, include a three line\n"
+	      "\tHarmSpringConsts.txt file in the execution directory.\n");
       painCave.severity = OOPSE_WARNING;
       painCave.isFatal = 0;
       simError();   
@@ -47,22 +55,49 @@ Restraints::Restraints(double lambdaVal, double lambda
     } else  {
       
       springs.getline(inLine,999,'\n');
+      // the file is blank!
+      if (springs.eof()){
+      sprintf(painCave.errMsg,
+	      "HarmSpringConsts.txt file is not valid.\n"
+	      "\tThe file should contain four rows, the last three containing\n"
+	      "\ta label and the spring constant value. They should be listed\n"
+	      "\tin the following order: kDist (positional restrant), kTheta\n"
+	      "\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
+	      "\trestraint: rotation about the z-axis).\n");
+      painCave.severity = OOPSE_ERROR;
+      painCave.isFatal = 1;
+      simError();   
+      }
+      // read in spring constants and check to make sure it is a valid file
       springs.getline(inLine,999,'\n');
-      token = strtok(inLine,jolt);
-      token = strtok(NULL,jolt);
-      strcpy(inValue,token);
-      kDist = (atof(inValue));
-      springs.getline(inLine,999,'\n');
-      token = strtok(inLine,jolt);
-      token = strtok(NULL,jolt);
-      strcpy(inValue,token);
-      kTheta = (atof(inValue));
-      springs.getline(inLine,999,'\n');
-      token = strtok(inLine,jolt);
-      token = strtok(NULL,jolt);
-      strcpy(inValue,token);
-      kOmega = (atof(inValue));
-      springs.close();
+      while (!springs.eof()){
+	if (NULL != inLine){
+	  token = strtok(inLine,jolt);
+	  token = strtok(NULL,jolt);
+	  if (NULL != token){
+	    strcpy(inValue,token);
+	    resConsts.push_back(atof(inValue));
+	  }
+ 	}
+	springs.getline(inLine,999,'\n');
+      }
+      if (resConsts.size() == 3){
+	kDist = resConsts[0];
+	kTheta = resConsts[1];
+	kOmega = resConsts[2];
+      }
+      else {
+	sprintf(painCave.errMsg,
+		"HarmSpringConsts.txt file is not valid.\n"
+		"\tThe file should contain four rows, the last three containing\n"
+		"\ta label and the spring constant value. They should be listed\n"
+		"\tin the following order: kDist (positional restrant), kTheta\n"
+		"\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
+		"\trestraint: rotation about the z-axis).\n");
+	painCave.severity = OOPSE_ERROR;
+	painCave.isFatal = 1;
+	simError();  
+      }
     }
 #ifdef IS_MPI
   }
@@ -90,21 +125,21 @@ void Restraints::Calc_rVal(double position[3], int cur
 Restraints::~Restraints(){
 }
 
-void Restraints::Calc_rVal(double position[3], int currentMol){
-  delRx = position[0] - cofmPosX[currentMol];
-  delRy = position[1] - cofmPosY[currentMol];
-  delRz = position[2] - cofmPosZ[currentMol];
+void Restraints::Calc_rVal(double position[3], double refPosition[3]){
+  delRx = position[0] - refPosition[0];
+  delRy = position[1] - refPosition[1];
+  delRz = position[2] - refPosition[2];
 
   return;
 }
 
-void Restraints::Calc_body_thetaVal(double matrix[3][3], int currentMol){
-  ub0x = matrix[0][0]*uX0[currentMol] + matrix[0][1]*uY0[currentMol]
-    + matrix[0][2]*uZ0[currentMol];
-  ub0y = matrix[1][0]*uX0[currentMol] + matrix[1][1]*uY0[currentMol]
-    + matrix[1][2]*uZ0[currentMol];
-  ub0z = matrix[2][0]*uX0[currentMol] + matrix[2][1]*uY0[currentMol]
-    + matrix[2][2]*uZ0[currentMol];
+void Restraints::Calc_body_thetaVal(double matrix[3][3], double refUnit[3]){
+  ub0x = matrix[0][0]*refUnit[0] + matrix[0][1]*refUnit[1]
+    + matrix[0][2]*refUnit[2];
+  ub0y = matrix[1][0]*refUnit[0] + matrix[1][1]*refUnit[1]
+    + matrix[1][2]*refUnit[2];
+  ub0z = matrix[2][0]*refUnit[0] + matrix[2][1]*refUnit[1]
+    + matrix[2][2]*refUnit[2];
 
   normalize = sqrt(ub0x*ub0x + ub0y*ub0y + ub0z*ub0z);
   ub0x = ub0x/normalize;
@@ -159,11 +194,15 @@ double Restraints::Calc_Restraint_Forces(vector<StuntD
 double Restraints::Calc_Restraint_Forces(vector<StuntDouble*> vecParticles){
   double pos[3];
   double A[3][3];
+  double refPos[3];
+  double refVec[3];
   double tolerance;
   double tempPotent;
   double factor;
   double spaceTrq[3];
   double omegaPass;
+  GenericData* data;
+  DoubleGenericData* doubleData;
 
   tolerance = 5.72957795131e-7;
 
@@ -172,22 +211,96 @@ double Restraints::Calc_Restraint_Forces(vector<StuntD
   factor = 1 - pow(lambdaValue, lambdaK);
 
   for (i=0; i<vecParticles.size(); i++){
+    // obtain the current and reference positions
+    vecParticles[i]->getPos(pos);
+
+    data = vecParticles[i]->getProperty("refPosX");
+    if (data){
+      doubleData = dynamic_cast<DoubleGenericData*>(data);
+      if (!doubleData){
+	cerr << "Can't obtain refPosX from StuntDouble\n";
+	return 0.0;
+      }
+      else refPos[0] = doubleData->getData();
+    }
+    data = vecParticles[i]->getProperty("refPosY");
+    if (data){
+      doubleData = dynamic_cast<DoubleGenericData*>(data);
+      if (!doubleData){
+	cerr << "Can't obtain refPosY from StuntDouble\n";
+	return 0.0;
+      }
+      else refPos[1] = doubleData->getData();
+    }
+    data = vecParticles[i]->getProperty("refPosZ");
+    if (data){
+      doubleData = dynamic_cast<DoubleGenericData*>(data);
+      if (!doubleData){
+	cerr << "Can't obtain refPosZ from StuntDouble\n";
+	return 0.0;
+      }
+      else refPos[2] = doubleData->getData();
+    }
+
+    // calculate the displacement
+    Calc_rVal( pos, refPos );
+
+    // calculate the derivatives
+    dVdrx = -kDist*delRx;
+    dVdry = -kDist*delRy;
+    dVdrz = -kDist*delRz;
+    
+    // next we calculate the restraint forces
+    restraintFrc[0] = dVdrx;
+    restraintFrc[1] = dVdry;
+    restraintFrc[2] = dVdrz;
+    tempPotent = 0.5*kDist*(delRx*delRx + delRy*delRy + delRz*delRz);
+
+    // apply the lambda scaling factor to the forces
+    for (j = 0; j < 3; j++) restraintFrc[j] *= factor;
+
+    // and add the temporary force to the total force
+    vecParticles[i]->addFrc(restraintFrc);
+
+    // if the particle is directional, we accumulate the rot. restraints
     if (vecParticles[i]->isDirectional()){
-      vecParticles[i]->getPos(pos);
+ 
+      // get the current rotation matrix and reference vector
       vecParticles[i]->getA(A);
-      Calc_rVal( pos, i );
-      Calc_body_thetaVal( A, i );
+      
+      data = vecParticles[i]->getProperty("refVectorX");
+      if (data){
+	doubleData = dynamic_cast<DoubleGenericData*>(data);
+	if (!doubleData){
+	  cerr << "Can't obtain refVectorX from StuntDouble\n";
+	  return 0.0;
+	}
+	else refVec[0] = doubleData->getData();
+      }
+      data = vecParticles[i]->getProperty("refVectorY");
+      if (data){
+	doubleData = dynamic_cast<DoubleGenericData*>(data);
+	if (!doubleData){
+	  cerr << "Can't obtain refVectorY from StuntDouble\n";
+	  return 0.0;
+	}
+	else refVec[1] = doubleData->getData();
+      }
+      data = vecParticles[i]->getProperty("refVectorZ");
+      if (data){
+	doubleData = dynamic_cast<DoubleGenericData*>(data);
+	if (!doubleData){
+	  cerr << "Can't obtain refVectorZ from StuntDouble\n";
+	  return 0.0;
+	}
+	else refVec[2] = doubleData->getData();
+      }
+      
+      // calculate the theta and omega displacements
+      Calc_body_thetaVal( A, refVec );
       omegaPass = vecParticles[i]->getZangle();
       Calc_body_omegaVal( A, omegaPass );
 
-      if (omega > PI || omega < -PI)
-	cout << "oops... " << omega << "\n";
-
-      // first we calculate the derivatives
-      dVdrx = -kDist*delRx;
-      dVdry = -kDist*delRy;
-      dVdrz = -kDist*delRz;
-
       // uTx... and vTx... are the body-fixed z and y unit vectors
       uTx = 0.0;
       uTy = 0.0;
@@ -196,12 +309,12 @@ double Restraints::Calc_Restraint_Forces(vector<StuntD
       vTy = 1.0;
       vTz = 0.0;
 
-      dVdux = 0;
-      dVduy = 0;
-      dVduz = 0;
-      dVdvx = 0;
-      dVdvy = 0;
-      dVdvz = 0;
+      dVdux = 0.0;
+      dVduy = 0.0;
+      dVduz = 0.0;
+      dVdvx = 0.0;
+      dVdvy = 0.0;
+      dVdvz = 0.0;
 
       if (fabs(theta) > tolerance) {
 	dVdux = -(kTheta*theta/sin(theta))*ub0x;
@@ -215,12 +328,7 @@ double Restraints::Calc_Restraint_Forces(vector<StuntD
 	dVdvz = -(kOmega*omega/sin(omega))*vb0z;
       }
 
-      // next we calculate the restraint forces and torques
-      restraintFrc[0] = dVdrx;
-      restraintFrc[1] = dVdry;
-      restraintFrc[2] = dVdrz;
-      tempPotent = 0.5*kDist*(delRx*delRx + delRy*delRy + delRz*delRz);
-
+      // next we calculate the restraint torques
       restraintTrq[0] = 0.0;
       restraintTrq[1] = 0.0;
       restraintTrq[2] = 0.0;
@@ -238,13 +346,9 @@ double Restraints::Calc_Restraint_Forces(vector<StuntD
 	tempPotent += 0.5*(kTheta*theta*theta);
       }
 
-      for (j = 0; j < 3; j++) {
-	restraintFrc[j] *= factor;
-	restraintTrq[j] *= factor;
-      }
+      // apply the lambda scaling factor to these torques
+      for (j = 0; j < 3; j++) restraintTrq[j] *= factor;
 
-      harmPotent += tempPotent;
-
       // now we need to convert from body-fixed torques to space-fixed torques
       spaceTrq[0] = A[0][0]*restraintTrq[0] + A[1][0]*restraintTrq[1] 
 	+ A[2][0]*restraintTrq[2];
@@ -253,141 +357,37 @@ double Restraints::Calc_Restraint_Forces(vector<StuntD
       spaceTrq[2] = A[0][2]*restraintTrq[0] + A[1][2]*restraintTrq[1] 
 	+ A[2][2]*restraintTrq[2];
 
-      // now it's time to pass these temporary forces and torques
-      // to the total forces and torques
-      vecParticles[i]->addFrc(restraintFrc);
+      // now pass this temporary torque vector to the total torque
       vecParticles[i]->addTrq(spaceTrq);
     }
-  }
 
-  // and we can return the appropriately scaled potential energy
+    // update the total harmonic potential with this object's contribution
+    harmPotent += tempPotent;
+  }
+  
+  // we can finish by returning the appropriately scaled potential energy
   tempPotent = harmPotent * factor;
   return tempPotent;
 }
 
-void Restraints::Store_Init_Info(vector<StuntDouble*> vecParticles){
-  double pos[3];
-  double A[3][3];
-  double RfromQ[3][3];
-  double quat0, quat1, quat2, quat3;
-  double dot;
-//   char *token;
-//   char fileName[200];
-//   char angleName[200];
-//   char inLine[1000];
-//   char inValue[200];
-  const char *delimit = " \t\n;,";
+void Restraints::Write_zAngle_File(vector<StuntDouble*> vecParticles,
+				   int currTime,
+				   int nIntObj){
 
-  //open the idealCrystal.in file and zAngle.ang file
-  strcpy(fileName, "idealCrystal.in");
-  strcpy(angleName, "zAngle.ang");
-  
-  ifstream crystalIn(fileName);
-  ifstream angleIn(angleName);
+  char zOutName[200];
 
-  if (!crystalIn) { 
-    sprintf(painCave.errMsg,
-	    "Restraints Error: Unable to open idealCrystal.in for reading.\n"
-	    "\tMake sure a reference crystal file is in the current directory.\n");
-    painCave.isFatal = 1;
-    simError();   
-    
-    return;
-  }
+  std::cerr << nIntObj << " is the number of integrable objects\n";
 
-  if (!angleIn) { 
-    sprintf(painCave.errMsg,
-	    "Restraints Warning: The lack of a zAngle.ang file is mildly\n"
-	    "\tunsettling... This means the simulation is starting from the\n"
-	    "\tidealCrystal.in reference configuration, so the omega values\n"
-	    "\twill all be set to zero. If this is not the case, you should\n"
-	    "\tquestion your results.\n");
-    painCave.isFatal = 0;
-    simError();   
-  }
-
-  // A rather specific reader for OOPSE .eor files...
-  // Let's read in the perfect crystal file
-  crystalIn.getline(inLine,999,'\n');
-  crystalIn.getline(inLine,999,'\n');
+  //#ifndef IS_MPI
   
-  for (i=0; i<vecParticles.size(); i++) {
-    crystalIn.getline(inLine,999,'\n');
-    token = strtok(inLine,delimit);
-    token = strtok(NULL,delimit);
-    strcpy(inValue,token);
-    cofmPosX.push_back(atof(inValue));
-    token = strtok(NULL,delimit);
-    strcpy(inValue,token);
-    cofmPosY.push_back(atof(inValue));
-    token = strtok(NULL,delimit);
-    strcpy(inValue,token);
-    cofmPosZ.push_back(atof(inValue));
-    token = strtok(NULL,delimit);
-    token = strtok(NULL,delimit);
-    token = strtok(NULL,delimit);
-    token = strtok(NULL,delimit);
-    strcpy(inValue,token);
-    quat0 = atof(inValue);
-    token = strtok(NULL,delimit);
-    strcpy(inValue,token);
-    quat1 = atof(inValue);
-    token = strtok(NULL,delimit);
-    strcpy(inValue,token);
-    quat2 = atof(inValue);
-    token = strtok(NULL,delimit);
-    strcpy(inValue,token);
-    quat3 = atof(inValue);
-
-    // now build the rotation matrix and find the unit vectors
-    RfromQ[0][0] = quat0*quat0 + quat1*quat1 - quat2*quat2 - quat3*quat3;
-    RfromQ[0][1] = 2*(quat1*quat2 + quat0*quat3);
-    RfromQ[0][2] = 2*(quat1*quat3 - quat0*quat2);
-    RfromQ[1][0] = 2*(quat1*quat2 - quat0*quat3);
-    RfromQ[1][1] = quat0*quat0 - quat1*quat1 + quat2*quat2 - quat3*quat3;
-    RfromQ[1][2] = 2*(quat2*quat3 + quat0*quat1);
-    RfromQ[2][0] = 2*(quat1*quat3 + quat0*quat2);
-    RfromQ[2][1] = 2*(quat2*quat3 - quat0*quat1);
-    RfromQ[2][2] = quat0*quat0 - quat1*quat1 - quat2*quat2 + quat3*quat3;
-    
-    normalize = sqrt(RfromQ[2][0]*RfromQ[2][0] + RfromQ[2][1]*RfromQ[2][1] 
-		     + RfromQ[2][2]*RfromQ[2][2]);
-    uX0.push_back(RfromQ[2][0]/normalize);
-    uY0.push_back(RfromQ[2][1]/normalize);
-    uZ0.push_back(RfromQ[2][2]/normalize);
-
-    normalize = sqrt(RfromQ[1][0]*RfromQ[1][0] + RfromQ[1][1]*RfromQ[1][1]
-		     + RfromQ[1][2]*RfromQ[1][2]);
-    vX0.push_back(RfromQ[1][0]/normalize);
-    vY0.push_back(RfromQ[1][1]/normalize);
-    vZ0.push_back(RfromQ[1][2]/normalize);
-  }
-
-  // now we can read in the zAngle.ang file
-  if (angleIn){
-    angleIn.getline(inLine,999,'\n');
-    for (i=0; i<vecParticles.size(); i++) {
-      angleIn.getline(inLine,999,'\n');
-      token = strtok(inLine,delimit);
-      strcpy(inValue,token);
-      vecParticles[i]->setZangle(atof(inValue));
-    }
-  }
-
-  return;
-}
-
-void Restraints::Write_zAngle_File(vector<StuntDouble*> vecParticles){
-
-  char zOutName[200];
-
   strcpy(zOutName,"zAngle.ang");
-
+  
   ofstream angleOut(zOutName);
-  angleOut << "This file contains the omega values for the .eor file\n";
+  angleOut << currTime << ": omega values at this time\n";
   for (i=0; i<vecParticles.size(); i++) {
     angleOut << vecParticles[i]->getZangle() << "\n";
   }
+
   return;
 }