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 1287 by chrisfen, Wed Jun 23 20:18:48 2004 UTC

# Line 1 | Line 1
1   #include "GridBuilder.hpp"
2 – #include "MatVec3.h"
2   #define PI 3.14159265359
3  
4  
5 < GridBuilder::GridBuilder(RigidBody* rb, int bandWidth) {
5 > GridBuilder::GridBuilder(RigidBody* rb, int gridWidth) {
6    rbMol = rb;
7 <  bandwidth = bandWidth;
8 <  thetaStep = PI / bandwidth;
7 >  gridwidth = gridWidth;
8 >  thetaStep = PI / gridwidth;
9    thetaMin = thetaStep / 2.0;
10    phiStep = thetaStep * 2.0;
12 –        
13 –  //zero out the rot mats
14 –  for (i=0; i<3; i++) {
15 –    for (j=0; j<3; j++) {
16 –      rotX[i][j] = 0.0;
17 –      rotZ[i][j] = 0.0;
18 –      rbMatrix[i][j] = 0.0;
19 –    }
20 –  }
11   }
12  
13   GridBuilder::~GridBuilder() {
14   }
15  
16 < void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid, vector<double> sGrid,
17 <                              vector<double> epsGrid){
16 > void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid,
17 >                              vector<double> sGrid, vector<double> epsGrid){
18 >  ofstream sigmaOut("sigma.grid");
19 >  ofstream sOut("s.grid");
20 >  ofstream epsOut("eps.grid");
21    double startDist;
22 +  double phiVal;
23 +  double thetaVal;
24 +  double sigTemp, sTemp, epsTemp, sigProbe;
25    double minDist = 10.0; //minimum start distance
26          
27 +  sigList = sigmaGrid;
28 +  sList = sGrid;
29 +  epsList = epsGrid;
30    forcefield = forceField;
31 +  
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      
56 <  //first determine the start distance - we always start at least minDist away
56 >  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    startDist = rbMol->findMaxExtent() + minDist;
63    if (startDist < minDist)
64      startDist = minDist;
65 <        
66 <  initBody();
67 <  for (i=0; i<bandwidth; i++){          
68 <    for (j=0; j<bandwidth; j++){
65 >
66 >  //set the initial orientation of the body and loop over theta values
67 >
68 >  for (k =0; k < gridwidth; k++) {
69 >    thetaVal = thetaMin + k*thetaStep;
70 >    printf("Theta step %i\n", k);
71 >    for (j=0; j < gridwidth; j++) {
72 >      phiVal = j*phiStep;
73 >
74 >      rbMol->setEuler(0.0, thetaVal, phiVal);
75 >
76        releaseProbe(startDist);
77  
78 <      sigmaGrid.push_back(sigDist);
79 <      sGrid.push_back(sDist);
80 <      epsGrid.push_back(epsVal);
81 <                        
82 <      stepPhi(phiStep);
78 >      //translate the values to sigma, s, and epsilon of the rigid body
79 >      sigTemp = 2*sigDist - sigProbe;
80 >      sTemp = (2*(sDist - sigDist))/(0.122462048309) - sigProbe;
81 >      epsTemp = pow(epsVal, 2)/fabs(epsHe);
82 >      
83 >      sigList.push_back(sigTemp);
84 >      sList.push_back(sTemp);
85 >      epsList.push_back(epsTemp);
86      }
49 –    stepTheta(thetaStep);
87    }            
88   }
89  
53 – void GridBuilder::initBody(){
54 –  //set up the rigid body in the starting configuration
55 –  stepTheta(thetaMin);
56 – }
57 –
90   void GridBuilder::releaseProbe(double farPos){
91    int tooClose;
92    double tempPotEnergy;
# Line 67 | Line 99 | void GridBuilder::releaseProbe(double farPos){
99    tooClose = 0;
100    epsVal = 0;
101    rhoStep = 0.1; //the distance the probe atom moves between steps
102 <        
71 <        
102 >                
103    while (!tooClose){
104      calcEnergy();
105      potProgress.push_back(potEnergy);
# Line 106 | Line 137 | void GridBuilder::calcEnergy(){
137   }
138  
139   void GridBuilder::calcEnergy(){
140 <        
141 < }
140 >  double rXij, rYij, rZij;
141 >  double rijSquared;
142 >  double rValSquared, rValPowerSix;
143 >  double atomRpar, atomEps;
144 >  double rbAtomPos[3];
145 >    
146 >  potEnergy = 0.0;
147  
148 < void GridBuilder::stepTheta(double increment){
149 <  //zero out the euler angles
150 <  for (i=0; i<3; i++)
151 <    angles[i] = 0.0;
152 <        
153 <  //the second euler angle is for rotation about the x-axis (we use the zxz convention)
154 <  angles[1] = increment;
155 <        
156 <  //obtain the rotation matrix through the rigid body class
157 <  rbMol->doEulerToRotMat(angles, rotX);
158 <        
159 <  //rotate the rigid body
160 <  rbMol->getA(rbMatrix);
161 <  matMul3(rotX, rbMatrix, rotatedMat);
162 <  rbMol->setA(rotatedMat);      
163 < }
148 >  for(i=0; i<rbMol->getNumAtoms(); i++){
149 >    rbMol->getAtomPos(rbAtomPos, i);
150 >    
151 >    rXij = rbAtomPos[0];
152 >    rYij = rbAtomPos[1];
153 >    rZij = rbAtomPos[2] - probeCoor;
154 >    
155 >    rijSquared = rXij * rXij + rYij * rYij + rZij * rZij;
156 >    
157 >    //in the interest of keeping the code more compact, we are being less
158 >    //efficient by placing a switch statement in the calculation loop
159 >    switch(forcefield){
160 >      case 1:{
161 >        //we are using the CHARMm force field
162 >        atomRpar = rbMol->getAtomRpar(i);
163 >        atomEps = rbMol->getAtomEps(i);
164 >        
165 >        rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared);
166 >        rValPowerSix = rValSquared * rValSquared * rValSquared;
167 >        potEnergy += sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0));
168 >      }; break;
169 >      
170 >      case 2:{
171 >        //we are using the AMBER 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 3:{
181 >        //we are using Allen-Tildesley LJ parameters
182 >        atomRpar = rbMol->getAtomRpar(i);
183 >        atomEps = rbMol->getAtomEps(i);
184 >        
185 >        rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (4*rijSquared);
186 >        rValPowerSix = rValSquared * rValSquared * rValSquared;
187 >        potEnergy += 4*sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 1.0));
188 >        
189 >      }; break;
190 >      
191 >      case 4:{
192 >        //we are using the OPLS force field
193 >        atomRpar = rbMol->getAtomRpar(i);
194 >        atomEps = rbMol->getAtomEps(i);
195 >        
196 >        rValSquared = (pow(sqrt(rparHe+atomRpar),2)) / (rijSquared);
197 >        rValPowerSix = rValSquared * rValSquared * rValSquared;
198 >        potEnergy += 4*sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 1.0));
199 >      }; break;
200 >      
201 >      case 5:{
202 >        //we are using the GAFF force field
203 >        atomRpar = rbMol->getAtomRpar(i);
204 >        atomEps = rbMol->getAtomEps(i);
205 >        
206 >        rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared);
207 >        rValPowerSix = rValSquared * rValSquared * rValSquared;
208 >        potEnergy += sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0));
209 >      }; break;
210 >    }    
211 >  }
212 > }
213  
214 < void GridBuilder::stepPhi(double increment){
215 <  //zero out the euler angles
216 <  for (i=0; i<3; i++)
217 <    angles[i] = 0.0;
218 <        
219 <  //the phi euler angle is for rotation about the z-axis (we use the zxz convention)
220 <  angles[0] = increment;
221 <        
222 <  //obtain the rotation matrix through the rigid body class
223 <  rbMol->doEulerToRotMat(angles, rotZ);
139 <        
140 <  //rotate the rigid body
141 <  rbMol->getA(rbMatrix);
142 <  matMul3(rotZ, rbMatrix, rotatedMat);
143 <  rbMol->setA(rotatedMat);      
214 > void GridBuilder::printGridFiles(){
215 >  ofstream sigmaOut("sigma.grid");
216 >  ofstream sOut("s.grid");
217 >  ofstream epsOut("eps.grid");
218 >  
219 >  for (k=0; k<sigList.size(); k++){
220 >    sigmaOut << sigList[k] << "\n0\n";
221 >    sOut << sList[k] << "\n0\n";    
222 >    epsOut << epsList[k] << "\n0\n";
223 >  }
224   }
225 +

Diff Legend

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