trunk/SHAPES/GridBuilder.cpp

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


GridBuilder::GridBuilder(RigidBody* rb, int gridWidth) {
  rbMol = rb;
  gridwidth = gridWidth;
  thetaStep = PI / gridwidth;
  thetaMin = thetaStep / 2.0;
  phiStep = thetaStep * 2.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 sigTemp, sTemp, epsTemp, sigProbe;
  double minDist = 10.0; //minimum start distance
        
  sList = sGrid;
  sigList = sigmaGrid;
  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;
  }
    
  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;

  //set the initial orientation of the body and loop over theta values

  for (k =0; k < gridwidth; k++) {
    thetaVal = thetaMin + k*thetaStep;
    printf("Theta step %i\n", k);
    for (j=0; j < gridwidth; j++) {
      phiVal = j*phiStep;

      rbMol->setEuler(0.0, thetaVal, phiVal);

      releaseProbe(startDist);

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

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