trunk/SHAPES/GridBuilder.cpp

#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() {
}

void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid, vector<double> sGrid,
                              vector<double> epsGrid){
  double startDist;
  double minDist = 10.0; //minimum start distance
        
  forcefield = forceField;
    
  //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; i<bandwidth; i++){          
    for (j=0; j<bandwidth; j++){
      releaseProbe(startDist);

      sigmaGrid.push_back(sigDist);
      sGrid.push_back(sDist);
      epsGrid.push_back(epsVal);
                        
      stepPhi(phiStep);
    }
    stepTheta(thetaStep);
  }             
}

void GridBuilder::initBody(){
  //set up the rigid body in the starting configuration
  stepTheta(thetaMin);
}

void GridBuilder::releaseProbe(double farPos){
  int tooClose;
  double tempPotEnergy;
  double interpRange;
  double interpFrac;
        
  probeCoor = farPos;
  potProgress.clear();
  distProgress.clear();
  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);
                
    //if we've reached a new minimum, save the value and position
    if (potEnergy < epsVal){
      epsVal = potEnergy;
      sDist = probeCoor;
    }
                
    //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;
    }
  }
}

void GridBuilder::calcEnergy(){
        
}

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);      
}

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);      
}
Revision:	1280
Committed:	Fri Jun 18 19:40:31 2004 UTC (20 years ago) by chrisfen
File size:	3454 byte(s)
Log Message:	Fixed XCode auto-indent to match XEmacs
#	User	Rev	Content
1	chrisfen	1277	#include "GridBuilder.hpp"
2			#include "MatVec3.h"
3			#define PI 3.14159265359
4
5
6			GridBuilder::GridBuilder(RigidBody* rb, int bandWidth) {
7	chrisfen	1280	rbMol = rb;
8			bandwidth = bandWidth;
9			thetaStep = PI / bandwidth;
10			thetaMin = thetaStep / 2.0;
11			phiStep = thetaStep * 2.0;
12	chrisfen	1277
13	chrisfen	1280	//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			}
21	chrisfen	1277	}
22
23			GridBuilder::~GridBuilder() {
24			}
25
26			void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid, vector<double> sGrid,
27	chrisfen	1280	vector<double> epsGrid){
28			double startDist;
29			double minDist = 10.0; //minimum start distance
30	chrisfen	1277
31	chrisfen	1280	forcefield = forceField;
32
33			//first determine the start distance - we always start at least minDist away
34			startDist = rbMol->findMaxExtent() + minDist;
35			if (startDist < minDist)
36			startDist = minDist;
37	chrisfen	1277
38	chrisfen	1280	initBody();
39			for (i=0; i<bandwidth; i++){
40			for (j=0; j<bandwidth; j++){
41			releaseProbe(startDist);
42	chrisfen	1279
43	chrisfen	1280	sigmaGrid.push_back(sigDist);
44			sGrid.push_back(sDist);
45			epsGrid.push_back(epsVal);
46	chrisfen	1278
47	chrisfen	1280	stepPhi(phiStep);
48			}
49			stepTheta(thetaStep);
50			}
51	chrisfen	1277	}
52
53			void GridBuilder::initBody(){
54	chrisfen	1280	//set up the rigid body in the starting configuration
55			stepTheta(thetaMin);
56	chrisfen	1277	}
57
58			void GridBuilder::releaseProbe(double farPos){
59	chrisfen	1280	int tooClose;
60			double tempPotEnergy;
61			double interpRange;
62			double interpFrac;
63	chrisfen	1277
64	chrisfen	1280	probeCoor = farPos;
65			potProgress.clear();
66			distProgress.clear();
67			tooClose = 0;
68			epsVal = 0;
69			rhoStep = 0.1; //the distance the probe atom moves between steps
70	chrisfen	1277
71	chrisfen	1279
72	chrisfen	1280	while (!tooClose){
73			calcEnergy();
74			potProgress.push_back(potEnergy);
75			distProgress.push_back(probeCoor);
76	chrisfen	1277
77	chrisfen	1280	//if we've reached a new minimum, save the value and position
78			if (potEnergy < epsVal){
79			epsVal = potEnergy;
80			sDist = probeCoor;
81			}
82	chrisfen	1277
83	chrisfen	1280	//test if the probe reached the origin - if so, stop stepping closer
84			if (probeCoor < 0){
85			sigDist = 0.0;
86			tooClose = 1;
87			}
88	chrisfen	1277
89	chrisfen	1280	//test if the probe beyond the contact point - if not, take a step closer
90			if (potEnergy < 0){
91			sigDist = probeCoor;
92			tempPotEnergy = potEnergy;
93			probeCoor -= rhoStep;
94			}
95			else {
96			//do a linear interpolation to obtain the sigDist
97			interpRange = potEnergy - tempPotEnergy;
98			interpFrac = potEnergy / interpRange;
99			interpFrac = interpFrac * rhoStep;
100			sigDist = probeCoor + interpFrac;
101	chrisfen	1277
102	chrisfen	1280	//end the loop
103			tooClose = 1;
104			}
105			}
106	chrisfen	1277	}
107
108			void GridBuilder::calcEnergy(){
109
110			}
111
112			void GridBuilder::stepTheta(double increment){
113	chrisfen	1280	//zero out the euler angles
114			for (i=0; i<3; i++)
115			angles[i] = 0.0;
116	chrisfen	1277
117	chrisfen	1280	//the second euler angle is for rotation about the x-axis (we use the zxz convention)
118			angles[1] = increment;
119	chrisfen	1277
120	chrisfen	1280	//obtain the rotation matrix through the rigid body class
121			rbMol->doEulerToRotMat(angles, rotX);
122	chrisfen	1277
123	chrisfen	1280	//rotate the rigid body
124			rbMol->getA(rbMatrix);
125			matMul3(rotX, rbMatrix, rotatedMat);
126			rbMol->setA(rotatedMat);
127	chrisfen	1277	}
128
129			void GridBuilder::stepPhi(double increment){
130	chrisfen	1280	//zero out the euler angles
131			for (i=0; i<3; i++)
132			angles[i] = 0.0;
133	chrisfen	1277
134	chrisfen	1280	//the phi euler angle is for rotation about the z-axis (we use the zxz convention)
135			angles[0] = increment;
136	chrisfen	1277
137	chrisfen	1280	//obtain the rotation matrix through the rigid body class
138			rbMol->doEulerToRotMat(angles, rotZ);
139	chrisfen	1277
140	chrisfen	1280	//rotate the rigid body
141			rbMol->getA(rbMatrix);
142			matMul3(rotZ, rbMatrix, rotatedMat);
143			rbMol->setA(rotatedMat);
144	chrisfen	1277	}