# | 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) { |
5 | > | GridBuilder::GridBuilder(RigidBody* rb, int gridWidth) { |
6 | rbMol = rb; | |
7 | < | bandwidth = bandWidth; |
8 | < | thetaStep = PI / bandwidth; |
7 | > | gridwidth = gridWidth; |
8 | > | thetaStep = PI / gridwidth; |
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 | – | } |
11 | } | |
12 | ||
13 | GridBuilder::~GridBuilder() { | |
14 | } | |
15 | ||
16 | < | void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid, vector<double> sGrid, |
17 | < | vector<double> epsGrid){ |
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 | + | 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 | < | //first determine the start distance - we always start at least minDist away |
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 | < | initBody(); |
67 | < | for (i=0; i<bandwidth; i++){ |
68 | < | for (j=0; j<bandwidth; j++){ |
65 | > | |
66 | > | //set the initial orientation of the body and loop over theta values |
67 | > | |
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 | > | |
73 | > | rbMol->setEuler(0.0, thetaVal, phiVal); |
74 | > | |
75 | releaseProbe(startDist); | |
76 | ||
77 | < | sigmaGrid.push_back(sigDist); |
78 | < | sGrid.push_back(sDist); |
79 | < | epsGrid.push_back(epsVal); |
80 | < | |
81 | < | stepPhi(phiStep); |
77 | > | //translate the values to sigma, s, and epsilon of the rigid body |
78 | > | sigTemp = 2*sigDist - sigProbe; |
79 | > | sTemp = (2*(sDist - sigDist))/(0.122462048309) - sigProbe; |
80 | > | epsTemp = pow(epsVal, 2)/fabs(epsHe); |
81 | > | |
82 | > | sigList.push_back(sigTemp); |
83 | > | sList.push_back(sTemp); |
84 | > | epsList.push_back(epsTemp); |
85 | } | |
49 | – | stepTheta(thetaStep); |
86 | } | |
87 | } | |
88 | ||
53 | – | void GridBuilder::initBody(){ |
54 | – | //set up the rigid body in the starting configuration |
55 | – | stepTheta(thetaMin); |
56 | – | } |
57 | – | |
89 | void GridBuilder::releaseProbe(double farPos){ | |
90 | int tooClose; | |
91 | double tempPotEnergy; | |
# | Line 67 | Line 98 | void GridBuilder::releaseProbe(double farPos){ | |
98 | tooClose = 0; | |
99 | epsVal = 0; | |
100 | rhoStep = 0.1; //the distance the probe atom moves between steps | |
101 | < | |
71 | < | |
101 | > | |
102 | while (!tooClose){ | |
103 | calcEnergy(); | |
104 | potProgress.push_back(potEnergy); | |
# | Line 106 | Line 136 | void GridBuilder::calcEnergy(){ | |
136 | } | |
137 | ||
138 | void GridBuilder::calcEnergy(){ | |
139 | < | |
140 | < | } |
139 | > | double rXij, rYij, rZij; |
140 | > | double rijSquared; |
141 | > | double rValSquared, rValPowerSix; |
142 | > | double atomRpar, atomEps; |
143 | > | double rbAtomPos[3]; |
144 | > | |
145 | > | potEnergy = 0.0; |
146 | ||
147 | < | void GridBuilder::stepTheta(double increment){ |
148 | < | //zero out the euler angles |
149 | < | for (i=0; i<3; i++) |
150 | < | angles[i] = 0.0; |
151 | < | |
152 | < | //the second euler angle is for rotation about the x-axis (we use the zxz convention) |
153 | < | angles[1] = increment; |
154 | < | |
155 | < | //obtain the rotation matrix through the rigid body class |
156 | < | rbMol->doEulerToRotMat(angles, rotX); |
157 | < | |
158 | < | //rotate the rigid body |
159 | < | rbMol->getA(rbMatrix); |
160 | < | matMul3(rotX, rbMatrix, rotatedMat); |
161 | < | rbMol->setA(rotatedMat); |
162 | < | } |
147 | > | for(i=0; i<rbMol->getNumAtoms(); i++){ |
148 | > | rbMol->getAtomPos(rbAtomPos, i); |
149 | > | |
150 | > | rXij = rbAtomPos[0]; |
151 | > | rYij = rbAtomPos[1]; |
152 | > | rZij = rbAtomPos[2] - probeCoor; |
153 | > | |
154 | > | rijSquared = rXij * rXij + rYij * rYij + rZij * rZij; |
155 | > | |
156 | > | //in the interest of keeping the code more compact, we are being less |
157 | > | //efficient by placing a switch statement in the calculation loop |
158 | > | switch(forcefield){ |
159 | > | case 1:{ |
160 | > | //we are using the CHARMm force field |
161 | > | atomRpar = rbMol->getAtomRpar(i); |
162 | > | atomEps = rbMol->getAtomEps(i); |
163 | > | |
164 | > | rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared); |
165 | > | rValPowerSix = rValSquared * rValSquared * rValSquared; |
166 | > | potEnergy += sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0)); |
167 | > | }; break; |
168 | > | |
169 | > | case 2:{ |
170 | > | //we are using the AMBER force field |
171 | > | atomRpar = rbMol->getAtomRpar(i); |
172 | > | atomEps = rbMol->getAtomEps(i); |
173 | > | |
174 | > | rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared); |
175 | > | rValPowerSix = rValSquared * rValSquared * rValSquared; |
176 | > | potEnergy += sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0)); |
177 | > | }; break; |
178 | > | |
179 | > | case 3:{ |
180 | > | //we are using Allen-Tildesley LJ parameters |
181 | > | atomRpar = rbMol->getAtomRpar(i); |
182 | > | atomEps = rbMol->getAtomEps(i); |
183 | > | |
184 | > | rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (4*rijSquared); |
185 | > | rValPowerSix = rValSquared * rValSquared * rValSquared; |
186 | > | potEnergy += 4*sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 1.0)); |
187 | > | |
188 | > | }; break; |
189 | > | |
190 | > | case 4:{ |
191 | > | //we are using the OPLS force field |
192 | > | atomRpar = rbMol->getAtomRpar(i); |
193 | > | atomEps = rbMol->getAtomEps(i); |
194 | > | |
195 | > | rValSquared = (pow(sqrt(rparHe+atomRpar),2)) / (rijSquared); |
196 | > | rValPowerSix = rValSquared * rValSquared * rValSquared; |
197 | > | potEnergy += 4*sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 1.0)); |
198 | > | }; break; |
199 | > | |
200 | > | case 5:{ |
201 | > | //we are using the GAFF force field |
202 | > | atomRpar = rbMol->getAtomRpar(i); |
203 | > | atomEps = rbMol->getAtomEps(i); |
204 | > | |
205 | > | rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared); |
206 | > | rValPowerSix = rValSquared * rValSquared * rValSquared; |
207 | > | potEnergy += sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0)); |
208 | > | }; break; |
209 | > | } |
210 | > | } |
211 | > | } |
212 | ||
213 | < | void GridBuilder::stepPhi(double increment){ |
214 | < | //zero out the euler angles |
215 | < | for (i=0; i<3; i++) |
216 | < | angles[i] = 0.0; |
217 | < | |
218 | < | //the phi euler angle is for rotation about the z-axis (we use the zxz convention) |
219 | < | angles[0] = increment; |
220 | < | |
221 | < | //obtain the rotation matrix through the rigid body class |
222 | < | rbMol->doEulerToRotMat(angles, rotZ); |
139 | < | |
140 | < | //rotate the rigid body |
141 | < | rbMol->getA(rbMatrix); |
142 | < | matMul3(rotZ, rbMatrix, rotatedMat); |
143 | < | rbMol->setA(rotatedMat); |
213 | > | void GridBuilder::printGridFiles(){ |
214 | > | ofstream sigmaOut("sigma.grid"); |
215 | > | ofstream sOut("s.grid"); |
216 | > | ofstream epsOut("eps.grid"); |
217 | > | |
218 | > | for (k=0; k<sigList.size(); k++){ |
219 | > | sigmaOut << sigList[k] << "\n0\n"; |
220 | > | sOut << sList[k] << "\n0\n"; |
221 | > | epsOut << epsList[k] << "\n0\n"; |
222 | > | } |
223 | } | |
224 | + |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |