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

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;

  printf("startDist = %lf\n", startDist);

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

  for (k =0; k < bandwidth; k++) {
    thetaVal = thetaMin + k*thetaStep;
    for (j=0; j < bandwidth; j++) {
      phiVal = j*phiStep;

      printf("setting Euler, phi = %lf\ttheta = %lf\n", phiVal, thetaVal);

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

      releaseProbe(startDist);

      printf("found sigDist = %lf\t sDist = %lf \t epsVal = %lf\n", 
             sigDist, sDist, epsVal);

      sigList.push_back(sigDist);
      sList.push_back(sDist);
      epsList.push_back(epsVal);

    }
  }             
}

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;

  rbMol->getAtomPos(rbAtomPos, 0);

  printf("atom0 pos = %lf\t%lf\t%lf\n", rbAtomPos[0], rbAtomPos[1], rbAtomPos[2]);


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