--- trunk/SHAPES/GridBuilder.cpp	2004/06/18 19:40:31	1280
+++ trunk/SHAPES/GridBuilder.cpp	2004/06/25 17:51:51	1305
@@ -1,60 +1,91 @@
 #include "GridBuilder.hpp"
-#include "MatVec3.h"
 #define PI 3.14159265359
 
 
-GridBuilder::GridBuilder(RigidBody* rb, int bandWidth) {
+GridBuilder::GridBuilder(RigidBody* rb, int gridWidth) {
   rbMol = rb;
-  bandwidth = bandWidth;
-  thetaStep = PI / bandwidth;
+  gridwidth = gridWidth;
+  thetaStep = PI / gridwidth;
   thetaMin = thetaStep / 2.0;
   phiStep = thetaStep * 2.0;
-	
-  //zero out the rot mats
-  for (i=0; i<3; i++) {
-    for (j=0; j<3; j++) {
-      rotX[i][j] = 0.0;
-      rotZ[i][j] = 0.0;
-      rbMatrix[i][j] = 0.0;
-    }
-  }
 }
 
 GridBuilder::~GridBuilder() {
 }
 
-void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid, vector<double> sGrid,
-                              vector<double> epsGrid){
+void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid, 
+                              vector<double> sGrid, vector<double> epsGrid){
+  ofstream sigmaOut("sigma.grid");
+  ofstream sOut("s.grid");
+  ofstream epsOut("eps.grid");
   double startDist;
+  double phiVal;
+  double thetaVal;
+  double sigTemp, sTemp, epsTemp, sigProbe;
   double minDist = 10.0; //minimum start distance
 	
+  sigList = sigmaGrid;
+  sList = sGrid;
+  epsList = epsGrid;
   forcefield = forceField;
+  
+  //load the probe atom parameters
+  switch(forcefield){
+    case 1:{
+      rparHe = 1.4800;
+      epsHe = -0.021270;
+    }; break;
+    case 2:{
+      rparHe = 1.14;
+      epsHe = 0.0203;
+    }; break;
+    case 3:{
+      rparHe = 2.28;
+      epsHe = 0.020269601874;
+    }; break;
+    case 4:{
+      rparHe = 2.5560;
+      epsHe = 0.0200;
+    }; break;
+    case 5:{
+      rparHe = 1.14;
+      epsHe = 0.0203;
+    }; break;
+  }
     
-  //first determine the start distance - we always start at least minDist away
+  if (rparHe < 2.2)
+    sigProbe = 2*rparHe/1.12246204831;
+  else
+    sigProbe = rparHe;
+  
+  //determine the start distance - we always start at least minDist away
   startDist = rbMol->findMaxExtent() + minDist;
   if (startDist < minDist)
     startDist = minDist;
-	
-  initBody();
-  for (i=0; i<bandwidth; i++){		
-    for (j=0; j<bandwidth; j++){
+
+  //set the initial orientation of the body and loop over theta values
+
+  for (k =0; k < gridwidth; k++) {
+    thetaVal = thetaMin + k*thetaStep;
+    for (j=0; j < gridwidth; j++) {
+      phiVal = j*phiStep + 0.5*PI;
+
+      rbMol->setEuler(0.0, thetaVal, phiVal);
+
       releaseProbe(startDist);
 
-      sigmaGrid.push_back(sigDist);
-      sGrid.push_back(sDist);
-      epsGrid.push_back(epsVal);
-			
-      stepPhi(phiStep);
+      //translate the values to sigma, s, and epsilon of the rigid body
+      sigTemp = 2*sigDist - sigProbe;
+      sTemp = (2*(sDist - sigDist))/(0.122462048309) - sigProbe;
+      epsTemp = pow(epsVal, 2)/fabs(epsHe);
+      
+      sigList.push_back(sigTemp);
+      sList.push_back(sTemp);
+      epsList.push_back(epsTemp);
     }
-    stepTheta(thetaStep);
   }		
 }
 
-void GridBuilder::initBody(){
-  //set up the rigid body in the starting configuration
-  stepTheta(thetaMin);
-}
-
 void GridBuilder::releaseProbe(double farPos){
   int tooClose;
   double tempPotEnergy;
@@ -67,8 +98,7 @@ void GridBuilder::releaseProbe(double farPos){
   tooClose = 0;
   epsVal = 0;
   rhoStep = 0.1; //the distance the probe atom moves between steps
-	
-	
+		
   while (!tooClose){
     calcEnergy();
     potProgress.push_back(potEnergy);
@@ -106,39 +136,89 @@ void GridBuilder::calcEnergy(){
 }
 
 void GridBuilder::calcEnergy(){
-	
-}
+  double rXij, rYij, rZij;
+  double rijSquared;
+  double rValSquared, rValPowerSix;
+  double atomRpar, atomEps;
+  double rbAtomPos[3];
+    
+  potEnergy = 0.0;
 
-void GridBuilder::stepTheta(double increment){
-  //zero out the euler angles
-  for (i=0; i<3; i++)
-    angles[i] = 0.0;
-	
-  //the second euler angle is for rotation about the x-axis (we use the zxz convention)
-  angles[1] = increment;
-	
-  //obtain the rotation matrix through the rigid body class
-  rbMol->doEulerToRotMat(angles, rotX);
-	
-  //rotate the rigid body
-  rbMol->getA(rbMatrix);
-  matMul3(rotX, rbMatrix, rotatedMat);
-  rbMol->setA(rotatedMat);	
-}
+  for(i=0; i<rbMol->getNumAtoms(); i++){
+    rbMol->getAtomPos(rbAtomPos, i);
+    
+    rXij = rbAtomPos[0];
+    rYij = rbAtomPos[1];
+    rZij = rbAtomPos[2] - probeCoor;
+    
+    rijSquared = rXij * rXij + rYij * rYij + rZij * rZij;
+    
+    //in the interest of keeping the code more compact, we are being less 
+    //efficient by placing a switch statement in the calculation loop
+    switch(forcefield){
+      case 1:{
+        //we are using the CHARMm force field
+        atomRpar = rbMol->getAtomRpar(i);
+        atomEps = rbMol->getAtomEps(i);
+        
+        rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared);
+        rValPowerSix = rValSquared * rValSquared * rValSquared;
+        potEnergy += sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0));
+      }; break;
+      
+      case 2:{
+        //we are using the AMBER force field
+        atomRpar = rbMol->getAtomRpar(i);
+        atomEps = rbMol->getAtomEps(i);
+        
+        rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared);
+        rValPowerSix = rValSquared * rValSquared * rValSquared;
+        potEnergy += sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0));
+      }; break;
+      
+      case 3:{
+        //we are using Allen-Tildesley LJ parameters
+        atomRpar = rbMol->getAtomRpar(i);
+        atomEps = rbMol->getAtomEps(i);
+        
+        rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (4*rijSquared);
+        rValPowerSix = rValSquared * rValSquared * rValSquared;
+        potEnergy += 4*sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 1.0));
+        
+      }; break;
+      
+      case 4:{
+        //we are using the OPLS force field
+        atomRpar = rbMol->getAtomRpar(i);
+        atomEps = rbMol->getAtomEps(i);
+        
+        rValSquared = (pow(sqrt(rparHe+atomRpar),2)) / (rijSquared);
+        rValPowerSix = rValSquared * rValSquared * rValSquared;
+        potEnergy += 4*sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 1.0));
+      }; break;
+      
+      case 5:{
+        //we are using the GAFF force field
+        atomRpar = rbMol->getAtomRpar(i);
+        atomEps = rbMol->getAtomEps(i);
+        
+        rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared);
+        rValPowerSix = rValSquared * rValSquared * rValSquared;
+        potEnergy += sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0));
+      }; break;
+    }    
+  }
+} 
 
-void GridBuilder::stepPhi(double increment){
-  //zero out the euler angles
-  for (i=0; i<3; i++)
-    angles[i] = 0.0;
-	
-  //the phi euler angle is for rotation about the z-axis (we use the zxz convention)
-  angles[0] = increment;
-	
-  //obtain the rotation matrix through the rigid body class
-  rbMol->doEulerToRotMat(angles, rotZ);
-	
-  //rotate the rigid body
-  rbMol->getA(rbMatrix);
-  matMul3(rotZ, rbMatrix, rotatedMat);
-  rbMol->setA(rotatedMat);	
+void GridBuilder::printGridFiles(){
+  ofstream sigmaOut("sigma.grid");
+  ofstream sOut("s.grid");
+  ofstream epsOut("eps.grid");
+  
+  for (k=0; k<sigList.size(); k++){
+    sigmaOut << sigList[k] << "\n0\n";
+    sOut << sList[k] << "\n0\n";    
+    epsOut << epsList[k] << "\n0\n";
+  }
 }
+