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