-
Notifications
You must be signed in to change notification settings - Fork 4
/
gpu_energy.cc
257 lines (215 loc) · 9.42 KB
/
gpu_energy.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <cfloat>
#include <vector>
#include "heap.h"
#include "gpu_energy.h"
//#include "network.h"
using namespace std;
/* Construct a unit cell with the provided vector components*/
UNIT_CELL::UNIT_CELL(double vax, double vbx, double vby, double vcx, double vcy, double vcz){
va_x = vax;
vb_x = vbx; vb_y = vby;
vc_x = vcx; vc_y = vcy; vc_z = vcz;
double invDet = 1/(va_x*vb_y*vc_z);
inv_va_x = 1/va_x;
inv_vb_x = invDet*-1*vc_z*vb_x; inv_vb_y = 1/vb_y;
inv_vc_x = invDet* (vc_y*vb_x-vb_y*vc_x); inv_vc_y = invDet*-1*vc_y*va_x; inv_vc_z = 1/vc_z;
}
/** Transform the coordinates relative to the uc vectors into ones relative to
* xyz coordinate system and store them using the provided pointers. */
void UNIT_CELL::abc_to_xyz(double a, double b, double c, double &x, double &y, double &z){
x = a*va_x + b*vb_x + c*vc_x;
y = b*vb_y + c*vc_y;
z = c*vc_z;
}
/** Transform the coordinates relative to the xyz coordinate system vectors into ones relative to
* the unit cell vectors and store them using the provided pointers. */
void UNIT_CELL::xyz_to_abc(double x, double y, double z, double &a, double &b, double &c){
a = x*inv_va_x + y*inv_vb_x + z*inv_vc_x;
b = y*inv_vb_y + z*inv_vc_y;
c = z*inv_vc_z;
}
// Constructor that creates a dummy container
DISP_INFO::DISP_INFO(){
isReal = false;
}
/** Create a container with the provided displacement components and energy.
* Displacements are stored as characters to conserve space. */
DISP_INFO::DISP_INFO(int myA, int myB, int myC, float myMaxEnergy){
isReal = true;
a = (char)myA;
b = (char)myB;
c = (char)myC;
maxEnergy = myMaxEnergy;
}
// Returns true iff the displacements are equal
bool DISP_INFO::equalDisplacement(DISP_INFO other){
return (a == other.a) && (b == other.b) && (c == other.c);
}
/* Construct a TRIPLET with the three provided components.*/
TRIPLET::TRIPLET(int myX, int myY, int myZ){
vals[0] = x = myX;
vals[1] = y = myY;
vals[2] = z = myZ;
}
/* Access one of the TRIPLETs three values.*/
int& TRIPLET::operator[](int index){
if(index < 0 || index > 2){
cerr << "Error: Invalid index to [] operator for TRIPLET instance" << "\n"
<< "Exiting..." << "\n";
exit(1);
}
return vals[index];
}
/** Add each component of the triplet to that of the provided TRIPLET
* and return the result. */
TRIPLET TRIPLET::add(TRIPLET other){
return TRIPLET(x + other.x, y + other.y, z + other.z);
}
/* Return the 1-d index for the provided grid indices*/
int transformIndex(int x, int y, int z, int numX, int numY){
return x + y*numX + z*numX*numY;
}
/* Returns the energy of the point located at the provided grid indices.
* Assumes the grid indices are within the appropriate range.*/
float getEnergy(int x, int y, int z, int numX, int numY, float *energyGrid){
return energyGrid[transformIndex(x, y, z, numX, numY)];
}
/* Returns the energy of the point located in the grid referred to by the TRIPLET of indices.
* Assumes the grid indices are within the appropriate range.*/
float getEnergy(TRIPLET indices, int numX, int numY, float *energyGrid){
return getEnergy(indices[0], indices[1], indices[2], numX, numY, energyGrid);
}
/* Returns true iff the grid point referred to by the provided grid indices
* is indeed accessible. Assumes the grid indices are within the appropriate range. */
bool isAccessible(int x, int y, int z, int numX, int numY, float *accessGrid){
cerr << "Error: Accessibility check not implemented" << "\n";
exit(1);
return true;
}
/* Returns the integer nearest to the provided double.*/
int nearestInt(double num){
return (int)(floor(num + 0.5));
}
/* Translate the coordinate by unit cell increments so that it lies within the 0 to 1 range.*/
double translate_to_original_uc(double x){
double newX = x - (x < 0.0 ?(-(int)(0.5-x)):((int)(0.5+x)));
if(newX < 0)
return newX+1;
else
return newX;
}
/* Adjusts the indices so that they lie within the appropriate range
* and such that they are adjusted according to periodic boundary conditions.
* Stores a TRIPLET representing to which unit cell the indices originally referred.*/
void adjustIndices(TRIPLET &gridIndices, TRIPLET &shift, UNIT_CELL &uc){
double a, b, c;
uc.xyz_to_abc(gridIndices[0], gridIndices[1], gridIndices[2], a, b, c);
double newA, newB, newC;
newA = translate_to_original_uc(a); newB = translate_to_original_uc(b); newC = translate_to_original_uc(c);
double x, y, z;
uc.abc_to_xyz(newA, newB, newC, x, y, z);
shift = TRIPLET(nearestInt(a - newA), nearestInt(b - newB), nearestInt(c - newC));
gridIndices = TRIPLET(nearestInt(x), nearestInt(y), nearestInt(z));
}
/* Returns true iff the grid point referred to by the provided triplet of grid indices
* is indeed accessible. Assumes the grid indices are within the appropriate range. */
bool isAccessible(TRIPLET indices, int numX, int numY, float *accessGrid){
return isAccessible(indices.x, indices.y, indices.z, numX, numY, accessGrid);
}
/* Returns true if the first pair has energy >= that of the second pair.*/
bool hasHigherEnergy(pair<TRIPLET, DISP_INFO > p1, pair<TRIPLET, DISP_INFO > p2){
return p1.second.maxEnergy >= p2.second.maxEnergy;
}
// Helper function used in calculateMinEnergyBarrier()
bool findMinEnergyBarrier(int startX, int startY, int startZ, double &barrier,
int numX, int numY, int numZ, float *energyGrid, float *accessGrid, UNIT_CELL unit_cell){
barrier = DBL_MAX;
vector<DISP_INFO> visitedNodes = vector<DISP_INFO> (numX*numY*numZ, DISP_INFO());
// Place starting node on heap with (0,0,0) displacement
HEAP<pair<TRIPLET, DISP_INFO > > heap (hasHigherEnergy); // Lowest energy paths are favored
heap.insert(pair<TRIPLET, DISP_INFO> (TRIPLET(startX, startY, startZ), DISP_INFO(0, 0, 0, getEnergy(startX, startY, startZ, numX, numY, energyGrid))));
while(heap.size() != 0){
// Remove best remaining path
pair<TRIPLET, DISP_INFO> pathInfo = heap.pop();
int arrayIndex = transformIndex(pathInfo.first[0], pathInfo.first[1], pathInfo.first[2], numX, numY);
DISP_INFO nodeInfo = visitedNodes[arrayIndex];
if(nodeInfo.isReal){
if(nodeInfo.equalDisplacement(pathInfo.second)){
// Circling back to previous node so terminate path
continue;
}
else {
// Found the optimal path
barrier = min(pathInfo.second.maxEnergy, nodeInfo.maxEnergy);
return true;
}
}
else{
// Add all connected neigbors if appropriate
visitedNodes[arrayIndex] = pathInfo.second;
for(int i = 0; i < NUM_DIRECTIONS; i++){
TRIPLET dir = DIRECTIONS[i];
TRIPLET newIndices = dir.add(pathInfo.first);
TRIPLET change_in_uc(0,0,0);
adjustIndices(newIndices, change_in_uc, unit_cell);
if(isAccessible(newIndices, numX, numY, accessGrid)){
int da = change_in_uc[0] + pathInfo.second.a;
int db = change_in_uc[1] + pathInfo.second.b;
int dc = change_in_uc[2] + pathInfo.second.c;
heap.insert(pair<TRIPLET, DISP_INFO> (newIndices, DISP_INFO(da, db, dc, max(pathInfo.second.maxEnergy, getEnergy(newIndices, numX, numY, energyGrid)))));
}
}
}
}
return false;
}
/* Calculates the minimum energy barrier of a path that must pass through the point referred to by the provided
* grid indices. The resulting barrier is stored using the provided pointer. If no path is possible,
* the barrier is set to DBL_MAX. The function returns true iff a path is found. */
bool calculateMinEnergyBarrier(double &minEnergy, double &barrierEnergy, float box_x, float box_y, float box_z, int numX, int numY, int numZ, float *energyGrid, float *accessGrid,
double vax, double vbx, double vby, double vcx, double vcy, double vcz){
cout << "Identifying global accessible energy minimum" << "\n";
// Find global energy minimum and its corresponding grid indices
minEnergy = DBL_MAX;
int minZ = -1, minY = -1, minX = -1;
for(int i = 0; i < numZ; i++){
for(int j = 0; j < numY; j++){
for(int k = 0; k < numX; k++){
if(isAccessible(k, j, i, numX, numY, accessGrid)){
float energy = getEnergy(k, j, i, numX, numY, energyGrid);
if(energy < minEnergy){
minEnergy = energy;
minZ = i; minY = j; minX = k;
}
}
}
}
}
// Ensure that a minimum was found
if((minX == -1) || (minY == -1) || (minZ == -1)){
cerr << "Error: Unable to identify global energy minimum when calculating energy barriers." << "\n"
<< "Aborting energy calculation" << "\n";
return false;
}
cout << "Global accessible energy minimum found" << "\n";
// Construct the unit cell vectors in terms of grid points
double grid_spacing_x = box_x/numX, grid_spacing_y = box_y/numY, grid_spacing_z = box_z/numZ;
vax /= grid_spacing_x;
vbx /= grid_spacing_x; vby /= grid_spacing_y;
vcx /= grid_spacing_x; vcy /= grid_spacing_y; vcz /= grid_spacing_z;
UNIT_CELL unit_cell(vax, vbx, vby, vcx, vcy, vcz);
cout << "Calculating energy barrier for box of dimensions:" << box_x << " by " << box_y << " by " << box_z << "\n"
<< "Grid size: " << numX << " by " << numY << " by " << numZ << "\n"
<< "Grid unitcell vectors:"
<< "\t v_a:(" << vax << ", 0, 0)" << "\n"
<< "\t v_b:(" << vbx << ", " << vby << ", 0)" << "\n"
<< "\t v_c:{" << vcx << ", " << vcy << ", " << vcz << ")" << "\n"
<< "\n" << "\n";
// Calculate the barrier
bool result = findMinEnergyBarrier(minX, minY, minZ, barrierEnergy, numX, numY, numZ, energyGrid, accessGrid, unit_cell);
cout << "Energy barrier calculation completed" << "\n" << "\n";
return result;
}