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 1282 by chrisfen, Mon Jun 21 15:10:29 2004 UTC vs.
Revision 1285 by chrisfen, Tue Jun 22 18:04:58 2004 UTC

# Line 3 | Line 3
3   #define PI 3.14159265359
4  
5  
6 < GridBuilder::GridBuilder(RigidBody* rb, int bandWidth) {
6 > GridBuilder::GridBuilder(RigidBody* rb, int gridWidth) {
7    rbMol = rb;
8 <  bandwidth = bandWidth;
9 <  thetaStep = PI / bandwidth;
8 >  gridwidth = gridWidth;
9 >  thetaStep = PI / gridwidth;
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 –  }
12   }
13  
14   GridBuilder::~GridBuilder() {
15   }
16  
17 < void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid, vector<double> sGrid,
18 <                              vector<double> epsGrid){
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    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 <  //first determine the start distance - we always start at least minDist away
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    //set the initial orientation of the body and loop over theta values
68 <  phiVal = 0.0;
69 <  thetaVal = thetaMin;
70 <  rotBody(phiVal, thetaVal);
71 <  for (k=0; k<bandwidth; k++){  
72 <        //loop over phi values starting with phi = 0.0
73 <    for (j=0; j<bandwidth; j++){
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 <      sigList.push_back(sigDist);
80 <      sList.push_back(sDist);
81 <      epsList.push_back(epsVal);
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 <      phiVal += phiStep;
85 <      rotBody(phiVal, thetaVal);
84 >      sigList.push_back(sigTemp);
85 >      sList.push_back(sTemp);
86 >      epsList.push_back(epsTemp);
87      }
62 –    phiVal = 0.0;
63 –    thetaVal += thetaStep;
64 –    rotBody(phiVal, thetaVal);
65 –    printf("step theta %i\n",k);
88    }            
89   }
90  
# Line 78 | Line 100 | void GridBuilder::releaseProbe(double farPos){
100    tooClose = 0;
101    epsVal = 0;
102    rhoStep = 0.1; //the distance the probe atom moves between steps
103 <        
82 <        
103 >                
104    while (!tooClose){
105      calcEnergy();
106      potProgress.push_back(potEnergy);
# Line 120 | Line 141 | void GridBuilder::calcEnergy(){
141    double rXij, rYij, rZij;
142    double rijSquared;
143    double rValSquared, rValPowerSix;
123 –  double rparHe, epsHe;
144    double atomRpar, atomEps;
145    double rbAtomPos[3];
146 <  
127 <  //first get the probe atom parameters
128 <  switch(forcefield){
129 <    case 1:{
130 <      rparHe = 1.4800;
131 <      epsHe = -0.021270;
132 <    }; break;
133 <    case 2:{
134 <      rparHe = 1.14;
135 <      epsHe = 0.0203;
136 <    }; break;
137 <    case 3:{
138 <      rparHe = 2.28;
139 <      epsHe = 0.020269601874;
140 <    }; break;
141 <    case 4:{
142 <      rparHe = 2.5560;
143 <      epsHe = 0.0200;
144 <    }; break;
145 <    case 5:{
146 <      rparHe = 1.14;
147 <      epsHe = 0.0203;
148 <    }; break;
149 <  }
150 <  
146 >    
147    potEnergy = 0.0;
148 <  
148 >
149    for(i=0; i<rbMol->getNumAtoms(); i++){
150      rbMol->getAtomPos(rbAtomPos, i);
151      
# Line 159 | Line 155 | void GridBuilder::calcEnergy(){
155      
156      rijSquared = rXij * rXij + rYij * rYij + rZij * rZij;
157      
158 <    //in the interest of keeping the code more compact, we are being less efficient by placing
159 <    //a switch statement in the calculation loop
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
# Line 216 | Line 212 | void GridBuilder::calcEnergy(){
212    }
213   }
214  
219 – void GridBuilder::rotBody(double pValue, double tValue){
220 –  //zero out the euler angles
221 –  for (l=0; l<3; l++)
222 –    angles[i] = 0.0;
223 –        
224 –  //the phi euler angle is for rotation about the z-axis (we use the zxz convention)
225 –  angles[0] = pValue;
226 –  //the second euler angle is for rotation about the x-axis (we use the zxz convention)
227 –  angles[1] = tValue;
228 –        
229 –  //obtain the rotation matrix through the rigid body class
230 –  rbMol->doEulerToRotMat(angles, rotX);
231 –  
232 –  //start from the reference position
233 –  identityMat3(rbMatrix);
234 –  rbMol->setA(rbMatrix);
235 –  
236 –  //rotate the rigid body
237 –  matMul3(rotX, rbMatrix, rotatedMat);
238 –  rbMol->setA(rotatedMat);      
239 – }
240 –
215   void GridBuilder::printGridFiles(){
216    ofstream sigmaOut("sigma.grid");
217    ofstream sOut("s.grid");
# Line 248 | Line 222 | void GridBuilder::printGridFiles(){
222      sOut << sList[k] << "\n0\n";    
223      epsOut << epsList[k] << "\n0\n";
224    }
225 < }
225 > }

Diff Legend

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