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