diff --git a/anharmonic/BTE_RTA.py b/anharmonic/BTE_RTA.py index ad00ed61..33875955 100644 --- a/anharmonic/BTE_RTA.py +++ b/anharmonic/BTE_RTA.py @@ -301,7 +301,7 @@ class BTE_RTA: self._mesh = self._pp.get_mesh_numbers() if mesh_divisors is None: - self._mesh_divisors = np.array([1, 1, 1], dtype=int) + self._mesh_divisors = np.array([1, 1, 1], dtype='int32') else: self._mesh_divisors = [] for m, n in zip(self._mesh, mesh_divisors): @@ -309,7 +309,7 @@ class BTE_RTA: self._mesh_divisors.append(n) else: self._mesh_divisors.append(1) - self._mesh_divisors = np.array(self._mesh_divisors, dtype=int) + self._mesh_divisors = np.array(self._mesh_divisors, dtype='int32') if (self._mesh_divisors != mesh_divisors).any(): print "Mesh numbers are not dividable by mesh divisors." @@ -348,16 +348,16 @@ class BTE_RTA: print return ((np.array(normal_a, dtype=float), - np.array(normal_w, dtype=int), + np.array(normal_w, dtype='int32'), np.array(normal_f, dtype=float)), (np.array(umklapp_a, dtype=float), - np.array(umklapp_w, dtype=int), + np.array(umklapp_w, dtype='int32'), np.array(umklapp_f, dtype=float))) def _set_pointgroup_operations(self): exist_r_inv = False for rot in self._pp.get_symmetry().get_pointgroup_operations(): - if (rot == -np.eye(3, dtype=int)).all(): + if (rot == -np.eye(3, dtype='int32')).all(): exist_r_inv = True break diff --git a/anharmonic/__init__.py b/anharmonic/__init__.py index f83581c3..55ad2f6e 100644 --- a/anharmonic/__init__.py +++ b/anharmonic/__init__.py @@ -70,7 +70,8 @@ class Phono3py: r2q_TI_index=None, is_Peierls=False, symprec=1e-5, - log_level=0): + log_level=0, + lapack_zheev_uplo='L'): self._supercell = supercell self._primitive = primitive @@ -100,7 +101,8 @@ class Phono3py: symmetrize_fc3_q=self._symmetrize_fc3_q, is_Peierls=self._is_Peierls, log_level=self._log_level, - is_nosym=self._is_nosym) + is_nosym=self._is_nosym, + lapack_zheev_uplo=lapack_zheev_uplo) def set_dynamical_matrix(self, fc2, diff --git a/anharmonic/fc_interpolate.py b/anharmonic/fc_interpolate.py index e5194d20..47c6a7dc 100644 --- a/anharmonic/fc_interpolate.py +++ b/anharmonic/fc_interpolate.py @@ -4,127 +4,141 @@ from phonopy.harmonic.dynamical_matrix import get_smallest_vectors from anharmonic.r2q import get_fc3_reciprocal, get_fc2_reciprocal from phonopy.structure.cells import get_supercell, Primitive, print_cell -def get_fc_interpolation( fc3, - fc2, - supercell, - primitive, - q_mesh, - symprec=1e-5 ): +def get_fc_interpolation(fc3, + fc2, + supercell, + primitive, + q_mesh, + symprec=1e-5): vecs, multi = get_smallest_vectors( supercell, primitive, symprec ) p2s = primitive.get_primitive_to_supercell_map() s2p = primitive.get_supercell_to_primitive_map() - q_mesh = np.array( q_mesh ) - mesh_points = get_mesh_points( q_mesh ) + q_mesh = np.array(q_mesh) + mesh_points = get_mesh_points(q_mesh) - num_atom_prim = len( p2s ) + num_atom_prim = len(p2s) - supercell_intpl = get_supercell( primitive, np.diag( q_mesh ), symprec ) - primitive_intpl = Primitive( supercell_intpl, np.diag(1.0/q_mesh), symprec ) - print_cell( supercell_intpl ) + supercell_intpl = get_supercell(primitive, np.diag(q_mesh), symprec) + primitive_intpl = Primitive(supercell_intpl, np.diag(1.0/q_mesh), symprec) + print_cell(supercell_intpl) - fc3_interpolation = get_fc3_interpolation( fc3, - supercell_intpl, - primitive_intpl, - mesh_points, - vecs, - multi, - p2s, - s2p, - symprec ) + fc3_interpolation = get_fc3_interpolation(fc3, + supercell_intpl, + primitive_intpl, + mesh_points, + vecs, + multi, + p2s, + s2p, + symprec) - fc2_interpolation = get_fc2_interpolation( fc2, - supercell_intpl, - primitive_intpl, - mesh_points, - vecs, - multi, - p2s, - s2p, - symprec ) + fc2_interpolation = get_fc2_interpolation(fc2, + supercell_intpl, + primitive_intpl, + mesh_points, + vecs, + multi, + p2s, + s2p, + symprec) return fc2_interpolation.real, fc3_interpolation.real -def get_fc2_interpolation( fc2, - supercell, primitive, mesh_points, - vecs_orig, multi_orig, p2s_orig, s2p_orig, - symprec ): +def get_fc2_interpolation(fc2, + supercell, + primitive, + mesh_points, + vecs_orig, + multi_orig, + p2s_orig, + s2p_orig, + symprec): num_atom_super = supercell.get_number_of_atoms() num_atom_prim = primitive.get_number_of_atoms() - fc2_intpl = np.zeros( ( num_atom_prim, - num_atom_super, 3, 3 ), dtype=complex ) - vecs, multi = get_smallest_vectors( supercell, primitive, symprec ) + fc2_intpl = np.zeros((num_atom_prim, + num_atom_super, 3, 3), dtype=complex) + vecs, multi = get_smallest_vectors(supercell, primitive, symprec) s2p = primitive.get_supercell_to_primitive_map() p2p = primitive.get_primitive_to_primitive_map() - s2p_special = np.array( [ p2p[x] for x in s2p ] ) for q in mesh_points: - fc2_rec = get_fc2_reciprocal( vecs_orig, - multi_orig, - q, - p2s_orig, - s2p_orig, - fc2, - symprec ) + fc2_rec = get_fc2_reciprocal(vecs_orig, + multi_orig, + q, + p2s_orig, + s2p_orig, + fc2, + symprec) - fc2_intpl += get_py_fc2_realspace( fc2_rec, - q, - s2p, p2p, - vecs, multi, - symprec ) + fc2_intpl += get_py_fc2_realspace(fc2_rec, + q, + s2p, + p2p, + vecs, + multi, + symprec) - return fc2_intpl / len( mesh_points ) + return fc2_intpl / len(mesh_points) -def get_fc3_interpolation( fc3, - supercell, primitive, mesh_points, - vecs_orig, multi_orig, p2s_orig, s2p_orig, - symprec ): +def get_fc3_interpolation(fc3, + supercell, + primitive, + mesh_points, + vecs_orig, + multi_orig, + p2s_orig, + s2p_orig, + symprec ): - q1 = np.array([0,0,0]) + q1 = np.array([0, 0, 0]) num_atom_super = supercell.get_number_of_atoms() num_atom_prim = primitive.get_number_of_atoms() - fc3_intpl = np.zeros( ( num_atom_prim, - num_atom_super, - num_atom_super, 3, 3, 3 ), dtype=complex ) - vecs, multi = get_smallest_vectors( supercell, primitive, symprec ) + fc3_intpl = np.zeros((num_atom_prim, + num_atom_super, + num_atom_super, 3, 3, 3), dtype=complex) + vecs, multi = get_smallest_vectors(supercell, primitive, symprec) s2p = primitive.get_supercell_to_primitive_map() p2p = primitive.get_primitive_to_primitive_map() - s2p_special = np.array( [ p2p[x] for x in s2p ] ) + s2p_special = np.int32([p2p[x] for x in s2p]) - for i, q2 in enumerate( mesh_points ): - for j, q3 in enumerate( mesh_points ): + for i, q2 in enumerate(mesh_points): + for j, q3 in enumerate(mesh_points): sys.stdout.flush() print q2, q3, j+1+i*len(mesh_points), "/", len(mesh_points)**2 - fc3_rec = get_fc3_reciprocal( vecs_orig, - multi_orig, - np.array([ q1, q2, q3 ]), - p2s_orig, - s2p_orig, - fc3, - symprec=symprec ) + fc3_rec = get_fc3_reciprocal(vecs_orig, + multi_orig, + np.array([q1, q2, q3]), + p2s_orig, + s2p_orig, + fc3, + symprec=symprec) try: import anharmonic._phono3py as phono3c - phono3c.fc3_realspace( fc3_intpl, - vecs, - multi, - np.array( [ q1, q2, q3 ] ), - s2p_special, - fc3_rec, - symprec ) + phono3c.fc3_realspace(fc3_intpl, + vecs, + multi, + np.array([q1, q2, q3]), + s2p_special, + fc3_rec, + symprec) except ImportError: - fc3_intpl += get_py_fc3_realspace( fc3_rec, - q2, q3, - s2p_intpl, p2p_intpl, - vecs, multi, - symprec ) + fc3_intpl += get_py_fc3_realspace(fc3_rec, + q2, + q3, + s2p_intpl, + p2p_intpl, + vecs, + multi, + symprec) - return fc3_intpl / len( mesh_points ) ** 2 + return fc3_intpl / len(mesh_points) ** 2 def get_py_fc2_realspace( fc2_rec, q, s2p, p2p, vecs, multi, symprec=1e-5 ): diff --git a/anharmonic/jointDOS.py b/anharmonic/jointDOS.py index 55fbcbc0..765498cf 100644 --- a/anharmonic/jointDOS.py +++ b/anharmonic/jointDOS.py @@ -57,14 +57,16 @@ def get_jointDOS(fixed_grid_points, print "Triplets at q without considering symmetry" sys.stdout.flush() - triplets_at_q, weights_at_q, grid_points = get_nosym_triplets(mesh, gp) + (triplets_at_q, + weights_at_q, + grid_points) = get_nosym_triplets(mesh, gp) else: - triplets_at_q, weights_at_q, grid_points = \ - get_triplets_at_q(gp, - mesh, - primitive.get_cell(), - symmetry.get_pointgroup_operations(), - True) + triplets_at_q, weights_at_q, grid_points = get_triplets_at_q( + gp, + mesh, + primitive.get_cell(), + symmetry.get_pointgroup_operations(), + True) if verbose: print "Grid point (%d):" % gp, grid_points[gp] diff --git a/anharmonic/q2v.py b/anharmonic/q2v.py index 23fa6e5d..d33cb21a 100644 --- a/anharmonic/q2v.py +++ b/anharmonic/q2v.py @@ -33,7 +33,8 @@ class PhononPhonon: symmetrize_fc3_q=False, is_Peierls=False, log_level=False, - is_nosym=False): + is_nosym=False, + lapack_zheev_uplo='L'): self._freq_factor = freq_factor self._cutoff_frequency = 0.01 * self._freq_factor @@ -51,12 +52,12 @@ class PhononPhonon: self._is_Peierls = is_Peierls self._symmetrize_fc3_q = symmetrize_fc3_q self._is_nosym = is_nosym - self._p2s_map = np.array(primitive.get_primitive_to_supercell_map()) - self._s2p_map = np.array(primitive.get_supercell_to_primitive_map()) + self._lapack_zheev_uplo = lapack_zheev_uplo + + self._p2s_map = np.int32(primitive.get_primitive_to_supercell_map()) + self._s2p_map = np.int32(primitive.get_supercell_to_primitive_map()) self._conversion_factor = get_unit_conversion_factor(freq_factor) - # self._shortest_vectors, self._multiplicity = get_shortest_vectors( - # supercell, primitive, symprec) self._shortest_vectors, self._multiplicity = get_smallest_vectors( supercell, primitive, symprec) self._symmetry = Symmetry(primitive, symprec) @@ -159,7 +160,6 @@ class PhononPhonon: "triplets_q-%d%d%d-%d.dat" % (tuple(mesh) + (gp,))) else: self._print_log("Finding ir-triplets at %d\n" % gp) - (triplets_at_q, weights_at_q, self._grid_address) = get_triplets_at_q( @@ -168,7 +168,7 @@ class PhononPhonon: self._primitive.get_cell(), self._symmetry.get_pointgroup_operations(), is_time_reversal=True) - + if self._log_level > 1: t_filename = "triplets_q-%d%d%d-%d.dat" % (tuple(mesh) + (gp,)) write_triplets(triplets_at_q, weights_at_q, mesh, t_filename) @@ -197,9 +197,9 @@ class PhononPhonon: self._print_log("----- phonon-phonon interaction strength ------\n") num_atom = self._primitive.get_number_of_atoms() if band_indices == None: - self._band_indices = np.arange(num_atom * 3, dtype=int) + self._band_indices = np.arange(num_atom * 3, dtype='int32') else: - self._band_indices = np.array(band_indices, dtype=int) + self._band_indices = np.int32(band_indices) self._print_log(("Band indices: [" + " %d" * len(self._band_indices) + " ]\n") % tuple(self._band_indices + 1)) self.set_harmonic_phonons() @@ -302,10 +302,10 @@ class PhononPhonon: phono3c.interaction_strength( self._amplitude_at_q, self._frequencies_at_q, + len(self._q_grid), q0, q1s, q2s, - self._weights_at_q, svecs_fc2, multi_fc2, self._shortest_vectors, @@ -327,7 +327,8 @@ class PhononPhonon: self._freq_factor * self._factor, self._cutoff_frequency, self._symmetrize_fc3_q, - self._r2q_TI_index) + self._r2q_TI_index, + self._lapack_zheev_uplo) def _set_py_interaction_strength(self): @@ -519,8 +520,8 @@ def get_c_triplet_interaction_strength(amplitude, r2q_TI_index): import anharmonic._phono3py as phono3c - p2s_map = primitive.get_primitive_to_supercell_map() - s2p_map = primitive.get_supercell_to_primitive_map() + p2s_map = np.int32(primitive.get_primitive_to_supercell_map()) + s2p_map = np.int32(primitive.get_supercell_to_primitive_map()) phono3c.triplet_interaction_strength(amplitude, freqs, @@ -528,8 +529,8 @@ def get_c_triplet_interaction_strength(amplitude, shortest_vectors, multiplicity, q_set, - np.array(p2s_map), - np.array(s2p_map), + p2s_map, + s2p_map, fc3, primitive.get_masses(), band_indices, @@ -583,7 +584,7 @@ def get_py_triplet_interaction_strength(amplitude, multiplicity, np.array([q_set[k[0]], q_set[k[1]], - q_set[k[2]]]), + q_set[k[2]]], dtype=float), p2s_map, s2p_map, fc3, @@ -739,8 +740,13 @@ def print_fc3_q(fc3, q_set, shortest_vectors, multiplicity, primitive, symprec=1 p2s_map = primitive.get_primitive_to_supercell_map() s2p_map = primitive.get_supercell_to_primitive_map() - fc3_q = get_fc3_reciprocal(shortest_vectors, multiplicity, - q_set, p2s_map, s2p_map, fc3, symprec=symprec) + fc3_q = get_fc3_reciprocal(shortest_vectors, + multiplicity, + q_set, + p2s_map, + s2p_map, + fc3, + symprec=symprec) num_atom = primitive.get_number_of_atoms() for i in range(num_atom): diff --git a/anharmonic/r2q.py b/anharmonic/r2q.py index e566935d..19363c44 100644 --- a/anharmonic/r2q.py +++ b/anharmonic/r2q.py @@ -118,8 +118,8 @@ def get_c_fc3_reciprocal(shortest_vectors, shortest_vectors, multiplicity, q_set, - np.array(p2s_map), - np.array(s2p_map), + np.int32(p2s_map), + np.int32(s2p_map), fc3, r2q_TI_index) diff --git a/anharmonic/triplets.py b/anharmonic/triplets.py index 71610ab2..db2ef19b 100644 --- a/anharmonic/triplets.py +++ b/anharmonic/triplets.py @@ -11,7 +11,7 @@ def get_triplets_at_q(gp, third_q, grid_address) = spg.get_triplets_reciprocal_mesh_at_q(gp, mesh, - rotations, + np.int64(rotations), is_time_reversal) weights_at_q = [] triplets_at_q = [] @@ -20,18 +20,19 @@ def get_triplets_at_q(gp, weights_at_q.append(w) triplets_at_q.append([gp, i, q]) - weights_at_q = np.array(weights_at_q) + weights_at_q = np.array(weights_at_q, dtype='int32') + triplets_at_q = np.array(triplets_at_q, dtype='int32') assert np.prod(mesh) == weights_at_q.sum(), \ "Num grid points %d, sum of weight %d" % ( np.prod(mesh), weights_at_q.sum()) - return np.array(triplets_at_q), weights_at_q, grid_address + return triplets_at_q, weights_at_q, grid_address def get_nosym_triplets(mesh, grid_point0): grid_address = get_grid_address(mesh) - triplets = np.zeros((len(grid_address), 3), dtype=int) - weights = np.ones(len(grid_address), dtype=int) + triplets = np.zeros((len(grid_address), 3), dtype='int32') + weights = np.ones(len(grid_address), dtype='int32') for i, g1 in enumerate(grid_address): g2 = - (grid_address[grid_point0] + g1) q = get_address(g2, mesh) diff --git a/c/_phono3py.c b/c/_phono3py.c index a6e6f3a1..24868c85 100644 --- a/c/_phono3py.c +++ b/c/_phono3py.c @@ -26,14 +26,14 @@ static int distribute_fc3(double *fc3, const int third_atom_rot, const double *positions, const int num_atom, - const int rot[3][3], + const int * rot, const double *rot_cartesian, const double *trans, const double symprec); static int get_atom_by_symmetry(const int atom_number, const int num_atom, const double *positions, - const int rot[3][3], + const int *rot, const double *trans, const double symprec); static void tensor3_roation(double *fc3, @@ -79,7 +79,6 @@ static PyObject * py_get_interaction_strength(PyObject *self, PyObject *args) PyArrayObject* qvec0; PyArrayObject* qvec1s; PyArrayObject* qvec2s; - PyArrayObject* weights; PyArrayObject* shortest_vectors_fc2; PyArrayObject* multiplicity_fc2; PyArrayObject* shortest_vectors_fc3; @@ -98,15 +97,16 @@ static PyObject * py_get_interaction_strength(PyObject *self, PyObject *args) PyArrayObject* q_direction; PyArrayObject* dielectric_constant; double cutoff_frequency, nac_factor, freq_unit_conversion_factor; - int r2q_TI_index, is_symmetrize_fc3_q; + int num_grid_points, r2q_TI_index, is_symmetrize_fc3_q; + char uplo; - if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOOOOOOOOOdddii", + if (!PyArg_ParseTuple(args, "OOiOOOOOOOOOOOOOOOOOOOOdddiic", &litude, &frequencies, + &num_grid_points, &qvec0, &qvec1s, &qvec2s, - &weights, &shortest_vectors_fc2, &multiplicity_fc2, &shortest_vectors_fc3, @@ -128,13 +128,14 @@ static PyObject * py_get_interaction_strength(PyObject *self, PyObject *args) &freq_unit_conversion_factor, &cutoff_frequency, &is_symmetrize_fc3_q, - &r2q_TI_index)) { + &r2q_TI_index, + &uplo)) { return NULL; } int i, j; int svecs_dims[4]; - Array1D * s2p_fc2, * p2s_fc2, * s2p_fc3, * p2s_fc3, *w, *bands; + Array1D * s2p_fc2, * p2s_fc2, * s2p_fc3, * p2s_fc3, *bands; Array2D * multi_fc2, * multi_fc3; ShortestVecs * svecs_fc2, * svecs_fc3; @@ -143,18 +144,17 @@ static PyObject * py_get_interaction_strength(PyObject *self, PyObject *args) const double *q0 = (double*)qvec0->data; const double *q1s = (double*)qvec1s->data; const double *q2s = (double*)qvec2s->data; - const long *weights_long = (long*)weights->data; const double* masses_fc2 = (double*)atomic_masses_fc2->data; const double* masses_fc3 = (double*)atomic_masses_fc3->data; - const long* multiplicity_long_fc2 = (long*)multiplicity_fc2->data; - const long* multiplicity_long_fc3 = (long*)multiplicity_fc3->data; - const long* p2s_map_fc2_long = (long*)p2s_map_fc2->data; - const long* s2p_map_fc2_long = (long*)s2p_map_fc2->data; - const long* p2s_map_fc3_long = (long*)p2s_map_fc3->data; - const long* s2p_map_fc3_long = (long*)s2p_map_fc3->data; + const int* multiplicity_fc2_int = (int*)multiplicity_fc2->data; + const int* multiplicity_fc3_int = (int*)multiplicity_fc3->data; + const int* p2s_map_fc2_int = (int*)p2s_map_fc2->data; + const int* s2p_map_fc2_int = (int*)s2p_map_fc2->data; + const int* p2s_map_fc3_int = (int*)p2s_map_fc3->data; + const int* s2p_map_fc3_int = (int*)s2p_map_fc3->data; const double* fc2 = (double*)force_constants_second->data; const double* fc3 = (double*)force_constants_third->data; - const long* bands_long = (long*)band_indicies->data; + const int* bands_int = (int*)band_indicies->data; const double* rec_lat = (double*)reciprocal_lattice->data; double* born; if ((PyObject*)born_effective_charge == Py_None) { @@ -178,12 +178,11 @@ static PyObject * py_get_interaction_strength(PyObject *self, PyObject *args) const int num_satom_fc2 = (int)multiplicity_fc2->dimensions[0]; const int num_satom_fc3 = (int)multiplicity_fc3->dimensions[0]; const int num_patom = (int)multiplicity_fc3->dimensions[1]; - const int num_triplets = (int)weights->dimensions[0]; - + const int num_triplets = (int)amplitude->dimensions[0]; bands = alloc_Array1D(band_indicies->dimensions[0]); for (i = 0; i < bands->d1; i++) { - bands->data[i] = (int)bands_long[i]; + bands->data[i] = (int)bands_int[i]; } for (i = 0; i < 4; i++) { @@ -200,43 +199,40 @@ static PyObject * py_get_interaction_strength(PyObject *self, PyObject *args) multi_fc2 = alloc_Array2D(num_satom_fc2, num_patom); for (i = 0; i < num_satom_fc2; i++) { for (j = 0; j < num_patom; j++) { - multi_fc2->data[i][j] = multiplicity_long_fc2[i * num_patom + j]; + multi_fc2->data[i][j] = multiplicity_fc2_int[i * num_patom + j]; } } multi_fc3 = alloc_Array2D(num_satom_fc3, num_patom); for (i = 0; i < num_satom_fc3; i++) { for (j = 0; j < num_patom; j++) { - multi_fc3->data[i][j] = multiplicity_long_fc3[i * num_patom + j]; + multi_fc3->data[i][j] = multiplicity_fc3_int[i * num_patom + j]; } } s2p_fc2 = alloc_Array1D(num_satom_fc2); for (i = 0; i < num_satom_fc2; i++) { - s2p_fc2->data[i] = s2p_map_fc2_long[i]; + s2p_fc2->data[i] = s2p_map_fc2_int[i]; } p2s_fc2 = alloc_Array1D(num_patom); for (i = 0; i < num_patom; i++) { - p2s_fc2->data[i] = p2s_map_fc2_long[i]; + p2s_fc2->data[i] = p2s_map_fc2_int[i]; } s2p_fc3 = alloc_Array1D(num_satom_fc3); for (i = 0; i < num_satom_fc3; i++) { - s2p_fc3->data[i] = s2p_map_fc3_long[i]; + s2p_fc3->data[i] = s2p_map_fc3_int[i]; } p2s_fc3 = alloc_Array1D(num_patom); for (i = 0; i < num_patom; i++) { - p2s_fc3->data[i] = p2s_map_fc3_long[i]; - } - w = alloc_Array1D(num_triplets); - for (i = 0; i < num_triplets; i++) { - w->data[i] = weights_long[i]; + p2s_fc3->data[i] = p2s_map_fc3_int[i]; } get_interaction_strength(amps, freqs, + num_triplets, + num_grid_points, q0, q1s, q2s, - w, fc2, fc3, masses_fc2, @@ -258,7 +254,8 @@ static PyObject * py_get_interaction_strength(PyObject *self, PyObject *args) freq_unit_conversion_factor, cutoff_frequency, is_symmetrize_fc3_q, - r2q_TI_index); + r2q_TI_index, + uplo); free_Array1D(p2s_fc2); free_Array1D(s2p_fc2); @@ -266,7 +263,6 @@ static PyObject * py_get_interaction_strength(PyObject *self, PyObject *args) free_Array1D(s2p_fc3); free_Array2D(multi_fc2); free_Array2D(multi_fc3); - free_Array1D(w); free_Array1D(bands); free_shortest_vecs(svecs_fc2); free_shortest_vecs(svecs_fc3); @@ -320,10 +316,10 @@ static PyObject * py_get_triplet_interaction_strength(PyObject *self, const double* freqs = (double*)omegas->data; const npy_cdouble* eigvecs = (npy_cdouble*)eigenvectors->data; const double* masses = (double*)atomic_masses->data; - const long* multiplicity_long = (long*)multiplicity->data; - const long* p2s_map_long = (long*)p2s_map->data; - const long* s2p_map_long = (long*)s2p_map->data; - const long* band_indices_long = (long*)set_of_band_indices->data; + const int* multiplicity_int = (int*)multiplicity->data; + const int* p2s_map_int = (int*)p2s_map->data; + const int* s2p_map_int = (int*)s2p_map->data; + const int* band_indices_int = (int*)set_of_band_indices->data; const double* fc3 = (double*)force_constants_third->data; const double* q_vecs = (double*)q_triplet->data; @@ -340,21 +336,21 @@ static PyObject * py_get_triplet_interaction_strength(PyObject *self, multi = alloc_Array2D(num_satom, num_patom); for (i = 0; i < num_satom; i++) { for (j = 0; j < num_patom; j++) { - multi->data[i][j] = multiplicity_long[i * num_patom + j]; + multi->data[i][j] = multiplicity_int[i * num_patom + j]; } } s2p = alloc_Array1D(num_satom); for (i = 0; i < num_satom; i++) { - s2p->data[i] = s2p_map_long[i]; + s2p->data[i] = s2p_map_int[i]; } p2s = alloc_Array1D(num_patom); for (i = 0; i < num_patom; i++) { - p2s->data[i] = p2s_map_long[i]; + p2s->data[i] = p2s_map_int[i]; } band_indices = alloc_Array1D(num_band0); for (i = 0; i < num_band0; i++) { - band_indices->data[i] = band_indices_long[i]; + band_indices->data[i] = band_indices_int[i]; } lpk_eigvecs = (lapack_complex_double*) @@ -363,7 +359,7 @@ static PyObject * py_get_triplet_interaction_strength(PyObject *self, lpk_eigvecs[i] = lapack_make_complex_double (eigvecs[i].real, eigvecs[i].imag); } - + get_triplet_interaction_strength(amps, fc3, q_vecs, @@ -485,9 +481,9 @@ static PyObject * py_get_fc3_reciprocal(PyObject *self, PyObject *args) npy_cdouble* fc3_q = (npy_cdouble*)force_constants_third_reciprocal->data; - const long* multiplicity_long = (long*)multiplicity->data; - const long* p2s_map_long = (long*)p2s_map->data; - const long* s2p_map_long = (long*)s2p_map->data; + const int* multiplicity_int = (int*)multiplicity->data; + const int* p2s_map_int = (int*)p2s_map->data; + const int* s2p_map_int = (int*)s2p_map->data; const double* fc3 = (double*)force_constants_third->data; const int num_satom = (int)multiplicity->dimensions[0]; @@ -503,17 +499,17 @@ static PyObject * py_get_fc3_reciprocal(PyObject *self, PyObject *args) multi = alloc_Array2D(num_satom, num_patom); for (i = 0; i < num_satom; i++) { for (j = 0; j < num_patom; j++) { - multi->data[i][j] = multiplicity_long[i * num_patom + j]; + multi->data[i][j] = multiplicity_int[i * num_patom + j]; } } s2p = alloc_Array1D(num_satom); p2s = alloc_Array1D(num_patom); for (i = 0; i < num_satom; i++) { - s2p->data[i] = s2p_map_long[i]; + s2p->data[i] = s2p_map_int[i]; } for (i = 0; i < num_patom; i++) { - p2s->data[i] = p2s_map_long[i]; + p2s->data[i] = p2s_map_int[i]; } q = alloc_DArray2D(3, 3); @@ -571,16 +567,15 @@ static PyObject * py_get_fc3_realspace(PyObject *self, PyObject *args) int i, j; int * svecs_dimensions; - Array1D * s2p; Array2D * multi; ShortestVecs * svecs; lapack_complex_double *lpk_fc3_real, *lpk_fc3_rec; npy_cdouble* fc3_real = (npy_cdouble*)fc3_realspace->data; const npy_cdouble* fc3_rec = (npy_cdouble*)fc3_reciprocal->data; - const long* multi_long = (long*)multiplicity->data; + const int* multi_int = (int*)multiplicity->data; const double* q_triplet = (double*)qpoints_triplet->data; - const long* s2p_map_long = (long*)s2p_map->data; + const int* s2p_map_int = (int*)s2p_map->data; const int num_satom = (int)multiplicity->dimensions[0]; const int num_patom = (int)multiplicity->dimensions[1]; @@ -595,15 +590,10 @@ static PyObject * py_get_fc3_realspace(PyObject *self, PyObject *args) multi = alloc_Array2D(num_satom, num_patom); for (i = 0; i < num_satom; i++) { for (j = 0; j < num_patom; j++) { - multi->data[i][j] = multi_long[i * num_patom + j]; + multi->data[i][j] = multi_int[i * num_patom + j]; } } - s2p = alloc_Array1D(num_satom); - for (i = 0; i < num_satom; i++) { - s2p->data[i] = s2p_map_long[i]; - } - lpk_fc3_real = (lapack_complex_double*) malloc(sizeof(lapack_complex_double) * num_satom * num_satom * num_satom * 27); @@ -618,7 +608,7 @@ static PyObject * py_get_fc3_realspace(PyObject *self, PyObject *args) svecs, multi, q_triplet, - s2p, + s2p_map_int, lpk_fc3_rec); for (i = 0; i < num_satom * num_satom * num_satom * 27; i++) { fc3_real[i].real = lapack_complex_double_real(lpk_fc3_real[i]); @@ -626,7 +616,6 @@ static PyObject * py_get_fc3_realspace(PyObject *self, PyObject *args) } free(lpk_fc3_real); free(lpk_fc3_rec); - free_Array1D(s2p); free_Array2D(multi); free_shortest_vecs(svecs); @@ -671,19 +660,13 @@ static PyObject * py_get_gamma(PyObject *self, PyObject *args) double* dfun = (double*)gammas->data; const double* o = (double*)omegas->data; const double* amp = (double*)amplitudes->data; - const long* w_long = (long*)weights->data; + const int* w = (int*)weights->data; const double* f = (double*)frequencies->data; const int num_band0 = (int)amplitudes->dimensions[1]; const int num_band = (int)amplitudes->dimensions[2]; const int num_omega = (int)omegas->dimensions[0]; const int num_triplet = (int)weights->dimensions[0]; - int i; - int w[num_triplet]; - for (i = 0; i < num_triplet; i++) { - w[i] = w_long[i]; - } - get_gamma(dfun, num_omega, num_triplet, @@ -722,16 +705,11 @@ static PyObject * py_get_jointDOS(PyObject *self, PyObject *args) double* jdos = (double*)jointdos->data; const double* o = (double*)omegas->data; - const long* w_long = (long*)weights->data; + const int* w = (int*)weights->data; const double* f = (double*)frequencies->data; const int num_band = (int)frequencies->dimensions[2]; const int num_omega = (int)omegas->dimensions[0]; const int num_triplet = (int)weights->dimensions[0]; - int i; - int w[num_triplet]; - for (i = 0; i < num_triplet; i++) { - w[i] = w_long[i]; - } get_jointDOS(jdos, num_omega, @@ -809,16 +787,9 @@ static PyObject * py_distribute_fc3(PyObject *self, PyObject *args) &symprec)) return NULL; - int i, j; double* fc3 = (double*)force_constants_third->data; const double* pos = (double*)positions->data; - const long* rot_long = (long*)rotation->data; - int rot_int[3][3]; - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - rot_int[i][j] = rot_long[i*3 + j]; - } - } + const int* rot = (int*)rotation->data; const double* rot_cart_inv = (double*)rotation_cart_inv->data; const double* trans = (double*)translation->data; const int num_atom = (int)positions->dimensions[0]; @@ -828,7 +799,7 @@ static PyObject * py_distribute_fc3(PyObject *self, PyObject *args) third_atom_rot, pos, num_atom, - rot_int, + rot, rot_cart_inv, trans, symprec)); @@ -858,7 +829,7 @@ static PyObject * py_phonopy_zheev(PyObject *self, PyObject *args) a[i] = lapack_make_complex_double(dynmat[i].real, dynmat[i].imag); } - info = phonopy_zheev(eigvals, a, dimension); + info = phonopy_zheev(eigvals, a, dimension, 'L'); for (i = 0; i < dimension * dimension; i++) { dynmat[i].real = lapack_complex_double_real(a[i]); @@ -876,7 +847,7 @@ static int distribute_fc3(double *fc3, const int third_atom_rot, const double *positions, const int num_atom, - const int rot[3][3], + const int *rot, const double *rot_cart_inv, const double *trans, const double symprec) @@ -919,7 +890,7 @@ static int distribute_fc3(double *fc3, static int get_atom_by_symmetry(const int atom_number, const int num_atom, const double *positions, - const int rot[3][3], + const int *rot, const double *trans, const double symprec) { @@ -929,7 +900,7 @@ static int get_atom_by_symmetry(const int atom_number, for (i = 0; i < 3; i++) { rot_pos[i] = trans[i]; for (j = 0; j < 3; j++) { - rot_pos[i] += rot[i][j] * positions[atom_number * 3 + j]; + rot_pos[i] += rot[i * 3 + j] * positions[atom_number * 3 + j]; } } diff --git a/c/_phonopy.c b/c/_phonopy.c index 9ab9e7bd..5991c9c8 100644 --- a/c/_phonopy.c +++ b/c/_phonopy.c @@ -59,7 +59,7 @@ static int distribute_fc2(double * fc2, const int atom_disp, const int map_atom_disp, const double * r_cart, - int r[3][3], + const int * r, const double * t, const double symprec); static int get_rotated_forces(double * rotated_forces, @@ -67,7 +67,7 @@ static int get_rotated_forces(double * rotated_forces, const int num_pos, const int atom_number, const double * f, - int (*r)[3][3], + const int * r, const int num_rot, const double symprec); static int nint(const double a); @@ -96,8 +96,8 @@ static PyObject * py_get_dynamical_matrix(PyObject *self, PyObject *args) PyArrayObject* q_vector; PyArrayObject* multiplicity; PyArrayObject* mass; - PyArrayObject* s2p_map; - PyArrayObject* p2s_map; + PyArrayObject* super2prim_map; + PyArrayObject* prim2super_map; if (!PyArg_ParseTuple(args, "OOOOOOOOO", &dynamical_matrix_real, @@ -107,39 +107,21 @@ static PyObject * py_get_dynamical_matrix(PyObject *self, PyObject *args) &r_vector, &multiplicity, &mass, - &s2p_map, - &p2s_map)) + &super2prim_map, + &prim2super_map)) return NULL; - int i; double* dm_r = (double*)dynamical_matrix_real->data; double* dm_i = (double*)dynamical_matrix_imag->data; const double* fc = (double*)force_constants->data; const double* q = (double*)q_vector->data; const double* r = (double*)r_vector->data; const double* m = (double*)mass->data; - const long* multi_long = (long*)multiplicity->data; - const long* s2p_map_long = (long*)s2p_map->data; - const long* p2s_map_long = (long*)p2s_map->data; - const int num_patom = p2s_map->dimensions[0]; - const int num_satom = s2p_map->dimensions[0]; - - int *multi, *s2p_map_int, *p2s_map_int; - - multi = (int*) malloc(num_patom * num_satom * sizeof(int)); - for (i = 0; i < num_patom*num_satom; i++) { - multi[i] = (int)multi_long[i]; - } - - s2p_map_int = (int*) malloc(num_satom * sizeof(int)); - for (i = 0; i < num_satom; i++) { - s2p_map_int[i] = (int)s2p_map_long[i]; - } - - p2s_map_int = (int*) malloc(num_patom * sizeof(int)); - for (i = 0; i < num_patom; i++) { - p2s_map_int[i] = (int)p2s_map_long[i]; - } + const int* multi = (int*)multiplicity->data; + const int* s2p_map = (int*)super2prim_map->data; + const int* p2s_map = (int*)prim2super_map->data; + const int num_patom = prim2super_map->dimensions[0]; + const int num_satom = super2prim_map->dimensions[0]; get_dynamical_matrix_at_q(dm_r, dm_i, @@ -150,14 +132,10 @@ static PyObject * py_get_dynamical_matrix(PyObject *self, PyObject *args) r, multi, m, - s2p_map_int, - p2s_map_int, + s2p_map, + p2s_map, NULL); - free(multi); - free(s2p_map_int); - free(p2s_map_int); - Py_RETURN_NONE; } @@ -172,8 +150,8 @@ static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args) PyArrayObject* q_vector; PyArrayObject* multiplicity; PyArrayObject* mass; - PyArrayObject* s2p_map; - PyArrayObject* p2s_map; + PyArrayObject* super2prim_map; + PyArrayObject* prim2super_map; PyArrayObject* born; double factor; @@ -185,8 +163,8 @@ static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args) &r_vector, &multiplicity, &mass, - &s2p_map, - &p2s_map, + &super2prim_map, + &prim2super_map, &q_cart_vector, &born, &factor)) @@ -200,34 +178,18 @@ static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args) const double* r = (double*)r_vector->data; const double* m = (double*)mass->data; const double* z = (double*)born->data; - const long* multi_long = (long*)multiplicity->data; - const long* s2p_map_long = (long*)s2p_map->data; - const long* p2s_map_long = (long*)p2s_map->data; - const int num_patom = p2s_map->dimensions[0]; - const int num_satom = s2p_map->dimensions[0]; + const int* multi = (int*)multiplicity->data; + const int* s2p_map = (int*)super2prim_map->data; + const int* p2s_map = (int*)prim2super_map->data; + const int num_patom = prim2super_map->dimensions[0]; + const int num_satom = super2prim_map->dimensions[0]; - int i, n; - int *multi, *s2p_map_int, *p2s_map_int; + int n; double *charge_sum; charge_sum = (double*) malloc(sizeof(double) * num_patom * num_patom * 9); n = num_satom / num_patom; - multi = (int*) malloc(num_patom * num_satom * sizeof(int)); - for (i = 0; i < num_patom * num_satom; i++) { - multi[i] = (int)multi_long[i]; - } - - s2p_map_int = (int*) malloc(num_satom * sizeof(int)); - for (i = 0; i < num_satom; i++) { - s2p_map_int[i] = (int)s2p_map_long[i]; - } - - p2s_map_int = (int*) malloc(num_patom * sizeof(int)); - for (i = 0; i < num_patom; i++) { - p2s_map_int[i] = (int)p2s_map_long[i]; - } - get_charge_sum(charge_sum, num_patom, factor / n, q_cart, z); get_dynamical_matrix_at_q(dm_r, dm_i, @@ -238,14 +200,11 @@ static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args) r, multi, m, - s2p_map_int, - p2s_map_int, + s2p_map, + p2s_map, charge_sum); free(charge_sum); - free(multi); - free(s2p_map_int); - free(p2s_map_int); Py_RETURN_NONE; } @@ -265,7 +224,7 @@ static PyObject * py_get_thermal_properties(PyObject *self, PyObject *args) } const double* freqs = (double*)frequencies->data; - const long* w = (long*)weights->data; + const int* w = (int*)weights->data; const int num_qpoints = frequencies->dimensions[0]; const int num_bands = frequencies->dimensions[1]; @@ -356,22 +315,13 @@ static PyObject * py_distribute_fc2(PyObject *self, PyObject *args) return NULL; } - const long* r_long = (long*)rotation->data; + const int* r = (int*)rotation->data; const double* r_cart = (double*)rotation_cart->data; double* fc2 = (double*)force_constants->data; const double* t = (double*)translation->data; const double* pos = (double*)positions->data; const int num_pos = positions->dimensions[0]; - - int i, j; - int r[3][3]; - - for (i = 0; i < 3; i++){ - for (j = 0; j < 3; j++){ - r[i][j] = (int) r_long[ i * 3 + j ]; - } - } distribute_fc2(fc2, pos, num_pos, @@ -391,7 +341,7 @@ static int distribute_fc2(double * fc2, const int atom_disp, const int map_atom_disp, const double * r_cart, - int r[3][3], + const int * r, const double * t, const double symprec) { @@ -402,18 +352,18 @@ static int distribute_fc2(double * fc2, #pragma omp parallel for private(j, k, l, m, rot_pos, diff, is_found, rot_atom, address_new, address) for (i = 0; i < num_pos; i++) { for (j = 0; j < 3; j++) { - rot_pos[ j ] = t[j]; + rot_pos[j] = t[j]; for (k = 0; k < 3; k++) { - rot_pos[ j ] += r[ j ][ k ] * pos[ i * 3 + k ]; + rot_pos[j] += r[j * 3 + k] * pos[i * 3 + k]; } } for (j = 0; j < num_pos; j++) { is_found = 1; for (k = 0; k < 3; k++) { - diff[ k ] = pos[ j * 3 + k ] - rot_pos[ k ]; - diff[ k ] -= nint(diff[k]); - if (fabs(diff[ k ]) > symprec) { + diff[k] = pos[j * 3 + k] - rot_pos[k]; + diff[k] -= nint(diff[k]); + if (fabs(diff[k]) > symprec) { is_found = 0; break; } @@ -436,9 +386,9 @@ static int distribute_fc2(double * fc2, for (k = 0; k < 3; k++) { for (l = 0; l < 3; l++) { for (m = 0; m < 3; m++) { - fc2[ address_new + j * 3 + k ] += - r_cart[ l * 3 + j ] * r_cart[ m * 3 + k ] * - fc2[ address + l * 3 + m ]; + fc2[address_new + j * 3 + k] += + r_cart[l * 3 + j] * r_cart[m * 3 + k] * + fc2[address + l * 3 + m]; } } } @@ -469,34 +419,21 @@ static PyObject * py_get_rotated_forces(PyObject *self, PyObject *args) return NULL; } - const long* rots_long = (long*)rotations->data; + const int* r = (int*)rotations->data; const int num_rot = rotations->dimensions[0]; const double* pos = (double*)positions->data; const int num_pos = positions->dimensions[0]; const double* f = (double*)forces->data; double* rot_forces = (double*)rotated_forces->data; - int i, j, k; - int (*r)[3][3]; - if (! ((r = (int (*)[3][3]) malloc(sizeof(int[3][3]) * num_rot)) == NULL)) { - for (i = 0; i < num_rot; i++){ - for (j = 0; j < 3; j++){ - for (k = 0; k < 3; k++){ - r[i][j][k] = (int) rots_long[ i * 9 + j * 3 + k ]; - } - } - } - get_rotated_forces(rot_forces, - pos, - num_pos, - atom_number, - f, - r, - num_rot, - symprec); - - free(r); - } + get_rotated_forces(rot_forces, + pos, + num_pos, + atom_number, + f, + r, + num_rot, + symprec); Py_RETURN_NONE; } @@ -507,7 +444,7 @@ static int get_rotated_forces(double * rotated_forces, const int num_pos, const int atom_number, const double * f, - int (*r)[3][3], + const int * r, const int num_rot, const double symprec) { @@ -517,25 +454,25 @@ static int get_rotated_forces(double * rotated_forces, #pragma omp parallel for private(j, k, rot_pos, diff, is_found) for (i = 0; i < num_rot; i++) { for (j = 0; j < 3; j++) { - rot_pos[ j ] = 0; + rot_pos[j] = 0; for (k = 0; k < 3; k++) { - rot_pos[ j ] += r[ i ][ j ][ k ] * pos[ atom_number * 3 + k ]; + rot_pos[j] += r[i * 9 + j * 3 + k] * pos[atom_number * 3 + k]; } } for (j = 0; j < num_pos; j++) { is_found = 1; for (k = 0; k < 3; k++) { - diff[ k ] = pos[ j * 3 + k ] - rot_pos[ k ]; - diff[ k ] -= nint(diff[k]); - if (fabs(diff[ k ]) > symprec) { + diff[k] = pos[j * 3 + k] - rot_pos[k]; + diff[k] -= nint(diff[k]); + if (fabs(diff[k]) > symprec) { is_found = 0; break; } } if (is_found) { for (k = 0; k < 3; k++) { - rotated_forces[ i * 3 + k ] = f[ j * 3 + k ]; + rotated_forces[i * 3 + k] = f[j * 3 + k]; } break; } diff --git a/c/anharmonic/interaction_strength.c b/c/anharmonic/interaction_strength.c index d75e37e6..ef32aed0 100644 --- a/c/anharmonic/interaction_strength.c +++ b/c/anharmonic/interaction_strength.c @@ -45,14 +45,16 @@ static int get_phonons(lapack_complex_double *a, const double *dielectric, const double *reciprocal_lattice, const double *q_direction, - const double nac_factor); + const double nac_factor, + const char uplo); int get_interaction_strength(double *amps, double *freqs, + const int num_triplets, + const int num_grid_points, const double *q0, const double *q1s, const double *q2s, - const Array1D *weights, const double *fc2, const double *fc3, const double *masses_fc2, @@ -74,19 +76,15 @@ int get_interaction_strength(double *amps, const double freq_unit_factor, const double cutoff_frequency, const int is_symmetrize_fc3_q, - const int r2q_TI_index) + const int r2q_TI_index, + const char uplo) { - int i, j, num_triplets, num_patom, num_grid_points; + int i, j, num_patom; double *q_vecs; double *w0 ,*w; lapack_complex_double *a0, *a; num_patom = p2s_fc2->d1; - num_triplets = weights->d1; - num_grid_points = 0; - for (i = 0; i < weights->d1; i++) { - num_grid_points += weights->data[i]; - } w0 = (double*)malloc(sizeof(double) * num_patom * 3); a0 = (lapack_complex_double*) @@ -105,7 +103,8 @@ int get_interaction_strength(double *amps, dielectric, reciprocal_lattice, q_direction, - nac_factor); + nac_factor, + uplo); #pragma omp parallel for private(j, w, a, q_vecs) for (i = 0; i < num_triplets; i++) { @@ -139,7 +138,8 @@ int get_interaction_strength(double *amps, dielectric, reciprocal_lattice, NULL, - nac_factor); + nac_factor, + uplo); get_phonons(a + num_patom * num_patom * 18, w + num_patom * 6, @@ -154,14 +154,14 @@ int get_interaction_strength(double *amps, dielectric, reciprocal_lattice, NULL, - nac_factor); + nac_factor, + uplo); for (j = 0; j < num_patom * 9; j++) { freqs[i * num_patom * 9 + j] = sqrt(fabs(w[j])) * freq_unit_factor * ((w[j] > 0) - (w[j] < 0)); } - get_triplet_interaction_strength (amps + i * band_indices->d1 * num_patom * num_patom * 9, fc3, @@ -187,160 +187,7 @@ int get_interaction_strength(double *amps, for (i = 0; i < num_triplets * band_indices->d1 * num_patom * num_patom * 9; i++) { - amps[i] = lapack_make_complex_double - (lapack_complex_double_real(amps[i]) / num_grid_points, - lapack_complex_double_imag(amps[i]) / num_grid_points); - } - - free(w0); - free(a0); - - return 1; -} - -int get_interaction_strength2(double *amps, - double *freqs, - const double *q0, - const double *q1s, - const double *q2s, - const Array1D *weights, - const double *fc2, - const double *fc3, - const double *masses_fc2, - const double *masses_fc3, - const Array1D *p2s_fc2, - const Array1D *s2p_fc2, - const Array1D *p2s_fc3, - const Array1D *s2p_fc3, - const Array2D *multi_fc2, - const ShortestVecs *svecs_fc2, - const Array2D *multi_fc3, - const ShortestVecs *svecs_fc3, - const Array1D *band_indices, - const double *born, - const double *dielectric, - const double *reciprocal_lattice, - const double *q_direction, - const double nac_factor, - const double freq_unit_factor, - const double cutoff_frequency, - const int is_symmetrize_fc3_q, - const int r2q_TI_index) -{ - int i, j, num_triplets, num_patom, num_grid_points; - double *q_vecs; - double *w0 ,*w; - lapack_complex_double *a0, *a; - - num_patom = p2s_fc2->d1; - num_triplets = weights->d1; - num_grid_points = 0; - for (i = 0; i < weights->d1; i++) { - num_grid_points += weights->data[i]; - } - - w0 = (double*)malloc(sizeof(double) * num_patom * 3); - a0 = (lapack_complex_double*) - malloc(sizeof(lapack_complex_double) * num_patom * num_patom * 9); - - get_phonons(a0, - w0, - q0, - fc2, - masses_fc2, - p2s_fc2, - s2p_fc2, - multi_fc2, - svecs_fc2, - born, - dielectric, - reciprocal_lattice, - q_direction, - nac_factor); - -#pragma omp parallel for private(j, w, a, q_vecs) - for (i = 0; i < num_triplets; i++) { - w = (double*)malloc(3 * sizeof(double) * num_patom * 3); - a = (lapack_complex_double*) - malloc(3 * sizeof(lapack_complex_double) * num_patom * num_patom * 9); - q_vecs = (double*)malloc(3 * sizeof(double) * 3); - - for (j = 0; j < num_patom * num_patom * 9; j++) { - a[j] = a0[j]; - } - for (j = 0; j < num_patom * 3; j++) { - w[j] = w0[j]; - } - for (j = 0; j < 3; j++) { - q_vecs[j] = q0[j]; - q_vecs[j + 3] = q1s[i * 3 + j]; - q_vecs[j + 6] = q2s[i * 3 + j]; - } - - get_phonons(a + num_patom * num_patom * 9, - w + num_patom * 3, - q_vecs + 3, - fc2, - masses_fc2, - p2s_fc2, - s2p_fc2, - multi_fc2, - svecs_fc2, - born, - dielectric, - reciprocal_lattice, - NULL, - nac_factor); - - get_phonons(a + num_patom * num_patom * 18, - w + num_patom * 6, - q_vecs + 6, - fc2, - masses_fc2, - p2s_fc2, - s2p_fc2, - multi_fc2, - svecs_fc2, - born, - dielectric, - reciprocal_lattice, - NULL, - nac_factor); - - for (j = 0; j < num_patom * 9; j++) { - freqs[i * num_patom * 9 + j] = - sqrt(fabs(w[j])) * freq_unit_factor * ((w[j] > 0) - (w[j] < 0)); - } - - - get_triplet_interaction_strength - (amps + i * band_indices->d1 * num_patom * num_patom * 9, - fc3, - q_vecs, - a, - freqs + i * num_patom * 9, - masses_fc3, - p2s_fc3, - s2p_fc3, - multi_fc3, - svecs_fc3, - band_indices, - cutoff_frequency, - is_symmetrize_fc3_q, - r2q_TI_index); - - free(q_vecs); - free(w); - free(a); - } - -#pragma omp parallel for - for (i = 0; - i < num_triplets * band_indices->d1 * num_patom * num_patom * 9; - i++) { - amps[i] = lapack_make_complex_double - (lapack_complex_double_real(amps[i]) / num_grid_points, - lapack_complex_double_imag(amps[i]) / num_grid_points); + amps[i] /= num_grid_points; } free(w0); @@ -412,9 +259,9 @@ int get_triplet_interaction_strength(double *amps, int get_fc3_realspace(lapack_complex_double* fc3_real, const ShortestVecs* svecs, - const Array2D * multi, + const Array2D* multi, const double* q_triplet, - const Array1D * s2p, + const int* s2p, const lapack_complex_double* fc3_rec) { int i, j, k, l, m, n; @@ -448,8 +295,8 @@ int get_fc3_realspace(lapack_complex_double* fc3_real, for (m = 0; m < 3; m++) { for (n = 0; n < 3; n++) { fc3_elem = prod(fc3_rec[i * num_patom * num_patom * 27 + - s2p->data[j] * num_patom * 27 + - s2p->data[k] * 27 + l * 9 + m * 3 + n], + s2p[j] * num_patom * 27 + + s2p[k] * 27 + l * 9 + m * 3 + n], prod(phase2, phase3)); fc3_real[i * num_satom * num_satom * 27 + @@ -793,7 +640,8 @@ static int get_phonons(lapack_complex_double *a, const double *dielectric, const double *reciprocal_lattice, const double *q_direction, - const double nac_factor) + const double nac_factor, + const char uplo) { int i, j, num_patom, num_satom; double q_cart[3]; @@ -845,7 +693,7 @@ static int get_phonons(lapack_complex_double *a, num_patom, dielectric_factor, q_cart, - born); + born); } } else { charge_sum = NULL; @@ -877,5 +725,5 @@ static int get_phonons(lapack_complex_double *a, free(dm_real); free(dm_imag); - return phonopy_zheev(w, a, num_patom * 3); + return phonopy_zheev(w, a, num_patom * 3, uplo); } diff --git a/c/anharmonic_h/interaction_strength.h b/c/anharmonic_h/interaction_strength.h index 7693eee8..6e0e5b11 100644 --- a/c/anharmonic_h/interaction_strength.h +++ b/c/anharmonic_h/interaction_strength.h @@ -4,40 +4,13 @@ #include #include "alloc_array.h" -int get_interaction_strength2(double *amps, - double *freqs, - const double *q0, - const double *q1s, - const double *q2s, - const Array1D *weights, - const double *fc2, - const double *fc3, - const double *masses_fc2, - const double *masses_fc3, - const Array1D *p2s_fc2, - const Array1D *s2p_fc2, - const Array1D *p2s_fc3, - const Array1D *s2p_fc3, - const Array2D *multi_fc2, - const ShortestVecs *svecs_fc2, - const Array2D *multi_fc3, - const ShortestVecs *svecs_fc3, - const Array1D *band_indices, - const double *born, - const double *dielectric, - const double *reciprocal_lattice, - const double *q_direction, - const double nac_factor, - const double freq_unit_factor, - const double cutoff_frequency, - const int is_symmetrize_fc3_q, - const int r2q_TI_index); int get_interaction_strength(double *amps, double *freqs, + const int num_triplets, + const int num_grid_points, const double *q0, const double *q1s, const double *q2s, - const Array1D *weights, const double *fc2, const double *fc3, const double *masses_fc2, @@ -59,7 +32,8 @@ int get_interaction_strength(double *amps, const double freq_unit_factor, const double cutoff_frequency, const int is_symmetrize_fc3_q, - const int r2q_TI_index); + const int r2q_TI_index, + const char uplo); int get_triplet_interaction_strength(double *amps, const double *fc3, const double *q_vecs, @@ -78,7 +52,7 @@ int get_fc3_realspace(lapack_complex_double* fc3_real, const ShortestVecs * svecs, const Array2D * multi, const double* q_triplet, - const Array1D * s2p, + const int * s2p, const lapack_complex_double* fc3_rec); int get_fc3_reciprocal(lapack_complex_double* fc3_q, const ShortestVecs * svecs, diff --git a/c/harmonic/lapack_wrapper.c b/c/harmonic/lapack_wrapper.c index 2228a364..21c24b8f 100644 --- a/c/harmonic/lapack_wrapper.c +++ b/c/harmonic/lapack_wrapper.c @@ -3,10 +3,11 @@ #include int phonopy_zheev(double *w, lapack_complex_double *a, - const int n) + const int n, + const char uplo) { lapack_int info; - info = LAPACKE_zheev(LAPACK_ROW_MAJOR,'V', 'U', + info = LAPACKE_zheev(LAPACK_ROW_MAJOR,'V', uplo, (lapack_int)n, a, (lapack_int)n, w); return (int)info; } diff --git a/c/harmonic_h/lapack_wrapper.h b/c/harmonic_h/lapack_wrapper.h index e95d6f9d..4ff7570a 100644 --- a/c/harmonic_h/lapack_wrapper.h +++ b/c/harmonic_h/lapack_wrapper.h @@ -4,6 +4,7 @@ #include int phonopy_zheev(double *w, lapack_complex_double *a, - const int n); + const int n, + const char uplo); #endif diff --git a/phonopy/harmonic/dynamical_matrix.py b/phonopy/harmonic/dynamical_matrix.py index b9b50305..11156758 100644 --- a/phonopy/harmonic/dynamical_matrix.py +++ b/phonopy/harmonic/dynamical_matrix.py @@ -67,11 +67,11 @@ class DynamicalMatrix: self._decimals = decimals self._symprec = symprec - self._p2s_map = np.array(primitive.get_primitive_to_supercell_map()) - self._s2p_map = np.array(primitive.get_supercell_to_primitive_map()) + self._p2s_map = np.int32(primitive.get_primitive_to_supercell_map()) + self._s2p_map = np.int32(primitive.get_supercell_to_primitive_map()) p2p_map = primitive.get_primitive_to_primitive_map() - self._p2p_map = np.array([p2p_map[self._s2p_map[i]] - for i in range(len(self._s2p_map))]) + self._p2p_map = [p2p_map[self._s2p_map[i]] + for i in range(len(self._s2p_map))] self._smallest_vectors, self._multiplicity = \ get_smallest_vectors(supercell, primitive, symprec) self._mass = self._pcell.get_masses() @@ -185,7 +185,7 @@ class DynamicalMatrix: for p_i, s_i in enumerate(self._p2s_map): # run in primitive for s_j in range(r.shape[0]): # run in supercell for tmp_p_j, tmp_s_j in enumerate(self._p2s_map): - if self._s2p_map[ s_j ] == tmp_s_j: + if self._s2p_map[s_j] == tmp_s_j: p_j = tmp_p_j for k in range(m[s_j][p_i]): print " %4d %4d %4d %4d %4d %10.5f" % \ @@ -204,8 +204,7 @@ class DynamicalMatrix: size_super = fc.shape[0] dynamical_matrix_real = np.zeros((size_prim * 3, size_prim * 3), dtype=float) - dynamical_matrix_image = np.zeros((size_prim * 3, size_prim * 3), - dtype=float) + dynamical_matrix_image = np.zeros_like(dynamical_matrix_real) phonoc.dynamical_matrix(dynamical_matrix_real, dynamical_matrix_image, fc, @@ -295,8 +294,8 @@ class DynamicalMatrixNAC(DynamicalMatrix): self._damping_factor ** 2) for i in range(num_atom): for j in range(num_atom): - nac_q[ i*3:(i+1)*3, j*3:(j+1)*3 ] = \ - charge_sum[ i, j ] * constant / np.sqrt(m[i] * m[j]) + nac_q[i*3:(i+1)*3, j*3:(j+1)*3] = \ + charge_sum[i, j] * constant / np.sqrt(m[i] * m[j]) DynamicalMatrix.set_dynamical_matrix(self, q_red, verbose) self._dynamical_matrix += nac_q @@ -313,8 +312,8 @@ class DynamicalMatrixNAC(DynamicalMatrix): nac_q = np.zeros((num_atom * 3, num_atom * 3), dtype=float) for i in range(num_atom): for j in range(num_atom): - nac_q[ i*3:(i+1)*3, j*3:(j+1)*3 ] = \ - charge_sum[ i, j ] * constant + nac_q[i*3:(i+1)*3, j*3:(j+1)*3] = \ + charge_sum[i, j] * constant self._set_NAC_force_constants(fc, nac_q) self._force_constants = fc DynamicalMatrix.set_dynamical_matrix(self, q_red, verbose) @@ -327,12 +326,12 @@ class DynamicalMatrixNAC(DynamicalMatrix): # In contructing dynamical matrix in phonopy # fc of left indices with s1 == self._s2p_map[ s1 ] are # only used. - if not (s1==self._s2p_map[ s1 ]): + if not (s1==self._s2p_map[s1]): continue p1 = self._p2p_map[s1] for s2 in range(self._scell.get_number_of_atoms()): p2 = self._p2p_map[s2] - fc[ s1, s2 ] += nac_q[ p1*3:(p1+1)*3, p2*3:(p2+1)*3 ] / N + fc[s1, s2] += nac_q[p1*3:(p1+1)*3, p2*3:(p2+1)*3] / N def _get_charge_sum(self, num_atom, q): charge_sum = np.zeros((num_atom, num_atom, 3, 3), dtype=float) @@ -340,7 +339,7 @@ class DynamicalMatrixNAC(DynamicalMatrix): for j in range(num_atom): for a in (0, 1, 2): for b in (0, 1, 2): - charge_sum[ i, j, a, b ] = \ + charge_sum[i, j, a, b] = \ np.dot(q, self._born[i, :, a]) * np.dot(q, self._born[j, :, b]) return charge_sum @@ -417,7 +416,7 @@ def get_smallest_vectors(supercell, primitive, symprec): size_super = supercell.get_number_of_atoms() size_prim = primitive.get_number_of_atoms() r = np.zeros((size_super, size_prim, 27, 3), dtype=float) - multiplicity = np.zeros((size_super, size_prim), dtype=int) + multiplicity = np.zeros((size_super, size_prim), dtype='int32') for i in range(size_super): # run in supercell for j, s_j in enumerate(p2s_map): # run in primitive diff --git a/phonopy/harmonic/force_constants.py b/phonopy/harmonic/force_constants.py index 8c2d3091..11b270ef 100644 --- a/phonopy/harmonic/force_constants.py +++ b/phonopy/harmonic/force_constants.py @@ -166,8 +166,8 @@ def distribute_force_constants(force_constants, map_atom_disp, i, rot_cartesian, - rotations[ map_sym ], - trans[ map_sym ], + rotations[map_sym], + trans[map_sym], symprec) force_constants[atom_disp, i] /= len(map_atom_disps) @@ -348,7 +348,8 @@ def solve_force_constants_disps(force_constants, try: import phonopy._phonopy as phonoc for forces in sets_of_forces: - rotated_forces = np.zeros(len(site_symmetry) * 3) + rotated_forces = np.zeros(len(site_symmetry) * 3, + dtype=float) phonoc.rotated_forces(rotated_forces, positions, i, diff --git a/phonopy/phonon/mesh.py b/phonopy/phonon/mesh.py index 218e1146..eae2d367 100644 --- a/phonopy/phonon/mesh.py +++ b/phonopy/phonon/mesh.py @@ -79,7 +79,7 @@ def _get_qpoint_symmetry(mesh, is_time_reversal, symprec) ir_list = np.unique(mapping) - weights = np.zeros(ir_list.shape[0], dtype=int) + weights = np.zeros(ir_list.shape[0], dtype='int32') qpoints = np.zeros((ir_list.shape[0], 3), dtype=float) for i, g in enumerate(ir_list): @@ -107,7 +107,7 @@ def _get_qpoint_no_symmetry(mesh, shift): q[2] - (q[2] > 0.5)])) qpoints = np.array(qpoints) - weights = np.ones(qpoints.shape[0], dtype=int) + weights = np.ones(qpoints.shape[0], dtype='int32') return qpoints, weights diff --git a/phonopy/phonon/thermal_properties.py b/phonopy/phonon/thermal_properties.py index fe47c1df..aa73b39d 100644 --- a/phonopy/phonon/thermal_properties.py +++ b/phonopy/phonon/thermal_properties.py @@ -40,7 +40,7 @@ class ThermalPropertiesBase: self._temperature = 0 self._frequencies = frequencies if weights == None: - self._weights = np.ones(frequencies.shape[0], dtype=int) + self._weights = np.ones(frequencies.shape[0], dtype='int32') else: self._weights = weights self._nqpoint = frequencies.shape[0] diff --git a/phonopy/structure/cells.py b/phonopy/structure/cells.py index b45e81e2..acdbbeb0 100644 --- a/phonopy/structure/cells.py +++ b/phonopy/structure/cells.py @@ -256,12 +256,12 @@ class Primitive(Atoms): inv_F = np.linalg.inv(self.frame) s2p_map = [] for i in range(pos.shape[0]): - s_pos = np.dot( pos[i], inv_F.T ) + s_pos = np.dot(pos[i], inv_F.T) for j in self.p2s_map: - p_pos = np.dot( pos[j], inv_F.T ) + p_pos = np.dot(pos[j], inv_F.T) diff = p_pos - s_pos diff -= diff.round() - if ( abs(diff) < self.symprec ).all(): + if (abs(diff) < self.symprec).all(): s2p_map.append(j) break self.s2p_map = s2p_map diff --git a/phonopy/structure/spglib.py b/phonopy/structure/spglib.py index d78b5615..aaf92350 100644 --- a/phonopy/structure/spglib.py +++ b/phonopy/structure/spglib.py @@ -45,8 +45,8 @@ def get_symmetry(bulk, symprec=1e-5, angle_tolerance=-1.0): symprec, angle_tolerance) - return {'rotations': rotation[:num_sym], - 'translations': translation[:num_sym]} + return {'rotations': np.int32(rotation[:num_sym].copy()), + 'translations': np.int32(translation[:num_sym].copy())} def get_symmetry_dataset(bulk, symprec=1e-5, angle_tolerance=-1.0): """ @@ -88,11 +88,11 @@ def get_symmetry_dataset(bulk, symprec=1e-5, angle_tolerance=-1.0): dataset['hall'] = dataset['hall'].strip() dataset['transformation_matrix'] = np.array(dataset['transformation_matrix']) dataset['origin_shift'] = np.array(dataset['origin_shift']) - dataset['rotations'] = np.array(dataset['rotations']) + dataset['rotations'] = np.int32(dataset['rotations']) dataset['translations'] = np.array(dataset['translations']) letters = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ" dataset['wyckoffs'] = [letters[x] for x in dataset['wyckoffs']] - dataset['equivalent_atoms'] = np.array(dataset['equivalent_atoms']) + dataset['equivalent_atoms'] = np.int32(dataset['equivalent_atoms']) return dataset diff --git a/phonopy/structure/symmetry.py b/phonopy/structure/symmetry.py index 82de815c..56b8af99 100644 --- a/phonopy/structure/symmetry.py +++ b/phonopy/structure/symmetry.py @@ -64,7 +64,7 @@ class Symmetry: self.map_atoms = None if not is_symmetry: self._set_nosym() - elif cell.get_magnetic_moments() == None: + elif cell.get_magnetic_moments() is None: self._symmetry_dataset() else: self._symmetry_operations() @@ -85,62 +85,12 @@ class Symmetry: return {'rotations': operation['rotations'][operation_number], 'translations': operation['translations'][operation_number]} - def _map_atoms(self): - rotations = self.symmetry_operations['rotations'] - translations = self.symmetry_operations['translations'] - positions = self.__cell.get_scaled_positions() - lattice = self.__cell.get_cell() - map_atoms = range(self.__cell.get_number_of_atoms()) - for i, p in enumerate(positions): - is_found = False - for j in range(i): - for r, t in zip(rotations, translations): - diff = np.dot(p, r.T) + t - positions[j] - diff -= diff.round() - dist = np.linalg.norm(np.dot(diff, lattice)) - if dist < self.symprec: - map_atoms[i] = j - is_found = True - break - if is_found: - break - self.map_atoms = np.array(map_atoms, dtype=int) - - def _symmetry_operations(self): - self.symmetry_operations = \ - spg.get_symmetry(self.__cell, self.symprec) - self._map_atoms() - - def _pointgroup_operations(self): - rotations = [] - for rot in self.symmetry_operations['rotations']: - is_same = False - for tmp_rot in rotations: - if (tmp_rot==rot).all(): - is_same = True - break - if not is_same: - rotations.append(rot) - - self.pointgroup_operations = np.array(rotations) - self.pointgroup = get_pointgroup(self.pointgroup_operations)[0] - def get_pointgroup_operations(self): return self.pointgroup_operations def get_pointgroup(self): return self.pointgroup - def _symmetry_dataset(self): - self.dataset = spg.get_symmetry_dataset(self.__cell, self.symprec) - self.symmetry_operations = \ - {'rotations': self.dataset['rotations'], - 'translations': self.dataset['translations']} - self.international_table = "%s (%d)" % (self.dataset['international'], - self.dataset['number']) - self.wyckoff_letters = self.dataset['wyckoffs'] - self.map_atoms = self.dataset['equivalent_atoms'] - def get_international_table(self): return self.international_table @@ -152,26 +102,6 @@ class Symmetry: """ return self.dataset - def _map_operations(self): - ops = self.symmetry_operations - pos = self.__cell.get_scaled_positions() - map_operations = np.zeros(len(pos), dtype=int) - independent_atoms = [] - - for i, eq_atom in enumerate(self.map_atoms): - if i == eq_atom: - independent_atoms.append(i) - for j, (r, t) in enumerate( - zip(ops['rotations'], ops['translations'])): - - diff = np.dot(pos[i], r.T) + t - pos[eq_atom] - if (abs(diff - diff.round()) < self.symprec).all(): - map_operations[i] = j - break - - self.independent_atoms = np.array(independent_atoms) - self.map_operations = map_operations - def get_independent_atoms(self): return self.independent_atoms @@ -194,11 +124,81 @@ class Symmetry: if (abs(diff - diff.round()) < symprec).all(): site_symmetries.append(r) - return np.array(site_symmetries) + return np.int32(site_symmetries) def get_symmetry_tolerance(self): return self.symprec + def _map_atoms(self): + rotations = self.symmetry_operations['rotations'] + translations = self.symmetry_operations['translations'] + positions = self.__cell.get_scaled_positions() + lattice = self.__cell.get_cell() + map_atoms = range(self.__cell.get_number_of_atoms()) + for i, p in enumerate(positions): + is_found = False + for j in range(i): + for r, t in zip(rotations, translations): + diff = np.dot(p, r.T) + t - positions[j] + diff -= diff.round() + dist = np.linalg.norm(np.dot(diff, lattice)) + if dist < self.symprec: + map_atoms[i] = j + is_found = True + break + if is_found: + break + self.map_atoms = np.array(map_atoms, dtype=int) + + def _symmetry_dataset(self): + self.dataset = spg.get_symmetry_dataset(self.__cell, self.symprec) + self.symmetry_operations = \ + {'rotations': self.dataset['rotations'], + 'translations': self.dataset['translations']} + self.international_table = "%s (%d)" % (self.dataset['international'], + self.dataset['number']) + self.wyckoff_letters = self.dataset['wyckoffs'] + self.map_atoms = self.dataset['equivalent_atoms'] + + def _symmetry_operations(self): + self.symmetry_operations = \ + spg.get_symmetry(self.__cell, self.symprec) + self._map_atoms() + + def _pointgroup_operations(self): + rotations = [] + for rot in self.symmetry_operations['rotations']: + is_same = False + for tmp_rot in rotations: + if (tmp_rot==rot).all(): + is_same = True + break + if not is_same: + rotations.append(rot) + + self.pointgroup_operations = np.array(rotations) + self.pointgroup = get_pointgroup(self.pointgroup_operations)[0] + + def _map_operations(self): + ops = self.symmetry_operations + pos = self.__cell.get_scaled_positions() + map_operations = np.zeros(len(pos), dtype=int) + independent_atoms = [] + + for i, eq_atom in enumerate(self.map_atoms): + if i == eq_atom: + independent_atoms.append(i) + for j, (r, t) in enumerate( + zip(ops['rotations'], ops['translations'])): + + diff = np.dot(pos[i], r.T) + t - pos[eq_atom] + if (abs(diff - diff.round()) < self.symprec).all(): + map_operations[i] = j + break + + self.independent_atoms = np.array(independent_atoms) + self.map_operations = map_operations + def _set_nosym(self): translations = [] rotations = [] @@ -217,7 +217,7 @@ class Symmetry: trans = p - positions[ipos0] trans -= np.floor(trans) translations.append(trans) - rotations.append(np.eye(3, dtype=int)) + rotations.append(np.eye(3, dtype='int32')) self.map_atoms = s2u_map else: @@ -225,8 +225,9 @@ class Symmetry: translations.append(np.zeros(3, dtype=float)) self.map_atoms = range(self.__cell.get_number_of_atoms()) - self.symmetry_operations = {'rotations': np.array(rotations), - 'translations': np.array(translations)} + self.symmetry_operations = {'rotations': np.int32(rotations), + 'translations': np.array(translations, + dtype=float)} self.international_table = 'P1 (1)' self.wyckoff_letters = ['a'] * self.__cell.get_number_of_atoms() diff --git a/scripts/phono3py b/scripts/phono3py index 6cc68153..a45640ad 100755 --- a/scripts/phono3py +++ b/scripts/phono3py @@ -139,6 +139,7 @@ parser.set_defaults(amplitude=None, tstep=None, temperatures=None, verbose=True, + uplo='L', write_amplitude=False, write_gamma=False, write_grid_points=False) @@ -307,6 +308,10 @@ parser.add_option("--loglevel", dest="log_level", type="int", help="Log level") parser.add_option("--ts", dest="temperatures", type="string", help="Temperatures for damping functions") +parser.add_option("--uplo", + dest="uplo", + type="string", + help="Lapack zheev UPLO") parser.add_option("--wgp", "--write_grid_points", dest="write_grid_points", action="store_true", help="Write grid address of irreducible grid points for specified mesh numbers to ir_grid_address.yaml") @@ -749,7 +754,8 @@ else: r2q_TI_index=options.r2q_TI_index, is_Peierls=options.is_Peierls, symprec=options.symprec, - log_level=log_level) + log_level=log_level, + lapack_zheev_uplo=options.uplo) phono3py.set_dynamical_matrix(fc2, supercell_dm,