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

Untablify

This commit is contained in:
Atsushi Togo 2018-01-03 21:45:04 +09:00
parent 6870fec4b8
commit 61cec3ca26
5 changed files with 493 additions and 493 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

334
c/_phonopy.c
View File

@ -58,15 +58,15 @@ static PyObject * py_compute_permutation(PyObject *self, PyObject *args);
static PyObject * py_gsv_copy_smallest_vectors(PyObject *self, PyObject *args); static PyObject * py_gsv_copy_smallest_vectors(PyObject *self, PyObject *args);
static int distribute_fc2(double *fc2, static int distribute_fc2(double *fc2,
PHPYCONST double lat[3][3], PHPYCONST double lat[3][3],
PHPYCONST double (*pos)[3], PHPYCONST double (*pos)[3],
const int num_pos, const int num_pos,
const int atom_disp, const int atom_disp,
const int map_atom_disp, const int map_atom_disp,
PHPYCONST double r_cart[3][3], PHPYCONST double r_cart[3][3],
PHPYCONST int r[3][3], PHPYCONST int r[3][3],
const double t[3], const double t[3],
const double symprec); const double symprec);
static void distribute_fc2_with_mappings(double (*fc2)[3][3], static void distribute_fc2_with_mappings(double (*fc2)[3][3],
const int * atom_list, const int * atom_list,
@ -79,11 +79,11 @@ static void distribute_fc2_with_mappings(double (*fc2)[3][3],
const int num_pos); const int num_pos);
static int compute_permutation(int * rot_atom, static int compute_permutation(int * rot_atom,
PHPYCONST double lat[3][3], PHPYCONST double lat[3][3],
PHPYCONST double (*pos)[3], PHPYCONST double (*pos)[3],
PHPYCONST double (*rot_pos)[3], PHPYCONST double (*rot_pos)[3],
const int num_pos, const int num_pos,
const double symprec); const double symprec);
static void gsv_copy_smallest_vectors(double (*shortest_vectors)[27][3], static void gsv_copy_smallest_vectors(double (*shortest_vectors)[27][3],
int * multiplicity, int * multiplicity,
@ -105,11 +105,11 @@ static PyObject * py_get_tetrahedra_frequenies(PyObject *self, PyObject *args);
static PyObject * py_tetrahedron_method_dos(PyObject *self, PyObject *args); static PyObject * py_tetrahedron_method_dos(PyObject *self, PyObject *args);
static double get_free_energy_omega(const double temperature, static double get_free_energy_omega(const double temperature,
const double omega); const double omega);
static double get_entropy_omega(const double temperature, static double get_entropy_omega(const double temperature,
const double omega); const double omega);
static double get_heat_capacity_omega(const double temperature, static double get_heat_capacity_omega(const double temperature,
const double omega); const double omega);
/* static double get_energy_omega(double temperature, double omega); */ /* static double get_energy_omega(double temperature, double omega); */
static int check_overlap(PHPYCONST double (*pos)[3], static int check_overlap(PHPYCONST double (*pos)[3],
const int num_pos, const int num_pos,
@ -334,14 +334,14 @@ static PyObject * py_get_dynamical_matrix(PyObject *self, PyObject *args)
int num_satom; int num_satom;
if (!PyArg_ParseTuple(args, "OOOOOOOO", if (!PyArg_ParseTuple(args, "OOOOOOOO",
&dynamical_matrix, &dynamical_matrix,
&force_constants, &force_constants,
&q_vector, &q_vector,
&r_vector, &r_vector,
&multiplicity, &multiplicity,
&mass, &mass,
&super2prim_map, &super2prim_map,
&prim2super_map)) &prim2super_map))
return NULL; return NULL;
dm = (double*)PyArray_DATA(dynamical_matrix); dm = (double*)PyArray_DATA(dynamical_matrix);
@ -356,17 +356,17 @@ static PyObject * py_get_dynamical_matrix(PyObject *self, PyObject *args)
num_satom = PyArray_DIMS(super2prim_map)[0]; num_satom = PyArray_DIMS(super2prim_map)[0];
get_dynamical_matrix_at_q(dm, get_dynamical_matrix_at_q(dm,
num_patom, num_patom,
num_satom, num_satom,
fc, fc,
q, q,
r, r,
multi, multi,
m, m,
s2p_map, s2p_map,
p2s_map, p2s_map,
NULL, NULL,
1); 1);
Py_RETURN_NONE; Py_RETURN_NONE;
} }
@ -403,17 +403,17 @@ static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args)
double *charge_sum; double *charge_sum;
if (!PyArg_ParseTuple(args, "OOOOOOOOOOd", if (!PyArg_ParseTuple(args, "OOOOOOOOOOd",
&dynamical_matrix, &dynamical_matrix,
&force_constants, &force_constants,
&q_vector, &q_vector,
&r_vector, &r_vector,
&multiplicity, &multiplicity,
&mass, &mass,
&super2prim_map, &super2prim_map,
&prim2super_map, &prim2super_map,
&q_cart_vector, &q_cart_vector,
&born, &born,
&factor)) &factor))
return NULL; return NULL;
dm = (double*)PyArray_DATA(dynamical_matrix); dm = (double*)PyArray_DATA(dynamical_matrix);
@ -434,17 +434,17 @@ static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args)
get_charge_sum(charge_sum, num_patom, factor / n, q_cart, z); get_charge_sum(charge_sum, num_patom, factor / n, q_cart, z);
get_dynamical_matrix_at_q(dm, get_dynamical_matrix_at_q(dm,
num_patom, num_patom,
num_satom, num_satom,
fc, fc,
q, q,
r, r,
multi, multi,
m, m,
s2p_map, s2p_map,
p2s_map, p2s_map,
charge_sum, charge_sum,
1); 1);
free(charge_sum); free(charge_sum);
@ -473,14 +473,14 @@ static PyObject * py_get_dipole_dipole(PyObject *self, PyObject *args)
int num_patom, num_K; int num_patom, num_K;
if (!PyArg_ParseTuple(args, "OOOOOOOdd", if (!PyArg_ParseTuple(args, "OOOOOOOdd",
&dd_py, &dd_py,
&K_list_py, &K_list_py,
&q_vector_py, &q_vector_py,
&q_direction_py, &q_direction_py,
&born_py, &born_py,
&dielectric_py, &dielectric_py,
&pos_py, &pos_py,
&factor, &factor,
&tolerance)) &tolerance))
return NULL; return NULL;
@ -549,19 +549,19 @@ static PyObject * py_get_derivative_dynmat(PyObject *self, PyObject *args)
double *q_dir; double *q_dir;
if (!PyArg_ParseTuple(args, "OOOOOOOOOdOOO", if (!PyArg_ParseTuple(args, "OOOOOOOOOdOOO",
&derivative_dynmat, &derivative_dynmat,
&force_constants, &force_constants,
&q_vector, &q_vector,
&lattice, /* column vectors */ &lattice, /* column vectors */
&r_vector, &r_vector,
&multiplicity, &multiplicity,
&mass, &mass,
&super2prim_map, &super2prim_map,
&prim2super_map, &prim2super_map,
&nac_factor, &nac_factor,
&born, &born,
&dielectric, &dielectric,
&q_direction)) { &q_direction)) {
return NULL; return NULL;
} }
@ -594,20 +594,20 @@ static PyObject * py_get_derivative_dynmat(PyObject *self, PyObject *args)
} }
get_derivative_dynmat_at_q(ddm, get_derivative_dynmat_at_q(ddm,
num_patom, num_patom,
num_satom, num_satom,
fc, fc,
q, q,
lat, lat,
r, r,
multi, multi,
m, m,
s2p_map, s2p_map,
p2s_map, p2s_map,
nac_factor, nac_factor,
z, z,
epsilon, epsilon,
q_dir); q_dir);
Py_RETURN_NONE; Py_RETURN_NONE;
} }
@ -635,9 +635,9 @@ static PyObject * py_get_thermal_properties(PyObject *self, PyObject *args)
if (!PyArg_ParseTuple(args, "OOOO", if (!PyArg_ParseTuple(args, "OOOO",
&thermal_props_py, &thermal_props_py,
&temperatures_py, &temperatures_py,
&frequencies_py, &frequencies_py,
&weights_py)) { &weights_py)) {
return NULL; return NULL;
} }
@ -698,7 +698,7 @@ static PyObject * py_get_thermal_properties(PyObject *self, PyObject *args)
} }
static double get_free_energy_omega(const double temperature, static double get_free_energy_omega(const double temperature,
const double omega) const double omega)
{ {
/* temperature is defined by T (K) */ /* temperature is defined by T (K) */
/* omega must be normalized to eV. */ /* omega must be normalized to eV. */
@ -706,7 +706,7 @@ static double get_free_energy_omega(const double temperature,
} }
static double get_entropy_omega(const double temperature, static double get_entropy_omega(const double temperature,
const double omega) const double omega)
{ {
/* temperature is defined by T (K) */ /* temperature is defined by T (K) */
/* omega must be normalized to eV. */ /* omega must be normalized to eV. */
@ -717,7 +717,7 @@ static double get_entropy_omega(const double temperature,
} }
static double get_heat_capacity_omega(const double temperature, static double get_heat_capacity_omega(const double temperature,
const double omega) const double omega)
{ {
/* temperature is defined by T (K) */ /* temperature is defined by T (K) */
/* omega must be normalized to eV. */ /* omega must be normalized to eV. */
@ -881,15 +881,15 @@ static PyObject * py_distribute_fc2(PyObject *self, PyObject *args)
num_pos = PyArray_DIMS(positions_py)[0]; num_pos = PyArray_DIMS(positions_py)[0];
distribute_fc2(fc2, distribute_fc2(fc2,
lat, lat,
pos, pos,
num_pos, num_pos,
atom_disp, atom_disp,
map_atom_disp, map_atom_disp,
r_cart, r_cart,
r, r,
t, t,
symprec); symprec);
Py_RETURN_NONE; Py_RETURN_NONE;
} }
@ -1070,12 +1070,12 @@ static PyObject *py_thm_neighboring_grid_points(PyObject *self, PyObject *args)
int *bz_map; int *bz_map;
if (!PyArg_ParseTuple(args, "OiOOOO", if (!PyArg_ParseTuple(args, "OiOOOO",
&relative_grid_points_py, &relative_grid_points_py,
&grid_point, &grid_point,
&relative_grid_address_py, &relative_grid_address_py,
&mesh_py, &mesh_py,
&bz_grid_address_py, &bz_grid_address_py,
&bz_map_py)) { &bz_map_py)) {
return NULL; return NULL;
} }
@ -1087,12 +1087,12 @@ static PyObject *py_thm_neighboring_grid_points(PyObject *self, PyObject *args)
bz_map = (int*)PyArray_DATA(bz_map_py); bz_map = (int*)PyArray_DATA(bz_map_py);
thm_get_neighboring_grid_points(relative_grid_points, thm_get_neighboring_grid_points(relative_grid_points,
grid_point, grid_point,
relative_grid_address, relative_grid_address,
num_relative_grid_address, num_relative_grid_address,
mesh, mesh,
bz_grid_address, bz_grid_address,
bz_map); bz_map);
Py_RETURN_NONE; Py_RETURN_NONE;
} }
@ -1106,8 +1106,8 @@ py_thm_relative_grid_address(PyObject *self, PyObject *args)
double (*reciprocal_lattice)[3]; double (*reciprocal_lattice)[3];
if (!PyArg_ParseTuple(args, "OO", if (!PyArg_ParseTuple(args, "OO",
&relative_grid_address_py, &relative_grid_address_py,
&reciprocal_lattice_py)) { &reciprocal_lattice_py)) {
return NULL; return NULL;
} }
@ -1127,7 +1127,7 @@ py_thm_all_relative_grid_address(PyObject *self, PyObject *args)
int (*relative_grid_address)[24][4][3]; int (*relative_grid_address)[24][4][3];
if (!PyArg_ParseTuple(args, "O", if (!PyArg_ParseTuple(args, "O",
&relative_grid_address_py)) { &relative_grid_address_py)) {
return NULL; return NULL;
} }
@ -1150,9 +1150,9 @@ py_thm_integration_weight(PyObject *self, PyObject *args)
double iw; double iw;
if (!PyArg_ParseTuple(args, "dOs", if (!PyArg_ParseTuple(args, "dOs",
&omega, &omega,
&tetrahedra_omegas_py, &tetrahedra_omegas_py,
&function)) { &function)) {
return NULL; return NULL;
} }
@ -1179,10 +1179,10 @@ py_thm_integration_weight_at_omegas(PyObject *self, PyObject *args)
double (*tetrahedra_omegas)[4]; double (*tetrahedra_omegas)[4];
if (!PyArg_ParseTuple(args, "OOOs", if (!PyArg_ParseTuple(args, "OOOs",
&integration_weights_py, &integration_weights_py,
&omegas_py, &omegas_py,
&tetrahedra_omegas_py, &tetrahedra_omegas_py,
&function)) { &function)) {
return NULL; return NULL;
} }
@ -1192,10 +1192,10 @@ py_thm_integration_weight_at_omegas(PyObject *self, PyObject *args)
tetrahedra_omegas = (double(*)[4])PyArray_DATA(tetrahedra_omegas_py); tetrahedra_omegas = (double(*)[4])PyArray_DATA(tetrahedra_omegas_py);
thm_get_integration_weight_at_omegas(iw, thm_get_integration_weight_at_omegas(iw,
num_omegas, num_omegas,
omegas, omegas,
tetrahedra_omegas, tetrahedra_omegas,
function[0]); function[0]);
Py_RETURN_NONE; Py_RETURN_NONE;
} }
@ -1226,13 +1226,13 @@ static PyObject * py_get_tetrahedra_frequenies(PyObject *self, PyObject *args)
int address_double[3]; int address_double[3];
if (!PyArg_ParseTuple(args, "OOOOOOO", if (!PyArg_ParseTuple(args, "OOOOOOO",
&freq_tetras_py, &freq_tetras_py,
&grid_points_py, &grid_points_py,
&mesh_py, &mesh_py,
&grid_address_py, &grid_address_py,
&gp_ir_index_py, &gp_ir_index_py,
&relative_grid_address_py, &relative_grid_address_py,
&frequencies_py)) { &frequencies_py)) {
return NULL; return NULL;
} }
@ -1250,16 +1250,16 @@ static PyObject * py_get_tetrahedra_frequenies(PyObject *self, PyObject *args)
#pragma omp parallel for private(k, g_addr, gp, address_double) #pragma omp parallel for private(k, g_addr, gp, address_double)
for (j = 0; j < num_band * 96; j++) { for (j = 0; j < num_band * 96; j++) {
for (k = 0; k < 3; k++) { for (k = 0; k < 3; k++) {
g_addr[k] = grid_address[grid_points[i]][k] + g_addr[k] = grid_address[grid_points[i]][k] +
relative_grid_address[j % 96][k]; relative_grid_address[j % 96][k];
} }
kgd_get_grid_address_double_mesh(address_double, kgd_get_grid_address_double_mesh(address_double,
g_addr, g_addr,
mesh, mesh,
is_shift); is_shift);
gp = kgd_get_grid_point_double_mesh(address_double, mesh); gp = kgd_get_grid_point_double_mesh(address_double, mesh);
freq_tetras[i * num_band * 96 + j] = freq_tetras[i * num_band * 96 + j] =
frequencies[gp_ir_index[gp] * num_band + j / 96]; frequencies[gp_ir_index[gp] * num_band + j / 96];
} }
} }
@ -1305,14 +1305,14 @@ static PyObject * py_tetrahedron_method_dos(PyObject *self, PyObject *args)
weights = NULL; weights = NULL;
if (!PyArg_ParseTuple(args, "OOOOOOOO", if (!PyArg_ParseTuple(args, "OOOOOOOO",
&dos_py, &dos_py,
&mesh_py, &mesh_py,
&freq_points_py, &freq_points_py,
&frequencies_py, &frequencies_py,
&coef_py, &coef_py,
&grid_address_py, &grid_address_py,
&grid_mapping_table_py, &grid_mapping_table_py,
&relative_grid_address_py)) { &relative_grid_address_py)) {
return NULL; return NULL;
} }
@ -1397,15 +1397,15 @@ static PyObject * py_tetrahedron_method_dos(PyObject *self, PyObject *args)
} }
static int distribute_fc2(double *fc2, static int distribute_fc2(double *fc2,
PHPYCONST double lat[3][3], PHPYCONST double lat[3][3],
PHPYCONST double (*pos)[3], PHPYCONST double (*pos)[3],
const int num_pos, const int num_pos,
const int atom_disp, const int atom_disp,
const int map_atom_disp, const int map_atom_disp,
PHPYCONST double r_cart[3][3], PHPYCONST double r_cart[3][3],
PHPYCONST int r[3][3], PHPYCONST int r[3][3],
const double t[3], const double t[3],
const double symprec) const double symprec)
{ {
int i, j, k, l, m, address_new, address; int i, j, k, l, m, address_new, address;
int is_found, rot_atom; int is_found, rot_atom;
@ -1431,13 +1431,13 @@ static int distribute_fc2(double *fc2,
address_new = atom_disp * num_pos * 9 + i * 9; address_new = atom_disp * num_pos * 9 + i * 9;
for (j = 0; j < 3; j++) { for (j = 0; j < 3; j++) {
for (k = 0; k < 3; k++) { for (k = 0; k < 3; k++) {
for (l = 0; l < 3; l++) { for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) { for (m = 0; m < 3; m++) {
fc2[address_new + j * 3 + k] += fc2[address_new + j * 3 + k] +=
r_cart[l][j] * r_cart[m][k] * r_cart[l][j] * r_cart[m][k] *
fc2[address + l * 3 + m]; fc2[address + l * 3 + m];
} }
} }
} }
} }
end: end:

344
c/harmonic/derivative_dynmat.c
View File

@ -37,46 +37,46 @@
#define PI 3.14159265358979323846 #define PI 3.14159265358979323846
static void get_derivative_nac(double *ddnac, static void get_derivative_nac(double *ddnac,
double *dnac, double *dnac,
const int num_patom, const int num_patom,
const double *lattice, const double *lattice,
const double *mass, const double *mass,
const double *q, const double *q,
const double *born, const double *born,
const double *dielectric, const double *dielectric,
const double *q_direction, const double *q_direction,
const double factor); const double factor);
static double get_A(const int atom_i, static double get_A(const int atom_i,
const int cart_i, const int cart_i,
const double q[3], const double q[3],
const double *born); const double *born);
static double get_C(const double q[3], static double get_C(const double q[3],
const double *dielectric); const double *dielectric);
static double get_dA(const int atom_i, static double get_dA(const int atom_i,
const int cart_i, const int cart_i,
const int cart_j, const int cart_j,
const double *born); const double *born);
static double get_dC(const int cart_i, static double get_dC(const int cart_i,
const int cart_j, const int cart_j,
const int cart_k, const int cart_k,
const double q[3], const double q[3],
const double *dielectric); const double *dielectric);
void get_derivative_dynmat_at_q(double *derivative_dynmat, void get_derivative_dynmat_at_q(double *derivative_dynmat,
const int num_patom, const int num_patom,
const int num_satom, const int num_satom,
const double *fc, const double *fc,
const double *q, const double *q,
const double *lattice, /* column vector */ const double *lattice, /* column vector */
const double *r, const double *r,
const int *multi, const int *multi,
const double *mass, const double *mass,
const int *s2p_map, const int *s2p_map,
const int *p2s_map, const int *p2s_map,
const double nac_factor, const double nac_factor,
const double *born, const double *born,
const double *dielectric, const double *dielectric,
const double *q_direction) const double *q_direction)
{ {
int i, j, k, l, m, n, adrs, adrsT, is_nac; int i, j, k, l, m, n, adrs, adrsT, is_nac;
double coef[3], real_coef[3], imag_coef[3]; double coef[3], real_coef[3], imag_coef[3];
@ -88,15 +88,15 @@ void get_derivative_dynmat_at_q(double *derivative_dynmat,
is_nac = 1; is_nac = 1;
if (q_direction) { if (q_direction) {
if (fabs(q_direction[0]) < 1e-5 && if (fabs(q_direction[0]) < 1e-5 &&
fabs(q_direction[1]) < 1e-5 && fabs(q_direction[1]) < 1e-5 &&
fabs(q_direction[2]) < 1e-5) { fabs(q_direction[2]) < 1e-5) {
is_nac = 0; is_nac = 0;
} }
} else { } else {
if (fabs(q[0]) < 1e-5 && if (fabs(q[0]) < 1e-5 &&
fabs(q[1]) < 1e-5 && fabs(q[1]) < 1e-5 &&
fabs(q[2]) < 1e-5) { fabs(q[2]) < 1e-5) {
is_nac = 0; is_nac = 0;
} }
} }
} else { } else {
@ -108,15 +108,15 @@ void get_derivative_dynmat_at_q(double *derivative_dynmat,
dnac = (double*) malloc(sizeof(double) * num_patom * num_patom * 9); dnac = (double*) malloc(sizeof(double) * num_patom * num_patom * 9);
factor = nac_factor * num_patom / num_satom; factor = nac_factor * num_patom / num_satom;
get_derivative_nac(ddnac, get_derivative_nac(ddnac,
dnac, dnac,
num_patom, num_patom,
lattice, lattice,
mass, mass,
q, q,
born, born,
dielectric, dielectric,
q_direction, q_direction,
factor); factor);
} }
for (i = 0; i < num_patom; i++) { for (i = 0; i < num_patom; i++) {
@ -124,90 +124,90 @@ void get_derivative_dynmat_at_q(double *derivative_dynmat,
mass_sqrt = sqrt(mass[i] * mass[j]); mass_sqrt = sqrt(mass[i] * mass[j]);
for (k = 0; k < 3; k++) { for (k = 0; k < 3; k++) {
for (l = 0; l < 3; l++) { for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) { for (m = 0; m < 3; m++) {
ddm_real[m][k][l] = 0; ddm_real[m][k][l] = 0;
ddm_imag[m][k][l] = 0; ddm_imag[m][k][l] = 0;
} }
} }
} }
for (k = 0; k < num_satom; k++) { /* Lattice points of right index of fc */ for (k = 0; k < num_satom; k++) { /* Lattice points of right index of fc */
if (s2p_map[k] != p2s_map[j]) { if (s2p_map[k] != p2s_map[j]) {
continue; continue;
} }
real_phase = 0; real_phase = 0;
imag_phase = 0; imag_phase = 0;
for (l = 0; l < 3; l++) { for (l = 0; l < 3; l++) {
real_coef[l] = 0; real_coef[l] = 0;
imag_coef[l] = 0; imag_coef[l] = 0;
} }
for (l = 0; l < multi[k * num_patom + i]; l++) { for (l = 0; l < multi[k * num_patom + i]; l++) {
phase = 0; phase = 0;
for (m = 0; m < 3; m++) { for (m = 0; m < 3; m++) {
phase += q[m] * r[k * num_patom * 81 + i * 81 + l * 3 + m]; phase += q[m] * r[k * num_patom * 81 + i * 81 + l * 3 + m];
} }
s = sin(phase * 2 * PI); s = sin(phase * 2 * PI);
c = cos(phase * 2 * PI); c = cos(phase * 2 * PI);
real_phase += c; real_phase += c;
imag_phase += s; imag_phase += s;
for (m = 0; m < 3; m++) { for (m = 0; m < 3; m++) {
coef[m] = 0; coef[m] = 0;
for (n = 0; n < 3; n++) { for (n = 0; n < 3; n++) {
coef[m] += 2 * PI * coef[m] += 2 * PI *
lattice[m * 3 + n] * r[k * num_patom * 81 + i * 81 + l * 3 + n]; lattice[m * 3 + n] * r[k * num_patom * 81 + i * 81 + l * 3 + n];
} }
} }
for (m =0; m < 3; m++) { for (m =0; m < 3; m++) {
real_coef[m] -= coef[m] * s; real_coef[m] -= coef[m] * s;
imag_coef[m] += coef[m] * c; imag_coef[m] += coef[m] * c;
} }
} }
real_phase /= multi[k * num_patom + i]; real_phase /= multi[k * num_patom + i];
imag_phase /= multi[k * num_patom + i]; imag_phase /= multi[k * num_patom + i];
for (l = 0; l < 3; l++) { for (l = 0; l < 3; l++) {
real_coef[l] /= multi[k * num_patom + i]; real_coef[l] /= multi[k * num_patom + i];
imag_coef[l] /= multi[k * num_patom + i]; imag_coef[l] /= multi[k * num_patom + i];
} }
for (l = 0; l < 3; l++) { for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) { for (m = 0; m < 3; m++) {
fc_elem = fc[p2s_map[i] * num_satom * 9 + k * 9 + l * 3 + m] fc_elem = fc[p2s_map[i] * num_satom * 9 + k * 9 + l * 3 + m]
/ mass_sqrt; / mass_sqrt;
if (is_nac) { if (is_nac) {
fc_elem += dnac[i * 9 * num_patom + j * 9 + l * 3 + m]; fc_elem += dnac[i * 9 * num_patom + j * 9 + l * 3 + m];
} }
for (n = 0; n < 3; n++) { for (n = 0; n < 3; n++) {
ddm_real[n][l][m] += fc_elem * real_coef[n]; ddm_real[n][l][m] += fc_elem * real_coef[n];
ddm_imag[n][l][m] += fc_elem * imag_coef[n]; ddm_imag[n][l][m] += fc_elem * imag_coef[n];
if (is_nac) { if (is_nac) {
ddm_real[n][l][m] += ddm_real[n][l][m] +=
ddnac[n * num_patom * num_patom * 9 + ddnac[n * num_patom * num_patom * 9 +
i * 9 * num_patom + j * 9 + l * 3 + m] * real_phase; i * 9 * num_patom + j * 9 + l * 3 + m] * real_phase;
ddm_imag[n][l][m] += ddm_imag[n][l][m] +=
ddnac[n * num_patom * num_patom * 9 + ddnac[n * num_patom * num_patom * 9 +
i * 9 * num_patom + j * 9 + l * 3 + m] * imag_phase; i * 9 * num_patom + j * 9 + l * 3 + m] * imag_phase;
} }
} }
} }
} }
} }
for (k = 0; k < 3; k++) { for (k = 0; k < 3; k++) {
for (l = 0; l < 3; l++) { for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) { for (m = 0; m < 3; m++) {
adrs = (k * num_patom * num_patom * 18 + adrs = (k * num_patom * num_patom * 18 +
(i * 3 + l) * num_patom * 6 + j * 6 + m * 2); (i * 3 + l) * num_patom * 6 + j * 6 + m * 2);
derivative_dynmat[adrs] += ddm_real[k][l][m]; derivative_dynmat[adrs] += ddm_real[k][l][m];
derivative_dynmat[adrs + 1] += ddm_imag[k][l][m]; derivative_dynmat[adrs + 1] += ddm_imag[k][l][m];
} }
} }
} }
} }
} }
@ -216,14 +216,14 @@ void get_derivative_dynmat_at_q(double *derivative_dynmat,
for (i = 0; i < 3; i++) { for (i = 0; i < 3; i++) {
for (j = i; j < num_patom * 3; j++) { for (j = i; j < num_patom * 3; j++) {
for (k = 0; k < num_patom * 3; k++) { for (k = 0; k < num_patom * 3; k++) {
adrs = i * num_patom * num_patom * 18 + j * num_patom * 6 + k * 2; adrs = i * num_patom * num_patom * 18 + j * num_patom * 6 + k * 2;
adrsT = i * num_patom * num_patom * 18 + k * num_patom * 6 + j * 2; adrsT = i * num_patom * num_patom * 18 + k * num_patom * 6 + j * 2;
derivative_dynmat[adrs] += derivative_dynmat[adrsT]; derivative_dynmat[adrs] += derivative_dynmat[adrsT];
derivative_dynmat[adrs] /= 2; derivative_dynmat[adrs] /= 2;
derivative_dynmat[adrs + 1] -= derivative_dynmat[adrsT+ 1]; derivative_dynmat[adrs + 1] -= derivative_dynmat[adrsT+ 1];
derivative_dynmat[adrs + 1] /= 2; derivative_dynmat[adrs + 1] /= 2;
derivative_dynmat[adrsT] = derivative_dynmat[adrs]; derivative_dynmat[adrsT] = derivative_dynmat[adrs];
derivative_dynmat[adrsT + 1] = -derivative_dynmat[adrs + 1]; derivative_dynmat[adrsT + 1] = -derivative_dynmat[adrs + 1];
} }
} }
} }
@ -238,15 +238,15 @@ void get_derivative_dynmat_at_q(double *derivative_dynmat,
/* D_nac = a * AB/C */ /* D_nac = a * AB/C */
/* dD_nac = a * D_nac * (A'/A + B'/B - C'/C) */ /* dD_nac = a * D_nac * (A'/A + B'/B - C'/C) */
static void get_derivative_nac(double *ddnac, static void get_derivative_nac(double *ddnac,
double *dnac, double *dnac,
const int num_patom, const int num_patom,
const double *lattice, const double *lattice,
const double *mass, const double *mass,
const double *q, const double *q,
const double *born, const double *born,
const double *dielectric, const double *dielectric,
const double *q_direction, const double *q_direction,
const double factor) const double factor)
{ {
int i, j, k, l, m; int i, j, k, l, m;
double a, b, c, da, db, dc, volume, mass_sqrt; double a, b, c, da, db, dc, volume, mass_sqrt;
@ -271,9 +271,9 @@ static void get_derivative_nac(double *ddnac,
q_cart[i] = 0; q_cart[i] = 0;
for (j = 0; j < 3; j++) { for (j = 0; j < 3; j++) {
if (q_direction) { if (q_direction) {
q_cart[i] += rec_lat[i * 3 + j] * q_direction[j]; q_cart[i] += rec_lat[i * 3 + j] * q_direction[j];
} else { } else {
q_cart[i] += rec_lat[i * 3 + j] * q[j]; q_cart[i] += rec_lat[i * 3 + j] * q[j];
} }
} }
} }
@ -284,31 +284,31 @@ static void get_derivative_nac(double *ddnac,
for (j = 0; j < num_patom; j++) { /* atom_j */ for (j = 0; j < num_patom; j++) { /* atom_j */
mass_sqrt = sqrt(mass[i] * mass[j]); mass_sqrt = sqrt(mass[i] * mass[j]);
for (k = 0; k < 3; k++) { /* derivative direction */ for (k = 0; k < 3; k++) { /* derivative direction */
for (l = 0; l < 3; l++) { /* alpha */ for (l = 0; l < 3; l++) { /* alpha */
a = get_A(i, l, q_cart, born); a = get_A(i, l, q_cart, born);
da = get_dA(i, l, k, born); da = get_dA(i, l, k, born);
for (m = 0; m < 3; m++) { /* beta */ for (m = 0; m < 3; m++) { /* beta */
b = get_A(j, m, q_cart, born); b = get_A(j, m, q_cart, born);
db = get_dA(j, m, k, born); db = get_dA(j, m, k, born);
dc = get_dC(l, m, k, q_cart, dielectric); dc = get_dC(l, m, k, q_cart, dielectric);
ddnac[k * num_patom * num_patom * 9 + ddnac[k * num_patom * num_patom * 9 +
i * 9 * num_patom + j * 9 + l * 3 + m] = i * 9 * num_patom + j * 9 + l * 3 + m] =
(da * b + db * a - a * b * dc / c) / (c * mass_sqrt) * factor; (da * b + db * a - a * b * dc / c) / (c * mass_sqrt) * factor;
if (k == 0) { if (k == 0) {
dnac[i * 9 * num_patom + j * 9 + l * 3 + m] = dnac[i * 9 * num_patom + j * 9 + l * 3 + m] =
a * b / (c * mass_sqrt) * factor; a * b / (c * mass_sqrt) * factor;
} }
} }
} }
} }
} }
} }
} }
static double get_A(const int atom_i, static double get_A(const int atom_i,
const int cart_i, const int cart_i,
const double q[3], const double q[3],
const double *born) const double *born)
{ {
int i; int i;
double sum; double sum;
@ -322,7 +322,7 @@ static double get_A(const int atom_i,
} }
static double get_C(const double q[3], static double get_C(const double q[3],
const double *dielectric) const double *dielectric)
{ {
int i, j; int i, j;
double sum; double sum;
@ -338,33 +338,33 @@ static double get_C(const double q[3],
} }
static double get_dA(const int atom_i, static double get_dA(const int atom_i,
const int cart_i, const int cart_i,
const int cart_j, const int cart_j,
const double *born) const double *born)
{ {
return born[atom_i * 9 + cart_j * 3 + cart_i]; return born[atom_i * 9 + cart_j * 3 + cart_i];
} }
static double get_dC(const int cart_i, static double get_dC(const int cart_i,
const int cart_j, const int cart_j,
const int cart_k, const int cart_k,
const double q[3], const double q[3],
const double *dielectric) const double *dielectric)
{ {
if (cart_k == 0) { if (cart_k == 0) {
return (2 * q[0] * dielectric[0] + return (2 * q[0] * dielectric[0] +
q[1] * (dielectric[1] + dielectric[3]) + q[1] * (dielectric[1] + dielectric[3]) +
q[2] * (dielectric[2] + dielectric[6])); q[2] * (dielectric[2] + dielectric[6]));
} }
if (cart_k == 1) { if (cart_k == 1) {
return (2 * q[1] * dielectric[4] + return (2 * q[1] * dielectric[4] +
q[2] * (dielectric[5] + dielectric[7]) + q[2] * (dielectric[5] + dielectric[7]) +
q[0] * (dielectric[1] + dielectric[3])); q[0] * (dielectric[1] + dielectric[3]));
} }
if (cart_k == 2) { if (cart_k == 2) {
return (2 * q[2] * dielectric[8] + return (2 * q[2] * dielectric[8] +
q[0] * (dielectric[2] + dielectric[6]) + q[0] * (dielectric[2] + dielectric[6]) +
q[1] * (dielectric[5] + dielectric[7])); q[1] * (dielectric[5] + dielectric[7]));
} }
return 0; return 0;
} }

226
c/harmonic/dynmat.c
View File

@ -38,47 +38,47 @@
#define PI 3.14159265358979323846 #define PI 3.14159265358979323846
static void get_dynmat_ij(double *dynamical_matrix, static void get_dynmat_ij(double *dynamical_matrix,
const int num_patom, const int num_patom,
const int num_satom, const int num_satom,
const double *fc, const double *fc,
const double *q, const double *q,
const double *r, const double *r,
const int *multi, const int *multi,
const double *mass, const double *mass,
const int *s2p_map, const int *s2p_map,
const int *p2s_map, const int *p2s_map,
const double *charge_sum, const double *charge_sum,
const int i, const int i,
const int j); const int j);
static void get_dm(double dm_real[3][3], static void get_dm(double dm_real[3][3],
double dm_imag[3][3], double dm_imag[3][3],
const int num_patom, const int num_patom,
const int num_satom, const int num_satom,
const double *fc, const double *fc,
const double *q, const double *q,
const double *r, const double *r,
const int *multi, const int *multi,
const double mass_sqrt, const double mass_sqrt,
const int *p2s_map, const int *p2s_map,
const double *charge_sum, const double *charge_sum,
const int i, const int i,
const int j, const int j,
const int k); const int k);
static double get_dielectric_part(const double q[3], static double get_dielectric_part(const double q[3],
const double *dielectric); const double *dielectric);
int get_dynamical_matrix_at_q(double *dynamical_matrix, int get_dynamical_matrix_at_q(double *dynamical_matrix,
const int num_patom, const int num_patom,
const int num_satom, const int num_satom,
const double *fc, const double *fc,
const double *q, const double *q,
const double *r, const double *r,
const int *multi, const int *multi,
const double *mass, const double *mass,
const int *s2p_map, const int *s2p_map,
const int *p2s_map, const int *p2s_map,
const double *charge_sum, const double *charge_sum,
const int with_openmp) const int with_openmp)
{ {
int i, j, ij, adrs, adrsT; int i, j, ij, adrs, adrsT;
@ -86,35 +86,35 @@ int get_dynamical_matrix_at_q(double *dynamical_matrix,
#pragma omp parallel for #pragma omp parallel for
for (ij = 0; ij < num_patom * num_patom ; ij++) { for (ij = 0; ij < num_patom * num_patom ; ij++) {
get_dynmat_ij(dynamical_matrix, get_dynmat_ij(dynamical_matrix,
num_patom, num_patom,
num_satom, num_satom,
fc, fc,
q, q,
r, r,
multi, multi,
mass, mass,
s2p_map, s2p_map,
p2s_map, p2s_map,
charge_sum, charge_sum,
ij / num_patom, /* i */ ij / num_patom, /* i */
ij % num_patom); /* j */ ij % num_patom); /* j */
} }
} else { } else {
for (i = 0; i < num_patom; i++) { for (i = 0; i < num_patom; i++) {
for (j = 0; j < num_patom; j++) { for (j = 0; j < num_patom; j++) {
get_dynmat_ij(dynamical_matrix, get_dynmat_ij(dynamical_matrix,
num_patom, num_patom,
num_satom, num_satom,
fc, fc,
q, q,
r, r,
multi, multi,
mass, mass,
s2p_map, s2p_map,
p2s_map, p2s_map,
charge_sum, charge_sum,
i, i,
j); j);
} }
} }
} }
@ -247,10 +247,10 @@ void get_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real, imag)] */
} }
void get_charge_sum(double *charge_sum, void get_charge_sum(double *charge_sum,
const int num_patom, const int num_patom,
const double factor, /* 4pi/V*unit-conv and denominator */ const double factor, /* 4pi/V*unit-conv and denominator */
const double q_vector[3], const double q_vector[3],
const double *born) const double *born)
{ {
int i, j, k, a, b; int i, j, k, a, b;
double (*q_born)[3]; double (*q_born)[3];
@ -265,7 +265,7 @@ void get_charge_sum(double *charge_sum,
for (i = 0; i < num_patom; i++) { for (i = 0; i < num_patom; i++) {
for (j = 0; j < 3; j++) { for (j = 0; j < 3; j++) {
for (k = 0; k < 3; k++) { for (k = 0; k < 3; k++) {
q_born[i][j] += q_vector[k] * born[i * 9 + k * 3 + j]; q_born[i][j] += q_vector[k] * born[i * 9 + k * 3 + j];
} }
} }
} }
@ -273,10 +273,10 @@ void get_charge_sum(double *charge_sum,
for (i = 0; i < num_patom; i++) { for (i = 0; i < num_patom; i++) {
for (j = 0; j < num_patom; j++) { for (j = 0; j < num_patom; j++) {
for (a = 0; a < 3; a++) { for (a = 0; a < 3; a++) {
for (b = 0; b < 3; b++) { for (b = 0; b < 3; b++) {
charge_sum[i * 9 * num_patom + j * 9 + a * 3 + b] = charge_sum[i * 9 * num_patom + j * 9 + a * 3 + b] =
q_born[i][a] * q_born[j][b] * factor; q_born[i][a] * q_born[j][b] * factor;
} }
} }
} }
} }
@ -287,18 +287,18 @@ void get_charge_sum(double *charge_sum,
static void get_dynmat_ij(double *dynamical_matrix, static void get_dynmat_ij(double *dynamical_matrix,
const int num_patom, const int num_patom,
const int num_satom, const int num_satom,
const double *fc, const double *fc,
const double *q, const double *q,
const double *r, const double *r,
const int *multi, const int *multi,
const double *mass, const double *mass,
const int *s2p_map, const int *s2p_map,
const int *p2s_map, const int *p2s_map,
const double *charge_sum, const double *charge_sum,
const int i, const int i,
const int j) const int j)
{ {
int k, l, adrs; int k, l, adrs;
double mass_sqrt; double mass_sqrt;
@ -318,19 +318,19 @@ static void get_dynmat_ij(double *dynamical_matrix,
continue; continue;
} }
get_dm(dm_real, get_dm(dm_real,
dm_imag, dm_imag,
num_patom, num_patom,
num_satom, num_satom,
fc, fc,
q, q,
r, r,
multi, multi,
mass_sqrt, mass_sqrt,
p2s_map, p2s_map,
charge_sum, charge_sum,
i, i,
j, j,
k); k);
} }
for (k = 0; k < 3; k++) { for (k = 0; k < 3; k++) {
@ -343,19 +343,19 @@ static void get_dynmat_ij(double *dynamical_matrix,
} }
static void get_dm(double dm_real[3][3], static void get_dm(double dm_real[3][3],
double dm_imag[3][3], double dm_imag[3][3],
const int num_patom, const int num_patom,
const int num_satom, const int num_satom,
const double *fc, const double *fc,
const double *q, const double *q,
const double *r, const double *r,
const int *multi, const int *multi,
const double mass_sqrt, const double mass_sqrt,
const int *p2s_map, const int *p2s_map,
const double *charge_sum, const double *charge_sum,
const int i, const int i,
const int j, const int j,
const int k) const int k)
{ {
int l, m; int l, m;
double phase, cos_phase, sin_phase, fc_elem; double phase, cos_phase, sin_phase, fc_elem;
@ -375,12 +375,12 @@ static void get_dm(double dm_real[3][3],
for (l = 0; l < 3; l++) { for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) { for (m = 0; m < 3; m++) {
if (charge_sum) { if (charge_sum) {
fc_elem = (fc[p2s_map[i] * num_satom * 9 + k * 9 + l * 3 + m] + fc_elem = (fc[p2s_map[i] * num_satom * 9 + k * 9 + l * 3 + m] +
charge_sum[i * num_patom * 9 + charge_sum[i * num_patom * 9 +
j * 9 + l * 3 + m]) / mass_sqrt; j * 9 + l * 3 + m]) / mass_sqrt;
} else { } else {
fc_elem = fc[p2s_map[i] * num_satom * 9 + fc_elem = fc[p2s_map[i] * num_satom * 9 +
k * 9 + l * 3 + m] / mass_sqrt; k * 9 + l * 3 + m] / mass_sqrt;
} }
dm_real[l][m] += fc_elem * cos_phase; dm_real[l][m] += fc_elem * cos_phase;
dm_imag[l][m] += fc_elem * sin_phase; dm_imag[l][m] += fc_elem * sin_phase;

28
c/harmonic_h/derivative_dynmat.h
View File

@ -36,19 +36,19 @@
#define __derivative_dynmat_H__ #define __derivative_dynmat_H__
void get_derivative_dynmat_at_q(double *derivative_dynmat, void get_derivative_dynmat_at_q(double *derivative_dynmat,
const int num_patom, const int num_patom,
const int num_satom, const int num_satom,
const double *fc, const double *fc,
const double *q, const double *q,
const double *lattice, /* column vector */ const double *lattice, /* column vector */
const double *r, const double *r,
const int *multi, const int *multi,
const double *mass, const double *mass,
const int *s2p_map, const int *s2p_map,
const int *p2s_map, const int *p2s_map,
const double nac_factor, const double nac_factor,
const double *born, const double *born,
const double *dielectric, const double *dielectric,
const double *q_direction); const double *q_direction);
#endif #endif

30
c/harmonic_h/dynmat.h
View File

@ -36,17 +36,17 @@
#define __dynmat_H__ #define __dynmat_H__
int get_dynamical_matrix_at_q(double *dynamical_matrix, int get_dynamical_matrix_at_q(double *dynamical_matrix,
const int num_patom, const int num_patom,
const int num_satom, const int num_satom,
const double *fc, const double *fc,
const double *q, const double *q,
const double *r, const double *r,
const int *multi, const int *multi,
const double *mass, const double *mass,
const int *s2p_map, const int *s2p_map,
const int *p2s_map, const int *p2s_map,
const double *charge_sum, const double *charge_sum,
const int with_openmp); const int with_openmp);
void get_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real, imag)] */ void get_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real, imag)] */
const double *K_list, /* [num_kvec, 3] */ const double *K_list, /* [num_kvec, 3] */
const int num_K, const int num_K,
@ -59,8 +59,8 @@ void get_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real, imag)] */
const double *pos, /* [natom, 3] */ const double *pos, /* [natom, 3] */
const double tolerance); const double tolerance);
void get_charge_sum(double *charge_sum, void get_charge_sum(double *charge_sum,
const int num_patom, const int num_patom,
const double factor, const double factor,
const double q_vector[3], const double q_vector[3],
const double *born); const double *born);
#endif #endif