trunk/SHAPES/GridBuilder.cpp

#include "GridBuilder.hpp"
#include "MatVec3.h"
#define PI 3.14159265359


GridBuilder::GridBuilder(RigidBody* rb, int bandWidth) {
  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;
    }
  }
}

GridBuilder::~GridBuilder() {
}

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 minDist = 10.0; //minimum start distance
        
  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;

  initBody();
  for (k=0; k<bandwidth; k++){          
    printf("step theta...\n");
    for (j=0; j<bandwidth; j++){
      releaseProbe(startDist);

      sigList.push_back(sigDist);
      sList.push_back(sDist);
      epsList.push_back(epsVal);
                        
      stepPhi(phiStep);
    }
    stepTheta(thetaStep);
  }             
  /*
  //write out the grid files
  printf("the grid size is %d\n",sigmaGrid.size());
  for (k=0; k<sigmaGrid.size(); k++){
    sigmaOut << sigmaGrid[k] << "\n0\n";
    sOut << sGrid[k] << "\n0\n";
    epsOut << epsGrid[k] << "\n0\n";
  }
   */
}

void GridBuilder::initBody(){
  //set up the rigid body in the starting configuration
  stepTheta(thetaMin);
}

void GridBuilder::releaseProbe(double farPos){
  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
        
        
  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;
    }
                
    //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;
                        
      //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 (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;
        
  //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);      
}

