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

Careful copy of numpy array; np.array(oldarray).copy() seems not safe for array with .T (transpose). So instead of .copy(), more explicit way is used: np.array(oldarray, dtype='some', order='C')

This commit is contained in:
Atsushi Togo 2014-01-26 00:22:43 +09:00
parent 2e7bb4f2cf
commit e587143ae0
12 changed files with 134 additions and 143 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
anharmonic/phonon3/__init__.py
View File

@ -16,7 +16,6 @@ class Phono3py:
primitive,
mesh,
fc3=None,
tetrahedron_method=False,
sigmas = [],
band_indices=None,
cutoff_frequency=1e-4,
@ -30,18 +29,12 @@ class Phono3py:
self._supercell = supercell
self._primitive = primitive
self._mesh = mesh
self._sigmas = sigmas
if band_indices is None:
self._band_indices = [
np.arange(primitive.get_number_of_atoms() * 3)]
else:
self._band_indices = band_indices
self._tetrahedron_method = tetrahedron_method
if tetrahedron_method:
self._sigmas = [None] + list(sigmas)
else:
self._sigmas = list(sigmas)
self._frequency_factor_to_THz = frequency_factor_to_THz
self._is_nosym = is_nosym
self._symmetrize_fc3_q = symmetrize_fc3_q
@ -270,7 +263,6 @@ class Phono3pyJointDos:
fc2,
nac_params=None,
sigmas=[],
tetrahedron_method=False,
frequency_step=None,
frequency_factor_to_THz=VaspToTHz,
frequency_scale_factor=None,
@ -284,7 +276,6 @@ class Phono3pyJointDos:
self._fc2 = fc2
self._nac_params = nac_params
self._sigmas = sigmas
self._tetrahedron_method = tetrahedron_method
self._frequency_step = frequency_step
self._frequency_factor_to_THz = frequency_factor_to_THz
self._frequency_scale_factor = frequency_scale_factor
@ -299,7 +290,6 @@ class Phono3pyJointDos:
self._supercell,
self._fc2,
nac_params=self._nac_params,
tetrahedron_method=self._tetrahedron_method,
frequency_step=self._frequency_step,
frequency_factor_to_THz=self._frequency_factor_to_THz,
frequency_scale_factor=self._frequency_scale_factor,
@ -325,17 +315,15 @@ class Phono3pyJointDos:
frequencies = self._jdos.get_phonons()[0]
print frequencies[gp]
if self._tetrahedron_method:
print "Tetrahedron method"
self._jdos.set_sigma(None)
self._jdos.run()
self._write(gp, sigma=None)
if self._sigmas:
for sigma in self._sigmas:
print "Sigma:", sigma
if sigma is None:
print "Tetrahedron method"
else:
print "Sigma:", sigma
self._jdos.set_sigma(sigma)
self._jdos.run()
self._write(gp, sigma)
self._write(gp, sigma=sigma)
else:
print "sigma or tetrahedron method has to be set."

7
anharmonic/phonon3/conductivity_RTA.py
View File

@ -168,7 +168,7 @@ def _set_gamma_from_file(br, filename=None):
print "Gamma at grid point %d doesn't exist." % gp
gamma_at_sigma.append(gamma_gp)
gamma.append(gamma_at_sigma)
br.set_gamma(np.array(gamma, dtype='double'))
br.set_gamma(np.array(gamma, dtype='double', order='C'))
class conductivity_RTA:
def __init__(self,
@ -293,7 +293,8 @@ class conductivity_RTA:
self._mesh_divisors)
self._qpoints = np.array(self._grid_address[self._grid_points] /
self._mesh.astype('double'), dtype='double')
self._mesh.astype('double'),
dtype='double', order='C')
self._sum_num_kstar = 0
self._grid_point_count = 0
@ -520,7 +521,7 @@ class conductivity_RTA:
self._sum_num_kstar += len(gv_by_gv)
return np.array(gv_by_gv, dtype='double')
return np.array(gv_by_gv, dtype='double', order='C')
def _get_gv_by_gv_on_star(self, group_velocity, rotation_map):
gv2_tensor = []

