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