C implementation of group velocity

This commit is contained in:
Atsushi Togo 2014-03-22 22:57:25 +01:00
parent c079bf38eb
commit f62a273559
8 changed files with 500 additions and 20 deletions

View File

@ -36,12 +36,14 @@
#include <stdio.h>
#include <numpy/arrayobject.h>
#include "dynmat.h"
#include "derivative_dynmat.h"
#define KB 8.6173382568083159E-05
/* Build dynamical matrix */
static PyObject * py_get_dynamical_matrix(PyObject *self, PyObject *args);
static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args);
static PyObject * py_get_derivative_dynmat(PyObject *self, PyObject *args);
static PyObject * py_get_thermal_properties(PyObject *self, PyObject *args);
static PyObject * py_distribute_fc2(PyObject *self, PyObject *args);
@ -66,6 +68,7 @@ static int nint(const double a);
static PyMethodDef functions[] = {
{"dynamical_matrix", py_get_dynamical_matrix, METH_VARARGS, "Dynamical matrix"},
{"nac_dynamical_matrix", py_get_nac_dynamical_matrix, METH_VARARGS, "NAC dynamical matrix"},
{"derivative_dynmat", py_get_derivative_dynmat, METH_VARARGS, "Q derivative of dynamical matrix"},
{"thermal_properties", py_get_thermal_properties, METH_VARARGS, "Thermal properties"},
{"distribute_fc2", py_distribute_fc2, METH_VARARGS, "Distribute force constants"},
{NULL, NULL, 0, NULL}
@ -199,6 +202,92 @@ static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_get_derivative_dynmat(PyObject *self, PyObject *args)
{
PyArrayObject* derivative_dynmat_real;
PyArrayObject* derivative_dynmat_imag;
PyArrayObject* force_constants;
PyArrayObject* r_vector;
PyArrayObject* lattice;
PyArrayObject* q_vector;
PyArrayObject* multiplicity;
PyArrayObject* mass;
PyArrayObject* super2prim_map;
PyArrayObject* prim2super_map;
PyArrayObject* born;
PyArrayObject* dielectric;
PyArrayObject* q_direction;
double nac_factor;
if (!PyArg_ParseTuple(args, "OOOOOOOOOOdOOO",
&derivative_dynmat_real,
&derivative_dynmat_imag,
&force_constants,
&q_vector,
&lattice, /* column vectors */
&r_vector,
&multiplicity,
&mass,
&super2prim_map,
&prim2super_map,
&nac_factor,
&born,
&dielectric,
&q_direction)) {
return NULL;
}
double* ddm_r = (double*)derivative_dynmat_real->data;
double* ddm_i = (double*)derivative_dynmat_imag->data;
const double* fc = (double*)force_constants->data;
const double* q = (double*)q_vector->data;
const double* lat = (double*)lattice->data;
const double* r = (double*)r_vector->data;
const double* m = (double*)mass->data;
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];
double *z;
double *epsilon;
double *q_dir;
if ((PyObject*)born == Py_None) {
z = NULL;
} else {
z = (double*)born->data;
}
if ((PyObject*)dielectric == Py_None) {
epsilon = NULL;
} else {
epsilon = (double*)dielectric->data;
}
if ((PyObject*)q_direction == Py_None) {
q_dir = NULL;
} else {
q_dir = (double*)q_direction->data;
}
get_derivative_dynmat_at_q(ddm_r,
ddm_i,
num_patom,
num_satom,
fc,
q,
lat,
r,
multi,
m,
s2p_map,
p2s_map,
nac_factor,
z,
epsilon,
q_dir);
Py_RETURN_NONE;
}
/* Thermal properties */
static PyObject * py_get_thermal_properties(PyObject *self, PyObject *args)
{

View File

@ -0,0 +1,323 @@
#include <math.h>
#include <stdlib.h>
static void get_derivative_nac(double *ddnac,
double *dnac,
const int num_patom,
const double *lattice,
const double *mass,
const double *q,
const double *born,
const double *dielectric,
const double *q_direction,
const double factor);
static double get_A(const int atom_i,
const int cart_i,
const double q[3],
const double *born);
static double get_C(const double q[3],
const double *dielectric);
static double get_dA(const int atom_i,
const int cart_i,
const int cart_j,
const double *born);
static double get_dC(const int cart_i,
const int cart_j,
const int cart_k,
const double q[3],
const double *dielectric);
void get_derivative_dynmat_at_q(double *derivative_dynmat_real,
double *derivative_dynmat_imag,
const int num_patom,
const int num_satom,
const double *fc,
const double *q,
const double *lattice, /* column vector */
const double *r,
const int *multi,
const double *mass,
const int *s2p_map,
const int *p2s_map,
const double nac_factor,
const double *born,
const double *dielectric,
const double *q_direction)
{
int i, j, k, l, m, n, is_nac;
double coef[3], real_coef[3], imag_coef[3];
double c, s, phase, mass_sqrt, fc_elem, factor, real_phase, imag_phase;
double ddm_real[3][3][3], ddm_imag[3][3][3];
double *ddnac, *dnac;
if (born) {
is_nac = 1;
if (q_direction) {
if (fabs(q_direction[0]) < 1e-5 &&
fabs(q_direction[1]) < 1e-5 &&
fabs(q_direction[2]) < 1e-5) {
is_nac = 0;
}
} else {
if (fabs(q[0]) < 1e-5 &&
fabs(q[1]) < 1e-5 &&
fabs(q[2]) < 1e-5) {
is_nac = 0;
}
}
} else {
is_nac = 0;
}
if (is_nac) {
ddnac = (double*) malloc(sizeof(double) * num_patom * num_patom * 27);
dnac = (double*) malloc(sizeof(double) * num_patom * num_patom * 9);
factor = nac_factor * num_patom / num_satom;
get_derivative_nac(ddnac,
dnac,
num_patom,
lattice,
mass,
q,
born,
dielectric,
q_direction,
factor);
}
for (i = 0; i < num_patom; i++) {
for (j = 0; j < num_patom; j++) {
mass_sqrt = sqrt(mass[i] * mass[j]);
for (k = 0; k < 3; k++) {
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
ddm_real[m][k][l] = 0;
ddm_imag[m][k][l] = 0;
}
}
}
for (k = 0; k < num_satom; k++) { /* Lattice points of right index of fc */
if (s2p_map[k] != p2s_map[j]) {
continue;
}
real_phase = 0;
imag_phase = 0;
for (l = 0; l < 3; l++) {
real_coef[l] = 0;
imag_coef[l] = 0;
}
for (l = 0; l < multi[k * num_patom + i]; l++) {
phase = 0;
for (m = 0; m < 3; m++) {
phase += q[m] * r[k * num_patom * 81 + i * 81 + l * 3 + m];
}
s = sin(phase * 2 * M_PI);
c = cos(phase * 2 * M_PI);
real_phase += c;
imag_phase += s;
for (m = 0; m < 3; m++) {
coef[m] = 0;
for (n = 0; n < 3; n++) {
coef[m] += 2 * M_PI *
lattice[m * 3 + n] * r[k * num_patom * 81 + i * 81 + l * 3 + n];
}
}
for (m =0; m < 3; m++) {
real_coef[m] -= coef[m] * s;
imag_coef[m] += coef[m] * c;
}
}
real_phase /= multi[k * num_patom + i];
imag_phase /= multi[k * num_patom + i];
for (l = 0; l < 3; l++) {
real_coef[l] /= multi[k * num_patom + i];
imag_coef[l] /= multi[k * num_patom + i];
}
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
fc_elem = fc[p2s_map[i] * num_satom * 9 + k * 9 + l * 3 + m]
/ mass_sqrt;
if (is_nac) {
fc_elem += dnac[i * 9 * num_patom + j * 9 + l * 3 + m];
}
for (n = 0; n < 3; n++) {
ddm_real[n][l][m] += fc_elem * real_coef[n];
ddm_imag[n][l][m] += fc_elem * imag_coef[n];
if (is_nac) {
ddm_real[n][l][m] +=
ddnac[n * num_patom * num_patom * 9 +
i * 9 * num_patom + j * 9 + l * 3 + m] * real_phase;
ddm_imag[n][l][m] +=
ddnac[n * num_patom * num_patom * 9 +
i * 9 * num_patom + j * 9 + l * 3 + m] * imag_phase;
}
}
}
}
}
for (k = 0; k < 3; k++) {
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
derivative_dynmat_real
[k * num_patom * num_patom * 9 +
(i * 3 + l) * num_patom * 3 + j * 3 + m] +=
ddm_real[k][l][m];
derivative_dynmat_imag
[k * num_patom * num_patom * 9 +
(i * 3 + l) * num_patom * 3 + j * 3 + m] +=
ddm_imag[k][l][m];
}
}
}
}
}
if (is_nac) {
free(ddnac);
free(dnac);
}
}
/* D_nac = a * AB/C */
/* dD_nac = a * D_nac * (A'/A + B'/B - C'/C) */
static void get_derivative_nac(double *ddnac,
double *dnac,
const int num_patom,
const double *lattice,
const double *mass,
const double *q,
const double *born,
const double *dielectric,
const double *q_direction,
const double factor)
{
int i, j, k, l, m;
double a, b, c, da, db, dc, volume, mass_sqrt;
double q_cart[3], rec_lat[9];
volume =
lattice[0] * (lattice[4] * lattice[8] - lattice[5] * lattice[7]) +
lattice[1] * (lattice[5] * lattice[6] - lattice[3] * lattice[8]) +
lattice[2] * (lattice[3] * lattice[7] - lattice[4] * lattice[6]);
rec_lat[0] = (lattice[4] * lattice[8] - lattice[5] * lattice[7]) / volume;
rec_lat[1] = (lattice[5] * lattice[6] - lattice[3] * lattice[8]) / volume;
rec_lat[2] = (lattice[3] * lattice[7] - lattice[4] * lattice[6]) / volume;
rec_lat[3] = (lattice[7] * lattice[2] - lattice[8] * lattice[1]) / volume;
rec_lat[4] = (lattice[8] * lattice[0] - lattice[6] * lattice[2]) / volume;
rec_lat[5] = (lattice[6] * lattice[1] - lattice[7] * lattice[0]) / volume;
rec_lat[6] = (lattice[1] * lattice[5] - lattice[2] * lattice[4]) / volume;
rec_lat[7] = (lattice[2] * lattice[3] - lattice[0] * lattice[5]) / volume;
rec_lat[8] = (lattice[0] * lattice[4] - lattice[1] * lattice[3]) / volume;
for (i = 0; i < 3; i++) {
q_cart[i] = 0;
for (j = 0; j < 3; j++) {
if (q_direction) {
q_cart[i] += rec_lat[i * 3 + j] * q_direction[j];
} else {
q_cart[i] += rec_lat[i * 3 + j] * q[j];
}
}
}
c = get_C(q_cart, dielectric);
for (i = 0; i < num_patom; i++) { /* atom_i */
for (j = 0; j < num_patom; j++) { /* atom_j */
mass_sqrt = sqrt(mass[i] * mass[j]);
for (k = 0; k < 3; k++) { /* derivative direction */
for (l = 0; l < 3; l++) { /* alpha */
a = get_A(i, l, q_cart, born);
da = get_dA(i, l, k, born);
for (m = 0; m < 3; m++) { /* beta */
b = get_A(j, m, q_cart, born);
db = get_dA(j, m, k, born);
dc = get_dC(l, m, k, q_cart, dielectric);
ddnac[k * num_patom * num_patom * 9 +
i * 9 * num_patom + j * 9 + l * 3 + m] =
(da * b + db * a - a * b * dc / c) / (c * mass_sqrt) * factor;
if (k == 0) {
dnac[i * 9 * num_patom + j * 9 + l * 3 + m] =
a * b / (c * mass_sqrt) * factor;
}
}
}
}
}
}
}
static double get_A(const int atom_i,
const int cart_i,
const double q[3],
const double *born)
{
int i;
double sum;
sum = 0;
for (i = 0; i < 3; i++) {
sum += q[i] * born[atom_i * 9 + i * 3 + cart_i];
}
return sum;
}
static double get_C(const double q[3],
const double *dielectric)
{
int i, j;
double sum;
sum = 0;
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
sum += q[i] * dielectric[i * 3 + j] * q[j];
}
}
return sum;
}
static double get_dA(const int atom_i,
const int cart_i,
const int cart_j,
const double *born)
{
return born[atom_i * 9 + cart_j * 3 + cart_i];
}
static double get_dC(const int cart_i,
const int cart_j,
const int cart_k,
const double q[3],
const double *dielectric)
{
if (cart_k == 0) {
return (2 * q[0] * dielectric[0] +
q[1] * (dielectric[1] + dielectric[3]) +
q[2] * (dielectric[2] + dielectric[6]));
}
if (cart_k == 1) {
return (2 * q[1] * dielectric[4] +
q[2] * (dielectric[5] + dielectric[7]) +
q[0] * (dielectric[1] + dielectric[3]));
}
if (cart_k == 2) {
return (2 * q[2] * dielectric[8] +
q[0] * (dielectric[2] + dielectric[6]) +
q[1] * (dielectric[5] + dielectric[7]));
}
return 0;
}

