Merge branch 'release-1.10.8'

anharmonic/file_IO.py

@@ -637,6 +637,7 @@ def write_kappa_to_hdf5(temperature,
  grid_point=None,
  band_index=None,
  sigma=None,
+ kappa_unit_conversion=None,
  filename=None,
  verbose=True):
  if band_index is None:
@@ -671,6 +672,9 @@ def write_kappa_to_hdf5(temperature,
  w.create_dataset('qpoint', data=qpoint)
  if weight is not None:
  w.create_dataset('weight', data=weight)
+ if kappa_unit_conversion is not None:
+ w.create_dataset('kappa_unit_conversion',
+ data=kappa_unit_conversion)
 
  if verbose:
  text = ""

anharmonic/phonon3/conductivity_RTA.py

@@ -4,7 +4,7 @@
 from phonopy.phonon.thermal_properties import mode_cv as get_mode_cv
 from anharmonic.file_IO import (write_kappa_to_hdf5, write_triplets,
  read_gamma_from_hdf5, write_grid_address)
-from anharmonic.phonon3.conductivity import Conductivity
+from anharmonic.phonon3.conductivity import Conductivity, unit_to_WmK
 from anharmonic.phonon3.imag_self_energy import ImagSelfEnergy
 from anharmonic.phonon3.triplets import get_grid_points_by_rotations
 
@@ -58,7 +58,7 @@ def get_thermal_conductivity_RTA(
  if not _set_gamma_from_file(br, filename=input_filename):
  print("Reading collisions failed.")
  return False
- 
+
  for i in br:
  if write_gamma:
  _write_gamma(br,
@@ -72,7 +72,10 @@ def get_thermal_conductivity_RTA(
  if write_kappa:
  if (grid_points is None and _all_bands_exist(interaction)):
  br.set_kappa_at_sigmas()
- _write_kappa(br, filename=output_filename, log_level=log_level)
+ _write_kappa(br,
+ interaction,
+ filename=output_filename,
+ log_level=log_level)
 
  return br
 
@@ -87,9 +90,14 @@ def _write_gamma(br, interaction, i, filename=None, verbose=True):
  gamma = br.get_gamma()
  gamma_isotope = br.get_gamma_isotope()
  sigmas = br.get_sigmas()
+ volume = interaction.get_primitive().get_volume()
 
  gp = grid_points[i]
  if _all_bands_exist(interaction):
+ if ave_pp is None:
+ ave_pp_i = None
+ else:
+ ave_pp_i = ave_pp[i]
  frequencies = interaction.get_phonons()[0][gp]
  for j, sigma in enumerate(sigmas):
  if gamma_isotope is not None:
@@ -101,18 +109,22 @@ def _write_gamma(br, interaction, i, filename=None, verbose=True):
  frequency=frequencies,
  group_velocity=group_velocities[i],
  heat_capacity=mode_heat_capacities[:, i],
- kappa=None,
  gamma=gamma[j, :, i],
  gamma_isotope=gamma_isotope_at_sigma,
- averaged_pp_interaction=ave_pp[i],
+ averaged_pp_interaction=ave_pp_i,
  mesh_divisors=mesh_divisors,
  grid_point=gp,
  sigma=sigma,
+ kappa_unit_conversion=unit_to_WmK / volume,
  filename=filename,
  verbose=verbose)
  else:
  for j, sigma in enumerate(sigmas):
  for k, bi in enumerate(interaction.get_band_indices()):
+ if ave_pp is None:
+ ave_pp_ik = None
+ else:
+ ave_pp_ik = ave_pp[i, k]
  frequencies = interaction.get_phonons()[0][gp, k]
  if gamma_isotope is not None:
  gamma_isotope_at_sigma = gamma_isotope[j, i, k]
@@ -124,19 +136,17 @@ def _write_gamma(br, interaction, i, filename=None, verbose=True):
  frequency=frequencies,
  group_velocity=group_velocities[i, k],
  heat_capacity=mode_heat_capacities[:, i, k],
- kappa=None,
  gamma=gamma[j, :, i, k],
  gamma_isotope=gamma_isotope_at_sigma,
- averaged_pp_interaction=ave_pp[i, k],
+ averaged_pp_interaction=ave_pp_ik,
  mesh_divisors=mesh_divisors,
  grid_point=gp,
  band_index=bi,
  sigma=sigma,
+ kappa_unit_conversion=unit_to_WmK / volume,
  filename=filename,
  verbose=verbose)
 
-
-
 def _all_bands_exist(interaction):
  band_indices = interaction.get_band_indices()
  num_band = interaction.get_primitive().get_number_of_atoms() * 3
@@ -157,7 +167,7 @@ def _write_triplets(interaction, filename=None):
  filename=filename)
  write_grid_address(grid_address, mesh, filename=filename)
 
-def _write_kappa(br, filename=None, log_level=0):
+def _write_kappa(br, interaction, filename=None, log_level=0):
  temperatures = br.get_temperatures()
  sigmas = br.get_sigmas()
  gamma = br.get_gamma()
@@ -175,7 +185,8 @@ def _write_kappa(br, filename=None, log_level=0):
  num_ignored_phonon_modes = br.get_number_of_ignored_phonon_modes()
  num_band = br.get_frequencies().shape[1]
  num_phonon_modes = br.get_number_of_sampling_grid_points() * num_band
-
+ volume = interaction.get_primitive().get_volume()
+
  for i, sigma in enumerate(sigmas):
  kappa_at_sigma = kappa[i]
  if gamma_isotope is not None:
@@ -194,7 +205,7 @@ def _write_kappa(br, filename=None, log_level=0):
  ("T(K)", "xx", "yy", "zz", "yz", "xz", "xy"))
  for j, (t, k) in enumerate(zip(temperatures, kappa_at_sigma)):
  print(("%7.1f" + " %10.3f" * 6 + " %d/%d") %
- ((t,) + tuple(k) + 
+ ((t,) + tuple(k) +
  (num_ignored_phonon_modes[i, j], num_phonon_modes)))
  else:
  print(("#%6s " + " %-10s" * 6) %
@@ -217,9 +228,10 @@ def _write_kappa(br, filename=None, log_level=0):
  weight=weights,
  mesh_divisors=mesh_divisors,
  sigma=sigma,
+ kappa_unit_conversion=unit_to_WmK / volume,
  filename=filename,
  verbose=log_level)
- 
+
 def _set_gamma_from_file(br, filename=None, verbose=True):
  sigmas = br.get_sigmas()
  mesh = br.get_mesh_numbers()
@@ -335,7 +347,7 @@ def __init__(self,
  self._symmetry = None
  self._point_operations = None
  self._rotations_cartesian = None
- 
+
  self._grid_points = None
  self._grid_weights = None
  self._grid_address = None
@@ -352,7 +364,7 @@ def __init__(self,
  self._averaged_pp_interaction = None
  self._num_ignored_phonon_modes = None
  self._num_sampling_grid_points = None
- 
+
  self._mesh = None
  self._mesh_divisors = None
  self._coarse_mesh = None
@@ -387,16 +399,16 @@ def __init__(self,
  def set_kappa_at_sigmas(self):
  num_band = self._primitive.get_number_of_atoms() * 3
  self._num_sampling_grid_points = 0
- 
+
  for i, grid_point in enumerate(self._grid_points):
  cv = self._cv[:, i, :]
  gp = self._grid_points[i]
  frequencies = self._frequencies[gp]
- 
+
  # Outer product of group velocities (v x v) [num_k*, num_freqs, 3, 3]
  gv_by_gv_tensor, order_kstar = self._get_gv_by_gv(i)
  self._num_sampling_grid_points += order_kstar
- 
+
  # Sum all vxv at k*
  gv_sum2 = np.zeros((6, num_band), dtype='double')
  for j, vxv in enumerate(
@@ -433,7 +445,7 @@ def get_averaged_pp_interaction(self):
 
  def set_averaged_pp_interaction(self, ave_pp):
  self._averaged_pp_interaction = ave_pp
- 
+
  def _run_at_grid_point(self):
  i = self._grid_point_count
  self._show_log_header(i)
@@ -453,14 +465,14 @@ def _run_at_grid_point(self):
  print("Calculating interaction...")
 
  self._set_gamma_at_sigmas(i)
- 
+
  if self._isotope is not None and not self._read_gamma_iso:
  self._set_gamma_isotope_at_sigmas(i)
 
  freqs = self._frequencies[grid_point][self._pp.get_band_indices()]
  self._cv[:, i, :] = self._get_cv(freqs)
  self._set_gv(i)
- 
+
  if self._log_level:
  self._show_log(self._qpoints[i], i)
 
@@ -495,7 +507,7 @@ def _allocate_values(self):
  self._collision = ImagSelfEnergy(
  self._pp,
  unit_conversion=self._gamma_unit_conversion)
- 
+
  def _set_gamma_at_sigmas(self, i):
  for j, sigma in enumerate(self._sigmas):
  if self._log_level:
@@ -522,14 +534,14 @@ def _set_gamma_at_sigmas(self, i):
  self._collision.set_temperature(t)
  self._collision.run()
  self._gamma[j, k, i] = self._collision.get_imag_self_energy()
- 
+
  def _get_gv_by_gv(self, i):
  rotation_map = get_grid_points_by_rotations(
  self._grid_address[self._grid_points[i]],
  self._point_operations,
  self._mesh)
  gv_by_gv = np.zeros((len(self._gv[i]), 3, 3), dtype='double')
- 
+
  for r in self._rotations_cartesian:
  gvs_rot = np.dot(self._gv[i], r.T)
  gv_by_gv += [np.outer(r_gv, r_gv) for r_gv in gvs_rot]
@@ -586,7 +598,7 @@ def _show_log(self, q, i):
  self._point_operations, self._rotations_cartesian)):
  if rotation_map[k] != j:
  continue
- 
+
  print(" k*%-2d (%5.2f %5.2f %5.2f)" %
  ((i + 1,) + tuple(np.dot(rot, q))))
  if self._is_full_pp or self._use_ave_pp:

c/_spglib.c

@@ -60,6 +60,7 @@ static PyObject * get_grid_points_by_rotations(PyObject *self, PyObject *args);
 static PyObject * get_BZ_grid_points_by_rotations(PyObject *self, PyObject *args);
 static PyObject * relocate_BZ_grid_address(PyObject *self, PyObject *args);
 static PyObject * get_symmetry_from_database(PyObject *self, PyObject *args);
+static PyObject * py_niggli_reduce(PyObject *self, PyObject *args);
 
 struct module_state {
  PyObject *error;
@@ -107,6 +108,7 @@ static PyMethodDef _spglib_methods[] = {
  "Rotated grid points in BZ are returned"},
  {"BZ_grid_address", relocate_BZ_grid_address, METH_VARARGS,
  "Relocate grid addresses inside Brillouin zone"},
+ {"niggli_reduce", py_niggli_reduce, METH_VARARGS, "Niggli reduction"},
  {NULL, NULL, 0, NULL}
 };
 
@@ -363,7 +365,7 @@ static PyObject * get_spacegroup_type(PyObject *self, PyObject *args)
  PyObject *array;
  SpglibSpacegroupType symbols;
 
- if (!PyArg_ParseTuple(args, "i",&hall_number)) {
+ if (!PyArg_ParseTuple(args, "i", &hall_number)) {
  return NULL;
  }
 
@@ -819,3 +821,18 @@ static PyObject * relocate_BZ_grid_address(PyObject *self, PyObject *args)
 
  return PyLong_FromLong((long) num_ir_gp);
 }
+
+static PyObject * py_niggli_reduce(PyObject *self, PyObject *args)
+{
+ PyArrayObject* lattice_py;
+ double eps;
+ if (!PyArg_ParseTuple(args, "Od", &lattice_py, &eps)) {
+ return NULL;
+ }
+
+ double (*lattice)[3] = (double(*)[3])PyArray_DATA(lattice_py);
+
+ int result = spg_niggli_reduce(lattice, eps);
+
+ return PyLong_FromLong((long) result);
+}

c/spglib/arithmetic.c

@@ -0,0 +1,65 @@
+/* Copyright (C) 2016 Atsushi Togo */
+/* All rights reserved. */
+
+/* This file is part of spglib. */
+
+/* Redistribution and use in source and binary forms, with or without */
+/* modification, are permitted provided that the following conditions */
+/* are met: */
+
+/* * Redistributions of source code must retain the above copyright */
+/* notice, this list of conditions and the following disclaimer. */
+
+/* * Redistributions in binary form must reproduce the above copyright */
+/* notice, this list of conditions and the following disclaimer in */
+/* the documentation and/or other materials provided with the */
+/* distribution. */
+
+/* * Neither the name of the phonopy project nor the names of its */
+/* contributors may be used to endorse or promote products derived */
+/* from this software without specific prior written permission. */
+
+/* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS */
+/* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT */
+/* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS */
+/* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE */
+/* COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */
+/* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, */
+/* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */
+/* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER */
+/* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT */
+/* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN */
+/* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
+/* POSSIBILITY OF SUCH DAMAGE. */
+
+#include <string.h>
+#include <stdio.h>
+#include "arithmetic.h"
+
+#include "debug.h"
+
+static int arithmetic_crystal_classes[231] = {
+ 0,
+ 1, 2, 3, 3, 4, 5, 5, 6, 6, 7,
+ 7, 8, 7, 7, 8, 9, 9, 9, 9, 10,
+ 10, 11, 12, 12, 13, 13, 13, 13, 13, 13,
+ 13, 13, 13, 13, 14, 14, 14, 15, 15, 15,
+ 15, 16, 16, 17, 17, 17, 18, 18, 18, 18,
+ 18, 18, 18, 18, 18, 18, 18, 18, 18, 18,
+ 18, 18, 19, 19, 19, 19, 19, 19, 20, 20,
+ 21, 21, 21, 21, 22, 22, 22, 22, 23, 23,
+ 24, 25, 26, 26, 26, 26, 27, 27, 28, 28,
+ 28, 28, 28, 28, 28, 28, 29, 29, 30, 30,
+ 30, 30, 30, 30, 30, 30, 31, 31, 31, 31,
+ 32, 32, 32, 32, 33, 33, 33, 33, 34, 34,
+ 35, 35, 36, 36, 36, 36, 36, 36, 36, 36,
+ 36, 36, 36, 36, 36, 36, 36, 36, 37, 37,
+ 37, 37, 38, 38, 38, 39, 40, 41, 42, 43,
+ 42, 43, 42, 43, 44, 45, 46, 45, 46, 47,
+ 47, 48, 48, 49, 49, 50, 50, 51, 51, 51,
+ 51, 51, 51, 52, 53, 53, 54, 54, 54, 54,
+ 54, 54, 55, 55, 55, 55, 56, 56, 57, 57,
+ 58, 58, 58, 58, 59, 60, 61, 59, 61, 62,
+ 62, 63, 63, 64, 62, 64, 65, 65, 66, 66,
+ 67, 65, 65, 67, 68, 69, 70, 68, 69, 70,
+ 71, 71, 71, 71, 72, 72, 72, 72, 73, 73};