Null node when useparallel=True

[rebeccafair]

When I run the following script:

import spglib
import numpy as np
from euphonic import ureg, ForceConstants
import brille as br

fc_file = 'castep_files/La2Zr2O7.castep_bin'
fc = ForceConstants.from_castep(fc_file)

crystal = fc.crystal
cell_vectors = crystal._cell_vectors
cell = crystal.to_spglib_cell()

dataset = spglib.get_symmetry_dataset(cell)
rotations = dataset['rotations']  # in fractional
translations = dataset['translations']  # in fractional

symmetry = br.Symmetry(rotations, translations)
direct = br.Direct(*cell)
direct.spacegroup = symmetry
bz = br.BrillouinZone(direct.star)

vol = bz.ir_polyhedron.volume
n_grid_points=2000
grid = br.BZTrellisQdc(bz, node_volume_fraction=vol/n_grid_points)

phonons = fc.calculate_qpoint_phonon_modes(grid.rlu)

evecs_basis = np.einsum('ba,ijkb->ijka', np.linalg.inv(cell_vectors),
                        phonons.eigenvectors)
n_atoms = crystal.n_atoms
frequencies = np.reshape(phonons._frequencies,
phonons._frequencies.shape + (1,))
freq_el = (1,)
freq_weight = (1., 0., 0.)
evecs = np.reshape(evecs_basis,
                   (evecs_basis.shape[0], 3*n_atoms, 3*n_atoms))
cost_function = 0
evecs_el = (0, 3*n_atoms, 0, 3, 0, cost_function)
evecs_weight = (0., 1., 0.)
grid.fill(frequencies, freq_el, freq_weight, evecs, evecs_el,
          evecs_weight)

N = 500
qpts = np.zeros((N, 3))
qpts[:, 0] = 1.5
qpts[:, 1] = np.linspace(-1.5, 2.0, N)
qpts[:, 2] = np.linspace(4.0, -3.0, N)
vals, vecs = grid.ir_interpolate_at(qpts, useparallel=True, threads=4)

I get the following error:

The node subscript The node subscript The node subscript  24  0  1  0 10  0 for the point [[-0.4701,  0.2109, -0.1826]]
 is either invalid or points to a null node! for the point [[ 0.4995, -0.0441,  0.0382]]
 is either invalid or points to a null node! 12  0 21 for the point [[ 0.2491, -0.0259,  0.4367]]
 is either invalid or points to a null node!

This only seems to happen when useparallel=True and threads>=1 and is a bit intermittent. I am using Brille 0.5.4 from PyPI on Windows.

[g5t]

After some investigation, it seems that the points listed in the error message are

1. variable (perhaps random)
2. always clearly outside of the irreducible Brillouin zone polyhedron

To try and work out what the problem could be, and confirm that the points are outside of the irreducible Brillouin zone, I resorted to using bz.ir_moveinto and brille.plotting:

q1, _,_,_ = bz.ir_moveinto(qpts, 1)
q4, _,_,_ = bz.ir_moveinto(qpts, 4)

ax =bp.plot(grid.BrillouinZone, show=False, units='rlu');
bp.plot(q1, axs=ax, show=False);
bp.plot(q4, axs=ax)

Plotted separately, for 1 thread:

and 4 threads (though this is run dependent)

There are clearly points outside of the irreducible Brillouin zone, which points to a problem with bz.ir_moveinto.

[g5t]

The implementation of BrillouinZone::ir_moveinto includes an OpenMP parallel for loop which called BrillouinZone::isinside_wedge which also contains a parallel for loop. The nested parallel for loops caused ir_moveinto to fail silently in some cases.

Since the call to isinside_wedge from ir_moveinto can be refactored to only-ever check one point the inner loop is redudant and a new private method has been written to perform the same isinside_wedge calculation for a single point. This change is now in a pull request which should alleviate the problem.

[g5t]

This was fixed by #62 and is in the release v0.5.7 and on PyPI