--- trunk/SHAPES/GridBuilder.cpp	2004/06/18 19:01:40	1279
+++ trunk/SHAPES/GridBuilder.cpp	2004/06/21 15:10:29	1282
@@ -4,142 +4,248 @@ GridBuilder::GridBuilder(RigidBody* rb, int bandWidth)
 
 
 GridBuilder::GridBuilder(RigidBody* rb, int bandWidth) {
-	rbMol = rb;
-	bandwidth = bandWidth;
-	thetaStep = PI / bandwidth;
-	thetaMin = thetaStep / 2.0;
-	phiStep = thetaStep * 2.0;
+  rbMol = rb;
+  bandwidth = bandWidth;
+  thetaStep = PI / bandwidth;
+  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;
-		}
-	}
+  //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){
-	double startDist;
-	double minDist = 10.0; //minimum start distance
+                              vector<double> epsGrid){
+  ofstream sigmaOut("sigma.grid");
+  ofstream sOut("s.grid");
+  ofstream epsOut("eps.grid");
+  double startDist;
+  double phiVal;
+  double thetaVal;
+  double minDist = 10.0; //minimum start distance
 	
-	//first 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++){
-			releaseProbe(startDist);
+  sList = sGrid;
+  sigList = sigmaGrid;
+  epsList = epsGrid;
+  forcefield = forceField;
+    
+  //first determine the start distance - we always start at least minDist away
+  startDist = rbMol->findMaxExtent() + minDist;
+  if (startDist < minDist)
+    startDist = minDist;
 
-			sigmaGrid.push_back(sigDist);
-			sGrid.push_back(sDist);
-			epsGrid.push_back(epsVal);
-			
-			stepPhi(phiStep);
-		}
-		stepTheta(thetaStep);
-	}
-		
-}
+  //set the initial orientation of the body and loop over theta values
+  phiVal = 0.0;
+  thetaVal = thetaMin;
+  rotBody(phiVal, thetaVal);
+  for (k=0; k<bandwidth; k++){	
+	//loop over phi values starting with phi = 0.0
+    for (j=0; j<bandwidth; j++){
+      releaseProbe(startDist);
 
-void GridBuilder::initBody(){
-	//set up the rigid body in the starting configuration
-	stepTheta(thetaMin);
+      sigList.push_back(sigDist);
+      sList.push_back(sDist);
+      epsList.push_back(epsVal);
+      
+      phiVal += phiStep;
+      rotBody(phiVal, thetaVal);
+    }
+    phiVal = 0.0;
+    thetaVal += thetaStep;
+    rotBody(phiVal, thetaVal);
+    printf("step theta %i\n",k);
+  }		
 }
 
 void GridBuilder::releaseProbe(double farPos){
-	int tooClose;
-	double tempPotEnergy;
-	double interpRange;
-	double interpFrac;
+  int tooClose;
+  double tempPotEnergy;
+  double interpRange;
+  double interpFrac;
 	
-	probeCoor = farPos;
-	potProgress.clear();
-	distProgress.clear();
-	tooClose = 0;
-	epsVal = 0;
-	rhoStep = 0.1; //the distance the probe atom moves between steps
+  probeCoor = farPos;
+  potProgress.clear();
+  distProgress.clear();
+  tooClose = 0;
+  epsVal = 0;
+  rhoStep = 0.1; //the distance the probe atom moves between steps
 	
 	
-	while (!tooClose){
-		calcEnergy();
-		potProgress.push_back(potEnergy);
-		distProgress.push_back(probeCoor);
+  while (!tooClose){
+    calcEnergy();
+    potProgress.push_back(potEnergy);
+    distProgress.push_back(probeCoor);
 		
-		//if we've reached a new minimum, save the value and position
-		if (potEnergy < epsVal){
-			epsVal = potEnergy;
-			sDist = probeCoor;
-		}
+    //if we've reached a new minimum, save the value and position
+    if (potEnergy < epsVal){
+      epsVal = potEnergy;
+      sDist = probeCoor;
+    }
 		
-		//test if the probe reached the origin - if so, stop stepping closer
-		if (probeCoor < 0){
-			sigDist = 0.0;
-			tooClose = 1;
-		}
+    //test if the probe reached the origin - if so, stop stepping closer
+    if (probeCoor < 0){
+      sigDist = 0.0;
+      tooClose = 1;
+    }
 		
-		//test if the probe beyond the contact point - if not, take a step closer
-		if (potEnergy < 0){
-			sigDist = probeCoor;
-			tempPotEnergy = potEnergy;
-			probeCoor -= rhoStep;
-		}
-		else {
-			//do a linear interpolation to obtain the sigDist 
-			interpRange = potEnergy - tempPotEnergy;
-			interpFrac = potEnergy / interpRange;
-			interpFrac = interpFrac * rhoStep;
-			sigDist = probeCoor + interpFrac;
+    //test if the probe beyond the contact point - if not, take a step closer
+    if (potEnergy < 0){
+      sigDist = probeCoor;
+      tempPotEnergy = potEnergy;
+      probeCoor -= rhoStep;
+    }
+    else {
+      //do a linear interpolation to obtain the sigDist 
+      interpRange = potEnergy - tempPotEnergy;
+      interpFrac = potEnergy / interpRange;
+      interpFrac = interpFrac * rhoStep;
+      sigDist = probeCoor + interpFrac;
 			
-			//end the loop
-			tooClose = 1;
-		}
-	}
+      //end the loop
+      tooClose = 1;
+    }
+  }
 }
 
 void GridBuilder::calcEnergy(){
-	
-}
+  double rXij, rYij, rZij;
+  double rijSquared;
+  double rValSquared, rValPowerSix;
+  double rparHe, epsHe;
+  double atomRpar, atomEps;
+  double rbAtomPos[3];
+  
+  //first get 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;
+  }
+  
+  potEnergy = 0.0;
+  
+  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::stepTheta(double increment){
-	//zero out the euler angles
-	for (i=0; i<3; i++)
-		angles[i] = 0.0;
+void GridBuilder::rotBody(double pValue, double tValue){
+  //zero out the euler angles
+  for (l=0; l<3; l++)
+    angles[i] = 0.0;
 	
-	//the second euler angle is for rotation about the x-axis (we use the zxz convention)
-	angles[1] = increment;
+  //the phi euler angle is for rotation about the z-axis (we use the zxz convention)
+  angles[0] = pValue;
+  //the second euler angle is for rotation about the x-axis (we use the zxz convention)
+  angles[1] = tValue;
 	
-	//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);
-	
+  //obtain the rotation matrix through the rigid body class
+  rbMol->doEulerToRotMat(angles, rotX);
+  
+  //start from the reference position
+  identityMat3(rbMatrix);
+  rbMol->setA(rbMatrix);
+  
+  //rotate the rigid body
+  matMul3(rotX, rbMatrix, rotatedMat);
+  rbMol->setA(rotatedMat);	
 }
 
-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";
+  }
+}
\ No newline at end of file