View File

@ -49,9 +49,9 @@ int get_dynamical_matrix_at_q(double *dynamical_matrix_real,
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
if (charge_sum) {
fc_elem = (fc[p2s_map[i] * num_satom * 9 + k * 9 + l * 3 + m] +
charge_sum[i * num_patom * 9 +
j * 9 + l * 3 + m]) / mass_sqrt;
fc_elem = (fc[p2s_map[i] * num_satom * 9 + k * 9 + l * 3 + m] +
charge_sum[i * num_patom * 9 +
j * 9 + l * 3 + m]) / mass_sqrt;
} else {
fc_elem = fc[p2s_map[i] * num_satom * 9 +
k * 9 + l * 3 + m] / mass_sqrt;

View File

@ -0,0 +1,21 @@
#ifndef __derivative_dynmat_H__
#define __derivative_dynmat_H__
void get_derivative_dynmat_at_q(double *derivative_dynmat_real,
double *derivative_dynmat_imag,
const int num_patom,
const int num_satom,
const double *fc,
const double *q,
const double *lattice, /* column vector */
const double *r,
const int *multi,
const double *mass,
const int *s2p_map,
const int *p2s_map,
const double nac_factor,
const double *born,
const double *dielectric,
const double *q_direction);
#endif

View File

@ -920,9 +920,7 @@ class Phonopy:
self._irreps.write_yaml(show_irreps=show_irreps)
# Group velocity
def set_group_velocity(self,
q_points=None,
q_length=1e-4):
def set_group_velocity(self, q_length=None):
self._set_dynamical_matrix()
self._group_velocity = GroupVelocity(
self._dynamical_matrix,
@ -930,10 +928,12 @@ class Phonopy:
symmetry=self._primitive_symmetry,
frequency_factor_to_THz=self._factor)
if q_points is not None:
self._group_velocity.set_q_points(q_points)
def get_group_velocity(self, q_point):
def get_group_velocity(self):
return self._group_velocity.get_group_velocity()
def get_group_velocity_at_q(self, q_point):
if self._group_velocity is None:
self.set_group_velocity()
self._group_velocity.set_q_points([q_point])
return self._group_velocity.get_group_velocity()[0]

View File

@ -51,13 +51,59 @@ class DerivativeOfDynamicalMatrix:
self._ddm = None
def run(self, q, q_direction=None):
self._run_py(q, q_direction=q_direction)
try:
import phonopy._phonopy as phonoc
self._run_c(q, q_direction=q_direction)
except ImportError:
self._run_py(q, q_direction=q_direction)
def set_derivative_order(self, order):
self._derivative_order = order
def get_derivative_of_dynamical_matrix(self):
return self._ddm
def _run_c(self, q, q_direction=None):
import phonopy._phonopy as phonoc
num_patom = len(self._p2s_map)
ddm_real = np.zeros((3, num_patom * 3, num_patom * 3), dtype='double')
ddm_imag = np.zeros((3, num_patom * 3, num_patom * 3), dtype='double')
mass = self._pcell.get_masses()
fc = self._force_constants
vectors = self._smallest_vectors
multiplicity = self._multiplicity
if self._dynmat.is_nac():
born = self._dynmat.get_born_effective_charges()
dielectric = self._dynmat.get_dielectric_constant()
nac_factor = self._dynmat.get_nac_factor()
if q_direction is None:
q_dir = None
else:
q_dir = np.array(q_direction, dtype='double', order='C')
else:
born = None
dielectric = None
nac_factor = 0
q_dir = None
phonoc.derivative_dynmat(ddm_real,
ddm_imag,
fc,
np.array(q, dtype='double'),
np.array(self._pcell.get_cell(),
dtype='double', order='C'),
vectors,
multiplicity,
mass,
self._s2p_map,
self._p2s_map,
nac_factor,
born,
dielectric,
q_dir)
self._ddm = np.array([ddm_real[i] + ddm_real[i].T +
1j * (ddm_imag[i] - ddm_imag[i].T)
for i in range(3)]) / 2
def _run_py(self, q, q_direction=None):
if self._dynmat.is_nac():
@ -94,17 +140,17 @@ class DerivativeOfDynamicalMatrix:
vecs_multi_cart = np.dot(vecs_multi, self._pcell.get_cell())
coef = (2j * np.pi * vecs_multi_cart) ** self._derivative_order
if self._dynmat.is_nac():
fc_elem = fc[s_i, k] + fc_nac[i, j]
else:
fc_elem = fc[s_i, k]
for l in range(3):
if self._dynmat.is_nac():
fc_elem = fc[s_i, k] + fc_nac[i, j]
else:
fc_elem = fc[s_i, k]
ddm_elem = fc_elem * (coef[:, l] * phase_multi).sum()
if self._dynmat.is_nac() and self._derivative_order == 1:
ddm_elem += d_nac[l, i, j] * phase_multi.sum()
ddm_local[l] += ddm_elem / multi / mass
ddm_local[l] += ddm_elem / mass / multi
ddm[:, (i * 3):(i * 3 + 3), (j * 3):(j * 3 + 3)] = ddm_local

View File

@ -32,4 +32,4 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
phonopy_version = "1.8.3.2"
phonopy_version = "1.8.4"

View File

@ -10,7 +10,8 @@ extension = Extension('phonopy._phonopy',
# extra_link_args=['-lgomp'],
include_dirs=['c/harmonic_h'] + include_dirs_numpy,
sources=['c/_phonopy.c',
'c/harmonic/dynmat.c'])
'c/harmonic/dynmat.c',
'c/harmonic/derivative_dynmat.c'])
extension_spglib = Extension(
'phonopy._spglib',
@ -38,7 +39,7 @@ extension_spglib = Extension(
setup(name='phonopy',
version='1.8.3.2',
version='1.8.4',
description='This is the phonopy module.',
author='Atsushi Togo',
author_email='atz.togo@gmail.com',