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 phiVal;
  double thetaVal;
  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;

  //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);

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