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

Relocate q-points in BZ

This commit is contained in:
Atsushi Togo 2013-09-23 22:45:15 +09:00
parent 98994780af
commit 2aacaef782
6 changed files with 104 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
anharmonic/phonon3/conductivity_RTA.py
View File

@ -134,15 +134,20 @@ class conductivity_RTA:
assert self._grid_weights.sum() == np.prod(self._mesh /
self._mesh_divisors)
rec_lat = np.linalg.inv(self._primitive.get_cell())
look_at_grid_address_in_BZ(rec_lat,
self._mesh,
self._grid_address,
self._grid_points)
primitive_lattice = np.linalg.inv(self._primitive.get_cell())
multiplicity = spg.relocate_BZ_grid_address(self._grid_address,
self._mesh,
primitive_lattice)
# rec_lat = np.linalg.inv(self._primitive.get_cell())
# look_at_grid_address_in_BZ(rec_lat,
# self._mesh,
# self._grid_address,
# self._grid_points)
def get_qpoints(self):
qpoints = np.array([self._grid_address[gp].astype(float) / self._mesh
for gp in self._grid_points], dtype='double')
qpoints = np.array(self._grid_address[self._grid_points] /
self._mesh.astype('double'), dtype='double')
return qpoints
def get_grid_points(self):

3
anharmonic/phonon3/interaction.py
View File

@ -164,7 +164,8 @@ class Interaction:
triplets_at_q, weights_at_q, grid_address = get_triplets_at_q(
grid_point,
self._mesh,
self._symmetry.get_pointgroup_operations())
self._symmetry.get_pointgroup_operations(),
np.linalg.inv(self._primitive.get_cell()))
self._triplets_at_q = triplets_at_q
self._weights_at_q = weights_at_q

11
anharmonic/phonon3/triplets.py
View File

