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, 10 months ago) by chrisfen
File size:	3454 byte(s)
Log Message:	Fixed XCode auto-indent to match XEmacs
#	Content
1	#include "GridBuilder.hpp"
2	#include "MatVec3.h"
3	#define PI 3.14159265359
4
5
6	GridBuilder::GridBuilder(RigidBody* rb, int bandWidth) {
7	rbMol = rb;
8	bandwidth = bandWidth;
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	}
21	}
22
23	GridBuilder::~GridBuilder() {
24	}
25
26	void GridBuilder::launchProbe(int forceField, vector<double> sigmaGrid, vector<double> sGrid,
27	vector<double> epsGrid){
28	double startDist;
29	double minDist = 10.0; //minimum start distance
30
31	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
38	initBody();
39	for (i=0; i<bandwidth; i++){
40	for (j=0; j<bandwidth; j++){
41	releaseProbe(startDist);
42
43	sigmaGrid.push_back(sigDist);
44	sGrid.push_back(sDist);
45	epsGrid.push_back(epsVal);
46
47	stepPhi(phiStep);
48	}
49	stepTheta(thetaStep);
50	}
51	}
52
53	void GridBuilder::initBody(){
54	//set up the rigid body in the starting configuration
55	stepTheta(thetaMin);
56	}
57
58	void GridBuilder::releaseProbe(double farPos){
59	int tooClose;
60	double tempPotEnergy;
61	double interpRange;
62	double interpFrac;
63
64	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
71
72	while (!tooClose){
73	calcEnergy();
74	potProgress.push_back(potEnergy);
75	distProgress.push_back(probeCoor);
76
77	//if we've reached a new minimum, save the value and position
78	if (potEnergy < epsVal){
79	epsVal = potEnergy;
80	sDist = probeCoor;
81	}
82
83	//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
89	//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
102	//end the loop
103	tooClose = 1;
104	}
105	}
106	}
107
108	void GridBuilder::calcEnergy(){
109
110	}
111
112	void GridBuilder::stepTheta(double increment){
113	//zero out the euler angles
114	for (i=0; i<3; i++)
115	angles[i] = 0.0;
116
117	//the second euler angle is for rotation about the x-axis (we use the zxz convention)
118	angles[1] = increment;
119
120	//obtain the rotation matrix through the rigid body class
121	rbMol->doEulerToRotMat(angles, rotX);
122
123	//rotate the rigid body
124	rbMol->getA(rbMatrix);
125	matMul3(rotX, rbMatrix, rotatedMat);
126	rbMol->setA(rotatedMat);
127	}
128
129	void GridBuilder::stepPhi(double increment){
130	//zero out the euler angles
131	for (i=0; i<3; i++)
132	angles[i] = 0.0;
133
134	//the phi euler angle is for rotation about the z-axis (we use the zxz convention)
135	angles[0] = increment;
136
137	//obtain the rotation matrix through the rigid body class
138	rbMol->doEulerToRotMat(angles, rotZ);
139
140	//rotate the rigid body
141	rbMol->getA(rbMatrix);
142	matMul3(rotZ, rbMatrix, rotatedMat);
143	rbMol->setA(rotatedMat);
144	}