void GridBuilder::stepPhi(double increment){
  //zero out the euler angles
  for (l=0; l<3; l++)
    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";
  }
}
Revision:	1281
Committed:	Mon Jun 21 13:38:55 2004 UTC (20 years, 2 months ago) by chrisfen
File size:	7193 byte(s)
Log Message:	Now calculates energies and builds the grid files. There are still bugs, but it compiles and runs...
#	User	Rev	Content
1	chrisfen	1277	#include "GridBuilder.hpp"
2			#include "MatVec3.h"
3			#define PI 3.14159265359
4
5
6			GridBuilder::GridBuilder(RigidBody* rb, int bandWidth) {
7	chrisfen	1280	rbMol = rb;
8			bandwidth = bandWidth;
9			thetaStep = PI / bandwidth;
10			thetaMin = thetaStep / 2.0;
11			phiStep = thetaStep * 2.0;
12	chrisfen	1277
13	chrisfen	1280	//zero out the rot mats
14			for (i=0; i<3; i++) {
15			for (j=0; j<3; j++) {
16			rotX[i][j] = 0.0;
17			rotZ[i][j] = 0.0;
18			rbMatrix[i][j] = 0.0;
19			}
20			}
21	chrisfen	1277	}
22
23			GridBuilder::~GridBuilder() {
24			}
25
26			void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid, vector<double> sGrid,
27	chrisfen	1280	vector<double> epsGrid){
28	chrisfen	1281	ofstream sigmaOut("sigma.grid");
29			ofstream sOut("s.grid");
30			ofstream epsOut("eps.grid");
31	chrisfen	1280	double startDist;
32			double minDist = 10.0; //minimum start distance
33	chrisfen	1277
34	chrisfen	1281	sList = sGrid;
35			sigList = sigmaGrid;
36			epsList = epsGrid;
37	chrisfen	1280	forcefield = forceField;
38
39			//first determine the start distance - we always start at least minDist away
40			startDist = rbMol->findMaxExtent() + minDist;
41			if (startDist < minDist)
42			startDist = minDist;
43	chrisfen	1281
44	chrisfen	1280	initBody();
45	chrisfen	1281	for (k=0; k<bandwidth; k++){
46			printf("step theta...\n");
47	chrisfen	1280	for (j=0; j<bandwidth; j++){
48			releaseProbe(startDist);
49	chrisfen	1279
50	chrisfen	1281	sigList.push_back(sigDist);
51			sList.push_back(sDist);
52			epsList.push_back(epsVal);
53	chrisfen	1278
54	chrisfen	1280	stepPhi(phiStep);
55			}
56			stepTheta(thetaStep);
57			}
58	chrisfen	1281	/*
59			//write out the grid files
60			printf("the grid size is %d\n",sigmaGrid.size());
61			for (k=0; k<sigmaGrid.size(); k++){
62			sigmaOut << sigmaGrid[k] << "\n0\n";
63			sOut << sGrid[k] << "\n0\n";
64			epsOut << epsGrid[k] << "\n0\n";
65			}
66			*/
67	chrisfen	1277	}
68
69			void GridBuilder::initBody(){
70	chrisfen	1280	//set up the rigid body in the starting configuration
71			stepTheta(thetaMin);
72	chrisfen	1277	}
73
74			void GridBuilder::releaseProbe(double farPos){
75	chrisfen	1280	int tooClose;
76			double tempPotEnergy;
77			double interpRange;
78			double interpFrac;
79	chrisfen	1277
80	chrisfen	1280	probeCoor = farPos;
81			potProgress.clear();
82			distProgress.clear();
83			tooClose = 0;
84			epsVal = 0;
85			rhoStep = 0.1; //the distance the probe atom moves between steps
86	chrisfen	1277
87	chrisfen	1279
88	chrisfen	1280	while (!tooClose){
89			calcEnergy();
90			potProgress.push_back(potEnergy);
91			distProgress.push_back(probeCoor);
92	chrisfen	1277
93	chrisfen	1280	//if we've reached a new minimum, save the value and position
94			if (potEnergy < epsVal){
95			epsVal = potEnergy;
96			sDist = probeCoor;
97			}
98	chrisfen	1277
99	chrisfen	1280	//test if the probe reached the origin - if so, stop stepping closer
100			if (probeCoor < 0){
101			sigDist = 0.0;
102			tooClose = 1;
103			}
104	chrisfen	1277
105	chrisfen	1280	//test if the probe beyond the contact point - if not, take a step closer
106			if (potEnergy < 0){
107			sigDist = probeCoor;
108			tempPotEnergy = potEnergy;
109			probeCoor -= rhoStep;
110			}
111			else {
112			//do a linear interpolation to obtain the sigDist
113			interpRange = potEnergy - tempPotEnergy;
114			interpFrac = potEnergy / interpRange;
115			interpFrac = interpFrac * rhoStep;
116			sigDist = probeCoor + interpFrac;
117	chrisfen	1277
118	chrisfen	1280	//end the loop
119			tooClose = 1;
120			}
121			}
122	chrisfen	1277	}
123
124			void GridBuilder::calcEnergy(){
125	chrisfen	1281	double rXij, rYij, rZij;
126			double rijSquared;
127			double rValSquared, rValPowerSix;
128			double rparHe, epsHe;
129			double atomRpar, atomEps;
130			double rbAtomPos[3];
131
132			//first get the probe atom parameters
133			switch(forcefield){
134			case 1:{
135			rparHe = 1.4800;
136			epsHe = -0.021270;
137			}; break;
138			case 2:{
139			rparHe = 1.14;
140			epsHe = 0.0203;
141			}; break;
142			case 3:{
143			rparHe = 2.28;
144			epsHe = 0.020269601874;
145			}; break;
146			case 4:{
147			rparHe = 2.5560;
148			epsHe = 0.0200;
149			}; break;
150			case 5:{
151			rparHe = 1.14;
152			epsHe = 0.0203;
153			}; break;
154			}
155
156			potEnergy = 0.0;
157
158			for(i=0; i<rbMol->getNumAtoms(); i++){
159			rbMol->getAtomPos(rbAtomPos, i);
160
161			rXij = rbAtomPos[0];
162			rYij = rbAtomPos[1];
163			rZij = rbAtomPos[2] - probeCoor;
164
165			rijSquared = rXij * rXij + rYij * rYij + rZij * rZij;
166
167			//in the interest of keeping the code more compact, we are being less efficient by placing
168			//a switch statement in the calculation loop
169			switch(forcefield){
170			case 1:{
171			//we are using the CHARMm force field
172			atomRpar = rbMol->getAtomRpar(i);
173			atomEps = rbMol->getAtomEps(i);
174
175			rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared);
176			rValPowerSix = rValSquared * rValSquared * rValSquared;
177			potEnergy += sqrt(epsHeatomEps)(rValPowerSix * (rValPowerSix - 2.0));
178			}; break;
179
180			case 2:{
181			//we are using the AMBER force field
182			atomRpar = rbMol->getAtomRpar(i);
183			atomEps = rbMol->getAtomEps(i);
184
185			rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared);
186			rValPowerSix = rValSquared * rValSquared * rValSquared;
187			potEnergy += sqrt(epsHeatomEps)(rValPowerSix * (rValPowerSix - 2.0));
188			}; break;
189
190			case 3:{
191			//we are using Allen-Tildesley LJ parameters
192			atomRpar = rbMol->getAtomRpar(i);
193			atomEps = rbMol->getAtomEps(i);
194
195			rValSquared = ((rparHe+atomRpar)(rparHe+atomRpar)) / (4rijSquared);
196			rValPowerSix = rValSquared * rValSquared * rValSquared;
197			potEnergy += 4sqrt(epsHeatomEps)(rValPowerSix (rValPowerSix - 1.0));
198
199			}; break;
200
201
202			case 4:{
203			//we are using the OPLS force field
204			atomRpar = rbMol->getAtomRpar(i);
205			atomEps = rbMol->getAtomEps(i);
206
207			rValSquared = (pow(sqrt(rparHe+atomRpar),2)) / (rijSquared);
208			rValPowerSix = rValSquared * rValSquared * rValSquared;
209			potEnergy += 4sqrt(epsHeatomEps)(rValPowerSix (rValPowerSix - 1.0));
210			}; break;
211
212			case 5:{
213			//we are using the GAFF force field
214			atomRpar = rbMol->getAtomRpar(i);
215			atomEps = rbMol->getAtomEps(i);
216
217			rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared);
218			rValPowerSix = rValSquared * rValSquared * rValSquared;
219			potEnergy += sqrt(epsHeatomEps)(rValPowerSix * (rValPowerSix - 2.0));
220			}; break;
221			}
222			}
223			}
224	chrisfen	1277
225			void GridBuilder::stepTheta(double increment){
226	chrisfen	1280	//zero out the euler angles
227	chrisfen	1281	for (l=0; l<3; l++)
228	chrisfen	1280	angles[i] = 0.0;
229	chrisfen	1277
230	chrisfen	1280	//the second euler angle is for rotation about the x-axis (we use the zxz convention)
231			angles[1] = increment;
232	chrisfen	1277
233	chrisfen	1280	//obtain the rotation matrix through the rigid body class
234			rbMol->doEulerToRotMat(angles, rotX);
235	chrisfen	1277
236	chrisfen	1280	//rotate the rigid body
237			rbMol->getA(rbMatrix);
238			matMul3(rotX, rbMatrix, rotatedMat);
239			rbMol->setA(rotatedMat);
240	chrisfen	1277	}
241
242			void GridBuilder::stepPhi(double increment){
243	chrisfen	1280	//zero out the euler angles
244	chrisfen	1281	for (l=0; l<3; l++)
245	chrisfen	1280	angles[i] = 0.0;
246	chrisfen	1277
247	chrisfen	1280	//the phi euler angle is for rotation about the z-axis (we use the zxz convention)
248			angles[0] = increment;
249	chrisfen	1277
250	chrisfen	1280	//obtain the rotation matrix through the rigid body class
251			rbMol->doEulerToRotMat(angles, rotZ);
252	chrisfen	1277
253	chrisfen	1280	//rotate the rigid body
254			rbMol->getA(rbMatrix);
255			matMul3(rotZ, rbMatrix, rotatedMat);
256			rbMol->setA(rotatedMat);
257	chrisfen	1277	}
258	chrisfen	1281
259			void GridBuilder::printGridFiles(){
260			ofstream sigmaOut("sigma.grid");
261			ofstream sOut("s.grid");
262			ofstream epsOut("eps.grid");
263
264			for (k=0; k<sigList.size(); k++){
265			sigmaOut << sigList[k] << "\n0\n";
266			sOut << sList[k] << "\n0\n";
267			epsOut << epsList[k] << "\n0\n";
268			}
269			}