--- trunk/SHAPES/GridBuilder.cpp 2004/06/17 21:35:22 1278 +++ trunk/SHAPES/GridBuilder.cpp 2004/06/25 17:54:51 1306 @@ -1,142 +1,226 @@ #include "GridBuilder.hpp" -#include "MatVec3.h" #define PI 3.14159265359 -GridBuilder::GridBuilder(RigidBody* rb, int bandWidth) { - rbMol = rb; - bandwidth = bandWidth; - thetaStep = PI / bandwidth; - thetaMin = thetaStep / 2.0; - phiStep = thetaStep * 2.0; - - //zero out the rot mats - for (i=0; i<3; i++) { - for (j=0; j<3; j++) { - rotX[i][j] = 0.0; - rotZ[i][j] = 0.0; - rbMatrix[i][j] = 0.0; - } - } +GridBuilder::GridBuilder(RigidBody* rb, int gridWidth) { + rbMol = rb; + gridwidth = gridWidth; + thetaStep = PI / gridwidth; + thetaMin = thetaStep / 2.0; + phiStep = thetaStep * 2.0; } GridBuilder::~GridBuilder() { } -void GridBuilder::launchProbe(int forceField, vector sigmaGrid, vector sGrid, - vector epsGrid){ - double startDist; - double minDist = 10.0; //minimum start distance +void GridBuilder::launchProbe(int forceField, vector sigmaGrid, + vector sGrid, vector epsGrid){ + ofstream sigmaOut("sigma.grid"); + ofstream sOut("s.grid"); + ofstream epsOut("eps.grid"); + double startDist; + double phiVal; + double thetaVal; + double sigTemp, sTemp, epsTemp, sigProbe; + double minDist = 10.0; //minimum start distance - //first determine the start distance - we always start at least minDist away - startDist = rbMol->findMaxExtent() + minDist; - if (startDist < minDist) - startDist = minDist; - - initBody(); - for (i=0; ifindMaxExtent() + minDist; + if (startDist < minDist) + startDist = minDist; -void GridBuilder::initBody(){ - //set up the rigid body in the starting configuration - stepTheta(thetaMin); + //set the initial orientation of the body and loop over theta values + + for (k =0; k < gridwidth; k++) { + thetaVal = thetaMin + k*thetaStep; + for (j=0; j < gridwidth; j++) { + phiVal = j*phiStep + 0.5*PI; + if (phiVal>=2*PI) + phiVal -= 2*PI; + + rbMol->setEuler(0.0, thetaVal, phiVal); + + releaseProbe(startDist); + + //translate the values to sigma, s, and epsilon of the rigid body + sigTemp = 2*sigDist - sigProbe; + sTemp = (2*(sDist - sigDist))/(0.122462048309) - sigProbe; + epsTemp = pow(epsVal, 2)/fabs(epsHe); + + sigList.push_back(sigTemp); + sList.push_back(sTemp); + epsList.push_back(epsTemp); + } + } } void GridBuilder::releaseProbe(double farPos){ - int tooClose; - double tempPotEnergy; - double interpRange; - double interpFrac; + int tooClose; + double tempPotEnergy; + double interpRange; + double interpFrac; - probeCoor = farPos; - tooClose = 0; - epsVal = 0; - rhoStep = 0.1; //the distance the probe atom moves between steps - - while (!tooClose){ - calcEnergy(); - potProgress.push_back(potEnergy); - distProgress.push_back(probeCoor); + probeCoor = farPos; + potProgress.clear(); + distProgress.clear(); + tooClose = 0; + epsVal = 0; + rhoStep = 0.1; //the distance the probe atom moves between steps - //if we've reached a new minimum, save the value and position - if (potEnergy < epsVal){ - epsVal = potEnergy; - sDist = probeCoor; - } + while (!tooClose){ + calcEnergy(); + potProgress.push_back(potEnergy); + distProgress.push_back(probeCoor); - //test if the probe reached the origin - if so, stop stepping closer - if (probeCoor < 0){ - sigDist = 0.0; - tooClose = 1; - } + //if we've reached a new minimum, save the value and position + if (potEnergy < epsVal){ + epsVal = potEnergy; + sDist = probeCoor; + } - //test if the probe beyond the contact point - if not, take a step closer - if (potEnergy < 0){ - sigDist = probeCoor; - tempPotEnergy = potEnergy; - probeCoor -= rhoStep; - } - else { - //do a linear interpolation to obtain the sigDist - interpRange = potEnergy - tempPotEnergy; - interpFrac = potEnergy / interpRange; - interpFrac = interpFrac * rhoStep; - sigDist = probeCoor + interpFrac; + //test if the probe reached the origin - if so, stop stepping closer + if (probeCoor < 0){ + sigDist = 0.0; + tooClose = 1; + } + + //test if the probe beyond the contact point - if not, take a step closer + if (potEnergy < 0){ + sigDist = probeCoor; + tempPotEnergy = potEnergy; + probeCoor -= rhoStep; + } + else { + //do a linear interpolation to obtain the sigDist + interpRange = potEnergy - tempPotEnergy; + interpFrac = potEnergy / interpRange; + interpFrac = interpFrac * rhoStep; + sigDist = probeCoor + interpFrac; - //end the loop - tooClose = 1; - } - } + //end the loop + tooClose = 1; + } + } } void GridBuilder::calcEnergy(){ - -} + double rXij, rYij, rZij; + double rijSquared; + double rValSquared, rValPowerSix; + double atomRpar, atomEps; + double rbAtomPos[3]; + + potEnergy = 0.0; -void GridBuilder::stepTheta(double increment){ - //zero out the euler angles - for (i=0; i<3; i++) - angles[i] = 0.0; - - //the second euler angle is for rotation about the x-axis (we use the zxz convention) - angles[1] = increment; - - //obtain the rotation matrix through the rigid body class - rbMol->doEulerToRotMat(angles, rotX); - - //rotate the rigid body - rbMol->getA(rbMatrix); - matMul3(rotX, rbMatrix, rotatedMat); - rbMol->setA(rotatedMat); - -} + for(i=0; igetNumAtoms(); i++){ + rbMol->getAtomPos(rbAtomPos, i); + + rXij = rbAtomPos[0]; + rYij = rbAtomPos[1]; + rZij = rbAtomPos[2] - probeCoor; + + rijSquared = rXij * rXij + rYij * rYij + rZij * rZij; + + //in the interest of keeping the code more compact, we are being less + //efficient by placing a switch statement in the calculation loop + switch(forcefield){ + case 1:{ + //we are using the CHARMm force field + atomRpar = rbMol->getAtomRpar(i); + atomEps = rbMol->getAtomEps(i); + + rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared); + rValPowerSix = rValSquared * rValSquared * rValSquared; + potEnergy += sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0)); + }; break; + + case 2:{ + //we are using the AMBER force field + atomRpar = rbMol->getAtomRpar(i); + atomEps = rbMol->getAtomEps(i); + + rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared); + rValPowerSix = rValSquared * rValSquared * rValSquared; + potEnergy += sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0)); + }; break; + + case 3:{ + //we are using Allen-Tildesley LJ parameters + atomRpar = rbMol->getAtomRpar(i); + atomEps = rbMol->getAtomEps(i); + + rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (4*rijSquared); + rValPowerSix = rValSquared * rValSquared * rValSquared; + potEnergy += 4*sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 1.0)); + + }; break; + + case 4:{ + //we are using the OPLS force field + atomRpar = rbMol->getAtomRpar(i); + atomEps = rbMol->getAtomEps(i); + + rValSquared = (pow(sqrt(rparHe+atomRpar),2)) / (rijSquared); + rValPowerSix = rValSquared * rValSquared * rValSquared; + potEnergy += 4*sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 1.0)); + }; break; + + case 5:{ + //we are using the GAFF force field + atomRpar = rbMol->getAtomRpar(i); + atomEps = rbMol->getAtomEps(i); + + rValSquared = ((rparHe+atomRpar)*(rparHe+atomRpar)) / (rijSquared); + rValPowerSix = rValSquared * rValSquared * rValSquared; + potEnergy += sqrt(epsHe*atomEps)*(rValPowerSix * (rValPowerSix - 2.0)); + }; break; + } + } +} -void GridBuilder::stepPhi(double increment){ - //zero out the euler angles - for (i=0; i<3; i++) - angles[i] = 0.0; - - //the phi euler angle is for rotation about the z-axis (we use the zxz convention) - angles[0] = increment; - - //obtain the rotation matrix through the rigid body class - rbMol->doEulerToRotMat(angles, rotZ); - - //rotate the rigid body - rbMol->getA(rbMatrix); - matMul3(rotZ, rbMatrix, rotatedMat); - rbMol->setA(rotatedMat); - +void GridBuilder::printGridFiles(){ + ofstream sigmaOut("sigma.grid"); + ofstream sOut("s.grid"); + ofstream epsOut("eps.grid"); + + for (k=0; k