| 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 |
|
sList = sGrid; |
| 43 |
|
if (startDist < minDist) |
| 44 |
|
startDist = minDist; |
| 45 |
|
|
| 46 |
< |
initBody(); |
| 47 |
< |
for (k=0; k<bandwidth; k++){ |
| 48 |
< |
printf("step theta...\n"); |
| 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 |
|
sigList.push_back(sigDist); |
| 56 |
|
sList.push_back(sDist); |
| 57 |
|
epsList.push_back(epsVal); |
| 58 |
< |
|
| 59 |
< |
stepPhi(phiStep); |
| 58 |
> |
|
| 59 |
> |
phiVal += phiStep; |
| 60 |
> |
rotBody(phiVal, thetaVal); |
| 61 |
|
} |
| 62 |
< |
stepTheta(thetaStep); |
| 62 |
> |
phiVal = 0.0; |
| 63 |
> |
thetaVal += thetaStep; |
| 64 |
> |
rotBody(phiVal, thetaVal); |
| 65 |
> |
printf("step theta %i\n",k); |
| 66 |
|
} |
| 58 |
– |
/* |
| 59 |
– |
//write out the grid files |
| 60 |
– |
printf("the grid size is %d\n",sigmaGrid.size()); |
| 61 |
– |
for (k=0; k<sigmaGrid.size(); k++){ |
| 62 |
– |
sigmaOut << sigmaGrid[k] << "\n0\n"; |
| 63 |
– |
sOut << sGrid[k] << "\n0\n"; |
| 64 |
– |
epsOut << epsGrid[k] << "\n0\n"; |
| 65 |
– |
} |
| 66 |
– |
*/ |
| 67 |
– |
} |
| 68 |
– |
|
| 69 |
– |
void GridBuilder::initBody(){ |
| 70 |
– |
//set up the rigid body in the starting configuration |
| 71 |
– |
stepTheta(thetaMin); |
| 67 |
|
} |
| 68 |
|
|
| 69 |
|
void GridBuilder::releaseProbe(double farPos){ |
| 192 |
|
potEnergy += 4*sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 1.0)); |
| 193 |
|
|
| 194 |
|
}; break; |
| 200 |
– |
|
| 195 |
|
|
| 196 |
|
case 4:{ |
| 197 |
|
//we are using the OPLS force field |
| 216 |
|
} |
| 217 |
|
} |
| 218 |
|
|
| 219 |
< |
void GridBuilder::stepTheta(double increment){ |
| 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] = increment; |
| 227 |
> |
angles[1] = tValue; |
| 228 |
|
|
| 229 |
|
//obtain the rotation matrix through the rigid body class |
| 230 |
|
rbMol->doEulerToRotMat(angles, rotX); |
| 231 |
< |
|
| 231 |
> |
|
| 232 |
> |
//start from the reference position |
| 233 |
> |
identityMat3(rbMatrix); |
| 234 |
> |
rbMol->setA(rbMatrix); |
| 235 |
> |
|
| 236 |
|
//rotate the rigid body |
| 237 |
– |
rbMol->getA(rbMatrix); |
| 237 |
|
matMul3(rotX, rbMatrix, rotatedMat); |
| 238 |
|
rbMol->setA(rotatedMat); |
| 239 |
|
} |
| 240 |
|
|
| 242 |
– |
void GridBuilder::stepPhi(double increment){ |
| 243 |
– |
//zero out the euler angles |
| 244 |
– |
for (l=0; l<3; l++) |
| 245 |
– |
angles[i] = 0.0; |
| 246 |
– |
|
| 247 |
– |
//the phi euler angle is for rotation about the z-axis (we use the zxz convention) |
| 248 |
– |
angles[0] = increment; |
| 249 |
– |
|
| 250 |
– |
//obtain the rotation matrix through the rigid body class |
| 251 |
– |
rbMol->doEulerToRotMat(angles, rotZ); |
| 252 |
– |
|
| 253 |
– |
//rotate the rigid body |
| 254 |
– |
rbMol->getA(rbMatrix); |
| 255 |
– |
matMul3(rotZ, rbMatrix, rotatedMat); |
| 256 |
– |
rbMol->setA(rotatedMat); |
| 257 |
– |
} |
| 258 |
– |
|
| 241 |
|
void GridBuilder::printGridFiles(){ |
| 242 |
|
ofstream sigmaOut("sigma.grid"); |
| 243 |
|
ofstream sOut("s.grid"); |
