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 |
– |
} |
12 |
|
} |
13 |
|
|
14 |
|
GridBuilder::~GridBuilder() { |
15 |
|
} |
16 |
|
|
17 |
< |
void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid, vector<double> sGrid, |
18 |
< |
vector<double> epsGrid){ |
17 |
> |
void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid, |
18 |
> |
vector<double> sGrid, vector<double> epsGrid){ |
19 |
|
ofstream sigmaOut("sigma.grid"); |
20 |
|
ofstream sOut("s.grid"); |
21 |
|
ofstream epsOut("eps.grid"); |
34 |
|
if (startDist < minDist) |
35 |
|
startDist = minDist; |
36 |
|
|
37 |
+ |
printf("startDist = %lf\n", startDist); |
38 |
+ |
|
39 |
|
//set the initial orientation of the body and loop over theta values |
40 |
< |
phiVal = 0.0; |
41 |
< |
thetaVal = thetaMin; |
42 |
< |
rotBody(phiVal, thetaVal); |
43 |
< |
for (k=0; k<bandwidth; k++){ |
44 |
< |
//loop over phi values starting with phi = 0.0 |
45 |
< |
for (j=0; j<bandwidth; j++){ |
40 |
> |
|
41 |
> |
for (k =0; k < bandwidth; k++) { |
42 |
> |
thetaVal = thetaMin + k*thetaStep; |
43 |
> |
for (j=0; j < bandwidth; j++) { |
44 |
> |
phiVal = j*phiStep; |
45 |
> |
|
46 |
> |
printf("setting Euler, phi = %lf\ttheta = %lf\n", phiVal, thetaVal); |
47 |
> |
|
48 |
> |
rbMol->setEuler(0.0, thetaVal, phiVal); |
49 |
> |
|
50 |
|
releaseProbe(startDist); |
51 |
|
|
52 |
+ |
printf("found sigDist = %lf\t sDist = %lf \t epsVal = %lf\n", |
53 |
+ |
sigDist, sDist, epsVal); |
54 |
+ |
|
55 |
|
sigList.push_back(sigDist); |
56 |
|
sList.push_back(sDist); |
57 |
|
epsList.push_back(epsVal); |
58 |
< |
|
59 |
< |
phiVal += phiStep; |
60 |
< |
rotBody(phiVal, thetaVal); |
58 |
> |
|
59 |
|
} |
62 |
– |
phiVal = 0.0; |
63 |
– |
thetaVal += thetaStep; |
64 |
– |
rotBody(phiVal, thetaVal); |
65 |
– |
printf("step theta %i\n",k); |
60 |
|
} |
61 |
|
} |
62 |
|
|
72 |
|
tooClose = 0; |
73 |
|
epsVal = 0; |
74 |
|
rhoStep = 0.1; //the distance the probe atom moves between steps |
75 |
< |
|
82 |
< |
|
75 |
> |
|
76 |
|
while (!tooClose){ |
77 |
|
calcEnergy(); |
78 |
|
potProgress.push_back(potEnergy); |
142 |
|
} |
143 |
|
|
144 |
|
potEnergy = 0.0; |
145 |
+ |
|
146 |
+ |
rbMol->getAtomPos(rbAtomPos, 0); |
147 |
+ |
|
148 |
+ |
printf("atom0 pos = %lf\t%lf\t%lf\n", rbAtomPos[0], rbAtomPos[1], rbAtomPos[2]); |
149 |
+ |
|
150 |
+ |
|
151 |
|
|
152 |
|
for(i=0; i<rbMol->getNumAtoms(); i++){ |
153 |
|
rbMol->getAtomPos(rbAtomPos, i); |
215 |
|
} |
216 |
|
} |
217 |
|
|
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 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 |
– |
//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 |
– |
|
218 |
|
void GridBuilder::printGridFiles(){ |
219 |
|
ofstream sigmaOut("sigma.grid"); |
220 |
|
ofstream sOut("s.grid"); |
225 |
|
sOut << sList[k] << "\n0\n"; |
226 |
|
epsOut << epsList[k] << "\n0\n"; |
227 |
|
} |
228 |
< |
} |
228 |
> |
} |