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 1277 by chrisfen, Thu Jun 17 21:29:17 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) {
6 <        rbMol = rb;
7 <        bandwidth = bandWidth;
8 <        thetaStep = PI / bandwidth;
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 <        }
5 > GridBuilder::GridBuilder(RigidBody* rb, int gridWidth) {
6 >  rbMol = rb;
7 >  gridwidth = gridWidth;
8 >  thetaStep = PI / gridwidth;
9 >  thetaMin = thetaStep / 2.0;
10 >  phiStep = thetaStep * 2.0;
11   }
12  
13   GridBuilder::~GridBuilder() {
14   }
15  
16 < void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid, vector<double> sGrid,
17 <        vector<double> epsGrid){
18 <        double startDist;
19 <        double minDist = 10.0; //minimum start distance
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 <        //first determine the start distance - we always start at least minDist away
28 <        startDist = rbMol->findMaxExtent() + minDist;
29 <        if (startDist < minDist)
30 <                startDist = minDist;
31 <        
32 <        initBody();
33 <        for (i=0; i<bandwidth; i++){            
34 <                for (j=0; j<bandwidth; j++){
35 <                        releaseProbe(startDist);
36 <                        stepPhi(phiStep);
37 <                }
38 <                stepTheta(thetaStep);
39 <        }
40 <                
41 < }
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 >  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 < void GridBuilder::initBody(){
67 <        //set up the rigid body in the starting configuration
68 <        stepTheta(thetaMin);
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 >      //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 >    }
89 >  }            
90   }
91  
92   void GridBuilder::releaseProbe(double farPos){
93 <        int tooClose;
94 <        double tempPotEnergy;
95 <        double interpRange;
96 <        double interpFrac;
93 >  int tooClose;
94 >  double tempPotEnergy;
95 >  double interpRange;
96 >  double interpFrac;
97          
98 <        probeCoor = farPos;
99 <        tooClose = 0;
100 <        epsVal = 0;
101 <        rhoStep = 0.1; //the distance the probe atom moves between steps
102 <        
103 <        while (!tooClose){
64 <                calcEnergy();
65 <                potProgress.push_back(potEnergy);
66 <                distProgress.push_back(probeCoor);
98 >  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                  
105 <                //if we've reached a new minimum, save the value and position
106 <                if (potEnergy < epsVal){
107 <                        epsVal = potEnergy;
108 <                        sDist = probeCoor;
72 <                }
105 >  while (!tooClose){
106 >    calcEnergy();
107 >    potProgress.push_back(potEnergy);
108 >    distProgress.push_back(probeCoor);
109                  
110 <                //test if the probe reached the origin - if so, stop stepping closer
111 <                if (probeCoor < 0){
112 <                        sigDist = 0.0;
113 <                        tooClose = 1;
114 <                }
110 >    //if we've reached a new minimum, save the value and position
111 >    if (potEnergy < epsVal){
112 >      epsVal = potEnergy;
113 >      sDist = probeCoor;
114 >    }
115                  
116 <                //test if the probe beyond the contact point - if not, take a step closer
117 <                if (potEnergy < 0){
118 <                        sigDist = probeCoor;
119 <                        tempPotEnergy = potEnergy;
120 <                        probeCoor -= rhoStep;
121 <                }
122 <                else {
123 <                        //do a linear interpolation to obtain the sigDist
124 <                        interpRange = potEnergy - tempPotEnergy;
125 <                        interpFrac = potEnergy / interpRange;
126 <                        interpFrac = interpFrac * rhoStep;
127 <                        sigDist = probeCoor + interpFrac;
116 >    //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 >                
122 >    //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                          
135 <                        //end the loop
136 <                        tooClose = 1;
137 <                }
138 <        }
135 >      //end the loop
136 >      tooClose = 1;
137 >    }
138 >  }
139   }
140  
141   void GridBuilder::calcEnergy(){
142 <        
143 < }
142 >  double rXij, rYij, rZij;
143 >  double rijSquared;
144 >  double rValSquared, rValPowerSix;
145 >  double atomRpar, atomEps;
146 >  double rbAtomPos[3];
147 >    
148 >  potEnergy = 0.0;
149  
150 < void GridBuilder::stepTheta(double increment){
151 <        //zero out the euler angles
152 <        for (i=0; i<3; i++)
153 <                angles[i] = 0.0;
154 <        
155 <        //the second euler angle is for rotation about the x-axis (we use the zxz convention)
156 <        angles[1] = increment;
157 <        
158 <        //obtain the rotation matrix through the rigid body class
159 <        rbMol->doEulerToRotMat(angles, rotX);
160 <        
161 <        //rotate the rigid body
162 <        rbMol->getA(rbMatrix);
163 <        matMul3(rotX, rbMatrix, rotatedMat);
164 <        rbMol->setA(rotatedMat);
165 <        
166 < }
150 >  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 >    //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
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(epsHe*atomEps)*(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(epsHe*atomEps)*(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)) / (4*rijSquared);
188 >        rValPowerSix = rValSquared * rValSquared * rValSquared;
189 >        potEnergy += 4*sqrt(epsHe*atomEps)*(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 += 4*sqrt(epsHe*atomEps)*(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(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0));
211 >      }; break;
212 >    }    
213 >  }
214 > }
215  
216 < void GridBuilder::stepPhi(double increment){
217 <        //zero out the euler angles
218 <        for (i=0; i<3; i++)
219 <                angles[i] = 0.0;
220 <        
221 <        //the phi euler angle is for rotation about the z-axis (we use the zxz convention)
222 <        angles[0] = increment;
223 <        
224 <        //obtain the rotation matrix through the rigid body class
225 <        rbMol->doEulerToRotMat(angles, rotZ);
131 <        
132 <        //rotate the rigid body
133 <        rbMol->getA(rbMatrix);
134 <        matMul3(rotZ, rbMatrix, rotatedMat);
135 <        rbMol->setA(rotatedMat);
136 <        
216 > 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   }
227 +

Diff Legend

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