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 1279 by chrisfen, Fri Jun 18 19:01:40 2004 UTC vs.
Revision 1280 by chrisfen, Fri Jun 18 19:40:31 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;
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 <        }
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 >  }
21   }
22  
23   GridBuilder::~GridBuilder() {
24   }
25  
26   void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid, vector<double> sGrid,
27 <        vector<double> epsGrid){
28 <        double startDist;
29 <        double minDist = 10.0; //minimum start distance
27 >                              vector<double> epsGrid){
28 >  double startDist;
29 >  double minDist = 10.0; //minimum start distance
30          
31 <        //first determine the start distance - we always start at least minDist away
32 <        startDist = rbMol->findMaxExtent() + minDist;
33 <        if (startDist < minDist)
34 <                startDist = minDist;
31 >  forcefield = forceField;
32 >    
33 >  //first determine the start distance - we always start at least minDist away
34 >  startDist = rbMol->findMaxExtent() + minDist;
35 >  if (startDist < minDist)
36 >    startDist = minDist;
37          
38 <        initBody();
39 <        for (i=0; i<bandwidth; i++){            
40 <                for (j=0; j<bandwidth; j++){
41 <                        releaseProbe(startDist);
38 >  initBody();
39 >  for (i=0; i<bandwidth; i++){          
40 >    for (j=0; j<bandwidth; j++){
41 >      releaseProbe(startDist);
42  
43 <                        sigmaGrid.push_back(sigDist);
44 <                        sGrid.push_back(sDist);
45 <                        epsGrid.push_back(epsVal);
43 >      sigmaGrid.push_back(sigDist);
44 >      sGrid.push_back(sDist);
45 >      epsGrid.push_back(epsVal);
46                          
47 <                        stepPhi(phiStep);
48 <                }
49 <                stepTheta(thetaStep);
50 <        }
49 <                
47 >      stepPhi(phiStep);
48 >    }
49 >    stepTheta(thetaStep);
50 >  }            
51   }
52  
53   void GridBuilder::initBody(){
54 <        //set up the rigid body in the starting configuration
55 <        stepTheta(thetaMin);
54 >  //set up the rigid body in the starting configuration
55 >  stepTheta(thetaMin);
56   }
57  
58   void GridBuilder::releaseProbe(double farPos){
59 <        int tooClose;
60 <        double tempPotEnergy;
61 <        double interpRange;
62 <        double interpFrac;
59 >  int tooClose;
60 >  double tempPotEnergy;
61 >  double interpRange;
62 >  double interpFrac;
63          
64 <        probeCoor = farPos;
65 <        potProgress.clear();
66 <        distProgress.clear();
67 <        tooClose = 0;
68 <        epsVal = 0;
69 <        rhoStep = 0.1; //the distance the probe atom moves between steps
64 >  probeCoor = farPos;
65 >  potProgress.clear();
66 >  distProgress.clear();
67 >  tooClose = 0;
68 >  epsVal = 0;
69 >  rhoStep = 0.1; //the distance the probe atom moves between steps
70          
71          
72 <        while (!tooClose){
73 <                calcEnergy();
74 <                potProgress.push_back(potEnergy);
75 <                distProgress.push_back(probeCoor);
72 >  while (!tooClose){
73 >    calcEnergy();
74 >    potProgress.push_back(potEnergy);
75 >    distProgress.push_back(probeCoor);
76                  
77 <                //if we've reached a new minimum, save the value and position
78 <                if (potEnergy < epsVal){
79 <                        epsVal = potEnergy;
80 <                        sDist = probeCoor;
81 <                }
77 >    //if we've reached a new minimum, save the value and position
78 >    if (potEnergy < epsVal){
79 >      epsVal = potEnergy;
80 >      sDist = probeCoor;
81 >    }
82                  
83 <                //test if the probe reached the origin - if so, stop stepping closer
84 <                if (probeCoor < 0){
85 <                        sigDist = 0.0;
86 <                        tooClose = 1;
87 <                }
83 >    //test if the probe reached the origin - if so, stop stepping closer
84 >    if (probeCoor < 0){
85 >      sigDist = 0.0;
86 >      tooClose = 1;
87 >    }
88                  
89 <                //test if the probe beyond the contact point - if not, take a step closer
90 <                if (potEnergy < 0){
91 <                        sigDist = probeCoor;
92 <                        tempPotEnergy = potEnergy;
93 <                        probeCoor -= rhoStep;
94 <                }
95 <                else {
96 <                        //do a linear interpolation to obtain the sigDist
97 <                        interpRange = potEnergy - tempPotEnergy;
98 <                        interpFrac = potEnergy / interpRange;
99 <                        interpFrac = interpFrac * rhoStep;
100 <                        sigDist = probeCoor + interpFrac;
89 >    //test if the probe beyond the contact point - if not, take a step closer
90 >    if (potEnergy < 0){
91 >      sigDist = probeCoor;
92 >      tempPotEnergy = potEnergy;
93 >      probeCoor -= rhoStep;
94 >    }
95 >    else {
96 >      //do a linear interpolation to obtain the sigDist
97 >      interpRange = potEnergy - tempPotEnergy;
98 >      interpFrac = potEnergy / interpRange;
99 >      interpFrac = interpFrac * rhoStep;
100 >      sigDist = probeCoor + interpFrac;
101                          
102 <                        //end the loop
103 <                        tooClose = 1;
104 <                }
105 <        }
102 >      //end the loop
103 >      tooClose = 1;
104 >    }
105 >  }
106   }
107  
108   void GridBuilder::calcEnergy(){
# Line 109 | Line 110 | void GridBuilder::stepTheta(double increment){
110   }
111  
112   void GridBuilder::stepTheta(double increment){
113 <        //zero out the euler angles
114 <        for (i=0; i<3; i++)
115 <                angles[i] = 0.0;
113 >  //zero out the euler angles
114 >  for (i=0; i<3; i++)
115 >    angles[i] = 0.0;
116          
117 <        //the second euler angle is for rotation about the x-axis (we use the zxz convention)
118 <        angles[1] = increment;
117 >  //the second euler angle is for rotation about the x-axis (we use the zxz convention)
118 >  angles[1] = increment;
119          
120 <        //obtain the rotation matrix through the rigid body class
121 <        rbMol->doEulerToRotMat(angles, rotX);
120 >  //obtain the rotation matrix through the rigid body class
121 >  rbMol->doEulerToRotMat(angles, rotX);
122          
123 <        //rotate the rigid body
124 <        rbMol->getA(rbMatrix);
125 <        matMul3(rotX, rbMatrix, rotatedMat);
126 <        rbMol->setA(rotatedMat);
126 <        
123 >  //rotate the rigid body
124 >  rbMol->getA(rbMatrix);
125 >  matMul3(rotX, rbMatrix, rotatedMat);
126 >  rbMol->setA(rotatedMat);      
127   }
128  
129   void GridBuilder::stepPhi(double increment){
130 <        //zero out the euler angles
131 <        for (i=0; i<3; i++)
132 <                angles[i] = 0.0;
130 >  //zero out the euler angles
131 >  for (i=0; i<3; i++)
132 >    angles[i] = 0.0;
133          
134 <        //the phi euler angle is for rotation about the z-axis (we use the zxz convention)
135 <        angles[0] = increment;
134 >  //the phi euler angle is for rotation about the z-axis (we use the zxz convention)
135 >  angles[0] = increment;
136          
137 <        //obtain the rotation matrix through the rigid body class
138 <        rbMol->doEulerToRotMat(angles, rotZ);
137 >  //obtain the rotation matrix through the rigid body class
138 >  rbMol->doEulerToRotMat(angles, rotZ);
139          
140 <        //rotate the rigid body
141 <        rbMol->getA(rbMatrix);
142 <        matMul3(rotZ, rbMatrix, rotatedMat);
143 <        rbMol->setA(rotatedMat);
144 <        
140 >  //rotate the rigid body
141 >  rbMol->getA(rbMatrix);
142 >  matMul3(rotZ, rbMatrix, rotatedMat);
143 >  rbMol->setA(rotatedMat);      
144   }

Diff Legend

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