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 1283 by gezelter, Mon Jun 21 15:54:27 2004 UTC

# Line 4 | Line 4 | GridBuilder::GridBuilder(RigidBody* rb, int bandWidth)
4  
5  
6   GridBuilder::GridBuilder(RigidBody* rb, int bandWidth) {
7 <        rbMol = rb;
8 <        bandwidth = bandWidth;
9 <        thetaStep = PI / bandwidth;
10 <        thetaMin = thetaStep / 2.0;
11 <        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 <        }
7 >  rbMol = rb;
8 >  bandwidth = bandWidth;
9 >  thetaStep = PI / bandwidth;
10 >  thetaMin = thetaStep / 2.0;
11 >  phiStep = thetaStep * 2.0;
12   }
13  
14   GridBuilder::~GridBuilder() {
15   }
16  
17 < void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid, vector<double> sGrid,
18 <        vector<double> epsGrid){
19 <        double startDist;
20 <        double minDist = 10.0; //minimum start distance
17 > void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid,
18 >                              vector<double> sGrid, vector<double> epsGrid){
19 >  ofstream sigmaOut("sigma.grid");
20 >  ofstream sOut("s.grid");
21 >  ofstream epsOut("eps.grid");
22 >  double startDist;
23 >  double phiVal;
24 >  double thetaVal;
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);
40 <                        stepPhi(phiStep);
41 <                }
42 <                stepTheta(thetaStep);
43 <        }
44 <                
45 < }
27 >  sList = sGrid;
28 >  sigList = sigmaGrid;
29 >  epsList = epsGrid;
30 >  forcefield = forceField;
31 >    
32 >  //first determine the start distance - we always start at least minDist away
33 >  startDist = rbMol->findMaxExtent() + minDist;
34 >  if (startDist < minDist)
35 >    startDist = minDist;
36  
37 < void GridBuilder::initBody(){
38 <        //set up the rigid body in the starting configuration
39 <        stepTheta(thetaMin);
37 >  printf("startDist = %lf\n", startDist);
38 >
39 >  //set the initial orientation of the body and loop over theta values
40 >
41 >  for (k =0; k < bandwidth; k++) {
42 >    thetaVal = thetaMin + k*thetaStep;
43 >    for (j=0; j < bandwidth; j++) {
44 >      phiVal = j*phiStep;
45 >
46 >      printf("setting Euler, phi = %lf\ttheta = %lf\n", phiVal, thetaVal);
47 >
48 >      rbMol->setEuler(0.0, thetaVal, phiVal);
49 >
50 >      releaseProbe(startDist);
51 >
52 >      printf("found sigDist = %lf\t sDist = %lf \t epsVal = %lf\n",
53 >             sigDist, sDist, epsVal);
54 >
55 >      sigList.push_back(sigDist);
56 >      sList.push_back(sDist);
57 >      epsList.push_back(epsVal);
58 >
59 >    }
60 >  }            
61   }
62  
63   void GridBuilder::releaseProbe(double farPos){
64 <        int tooClose;
65 <        double tempPotEnergy;
66 <        double interpRange;
67 <        double interpFrac;
64 >  int tooClose;
65 >  double tempPotEnergy;
66 >  double interpRange;
67 >  double interpFrac;
68          
69 <        probeCoor = farPos;
70 <        tooClose = 0;
71 <        epsVal = 0;
72 <        rhoStep = 0.1; //the distance the probe atom moves between steps
73 <        
74 <        while (!tooClose){
64 <                calcEnergy();
65 <                potProgress.push_back(potEnergy);
66 <                distProgress.push_back(probeCoor);
69 >  probeCoor = farPos;
70 >  potProgress.clear();
71 >  distProgress.clear();
72 >  tooClose = 0;
73 >  epsVal = 0;
74 >  rhoStep = 0.1; //the distance the probe atom moves between steps
75                  
76 <                //if we've reached a new minimum, save the value and position
77 <                if (potEnergy < epsVal){
78 <                        epsVal = potEnergy;
79 <                        sDist = probeCoor;
72 <                }
76 >  while (!tooClose){
77 >    calcEnergy();
78 >    potProgress.push_back(potEnergy);
79 >    distProgress.push_back(probeCoor);
80                  
81 <                //test if the probe reached the origin - if so, stop stepping closer
82 <                if (probeCoor < 0){
83 <                        sigDist = 0.0;
84 <                        tooClose = 1;
85 <                }
81 >    //if we've reached a new minimum, save the value and position
82 >    if (potEnergy < epsVal){
83 >      epsVal = potEnergy;
84 >      sDist = probeCoor;
85 >    }
86                  
87 <                //test if the probe beyond the contact point - if not, take a step closer
88 <                if (potEnergy < 0){
89 <                        sigDist = probeCoor;
90 <                        tempPotEnergy = potEnergy;
91 <                        probeCoor -= rhoStep;
92 <                }
93 <                else {
94 <                        //do a linear interpolation to obtain the sigDist
95 <                        interpRange = potEnergy - tempPotEnergy;
96 <                        interpFrac = potEnergy / interpRange;
97 <                        interpFrac = interpFrac * rhoStep;
98 <                        sigDist = probeCoor + interpFrac;
87 >    //test if the probe reached the origin - if so, stop stepping closer
88 >    if (probeCoor < 0){
89 >      sigDist = 0.0;
90 >      tooClose = 1;
91 >    }
92 >                
93 >    //test if the probe beyond the contact point - if not, take a step closer
94 >    if (potEnergy < 0){
95 >      sigDist = probeCoor;
96 >      tempPotEnergy = potEnergy;
97 >      probeCoor -= rhoStep;
98 >    }
99 >    else {
100 >      //do a linear interpolation to obtain the sigDist
101 >      interpRange = potEnergy - tempPotEnergy;
102 >      interpFrac = potEnergy / interpRange;
103 >      interpFrac = interpFrac * rhoStep;
104 >      sigDist = probeCoor + interpFrac;
105                          
106 <                        //end the loop
107 <                        tooClose = 1;
108 <                }
109 <        }
106 >      //end the loop
107 >      tooClose = 1;
108 >    }
109 >  }
110   }
111  
112   void GridBuilder::calcEnergy(){
113 <        
114 < }
113 >  double rXij, rYij, rZij;
114 >  double rijSquared;
115 >  double rValSquared, rValPowerSix;
116 >  double rparHe, epsHe;
117 >  double atomRpar, atomEps;
118 >  double rbAtomPos[3];
119 >  
120 >  //first get the probe atom parameters
121 >  switch(forcefield){
122 >    case 1:{
123 >      rparHe = 1.4800;
124 >      epsHe = -0.021270;
125 >    }; break;
126 >    case 2:{
127 >      rparHe = 1.14;
128 >      epsHe = 0.0203;
129 >    }; break;
130 >    case 3:{
131 >      rparHe = 2.28;
132 >      epsHe = 0.020269601874;
133 >    }; break;
134 >    case 4:{
135 >      rparHe = 2.5560;
136 >      epsHe = 0.0200;
137 >    }; break;
138 >    case 5:{
139 >      rparHe = 1.14;
140 >      epsHe = 0.0203;
141 >    }; break;
142 >  }
143 >  
144 >  potEnergy = 0.0;
145  
146 < void GridBuilder::stepTheta(double increment){
104 <        //zero out the euler angles
105 <        for (i=0; i<3; i++)
106 <                angles[i] = 0.0;
107 <        
108 <        //the second euler angle is for rotation about the x-axis (we use the zxz convention)
109 <        angles[1] = increment;
110 <        
111 <        //obtain the rotation matrix through the rigid body class
112 <        rbMol->doEulerToRotMat(angles, rotX);
113 <        
114 <        //rotate the rigid body
115 <        rbMol->getA(rbMatrix);
116 <        matMul3(rotX, rbMatrix, rotatedMat);
117 <        rbMol->setA(rotatedMat);
118 <        
119 < }
146 >  rbMol->getAtomPos(rbAtomPos, 0);
147  
148 < void GridBuilder::stepPhi(double increment){
149 <        //zero out the euler angles
150 <        for (i=0; i<3; i++)
151 <                angles[i] = 0.0;
152 <        
153 <        //the phi euler angle is for rotation about the z-axis (we use the zxz convention)
154 <        angles[0] = increment;
155 <        
156 <        //obtain the rotation matrix through the rigid body class
157 <        rbMol->doEulerToRotMat(angles, rotZ);
158 <        
159 <        //rotate the rigid body
160 <        rbMol->getA(rbMatrix);
161 <        matMul3(rotZ, rbMatrix, rotatedMat);
162 <        rbMol->setA(rotatedMat);
163 <        
148 >  printf("atom0 pos = %lf\t%lf\t%lf\n", rbAtomPos[0], rbAtomPos[1], rbAtomPos[2]);
149 >
150 >
151 >  
152 >  for(i=0; i<rbMol->getNumAtoms(); i++){
153 >    rbMol->getAtomPos(rbAtomPos, i);
154 >    
155 >    rXij = rbAtomPos[0];
156 >    rYij = rbAtomPos[1];
157 >    rZij = rbAtomPos[2] - probeCoor;
158 >    
159 >    rijSquared = rXij * rXij + rYij * rYij + rZij * rZij;
160 >    
161 >    //in the interest of keeping the code more compact, we are being less efficient by placing
162 >    //a switch statement in the calculation loop
163 >    switch(forcefield){
164 >      case 1:{
165 >        //we are using the CHARMm force field
166 >        atomRpar = rbMol->getAtomRpar(i);
167 >        atomEps = rbMol->getAtomEps(i);
168 >        
169 >        rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared);
170 >        rValPowerSix = rValSquared * rValSquared * rValSquared;
171 >        potEnergy += sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0));
172 >      }; break;
173 >      
174 >      case 2:{
175 >        //we are using the AMBER force field
176 >        atomRpar = rbMol->getAtomRpar(i);
177 >        atomEps = rbMol->getAtomEps(i);
178 >        
179 >        rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared);
180 >        rValPowerSix = rValSquared * rValSquared * rValSquared;
181 >        potEnergy += sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0));
182 >      }; break;
183 >      
184 >      case 3:{
185 >        //we are using Allen-Tildesley LJ parameters
186 >        atomRpar = rbMol->getAtomRpar(i);
187 >        atomEps = rbMol->getAtomEps(i);
188 >        
189 >        rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (4*rijSquared);
190 >        rValPowerSix = rValSquared * rValSquared * rValSquared;
191 >        potEnergy += 4*sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 1.0));
192 >        
193 >      }; break;
194 >      
195 >      case 4:{
196 >        //we are using the OPLS force field
197 >        atomRpar = rbMol->getAtomRpar(i);
198 >        atomEps = rbMol->getAtomEps(i);
199 >        
200 >        rValSquared = (pow(sqrt(rparHe+atomRpar),2)) / (rijSquared);
201 >        rValPowerSix = rValSquared * rValSquared * rValSquared;
202 >        potEnergy += 4*sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 1.0));
203 >      }; break;
204 >      
205 >      case 5:{
206 >        //we are using the GAFF force field
207 >        atomRpar = rbMol->getAtomRpar(i);
208 >        atomEps = rbMol->getAtomEps(i);
209 >        
210 >        rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared);
211 >        rValPowerSix = rValSquared * rValSquared * rValSquared;
212 >        potEnergy += sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0));
213 >      }; break;
214 >    }    
215 >  }
216 > }
217 >
218 > void GridBuilder::printGridFiles(){
219 >  ofstream sigmaOut("sigma.grid");
220 >  ofstream sOut("s.grid");
221 >  ofstream epsOut("eps.grid");
222 >  
223 >  for (k=0; k<sigList.size(); k++){
224 >    sigmaOut << sigList[k] << "\n0\n";
225 >    sOut << sList[k] << "\n0\n";    
226 >    epsOut << epsList[k] << "\n0\n";
227 >  }
228   }

Diff Legend

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