| 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"); |
| 22 |
|
double startDist; |
| 23 |
+ |
double phiVal; |
| 24 |
+ |
double thetaVal; |
| 25 |
|
double minDist = 10.0; //minimum start distance |
| 26 |
|
|
| 27 |
|
sList = sGrid; |
| 34 |
|
if (startDist < minDist) |
| 35 |
|
startDist = minDist; |
| 36 |
|
|
| 37 |
< |
initBody(); |
| 38 |
< |
for (k=0; k<bandwidth; k++){ |
| 39 |
< |
printf("step theta...\n"); |
| 40 |
< |
for (j=0; j<bandwidth; j++){ |
| 37 |
> |
printf("startDist = %lf\n", startDist); |
| 38 |
> |
|
| 39 |
> |
//set the initial orientation of the body and loop over theta values |
| 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 |
< |
|
| 54 |
< |
stepPhi(phiStep); |
| 58 |
> |
|
| 59 |
|
} |
| 56 |
– |
stepTheta(thetaStep); |
| 60 |
|
} |
| 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 |
– |
*/ |
| 61 |
|
} |
| 62 |
|
|
| 69 |
– |
void GridBuilder::initBody(){ |
| 70 |
– |
//set up the rigid body in the starting configuration |
| 71 |
– |
stepTheta(thetaMin); |
| 72 |
– |
} |
| 73 |
– |
|
| 63 |
|
void GridBuilder::releaseProbe(double farPos){ |
| 64 |
|
int tooClose; |
| 65 |
|
double tempPotEnergy; |
| 72 |
|
tooClose = 0; |
| 73 |
|
epsVal = 0; |
| 74 |
|
rhoStep = 0.1; //the distance the probe atom moves between steps |
| 75 |
< |
|
| 87 |
< |
|
| 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); |
| 191 |
|
potEnergy += 4*sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 1.0)); |
| 192 |
|
|
| 193 |
|
}; break; |
| 200 |
– |
|
| 194 |
|
|
| 195 |
|
case 4:{ |
| 196 |
|
//we are using the OPLS force field |
| 215 |
|
} |
| 216 |
|
} |
| 217 |
|
|
| 225 |
– |
void GridBuilder::stepTheta(double increment){ |
| 226 |
– |
//zero out the euler angles |
| 227 |
– |
for (l=0; l<3; l++) |
| 228 |
– |
angles[i] = 0.0; |
| 229 |
– |
|
| 230 |
– |
//the second euler angle is for rotation about the x-axis (we use the zxz convention) |
| 231 |
– |
angles[1] = increment; |
| 232 |
– |
|
| 233 |
– |
//obtain the rotation matrix through the rigid body class |
| 234 |
– |
rbMol->doEulerToRotMat(angles, rotX); |
| 235 |
– |
|
| 236 |
– |
//rotate the rigid body |
| 237 |
– |
rbMol->getA(rbMatrix); |
| 238 |
– |
matMul3(rotX, rbMatrix, rotatedMat); |
| 239 |
– |
rbMol->setA(rotatedMat); |
| 240 |
– |
} |
| 241 |
– |
|
| 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 |
– |
|
| 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 |
> |
} |
