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

Update Gonze NAC. Even faster.

This commit is contained in:
Atsushi Togo 2018-02-08 15:12:26 +09:00
parent cad967dcef
commit 7e87b61bbd
4 changed files with 181 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

70
c/_phonopy.c
View File

@ -52,6 +52,7 @@ py_perm_trans_symmetrize_compact_fc(PyObject *self, PyObject *args);
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_dipole_dipole(PyObject *self, PyObject *args);
static PyObject * py_get_dipole_dipole_q0(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_with_mappings(PyObject *self,
@ -133,12 +134,20 @@ static PyMethodDef _phonopy_methods[] = {
{"perm_trans_symmetrize_compact_fc", py_perm_trans_symmetrize_compact_fc,
METH_VARARGS,
"Enforce permutation and translational symmetry of compact force constants"},
{"dynamical_matrix", py_get_dynamical_matrix, METH_VARARGS, "Dynamical matrix"},
{"nac_dynamical_matrix", py_get_nac_dynamical_matrix, METH_VARARGS, "NAC dynamical matrix"},
{"dipole_dipole", py_get_dipole_dipole, METH_VARARGS, "Dipole-dipole interaction"},
{"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_with_mappings", py_distribute_fc2_with_mappings, METH_VARARGS,
{"dynamical_matrix", py_get_dynamical_matrix, METH_VARARGS,
"Dynamical matrix"},
{"nac_dynamical_matrix", py_get_nac_dynamical_matrix, METH_VARARGS,
"NAC dynamical matrix"},
{"dipole_dipole", py_get_dipole_dipole, METH_VARARGS,
"Dipole-dipole interaction"},
{"dipole_dipole_q0", py_get_dipole_dipole_q0, METH_VARARGS,
"q=0 terms of Dipole-dipole interaction"},
{"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_with_mappings", py_distribute_fc2_with_mappings,
METH_VARARGS,
"Distribute force constants for all atoms in atom_list using precomputed symmetry mappings."},
{"compute_permutation", py_compute_permutation, METH_VARARGS,
"Compute indices of original points in a set of rotated points."},
@ -591,6 +600,7 @@ static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args)
static PyObject * py_get_dipole_dipole(PyObject *self, PyObject *args)
{
PyArrayObject* dd_py;
PyArrayObject* dd_q0_py;
PyArrayObject* G_list_py;
PyArrayObject* q_vector_py;
PyArrayObject* q_direction_py;
@ -602,6 +612,7 @@ static PyObject * py_get_dipole_dipole(PyObject *self, PyObject *args)
double tolerance;
double* dd;
double* dd_q0;
double* G_list;
double* q_vector;
double* q_direction;
@ -610,8 +621,9 @@ static PyObject * py_get_dipole_dipole(PyObject *self, PyObject *args)
double *pos;
int num_patom, num_G;
if (!PyArg_ParseTuple(args, "OOOOOOOddd",
if (!PyArg_ParseTuple(args, "OOOOOOOOddd",
&dd_py,
&dd_q0_py,
&G_list_py,
&q_vector_py,
&q_direction_py,
@ -625,6 +637,7 @@ static PyObject * py_get_dipole_dipole(PyObject *self, PyObject *args)
dd = (double*)PyArray_DATA(dd_py);
dd_q0 = (double*)PyArray_DATA(dd_q0_py);
G_list = (double*)PyArray_DATA(G_list_py);
if ((PyObject*)q_direction_py == Py_None) {
q_direction = NULL;
@ -639,6 +652,7 @@ static PyObject * py_get_dipole_dipole(PyObject *self, PyObject *args)
num_patom = PyArray_DIMS(pos_py)[0];
get_dipole_dipole(dd, /* [natom, 3, natom, 3, (real, imag)] */
dd_q0, /* [natom, 3, 3, (real, imag)] */
G_list, /* [num_kvec, 3] */
num_G,
num_patom,
@ -654,7 +668,49 @@ static PyObject * py_get_dipole_dipole(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_get_dipole_dipole_q0(PyObject *self, PyObject *args)
{
PyArrayObject* dd_q0_py;
PyArrayObject* G_list_py;
PyArrayObject* dielectric_py;
PyArrayObject* pos_py;
double lambda;
double tolerance;
double* dd_q0;
double* G_list;
double* dielectric;
double *pos;
int num_patom, num_G;
if (!PyArg_ParseTuple(args, "OOOOdd",
&dd_q0_py,
&G_list_py,
&dielectric_py,
&pos_py,
&lambda,
&tolerance))
return NULL;
dd_q0 = (double*)PyArray_DATA(dd_q0_py);
G_list = (double*)PyArray_DATA(G_list_py);
dielectric = (double*)PyArray_DATA(dielectric_py);
pos = (double*)PyArray_DATA(pos_py);
num_G = PyArray_DIMS(G_list_py)[0];
num_patom = PyArray_DIMS(pos_py)[0];
get_dipole_dipole_q0(dd_q0, /* [natom, 3, 3, (real, imag)] */
G_list, /* [num_kvec, 3] */
num_G,
num_patom,
dielectric,
pos, /* [natom, 3] */
lambda, /* 4 * Lambda^2 */
tolerance);
Py_RETURN_NONE;
}
static PyObject * py_get_derivative_dynmat(PyObject *self, PyObject *args)
{

119
c/harmonic/dynmat.c
View File

@ -147,6 +147,7 @@ int get_dynamical_matrix_at_q(double *dynamical_matrix,
}
void get_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real, imag)] */
const double *dd_q0, /* [natom, 3, 3, (real, imag)] */
const double *G_list, /* [num_G, 3] */
const int num_G,
const int num_patom,
@ -160,16 +161,14 @@ void get_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real, imag)] */
const double tolerance)
{
int i, j, k, l, m, n, adrs, adrs_tmp, adrs_sum;
double zero_vec[3];
double *dd_tmp, *dd_sum;
double zz;
double *dd_tmp;
double zz;
dd_tmp = NULL;
dd_sum = NULL;
dd_tmp = (double*) malloc(sizeof(double) * num_patom * num_patom * 18);
dd_sum = (double*) malloc(sizeof(double) * num_patom * 18);
for (i = 0; i < num_patom * num_patom * 18; i++) {
dd[i] = 0;
dd_tmp[i] = 0;
}
@ -184,59 +183,19 @@ void get_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real, imag)] */
lambda,
tolerance);
for (i = 0; i < num_patom * num_patom * 18; i++) {
dd[i] = dd_tmp[i];
dd_tmp[i] = 0;
}
zero_vec[0] = 0;
zero_vec[1] = 0;
zero_vec[2] = 0;
get_KK(dd_tmp,
G_list,
num_G,
num_patom,
zero_vec,
NULL,
dielectric,
pos,
lambda,
tolerance);
for (i = 0; i < num_patom * 18; i++) {
dd_sum[i] = 0;
}
for (i = 0; i < num_patom; i++) {
for (j = 0; j < num_patom; j++) {
for (k = 0; k < 3; k++) { /* alpha */
for (l = 0; l < 3; l++) { /* beta */
adrs_tmp = i * num_patom * 18 + k * num_patom * 6 + j * 6 + l * 2;
adrs_sum = i * 18 + k * 6 + l * 2;
dd_sum[adrs_sum] += dd_tmp[adrs_tmp];
}
}
}
}
for (i = 0; i < num_patom * num_patom * 18; i++) {
dd_tmp[i] = dd[i];
}
for (i = 0; i < num_patom; i++) {
for (k = 0; k < 3; k++) { /* alpha */
for (l = 0; l < 3; l++) { /* beta */
adrs = i * num_patom * 18 + k * num_patom * 6 + i * 6 + l * 2;
adrs_sum = i * 18 + k * 6 + l * 2;
dd_tmp[adrs] -= dd_sum[adrs_sum];
dd_tmp[adrs] -= dd_q0[adrs_sum];
dd_tmp[adrs + 1] -= dd_q0[adrs_sum + 1];
}
}
}
for (i = 0; i < num_patom * num_patom * 18; i++) {
dd_tmp[i] *= factor;
dd[i] = 0;
}
for (i = 0; i < num_patom; i++) {
@ -259,11 +218,66 @@ void get_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real, imag)] */
free(dd_tmp);
dd_tmp = NULL;
free(dd_sum);
dd_sum = NULL;
}
void get_dipole_dipole_q0(double *dd_q0, /* [natom, 3, 3, (real, imag)] */
const double *G_list, /* [num_G, 3] */
const int num_G,
const int num_patom,
const double *dielectric,
const double *pos, /* [natom, 3] */
const double lambda,
const double tolerance)
{
int i, j, k, l, adrs_tmp, adrs_sum;
double zero_vec[3];
double *dd_tmp;
dd_tmp = NULL;
dd_tmp = (double*) malloc(sizeof(double) * num_patom * num_patom * 18);
for (i = 0; i < num_patom * num_patom * 18; i++) {
dd_tmp[i] = 0;
}
zero_vec[0] = 0;
zero_vec[1] = 0;
zero_vec[2] = 0;
get_KK(dd_tmp,
G_list,
num_G,
num_patom,
zero_vec,
NULL,
dielectric,
pos,
lambda,
tolerance);
for (i = 0; i < num_patom * 18; i++) {
dd_q0[i] = 0;
}
for (i = 0; i < num_patom; i++) {
for (j = 0; j < num_patom; j++) {
for (k = 0; k < 3; k++) { /* alpha */
for (l = 0; l < 3; l++) { /* beta */
adrs_tmp = i * num_patom * 18 + k * num_patom * 6 + j * 6 + l * 2;
adrs_sum = i * 18 + k * 6 + l * 2;
dd_q0[adrs_sum] += dd_tmp[adrs_tmp];
/* expected to be real, though */
dd_q0[adrs_sum +1] += dd_tmp[adrs_tmp + 1];
}
}
}
}
free(dd_tmp);
dd_tmp = NULL;
}
void get_charge_sum(double *charge_sum,
const int num_patom,
const double factor, /* 4pi/V*unit-conv and denominator */
@ -481,7 +495,10 @@ static void get_KK(double *dd_part, /* [natom, 3, natom, 3, (real, imag)] */
for (j = 0; j < num_patom; j++) {
phase = 0;
for (k = 0; k < 3; k++) {
phase += (pos[i * 3 + k] - pos[j * 3 + k]) * q_K[k];
/* For D-type dynamical matrix */
/* phase += (pos[i * 3 + k] - pos[j * 3 + k]) * q_K[k]; */
/* For C-type dynamical matrix */
phase += (pos[i * 3 + k] - pos[j * 3 + k]) * G_list[g * 3 + k];
}
phase *= 2 * PI;
cos_phase = cos(phase);

9
c/harmonic_h/dynmat.h
View File

@ -48,6 +48,7 @@ int get_dynamical_matrix_at_q(double *dynamical_matrix,
const double *charge_sum,
const int with_openmp);
void get_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real, imag)] */
const double *dd_q0, /* [natom, 3, 3, (real, imag)] */
const double *G_list, /* [num_G, 3] */
const int num_G,
const int num_patom,
@ -59,6 +60,14 @@ void get_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real, imag)] */
const double factor, /* 4pi/V*unit-conv */
const double lambda,
const double tolerance);
void get_dipole_dipole_q0(double *dd_q0, /* [natom, 3, 3, (real, imag)] */
const double *G_list, /* [num_G, 3] */
const int num_G,
const int num_patom,
const double *dielectric,
const double *pos, /* [natom, 3] */
const double lambda,
const double tolerance);
void get_charge_sum(double *charge_sum,
const int num_patom,
const double factor,

47
phonopy/harmonic/dynamical_matrix.py
View File

@ -249,6 +249,7 @@ class DynamicalMatrixNAC(DynamicalMatrix):
self._G_vec_list = None
self._G_cutoff = None
self._Lambda = None # 4*Lambda**2 is stored.
self._dd_q0 = None
self._nac = True
if nac_params is not None:
@ -287,6 +288,16 @@ class DynamicalMatrixNAC(DynamicalMatrix):
print("G-cutoff distance: %6.2f" % self._G_cutoff)
print("Number of G-points: %d" % len(self._G_list))
print("Lambda: %6.2f" % self._Lambda)
try:
import phonopy._phonopy as phonoc
self._set_c_dipole_dipole_q0()
except ImportError:
print("Python version of dipole-dipole calculation is not well "
"implemented.")
sys.exit(1)
self._set_Gonze_force_constants()
self._Gonze_count = 0
@ -423,17 +434,26 @@ class DynamicalMatrixNAC(DynamicalMatrix):
import phonopy._phonopy as phonoc
C = self._get_c_dipole_dipole(q, q_dir_cart)
except ImportError:
C = self._get_py_dipole_dipole(q, q_dir_cart)
print("Python version of dipole-dipole calculation is not well "
"implemented.")
sys.exit(1)
# D-type to C-type conversion and mass weighted
# mass = self._pcell.get_masses()
# num_atom = self._pcell.get_number_of_atoms()
# pos = self._pcell.get_positions()
# for i in range(num_atom):
# dpos = pos - pos[i]
# phase_factor = np.exp(2j * np.pi * np.dot(dpos, q))
# for j in range(num_atom):
# C[i, :, j, :] *= phase_factor[j] / np.sqrt(mass[i] * mass[j])
# Mass weighted
mass = self._pcell.get_masses()
num_atom = self._pcell.get_number_of_atoms()
pos = self._pcell.get_positions()
for i in range(num_atom):
dpos = pos - pos[i]
phase_factor = np.exp(2j * np.pi * np.dot(dpos, q))
for j in range(num_atom):
C[i, :, j, :] *= phase_factor[j] / np.sqrt(mass[i] * mass[j])
C[i, :, j, :] *= 1.0 / np.sqrt(mass[i] * mass[j])
C_dd = C.reshape(num_atom * 3, num_atom * 3)
@ -443,12 +463,13 @@ class DynamicalMatrixNAC(DynamicalMatrix):
import phonopy._phonopy as phonoc
pos = self._pcell.get_positions()
num_atom = self._pcell.get_number_of_atoms()
num_atom = len(pos)
volume = self._pcell.get_volume()
C = np.zeros((num_atom, 3, num_atom, 3),
dtype=self._dtype_complex, order='C')
phonoc.dipole_dipole(C.view(dtype='double'),
self._dd_q0.view(dtype='double'),
self._G_list,
q,
q_dir_cart,
@ -460,6 +481,20 @@ class DynamicalMatrixNAC(DynamicalMatrix):
self._symprec)
return C
def _set_c_dipole_dipole_q0(self):
import phonopy._phonopy as phonoc
pos = self._pcell.get_positions()
self._dd_q0 = np.zeros((len(pos), 3, 3),
dtype=self._dtype_complex, order='C')
phonoc.dipole_dipole_q0(self._dd_q0.view(dtype='double'),
self._G_list,
self._dielectric,
np.array(pos, dtype='double', order='C'),
self._Lambda,
self._symprec)
def _get_py_dipole_dipole(self, K_list, q, q_dir_cart):
pos = self._pcell.get_positions()
num_atom = self._pcell.get_number_of_atoms()