@ -4,6 +4,7 @@ import phonopy.structure.spglib as spg
def get_triplets_at_q(grid_point,
mesh,
point_group, # real space point group of space group
primitive_lattice, # column vectors
is_time_reversal=True):
weights, third_q, grid_address = spg.get_triplets_reciprocal_mesh_at_q(
grid_point,
@ -11,6 +12,10 @@ def get_triplets_at_q(grid_point,
point_group,
is_time_reversal)
multiplicity = spg.relocate_BZ_grid_address(grid_address,
mesh,
primitive_lattice)
triplets_at_q = spg.get_grid_triplets_at_q(
grid_point,
grid_address,
@ -18,7 +23,11 @@ def get_triplets_at_q(grid_point,
weights,
mesh)
weights_at_q = np.extract(weights > 0, weights)
grid_address = get_grid_address([x + (x % 2 == 0) for x in mesh])
extended_mesh = [x + (x % 2 == 0) for x in mesh]
grid_address = get_grid_address(extended_mesh)
multiplicity = spg.relocate_BZ_grid_address(grid_address,
extended_mesh,
primitive_lattice)
assert np.prod(mesh) == weights_at_q.sum(), \
"Num grid points %d, sum of weight %d" % (

49
c/spglib/kpoint.c
View File

@ -568,8 +568,8 @@ static void relocate_BZ_grid_address(int grid_address[][3],
for (i = 0; i < mesh[0] * mesh[1] * mesh[2]; i++) {
for (j = 0; j < 27; j++) {
for (k = 0; k < 3; k++) {
address[k] = grid_address[i][k];
address[k] += search_space[j][k] * mesh[k];
address[k] = grid_address[i][k] * 2;
address[k] += search_space[j][k] * mesh[k] * 2 + is_shift[k];
}
mat_multiply_matrix_vector_di3(vector, rec_lattice, address);
distance[j] = mat_norm_squared_d3(vector);
@ -718,16 +718,17 @@ static int get_ir_triplets_at_q(int weights[],
return num_ir_triplets;
}
/* Triplets: [grid_point, i, third_q[i]], weights[i] > 0 */
static void set_grid_triplets_at_q(int triplets[][3],
const int q_grid_point,
const int grid_point,
SPGCONST int grid_address[][3],
const int third_q[],
const int weights[],
const int mesh[3])
{
const int is_shift[3] = {0, 0, 0};
int i, j, k, num_edge, edge_pos, num_ir;
int grid_double[3][3], ex_mesh[3], ex_mesh_double[3];
int i, j, k, num_ir, smallest_index, smallest_g;
int grid_double[3][3], ex_mesh[3], ex_mesh_double[3], sum_g[27];
for (i = 0; i < 3; i++) {
ex_mesh[i] = mesh[i] + (mesh[i] % 2 == 0);
@ -735,7 +736,7 @@ static void set_grid_triplets_at_q(int triplets[][3],
}
for (i = 0; i < 3; i++) {
grid_double[0][i] = grid_address[q_grid_point][i] * 2;
grid_double[0][i] = grid_address[grid_point][i] * 2;
}
num_ir = 0;
@ -750,29 +751,27 @@ static void set_grid_triplets_at_q(int triplets[][3],
grid_double[2][j] = grid_address[third_q[i]][j] * 2;
}
for (j = 0; j < 3; j++) {
num_edge = 0;
edge_pos = -1;
for (j = 0; j < 27; j++) {
sum_g[j] = 0;
for (k = 0; k < 3; k++) {
if (abs(grid_double[k][j]) == mesh[j]) {
num_edge++;
edge_pos = k;
}
}
if (num_edge == 1) {
grid_double[edge_pos][j] = 0;
for (k = 0; k < 3; k++) {
if (k != edge_pos) {
grid_double[edge_pos][j] -= grid_double[k][j];
}
}
}
if (num_edge == 2) {
grid_double[edge_pos][j] = -grid_double[edge_pos][j];
sum_g[j] += abs(grid_double[0][k] + grid_double[1][k] +
grid_double[2][k] + search_space[j][k] * mesh[k] * 2);
}
}
smallest_index = 0;
smallest_g = sum_g[0];
for (j = 1; j < 27; j++) {
if (smallest_g > sum_g[j]) {
smallest_index = j;
smallest_g = sum_g[j];
}
}
for (j = 0; j < 3; j++) {
grid_double[2][j] += search_space[smallest_index][j] * mesh[j] * 2;
}
for (j = 0; j < 3; j++) {
get_vector_modulo(grid_double[j], ex_mesh_double);
triplets[num_ir][j] = get_grid_point(grid_double[j], ex_mesh, is_shift);

65
phonopy/structure/brillouin_zone.py
View File

@ -36,19 +36,6 @@ import numpy as np
from phonopy.structure.cells import get_reduced_bases
search_space = np.array([
[-1, -1, -1],
[-1, -1, 0],
[-1, -1, 1],
[-1, 0, -1],
[-1, 0, 0],
[-1, 0, 1],
[-1, 1, -1],
[-1, 1, 0],
[-1, 1, 1],
[0, -1, -1],
[0, -1, 0],
[0, -1, 1],
[0, 0, -1],
[0, 0, 0],
[0, 0, 1],
[0, 1, -1],
@ -62,7 +49,20 @@ search_space = np.array([
[1, 0, 1],
[1, 1, -1],
[1, 1, 0],
[1, 1, 1]], dtype='intc')
[1, 1, 1],
[-1, -1, -1],
[-1, -1, 0],
[-1, -1, 1],
[-1, 0, -1],
[-1, 0, 0],
[-1, 0, 1],
[-1, 1, -1],
[-1, 1, 0],
[-1, 1, 1],
[0, -1, -1],
[0, -1, 0],
[0, -1, 1],
[0, 0, -1]], dtype='intc')
def get_qpoints_in_Brillouin_zone(primitive_vectors,
qpoints,
@ -103,26 +103,47 @@ class BrillouinZone:
if __name__ == '__main__':
from phonopy.interface.vasp import read_vasp
from phonopy.structure.symmetry import Symmetry, get_lattice_vector_equivalence
from phonopy.structure.spglib import get_ir_reciprocal_mesh
from phonopy.structure.spglib import get_ir_reciprocal_mesh, relocate_BZ_grid_address
import sys
cell = read_vasp(sys.argv[1])
symmetry = Symmetry(cell)
mesh = [10, 10, 10]
mapping_table, grid_addrees = get_ir_reciprocal_mesh(
mesh = [4, 4, 4]
is_shift = np.array([1, 1, 1], dtype='intc')
mapping_table, grid_address = get_ir_reciprocal_mesh(
mesh,
cell,
is_shift=[0, 0, 0])
is_shift=is_shift)
ir_grid_points = np.unique(mapping_table)
qpoints = grid_addrees[ir_grid_points] / np.array(mesh, dtype='double')
grid_address_copy = grid_address.copy()
primitive_vectors = np.linalg.inv(cell.get_cell())
multiplicity = relocate_BZ_grid_address(grid_address_copy,
mesh,
np.linalg.inv(cell.get_cell()),
is_shift=is_shift)
# qpoints = grid_address[ir_grid_points] / np.array(mesh, dtype='double')
qpoints = (grid_address + is_shift.astype('double') / 2) / mesh
qpoints -= (qpoints > 0.5001) * 1
print multiplicity
for g, gc in zip(grid_address, grid_address_copy):
print g, gc,
if np.all(g == gc):
print
else:
print "*"
bz = BrillouinZone(primitive_vectors)
bz.run(qpoints)
sv = bz.get_qpoints()
sv = bz.get_shortest_qpoints()
for q, vs in zip(qpoints, sv):
print q
print q,
if np.allclose(q, vs[0]):
print
else:
print "*", np.linalg.norm(np.dot(primitive_vectors, q))
for v in vs:
print v, np.linalg.norm(np.dot(primitive_vectors, v))

13
phonopy/structure/spglib.py
View File

@ -244,6 +244,19 @@ def get_ir_reciprocal_mesh(mesh,
symprec)
return mapping, mesh_points
def relocate_BZ_grid_address(grid_address,
mesh,
reciprocal_lattice, # column vectors
is_shift=np.zeros(3, dtype='intc')):
multiplicity = np.zeros(len(grid_address), dtype='intc')
spg.BZ_grid_address(grid_address,
multiplicity,
np.array(mesh, dtype='intc').copy(),
np.array(reciprocal_lattice, dtype='double').copy(),
np.array(is_shift, dtype='intc').copy())
return multiplicity
def get_stabilized_reciprocal_mesh(mesh,
rotations,