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 1278 by chrisfen, Thu Jun 17 21:35:22 2004 UTC vs.
Revision 1285 by chrisfen, Tue Jun 22 18:04:58 2004 UTC

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

Diff Legend

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