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;
        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:	1279
Committed:	Fri Jun 18 19:01:40 2004 UTC (20 years, 3 months ago) by chrisfen
File size:	3236 byte(s)
Log Message:	Built principle axis determination and molecule rotation into calcRefCoords in RigidBody.cpp
#	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			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	chrisfen	1279
41	chrisfen	1278	sigmaGrid.push_back(sigDist);
42			sGrid.push_back(sDist);
43			epsGrid.push_back(epsVal);
44
45	chrisfen	1277	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	chrisfen	1279	potProgress.clear();
65			distProgress.clear();
66	chrisfen	1277	tooClose = 0;
67			epsVal = 0;
68			rhoStep = 0.1; //the distance the probe atom moves between steps
69
70	chrisfen	1279
71	chrisfen	1277	while (!tooClose){
72			calcEnergy();
73			potProgress.push_back(potEnergy);
74			distProgress.push_back(probeCoor);
75
76			//if we've reached a new minimum, save the value and position
77			if (potEnergy < epsVal){
78			epsVal = potEnergy;
79			sDist = probeCoor;
80			}
81
82			//test if the probe reached the origin - if so, stop stepping closer
83			if (probeCoor < 0){
84			sigDist = 0.0;
85			tooClose = 1;
86			}
87
88			//test if the probe beyond the contact point - if not, take a step closer
89			if (potEnergy < 0){
90			sigDist = probeCoor;
91			tempPotEnergy = potEnergy;
92			probeCoor -= rhoStep;
93			}
94			else {
95			//do a linear interpolation to obtain the sigDist
96			interpRange = potEnergy - tempPotEnergy;
97			interpFrac = potEnergy / interpRange;
98			interpFrac = interpFrac * rhoStep;
99			sigDist = probeCoor + interpFrac;
100
101			//end the loop
102			tooClose = 1;
103			}
104			}
105			}
106
107			void GridBuilder::calcEnergy(){
108
109			}
110
111			void GridBuilder::stepTheta(double increment){
112			//zero out the euler angles
113			for (i=0; i<3; i++)
114			angles[i] = 0.0;
115
116			//the second euler angle is for rotation about the x-axis (we use the zxz convention)
117			angles[1] = increment;
118
119			//obtain the rotation matrix through the rigid body class
120			rbMol->doEulerToRotMat(angles, rotX);
121
122			//rotate the rigid body
123			rbMol->getA(rbMatrix);
124			matMul3(rotX, rbMatrix, rotatedMat);
125			rbMol->setA(rotatedMat);
126
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
145			}