14
anharmonic/phonon3/imag_self_energy.py
View File

@ -532,14 +532,9 @@ class ImagSelfEnergy:
thm = TetrahedronMethod(reciprocal_lattice, mesh=self._mesh)
grid_address = self._interaction.get_grid_address()
bz_map = self._interaction.get_bz_map()
self._set_triplets_integration_weights_c(thm,
grid_address,
bz_map)
self._set_triplets_integration_weights_c(thm, grid_address, bz_map)
def _set_triplets_integration_weights_c(self,
thm,
grid_address,
bz_map):
def _set_triplets_integration_weights_c(self, thm, grid_address, bz_map):
import anharmonic._phono3py as phono3c
unique_vertices = thm.get_unique_tetrahedra_vertices()
for i, j in zip((1, 2), (1, -1)):
@ -574,10 +569,7 @@ class ImagSelfEnergy:
grid_address,
bz_map)
def _set_triplets_integration_weights_py(self,
thm,
grid_address,
bz_map):
def _set_triplets_integration_weights_py(self, thm, grid_address, bz_map):
tetrahedra_vertices = get_tetrahedra_vertices(thm.get_tetrahedra(),
self._mesh,
self._triplets_at_q,

8
anharmonic/phonon3/interaction.py
View File

@ -45,7 +45,7 @@ def set_phonon_c(dm,
svecs, multiplicity = dm.get_shortest_vectors()
masses = np.array(dm.get_primitive().get_masses(), dtype='double')
rec_lattice = np.array(
np.linalg.inv(dm.get_primitive().get_cell()), dtype='double')
np.linalg.inv(dm.get_primitive().get_cell()), dtype='double', order='C')
if dm.is_nac():
born = dm.get_born_effective_charges()
nac_factor = dm.get_nac_factor()
@ -113,13 +113,13 @@ class Interaction:
self._fc3 = fc3
self._supercell = supercell
self._primitive = primitive
self._mesh = np.array(mesh, dtype='intc').copy()
self._mesh = np.array(mesh, dtype='intc')
self._symmetry = symmetry
num_band = primitive.get_number_of_atoms() * 3
if band_indices is None:
self._band_indices = np.arange(num_band, dtype='intc')
else:
self._band_indices = np.array(band_indices, dtype='intc').copy()
self._band_indices = np.array(band_indices, dtype='intc')
self._frequency_factor_to_THz = frequency_factor_to_THz
if cutoff_frequency is None:
@ -281,7 +281,7 @@ class Interaction:
self._fc3,
svecs,
multiplicity,
np.array(masses, dtype='double'),
masses,
p2s,
s2p,
self._band_indices,

12
anharmonic/phonon3/joint_dos.py
View File

@ -15,8 +15,7 @@ class JointDos:
supercell,
fc2,
nac_params=None,
sigma=0.1,
tetrahedron_method=False,
sigma=None,
frequency_step=0.1,
frequency_factor_to_THz=VaspToTHz,
frequency_scale_factor=1.0,
@ -48,12 +47,7 @@ class JointDos:
self._set_dynamical_matrix()
self._symmetry = Symmetry(primitive, symprec)
if tetrahedron_method:
self._tetrahedron_method = TetrahedronMethod(
self._reciprocal_lattice, mesh=self._mesh)
else:
self._tetrahedron_method = None
self._tetrahedron_method = None
self._phonon_done = None
self._frequencies = None
self._eigenvectors = None
@ -112,6 +106,8 @@ class JointDos:
def _run_c(self, lang='Py'):
if self._sigma is None:
self._tetrahedron_method = TetrahedronMethod(
self._reciprocal_lattice, mesh=self._mesh)
if lang == 'C':
self._run_c_tetrahedron_method()
else:

24
anharmonic/phonon3/triplets.py
View File

@ -1,5 +1,6 @@
import numpy as np
import phonopy.structure.spglib as spg
from phonopy.structure.symmetry import Symmetry
def get_triplets_at_q(grid_point,
mesh,
@ -110,11 +111,11 @@ def reduce_grid_points(mesh_divisors,
dense_grid_points,
dense_grid_weights=None,
coarse_mesh_shifts=None):
divisors = np.array(mesh_divisors, dtype=int)
divisors = np.array(mesh_divisors, dtype='intc')
if (divisors == 1).all():
coarse_grid_points = np.array(dense_grid_points, dtype=int)
coarse_grid_points = np.array(dense_grid_points, dtype='intc')
if dense_grid_weights is not None:
coarse_grid_weights = np.array(dense_grid_weights, dtype=int)
coarse_grid_weights = np.array(dense_grid_weights, dtype='intc')
else:
grid_weights = []
if coarse_mesh_shifts is None:
@ -147,11 +148,14 @@ def from_coarse_to_dense_grid_points(dense_mesh,
def get_coarse_ir_grid_points(primitive,
mesh,
mesh_divs,
mesh_divisors,
coarse_mesh_shifts,
is_nosym=False):
if mesh_divs is None:
is_nosym=False,
symprec=1e-5):
if mesh_divisors is None:
mesh_divs = [1, 1, 1]
else:
mesh_divs = mesh_divisors
mesh = np.array(mesh, dtype='intc')
mesh_divs = np.array(mesh_divs, dtype='intc')
coarse_mesh = mesh / mesh_divs
@ -160,15 +164,15 @@ def get_coarse_ir_grid_points(primitive,
if is_nosym:
coarse_grid_address = get_grid_address(coarse_mesh)
coarse_grid_points = np.arange(np.prod(coarse_mesh),
dtype='intc')
coarse_grid_points = np.arange(np.prod(coarse_mesh), dtype='intc')
coarse_grid_weights = np.ones(len(coarse_grid_points), dtype='intc')
else:
symmetry = Symmetry(primitive, symprec)
(coarse_grid_points,
coarse_grid_weights,
coarse_grid_address) = get_ir_grid_points(
coarse_mesh,
primitive,
symmetry.get_pointgroup_operations(),
mesh_shifts=coarse_mesh_shifts)
grid_points = from_coarse_to_dense_grid_points(
mesh,
@ -211,7 +215,7 @@ search_space = np.array([
[0, -1, -1],
[0, -1, 0],
[0, -1, 1],
[0, 0, -1]], dtype='intc')
[0, 0, -1]], dtype='intc', order='C')
def get_grid_points_in_Brillouin_zone(primitive_vectors, # column vectors
mesh,

5
phonopy/harmonic/dynamical_matrix.py
View File

@ -247,7 +247,8 @@ class DynamicalMatrixNAC(DynamicalMatrix):
self._nac_params = nac_params
self._method = method
self._born = np.array(self._nac_params['born'], dtype='double')
self._born = np.array(self._nac_params['born'],
dtype='double', order='C')
factor = self._nac_params['factor']
if (isinstance(factor, list) or
isinstance(factor, tuple)):
@ -257,7 +258,7 @@ class DynamicalMatrixNAC(DynamicalMatrix):
self._unit_conversion = factor
self._damping_factor = DAMPING_FACTOR
self._dielectric = np.array(self._nac_params['dielectric'],
dtype='double')
dtype='double', order='C')
def set_dynamical_matrix(self, q_red, q_direction=None, verbose=False):
num_atom = self._pcell.get_number_of_atoms()

3
phonopy/phonon/thermal_properties.py
View File

@ -106,7 +106,8 @@ class ThermalProperties(ThermalPropertiesBase):
self._frequencies = np.where(frequencies > cutoff_frequency,
frequencies, -1)
else:
self._frequencies = np.array(frequencies, dtype='double') * THzToEv
self._frequencies = np.array(frequencies,
dtype='double', order='C') * THzToEv
self._set_high_T_entropy_and_zero_point_energy()
self._is_projection = is_projection

13
phonopy/structure/atoms.py
View File

@ -52,7 +52,7 @@ class Atoms:
if cell == None:
self.cell=None
else:
self.cell = np.array(cell, dtype=float)
self.cell = np.array(cell, dtype='double', order='C')
# position
self.scaled_positions = None
@ -68,7 +68,7 @@ class Atoms:
if numbers==None:
self.numbers = None
else:
self.numbers = np.array(numbers, dtype=int)
self.numbers = np.array(numbers, dtype='intc')
# masses
self.set_masses(masses)
@ -90,7 +90,7 @@ class Atoms:
def set_cell(self, cell):
self.cell = np.array(cell, dtype=float)
self.cell = np.array(cell, dtype='double', order='C')
def get_cell(self):
return self.cell.copy()
@ -103,7 +103,8 @@ class Atoms:
return np.dot(self.scaled_positions, self.cell)
def set_scaled_positions(self, scaled_positions):
self.scaled_positions = np.array(scaled_positions, dtype=float)
self.scaled_positions = np.array(scaled_positions,
dtype='double', order='C')
def get_scaled_positions(self):
return self.scaled_positions.copy()
@ -112,7 +113,7 @@ class Atoms:
if masses == None:
self.masses = None
else:
self.masses = np.array(masses, dtype=float)
self.masses = np.array(masses, dtype='double')
def get_masses(self):
return self.masses.copy()
@ -121,7 +122,7 @@ class Atoms:
if magmoms == None:
self.magmoms = None
else:
self.magmoms = np.array(magmoms, dtype=float)
self.magmoms = np.array(magmoms, dtype='double')
def get_magnetic_moments(self):
if self.magmoms == None:

140
phonopy/structure/spglib.py
View File

@ -15,15 +15,15 @@ def get_symmetry(bulk, use_magmoms=False, symprec=1e-5, angle_tolerance=-1.0):
"""
# Atomic positions have to be specified by scaled positions for spglib.
positions = bulk.get_scaled_positions().copy()
lattice = bulk.get_cell().T.copy()
numbers = np.intc(bulk.get_atomic_numbers()).copy()
positions = np.array(bulk.get_scaled_positions(), dtype='double', order='C')
lattice = np.array(bulk.get_cell().T, dtype='double', order='C')
numbers = np.array(bulk.get_atomic_numbers(), dtype='intc')
# Get number of symmetry operations and allocate symmetry operations
# multi = spg.multiplicity(cell, positions, numbers, symprec)
multi = 48 * bulk.get_number_of_atoms()
rotation = np.zeros((multi, 3, 3), dtype='intc')
translation = np.zeros((multi, 3))
translation = np.zeros((multi, 3), dtype='double')
# Get symmetry operations
if use_magmoms:
@ -45,8 +45,9 @@ def get_symmetry(bulk, use_magmoms=False, symprec=1e-5, angle_tolerance=-1.0):
symprec,
angle_tolerance)
return {'rotations': rotation[:num_sym].copy(),
'translations': translation[:num_sym].copy()}
return {'rotations': np.array(rotation[:num_sym], dtype='intc', order='C'),
'translations': np.array(translation[:num_sym],
dtype='double', order='C')}
def get_symmetry_dataset(bulk, symprec=1e-5, angle_tolerance=-1.0):
"""
@ -64,9 +65,10 @@ def get_symmetry_dataset(bulk, symprec=1e-5, angle_tolerance=-1.0):
wyckoffs:
Wyckoff letters
"""
positions = bulk.get_scaled_positions().copy()
lattice = bulk.get_cell().T.copy()
numbers = np.array(bulk.get_atomic_numbers(), dtype='intc').copy()
positions = np.array(bulk.get_scaled_positions(), dtype='double', order='C')
lattice = np.array(bulk.get_cell().T, dtype='double', order='C')
numbers = np.array(bulk.get_atomic_numbers(), dtype='intc')
keys = ('number',
'international',
'hall',
@ -87,12 +89,12 @@ def get_symmetry_dataset(bulk, symprec=1e-5, angle_tolerance=-1.0):
dataset['international'] = dataset['international'].strip()
dataset['hall'] = dataset['hall'].strip()
dataset['transformation_matrix'] = np.array(
dataset['transformation_matrix'], dtype='double')
dataset['origin_shift'] = np.array(dataset['origin_shift'],
dtype='double')
dataset['rotations'] = np.array(dataset['rotations'], dtype='intc')
dataset['transformation_matrix'], dtype='double', order='C')
dataset['origin_shift'] = np.array(dataset['origin_shift'], dtype='double')
dataset['rotations'] = np.array(dataset['rotations'],
dtype='intc', order='C')
dataset['translations'] = np.array(dataset['translations'],
dtype='double')
dtype='double', order='C')
letters = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ"
dataset['wyckoffs'] = [letters[x] for x in dataset['wyckoffs']]
dataset['equivalent_atoms'] = np.array(dataset['equivalent_atoms'],
@ -106,12 +108,12 @@ def get_spacegroup(bulk, symprec=1e-5, angle_tolerance=-1.0):
as a string.
"""
# Atomic positions have to be specified by scaled positions for spglib.
return spg.spacegroup(bulk.get_cell().T.copy(),
bulk.get_scaled_positions().copy(),
np.array(bulk.get_atomic_numbers(),
dtype='intc').copy(),
symprec,
angle_tolerance)
return spg.spacegroup(
np.array(bulk.get_cell().T, dtype='double', order='C'),
np.array(bulk.get_scaled_positions(), dtype='double', order='C'),
np.array(bulk.get_atomic_numbers(), dtype='intc'),
symprec,
angle_tolerance)
def get_pointgroup(rotations):
"""
@ -152,7 +154,7 @@ def get_pointgroup(rotations):
"""
# (symbol, pointgroup_number, transformation_matrix)
return spg.pointgroup(np.array(rotations, dtype='intc').copy())
return spg.pointgroup(np.array(rotations, dtype='intc', order='C'))
def refine_cell(bulk, symprec=1e-5, angle_tolerance=-1.0):
"""
@ -160,7 +162,7 @@ def refine_cell(bulk, symprec=1e-5, angle_tolerance=-1.0):
"""
# Atomic positions have to be specified by scaled positions for spglib.
num_atom = bulk.get_number_of_atoms()
lattice = bulk.get_cell().T.copy()
lattice = np.array(bulk.get_cell().T, dtype='double', order='C')
pos = np.zeros((num_atom * 4, 3), dtype='double')
pos[:num_atom] = bulk.get_scaled_positions()
@ -173,9 +175,9 @@ def refine_cell(bulk, symprec=1e-5, angle_tolerance=-1.0):
symprec,
angle_tolerance)
return (lattice.T.copy(),
pos[:num_atom_bravais].copy(),
numbers[:num_atom_bravais].copy())
return (np.array(lattice.T, dtype='double', order='C'),
np.array(pos[:num_atom_bravais], dtype='double', order='C'),
np.array(numbers[:num_atom_bravais], dtype='intc'))
def find_primitive(bulk, symprec=1e-5, angle_tolerance=-1.0):
"""
@ -185,9 +187,9 @@ def find_primitive(bulk, symprec=1e-5, angle_tolerance=-1.0):
"""
# Atomic positions have to be specified by scaled positions for spglib.
positions = bulk.get_scaled_positions().copy()
lattice = bulk.get_cell().T.copy()
numbers = np.array(bulk.get_atomic_numbers(), dtype='intc').copy()
positions = np.array(bulk.get_scaled_positions(), dtype='double', order='C')
lattice = np.array(bulk.get_cell().T, dtype='double', order='C')
numbers = np.array(bulk.get_atomic_numbers(), dtype='intc')
# lattice is transposed with respect to the definition of Atoms class
num_atom_prim = spg.primitive(lattice,
@ -196,9 +198,9 @@ def find_primitive(bulk, symprec=1e-5, angle_tolerance=-1.0):
symprec,
angle_tolerance)
if num_atom_prim > 0:
return (lattice.T.copy(),
positions[:num_atom_prim].copy(),
numbers[:num_atom_prim].copy())
return (np.array(lattice.T, dtype='double', order='C'),
np.array(positions[:num_atom_prim], dtype='double', order='C'),
np.array(numbers[:num_atom_prim], dtype='intc'))
else:
return None, None, None
@ -215,16 +217,16 @@ def get_ir_reciprocal_mesh(mesh,
mapping = np.zeros(np.prod(mesh), dtype='intc')
mesh_points = np.zeros((np.prod(mesh), 3), dtype='intc')
spg.ir_reciprocal_mesh(mesh_points,
mapping,
np.array(mesh, dtype='intc').copy(),
np.array(is_shift, dtype='intc').copy(),
is_time_reversal * 1,
bulk.get_cell().T.copy(),
bulk.get_scaled_positions().copy(),
np.array(bulk.get_atomic_numbers(),
dtype='intc').copy(),
symprec)
spg.ir_reciprocal_mesh(
mesh_points,
mapping,
np.array(mesh, dtype='intc'),
np.array(is_shift, dtype='intc'),
is_time_reversal * 1,
np.array(bulk.get_cell().T, dtype='double', order='C'),
np.array(bulk.get_scaled_positions(), dtype='double', order='C'),
np.array(bulk.get_atomic_numbers(), dtype='intc'),
symprec)
return mapping, mesh_points
@ -265,9 +267,9 @@ def relocate_BZ_grid_address(grid_address,
bz_grid_address,
bz_map,
grid_address,
np.array(mesh, dtype='intc').copy(),
np.array(reciprocal_lattice, dtype='double').copy(),
np.array(is_shift, dtype='intc').copy())
np.array(mesh, dtype='intc'),
np.array(reciprocal_lattice, dtype='double', order='C'),
np.array(is_shift, dtype='intc'))
return bz_grid_address[:num_bz_ir], bz_map
@ -287,18 +289,20 @@ def get_stabilized_reciprocal_mesh(mesh,
mapping = np.zeros(np.prod(mesh), dtype='intc')
mesh_points = np.zeros((np.prod(mesh), 3), dtype='intc')
qpoints = np.array(qpoints, dtype='double').copy()
qpoints = np.array(qpoints, dtype='double', order='C')
if qpoints.shape == (3,):
qpoints = np.array([qpoints], dtype='double')
qpoints = np.array([qpoints], dtype='double', order='C')
if qpoints.shape == (0,):
qpoints = np.array([[0, 0, 0]], dtype='double')
spg.stabilized_reciprocal_mesh(mesh_points,
mapping,
np.array(mesh, dtype='intc').copy(),
np.array(is_shift, dtype='intc'),
is_time_reversal * 1,
np.array(rotations, dtype='intc').copy(),
np.array(qpoints, dtype='double'))
qpoints = np.array([[0, 0, 0]], dtype='double', order='C')
spg.stabilized_reciprocal_mesh(
mesh_points,
mapping,
np.array(mesh, dtype='intc'),
np.array(is_shift, dtype='intc'),
is_time_reversal * 1,
np.array(rotations, dtype='intc', order='C'),
qpoints)
return mapping, mesh_points
@ -311,13 +315,14 @@ def get_triplets_reciprocal_mesh_at_q(fixed_grid_number,
third_q = np.zeros(np.prod(mesh), dtype='intc')
mesh_points = np.zeros((np.prod(mesh), 3), dtype='intc')
spg.triplets_reciprocal_mesh_at_q(weights,
mesh_points,
third_q,
fixed_grid_number,
np.array(mesh, dtype='intc').copy(),
is_time_reversal * 1,
np.array(rotations, dtype='intc').copy())
spg.triplets_reciprocal_mesh_at_q(
weights,
mesh_points,
third_q,
fixed_grid_number,
np.array(mesh, dtype='intc'),
is_time_reversal * 1,
np.array(rotations, dtype='intc', order='C'))
return weights, third_q, mesh_points
@ -334,7 +339,7 @@ def get_BZ_triplets_at_q(grid_point,
bz_grid_address,
bz_map,
weights,
np.array(mesh, dtype='intc').copy())
np.array(mesh, dtype='intc'))
return triplets
def get_neighboring_grid_points(grid_point,
@ -364,7 +369,7 @@ def get_triplets_tetrahedra_vertices(relative_grid_address,
spg.triplet_tetrahedra_vertices(
vertices_at_tp,
relative_grid_address,
np.array(mesh, dtype='intc').copy(),
np.array(mesh, dtype='intc'),
tp,
bz_grid_address,
bz_map)
@ -382,8 +387,9 @@ def get_tetrahedra_relative_grid_address(microzone_lattice):
relative_grid_address = np.zeros((24, 4, 3), dtype='intc')
spg.tetrahedra_relative_grid_address(
relative_grid_address, np.array(microzone_lattice,
dtype='double').copy())
relative_grid_address,
np.array(microzone_lattice, dtype='double', order='C'))
return relative_grid_address
@ -393,14 +399,14 @@ def get_tetrahedra_integration_weight(omegas,
if isinstance(omegas, float):
return spg.tetrahedra_integration_weight(
omegas,
np.array(tetrahedra_omegas, dtype='double').copy(),
np.array(tetrahedra_omegas, dtype='double', order='C'),
function)
else:
integration_weights = np.zeros(len(omegas), dtype='double')
spg.tetrahedra_integration_weight_at_omegas(
integration_weights,
np.array(omegas, dtype='double'),
np.array(tetrahedra_omegas, dtype='double').copy(),
np.array(tetrahedra_omegas, dtype='double', order='C'),
function)
return integration_weights

6
phonopy/structure/tetrahedron_method.py
View File

@ -42,14 +42,14 @@ parallelepiped_vertices = np.array([[0, 0, 0],
[0, 0, 1],
[1, 0, 1],
[0, 1, 1],
[1, 1, 1]], dtype='intc')
[1, 1, 1]], dtype='intc', order='C')
class TetrahedronMethod:
def __init__(self,
primitive_vectors,
mesh=[1, 1, 1]):
self._primitive_vectors = np.array(
primitive_vectors, dtype='double') / mesh # column vectors
primitive_vectors, dtype='double', order='C') / mesh # column vectors
self._vertices = None
self._relative_grid_addresses = None
self._central_indices = None
@ -83,7 +83,7 @@ class TetrahedronMethod:
break
if not found:
unique_vertices.append(adrs)
return np.array(unique_vertices, dtype='intc')
return np.array(unique_vertices, dtype='intc', order='C')
def set_tetrahedra_omegas(self, tetrahedra_omegas):
"""

21
scripts/phono3py
View File

@ -565,11 +565,12 @@ if options.write_grid_points:
(grid_points,
coarse_grid_weights,
grid_address) = get_coarse_ir_grid_points(
primitive,
mesh,
mesh_divs,
settings.get_coarse_mesh_shifts(),
is_nosym=options.no_kappa_stars)
primitive,
mesh,
mesh_divs,
settings.get_coarse_mesh_shifts(),
is_nosym=options.no_kappa_stars,
symprec=options.symprec)
write_ir_grid_points(mesh,
mesh_divs,
grid_points,
@ -845,6 +846,9 @@ elif settings.get_mesh_numbers() is not None:
else:
multiple_sigmas = settings.get_multiple_sigmas()
if settings.get_is_tetrahedron_method():
multiple_sigmas = [None] + multiple_sigmas
if settings.get_frequency_pitch() is None:
freq_step = 0.1
else:
@ -889,7 +893,6 @@ elif settings.get_mesh_numbers() is not None:
fc2,
nac_params=nac_params,
sigmas=multiple_sigmas,
tetrahedron_method=settings.get_is_tetrahedron_method(),
frequency_step=freq_step,
frequency_factor_to_THz=factor,
frequency_scale_factor=freq_scale,
@ -918,9 +921,8 @@ elif settings.get_mesh_numbers() is not None:
if settings.get_band_indices() is None:
band_indices = None
else:
band_indices = np.intc([x for bi in settings.get_band_indices()
for x in bi])
for gp in grid_points:
band_indices = np.hstack(self._band_indices).astype('intc')
for sigma in multiple_sigmas:
if log_level:
print "Sigma:", sigma
@ -935,7 +937,6 @@ elif settings.get_mesh_numbers() is not None:
primitive,
mesh,
fc3=fc3,
tetrahedron_method=settings.get_is_tetrahedron_method(),
sigmas=multiple_sigmas,
band_indices=settings.get_band_indices(),
cutoff_frequency=settings.get_cutoff_frequency(),