trunk/SHAPES/GridBuilder.cpp

#include "GridBuilder.hpp"
#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
        
  sigList = sigmaGrid;
  sList = sGrid;
  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;
    for (j=0; j < gridwidth; j++) {
      //s2kit10 is actually taking a grid from -pi/2 to 3pi/3, not 0 to 2pi...
      phiVal = j*phiStep - 0.5*PI;
       if (phiVal<0.0)
         phiVal += 2*PI;

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