scnc
/
phonopy
mirror of https://github.com/phonopy/phonopy.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Resolve some bottle necks for finding triplets

This commit is contained in:
Atsushi Togo 2013-04-24 23:57:27 +09:00
parent 2e749b0e7b
commit 0af157404b
3 changed files with 27 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
anharmonic/BTE_RTA.py
View File

@ -82,8 +82,10 @@ class BTE_RTA:
self._grid_address) = spg.get_ir_reciprocal_mesh(self._mesh,
self._primitive)
dense_grid_points = np.unique(grid_mapping_table)
dense_grid_weights = [np.sum(grid_mapping_table == g)
for g in dense_grid_points]
weights = np.zeros_like(grid_mapping_table)
for g in grid_mapping_table:
weights[g] += 1
dense_grid_weights = weights[dense_grid_points]
self._grid_points, self._grid_weights = reduce_grid_points(
self._mesh_divisors,

12
anharmonic/triplets.py
View File

@ -44,9 +44,9 @@ def get_grid_address(mesh):
for i in range(mesh[2]):
for j in range(mesh[1]):
for k in range(mesh[0]):
grid_address[count] = [k - (k > (mesh[2] // 2)) * mesh[0],
grid_address[count] = [k - (k > (mesh[0] // 2)) * mesh[0],
j - (j > (mesh[1] // 2)) * mesh[1],
i - (i > (mesh[0] // 2)) * mesh[2]]
i - (i > (mesh[2] // 2)) * mesh[2]]
count += 1
@ -77,9 +77,11 @@ def reduce_grid_points(mesh_divisors,
else:
grid_points = []
grid_weights = []
for i, dgp in enumerate(dense_grid_points):
if (grid_address[dgp] % divisors == 0).all():
grid_points.append(dgp)
for i, modulo in enumerate(
grid_address[dense_grid_points] % divisors):
if (modulo == 0).all():
grid_points.append(dense_grid_points[i])
if dense_grid_weights is not None:
grid_weights.append(dense_grid_weights[i])

31
c/spglib/kpoint.c
View File

@ -780,10 +780,10 @@ static int get_ir_triplets_at_q(int weights[],
const int mesh[3],
SPGCONST PointSymmetry * pointgroup)
{
int i, j, num_grid, weight_q, q_2, num_ir_q, num_ir_triplets;
int i, j, num_grid, q_2, num_ir_q, num_ir_triplets, ir_address;
int mesh_double[3], is_shift[3];
int grid_double0[3], grid_double1[3], grid_double2[3];
int *map_q, *ir_addresses;
int *map_q, *ir_addresses, *weight_q;
double tolerance;
double stabilizer_q[1][3];
PointSymmetry pointgroup_q;
@ -825,14 +825,21 @@ static int get_ir_triplets_at_q(int weights[],
#endif
ir_addresses = (int*) malloc(sizeof(int) * num_ir_q);
weight_q = (int*) malloc(sizeof(int) * num_grid);
num_ir_q = 0;
for (i = 0; i < num_grid; i++) {
if (map_q[i] == i) {
ir_addresses[num_ir_q] = i;
num_ir_q++;
}
weights[i] = 0;
weight_q[i] = 0;
}
#pragma omp parallel for
for (i = 0; i < num_grid; i++) {
weight_q[map_q[i]]++;
third_q[i] = -1;
weights[i] = 0;
}
#pragma omp parallel for private(j, grid_double1, grid_double2)
@ -847,26 +854,20 @@ static int get_ir_triplets_at_q(int weights[],
num_ir_triplets = 0;
for (i = 0; i < num_ir_q; i++) {
q_2 = third_q[ir_addresses[i]];
weight_q = 0;
#pragma omp parallel for reduction(+:weight_q)
for (j = 0; j < num_grid; j++) {
if (ir_addresses[i] == map_q[j]) {
weight_q++;
}
}
ir_address = ir_addresses[i];
q_2 = third_q[ir_address];
if (weights[map_q[q_2]]) {
weights[map_q[q_2]] += weight_q;
weights[map_q[q_2]] += weight_q[ir_address];
} else {
weights[ir_addresses[i]] = weight_q;
weights[ir_address] = weight_q[ir_address];
num_ir_triplets++;
}
}
free(map_q);
map_q = NULL;
free(weight_q);
weight_q = NULL;
free(ir_addresses);
ir_addresses = NULL;