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 1282 by chrisfen, Mon Jun 21 15:10:29 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 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 <  initBody();
47 <  for (i=0; i<bandwidth; i++){          
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 <      sigmaGrid.push_back(sigDist);
56 <      sGrid.push_back(sDist);
57 <      epsGrid.push_back(epsVal);
58 <                        
59 <      stepPhi(phiStep);
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 <    stepTheta(thetaStep);
62 >    phiVal = 0.0;
63 >    thetaVal += thetaStep;
64 >    rotBody(phiVal, thetaVal);
65 >    printf("step theta %i\n",k);
66    }            
67   }
68  
53 – void GridBuilder::initBody(){
54 –  //set up the rigid body in the starting configuration
55 –  stepTheta(thetaMin);
56 – }
57 –
69   void GridBuilder::releaseProbe(double farPos){
70    int tooClose;
71    double tempPotEnergy;
# Line 106 | Line 117 | void GridBuilder::calcEnergy(){
117   }
118  
119   void GridBuilder::calcEnergy(){
120 <        
121 < }
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(epsHe*atomEps)*(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(epsHe*atomEps)*(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)) / (4*rijSquared);
191 >        rValPowerSix = rValSquared * rValSquared * rValSquared;
192 >        potEnergy += 4*sqrt(epsHe*atomEps)*(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 += 4*sqrt(epsHe*atomEps)*(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(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0));
214 >      }; break;
215 >    }    
216 >  }
217 > }
218  
219 < void GridBuilder::stepTheta(double increment){
219 > void GridBuilder::rotBody(double pValue, double tValue){
220    //zero out the euler angles
221 <  for (i=0; i<3; i++)
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] = increment;
227 >  angles[1] = tValue;
228          
229    //obtain the rotation matrix through the rigid body class
230    rbMol->doEulerToRotMat(angles, rotX);
231 <        
231 >  
232 >  //start from the reference position
233 >  identityMat3(rbMatrix);
234 >  rbMol->setA(rbMatrix);
235 >  
236    //rotate the rigid body
124 –  rbMol->getA(rbMatrix);
237    matMul3(rotX, rbMatrix, rotatedMat);
238    rbMol->setA(rotatedMat);      
239   }
240  
241 < void GridBuilder::stepPhi(double increment){
242 <  //zero out the euler angles
243 <  for (i=0; i<3; i++)
244 <    angles[i] = 0.0;
245 <        
246 <  //the phi euler angle is for rotation about the z-axis (we use the zxz convention)
247 <  angles[0] = increment;
248 <        
249 <  //obtain the rotation matrix through the rigid body class
250 <  rbMol->doEulerToRotMat(angles, rotZ);
251 <        
140 <  //rotate the rigid body
141 <  rbMol->getA(rbMatrix);
142 <  matMul3(rotZ, rbMatrix, rotatedMat);
143 <  rbMol->setA(rotatedMat);      
144 < }
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 > }

Diff Legend

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