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 1281 by chrisfen, Mon Jun 21 13:38:55 2004 UTC vs.
Revision 1308 by chrisfen, Fri Jun 25 18:58:12 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          
34 –  sList = sGrid;
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 (k=0; k<bandwidth; k++){          
68 <    printf("step theta...\n");
69 <    for (j=0; j<bandwidth; j++){
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 >    for (j=0; j < gridwidth; j++) {
71 >      //s2kit10 is actually taking a grid from -pi/2 to 3pi/3, not 0 to 2pi...
72 >      phiVal = j*phiStep - 0.5*PI;
73 >       if (phiVal<0.0)
74 >         phiVal += 2*PI;
75 >
76 >      rbMol->setEuler(0.0, thetaVal, phiVal);
77 >
78        releaseProbe(startDist);
79  
80 <      sigList.push_back(sigDist);
81 <      sList.push_back(sDist);
82 <      epsList.push_back(epsVal);
83 <                        
84 <      stepPhi(phiStep);
80 >      //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      }
56 –    stepTheta(thetaStep);
89    }            
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 –   */
90   }
91  
69 – void GridBuilder::initBody(){
70 –  //set up the rigid body in the starting configuration
71 –  stepTheta(thetaMin);
72 – }
73 –
92   void GridBuilder::releaseProbe(double farPos){
93    int tooClose;
94    double tempPotEnergy;
# Line 83 | Line 101 | void GridBuilder::releaseProbe(double farPos){
101    tooClose = 0;
102    epsVal = 0;
103    rhoStep = 0.1; //the distance the probe atom moves between steps
104 <        
87 <        
104 >                
105    while (!tooClose){
106      calcEnergy();
107      potProgress.push_back(potEnergy);
# Line 125 | Line 142 | void GridBuilder::calcEnergy(){
142    double rXij, rYij, rZij;
143    double rijSquared;
144    double rValSquared, rValPowerSix;
128 –  double rparHe, epsHe;
145    double atomRpar, atomEps;
146    double rbAtomPos[3];
147 <  
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 <  
147 >    
148    potEnergy = 0.0;
149 <  
149 >
150    for(i=0; i<rbMol->getNumAtoms(); i++){
151      rbMol->getAtomPos(rbAtomPos, i);
152      
# Line 164 | Line 156 | void GridBuilder::calcEnergy(){
156      
157      rijSquared = rXij * rXij + rYij * rYij + rZij * rZij;
158      
159 <    //in the interest of keeping the code more compact, we are being less efficient by placing
160 <    //a switch statement in the calculation loop
159 >    //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      switch(forcefield){
162        case 1:{
163          //we are using the CHARMm force field
# Line 197 | Line 189 | void GridBuilder::calcEnergy(){
189          potEnergy += 4*sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 1.0));
190          
191        }; break;
200 –        
192        
193        case 4:{
194          //we are using the OPLS force field
# Line 222 | Line 213 | void GridBuilder::calcEnergy(){
213    }
214   }
215  
225 – void GridBuilder::stepTheta(double increment){
226 –  //zero out the euler angles
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)
231 –  angles[1] = increment;
232 –        
233 –  //obtain the rotation matrix through the rigid body class
234 –  rbMol->doEulerToRotMat(angles, rotX);
235 –        
236 –  //rotate the rigid body
237 –  rbMol->getA(rbMatrix);
238 –  matMul3(rotX, rbMatrix, rotatedMat);
239 –  rbMol->setA(rotatedMat);      
240 – }
241 –
242 – void GridBuilder::stepPhi(double increment){
243 –  //zero out the euler angles
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)
248 –  angles[0] = increment;
249 –        
250 –  //obtain the rotation matrix through the rigid body class
251 –  rbMol->doEulerToRotMat(angles, rotZ);
252 –        
253 –  //rotate the rigid body
254 –  rbMol->getA(rbMatrix);
255 –  matMul3(rotZ, rbMatrix, rotatedMat);
256 –  rbMol->setA(rotatedMat);      
257 – }
258 –
216   void GridBuilder::printGridFiles(){
217    ofstream sigmaOut("sigma.grid");
218    ofstream sOut("s.grid");
# Line 266 | Line 223 | void GridBuilder::printGridFiles(){
223      sOut << sList[k] << "\n0\n";    
224      epsOut << epsList[k] << "\n0\n";
225    }
226 < }
226 > }
227 >

Diff Legend

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