ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/SHAPES/GridBuilder.cpp
(Generate patch)

Comparing trunk/SHAPES/GridBuilder.cpp (file contents):
Revision 1280 by chrisfen, Fri Jun 18 19:40:31 2004 UTC vs.
Revision 1281 by chrisfen, Mon Jun 21 13:38:55 2004 UTC

# Line 25 | Line 25 | void GridBuilder::launchProbe(int forceField, vector<d
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 minDist = 10.0; //minimum start distance
33          
34 +  sList = sGrid;
35 +  sigList = sigmaGrid;
36 +  epsList = epsGrid;
37    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 <        
43 >
44    initBody();
45 <  for (i=0; i<bandwidth; i++){          
45 >  for (k=0; k<bandwidth; k++){          
46 >    printf("step theta...\n");
47      for (j=0; j<bandwidth; j++){
48        releaseProbe(startDist);
49  
50 <      sigmaGrid.push_back(sigDist);
51 <      sGrid.push_back(sDist);
52 <      epsGrid.push_back(epsVal);
50 >      sigList.push_back(sigDist);
51 >      sList.push_back(sDist);
52 >      epsList.push_back(epsVal);
53                          
54        stepPhi(phiStep);
55      }
56      stepTheta(thetaStep);
57    }            
58 +  /*
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   }
68  
69   void GridBuilder::initBody(){
# Line 106 | Line 122 | void GridBuilder::calcEnergy(){
122   }
123  
124   void GridBuilder::calcEnergy(){
125 <        
126 < }
125 >  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(epsHe*atomEps)*(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(epsHe*atomEps)*(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)) / (4*rijSquared);
196 >        rValPowerSix = rValSquared * rValSquared * rValSquared;
197 >        potEnergy += 4*sqrt(epsHe*atomEps)*(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 += 4*sqrt(epsHe*atomEps)*(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(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0));
220 >      }; break;
221 >    }    
222 >  }
223 > }
224  
225   void GridBuilder::stepTheta(double increment){
226    //zero out the euler angles
227 <  for (i=0; i<3; i++)
227 >  for (l=0; l<3; l++)
228      angles[i] = 0.0;
229          
230    //the second euler angle is for rotation about the x-axis (we use the zxz convention)
# Line 128 | Line 241 | void GridBuilder::stepPhi(double increment){
241  
242   void GridBuilder::stepPhi(double increment){
243    //zero out the euler angles
244 <  for (i=0; i<3; i++)
244 >  for (l=0; l<3; l++)
245      angles[i] = 0.0;
246          
247    //the phi euler angle is for rotation about the z-axis (we use the zxz convention)
# Line 142 | Line 255 | void GridBuilder::stepPhi(double increment){
255    matMul3(rotZ, rbMatrix, rotatedMat);
256    rbMol->setA(rotatedMat);      
257   }
258 +
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 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines