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

Rewriting numpy.int64 to numpy.int32

This commit is contained in:
Atsushi Togo 2013-05-02 19:04:19 +09:00
parent 165aec2d6e
commit c7fce6884a
21 changed files with 397 additions and 633 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
anharmonic/BTE_RTA.py
View File

@ -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

6
anharmonic/__init__.py
View File

@ -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,

180
anharmonic/fc_interpolate.py
View File

@ -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 ):

16
anharmonic/jointDOS.py
View File

@ -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]

40
anharmonic/q2v.py
View File

@ -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(
@ -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):

4
anharmonic/r2q.py
View File

@ -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)

11
anharmonic/triplets.py
View File

@ -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)

137
c/_phono3py.c
View File

@ -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",
&amplitude,
&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*)
@ -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];
}
}

167
c/_phonopy.c
View File

@ -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;
}

196
c/anharmonic/interaction_strength.c
View File

@ -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);
}

36
c/anharmonic_h/interaction_strength.h
View File

@ -4,40 +4,13 @@
#include <lapacke.h>
#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,

5
c/harmonic/lapack_wrapper.c
View File

@ -3,10 +3,11 @@
#include <lapacke.h>
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;
}

3
c/harmonic_h/lapack_wrapper.h
View File

@ -4,6 +4,7 @@
#include <lapacke.h>
int phonopy_zheev(double *w,
lapack_complex_double *a,
const int n);
const int n,
const char uplo);
#endif

29
phonopy/harmonic/dynamical_matrix.py
View File

@ -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

7
phonopy/harmonic/force_constants.py
View File

@ -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,

4
phonopy/phonon/mesh.py
View File

@ -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

2
phonopy/phonon/thermal_properties.py
View File

@ -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]

6
phonopy/structure/cells.py
View File

@ -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

8
phonopy/structure/spglib.py
View File

@ -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

151
phonopy/structure/symmetry.py
View File

@ -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()

8
scripts/phono3py
View File

@ -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,