# | Line 1 | Line 1 | |
---|---|---|
1 | #include "GridBuilder.hpp" | |
2 | – | #include "MatVec3.h" |
2 | #define PI 3.14159265359 | |
3 | ||
4 | ||
5 | < | GridBuilder::GridBuilder(RigidBody* rb, int bandWidth) { |
6 | < | rbMol = rb; |
7 | < | bandwidth = bandWidth; |
8 | < | thetaStep = PI / bandwidth; |
9 | < | thetaMin = thetaStep / 2.0; |
10 | < | 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 | < | } |
5 | > | GridBuilder::GridBuilder(RigidBody* rb, int gridWidth) { |
6 | > | rbMol = rb; |
7 | > | gridwidth = gridWidth; |
8 | > | thetaStep = PI / gridwidth; |
9 | > | thetaMin = thetaStep / 2.0; |
10 | > | phiStep = thetaStep * 2.0; |
11 | } | |
12 | ||
13 | GridBuilder::~GridBuilder() { | |
14 | } | |
15 | ||
16 | < | void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid, vector<double> sGrid, |
17 | < | vector<double> epsGrid){ |
18 | < | double startDist; |
19 | < | double minDist = 10.0; //minimum start distance |
16 | > | void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid, |
17 | > | vector<double> sGrid, vector<double> epsGrid){ |
18 | > | ofstream sigmaOut("sigma.grid"); |
19 | > | ofstream sOut("s.grid"); |
20 | > | ofstream epsOut("eps.grid"); |
21 | > | double startDist; |
22 | > | double phiVal; |
23 | > | double thetaVal; |
24 | > | double sigTemp, sTemp, epsTemp, sigProbe; |
25 | > | double minDist = 10.0; //minimum start distance |
26 | ||
27 | < | //first determine the start distance - we always start at least minDist away |
28 | < | startDist = rbMol->findMaxExtent() + minDist; |
29 | < | if (startDist < minDist) |
30 | < | startDist = minDist; |
31 | < | |
32 | < | initBody(); |
33 | < | for (i=0; i<bandwidth; i++){ |
34 | < | for (j=0; j<bandwidth; j++){ |
35 | < | releaseProbe(startDist); |
27 | > | sigList = sigmaGrid; |
28 | > | sList = sGrid; |
29 | > | epsList = epsGrid; |
30 | > | forcefield = forceField; |
31 | > | |
32 | > | //load the probe atom parameters |
33 | > | switch(forcefield){ |
34 | > | case 1:{ |
35 | > | rparHe = 1.4800; |
36 | > | epsHe = -0.021270; |
37 | > | }; break; |
38 | > | case 2:{ |
39 | > | rparHe = 1.14; |
40 | > | epsHe = 0.0203; |
41 | > | }; break; |
42 | > | case 3:{ |
43 | > | rparHe = 2.28; |
44 | > | epsHe = 0.020269601874; |
45 | > | }; break; |
46 | > | case 4:{ |
47 | > | rparHe = 2.5560; |
48 | > | epsHe = 0.0200; |
49 | > | }; break; |
50 | > | case 5:{ |
51 | > | rparHe = 1.14; |
52 | > | epsHe = 0.0203; |
53 | > | }; break; |
54 | > | } |
55 | > | |
56 | > | if (rparHe < 2.2) |
57 | > | sigProbe = 2*rparHe/1.12246204831; |
58 | > | else |
59 | > | sigProbe = rparHe; |
60 | > | |
61 | > | //determine the start distance - we always start at least minDist away |
62 | > | startDist = rbMol->findMaxExtent() + minDist; |
63 | > | if (startDist < minDist) |
64 | > | startDist = minDist; |
65 | ||
66 | < | sigmaGrid.push_back(sigDist); |
42 | < | sGrid.push_back(sDist); |
43 | < | epsGrid.push_back(epsVal); |
44 | < | |
45 | < | stepPhi(phiStep); |
46 | < | } |
47 | < | stepTheta(thetaStep); |
48 | < | } |
49 | < | |
50 | < | } |
66 | > | //set the initial orientation of the body and loop over theta values |
67 | ||
68 | < | void GridBuilder::initBody(){ |
69 | < | //set up the rigid body in the starting configuration |
70 | < | stepTheta(thetaMin); |
68 | > | for (k =0; k < gridwidth; k++) { |
69 | > | thetaVal = thetaMin + k*thetaStep; |
70 | > | for (j=0; j < gridwidth; j++) { |
71 | > | phiVal = j*phiStep + 0.5*PI; |
72 | > | if (phiVal>=2*PI) |
73 | > | phiVal -= 2*PI; |
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 | < | potProgress.clear(); |
99 | < | distProgress.clear(); |
100 | < | tooClose = 0; |
101 | < | epsVal = 0; |
102 | < | rhoStep = 0.1; //the distance the probe atom moves between steps |
69 | < | |
70 | < | |
71 | < | while (!tooClose){ |
72 | < | calcEnergy(); |
73 | < | potProgress.push_back(potEnergy); |
74 | < | 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; |
80 | < | } |
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); |
139 | < | |
140 | < | //rotate the rigid body |
141 | < | rbMol->getA(rbMatrix); |
142 | < | matMul3(rotZ, rbMatrix, rotatedMat); |
143 | < | rbMol->setA(rotatedMat); |
144 | < | |
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 | } | |
226 | + |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |